From 6c32788322609fb1196145e04ebf61edd11738a5 Mon Sep 17 00:00:00 2001 From: Fredrik Bagge Carlson Date: Wed, 20 Oct 2021 12:57:11 +0200 Subject: [PATCH 1/4] lqr -> lqrc, dlqr -> lqrd and similar lyap is not called lyapc since it's not our function --- README.md | 2 +- docs/src/examples/example.md | 2 +- example/dc_motor_lqg_design.jl | 2 +- src/ControlSystems.jl | 4 ++-- src/matrix_comps.jl | 8 ++++---- src/synthesis.jl | 28 ++++++++++++++-------------- src/types/lqg.jl | 6 +++--- test/test_plots.jl | 2 +- test/test_synthesis.jl | 10 +++++----- 9 files changed, 32 insertions(+), 32 deletions(-) diff --git a/README.md b/README.md index 03f446e4e..7beabf6ae 100644 --- a/README.md +++ b/README.md @@ -68,7 +68,7 @@ ss, tf, zpk ##### Analysis pole, tzero, norm, hinfnorm, linfnorm, ctrb, obsv, gangoffour, margin, markovparam, damp, dampreport, zpkdata, dcgain, covar, gram, sigma, sisomargin ##### Synthesis -care, dare, dlyap, lqr, dlqr, place, leadlink, laglink, leadlinkat, rstd, rstc, dab, balreal, baltrunc +care, dare, lyapd, lqr, lqrc, lqrd, place, leadlink, laglink, leadlinkat, rstd, rstc, dab, balreal, baltrunc ###### PID design pid, stabregionPID, loopshapingPI, pidplots ##### Time and Frequency response diff --git a/docs/src/examples/example.md b/docs/src/examples/example.md index 18b02122d..3958fa2d3 100644 --- a/docs/src/examples/example.md +++ b/docs/src/examples/example.md @@ -19,7 +19,7 @@ C = [1 0] sys = ss(A,B,C,0, Ts) Q = I R = I -L = dlqr(A,B,Q,R) # lqr(sys,Q,R) can also be used +L = lqrd(A,B,Q,R) # lqr(sys,Q,R) can also be used u(x,t) = -L*x .+ 1.5(t>=2.5)# Form control law (u is a function of t and x), a constant input disturbance is affecting the system from t≧2.5 t =0:Ts:5 diff --git a/example/dc_motor_lqg_design.jl b/example/dc_motor_lqg_design.jl index 53857f9c1..8e6c9a0f7 100644 --- a/example/dc_motor_lqg_design.jl +++ b/example/dc_motor_lqg_design.jl @@ -37,7 +37,7 @@ Q = [1. 0; 0 1] R = 20. -K = lqr(p60.A, p60.B, Q, R) +K = lqrc(p60.A, p60.B, Q, R) # needs to be modified if Nbar is not a scalar Nbar = 1. ./ (p60.D - (p60.C - p60.D*K) * inv(p60.A - p60.B*K) * p60.B) diff --git a/src/ControlSystems.jl b/src/ControlSystems.jl index 949500aa1..979d23309 100644 --- a/src/ControlSystems.jl +++ b/src/ControlSystems.jl @@ -17,9 +17,9 @@ export LTISystem, balance, care, dare, - dlyap, + lyapd, lqr, - dlqr, + lqrd, kalman, dkalman, lqg, diff --git a/src/matrix_comps.jl b/src/matrix_comps.jl index f62610361..37a7300b9 100644 --- a/src/matrix_comps.jl +++ b/src/matrix_comps.jl @@ -65,12 +65,12 @@ function dare(A, B, Q, R) return QZ.Z[(n+1):end, 1:n]/QZ.Z[1:n, 1:n]; end -"""`dlyap(A, Q)` +"""`lyapd(A, Q)` Compute the solution `X` to the discrete Lyapunov equation `AXA' - X + Q = 0`. """ -function dlyap(A::AbstractMatrix, Q) +function lyapd(A::AbstractMatrix, Q) lhs = kron(A, conj(A)) lhs = I - lhs x = lhs\reshape(Q, prod(size(Q)), 1) @@ -86,7 +86,7 @@ function gram(sys::AbstractStateSpace, opt::Symbol) if !isstable(sys) error("gram only valid for stable A") end - func = iscontinuous(sys) ? lyap : dlyap + func = iscontinuous(sys) ? lyap : lyapd if opt === :c # TODO probably remove type check in julia 0.7.0 return func(sys.A, sys.B*sys.B')#::Array{numeric_type(sys),2} # lyap is type-unstable @@ -165,7 +165,7 @@ function covar(sys::AbstractStateSpace, W) if !isstable(sys) return fill(Inf,(size(C,1),size(C,1))) end - func = iscontinuous(sys) ? lyap : dlyap + func = iscontinuous(sys) ? lyap : lyapd Q = try func(A, B*W*B') catch e diff --git a/src/synthesis.jl b/src/synthesis.jl index 0394bb08e..4984f388c 100644 --- a/src/synthesis.jl +++ b/src/synthesis.jl @@ -1,5 +1,5 @@ """ - lqr(A, B, Q, R) + lqrc(A, B, Q, R) Calculate the optimal gain matrix `K` for the state-feedback law `u = -K*x` that minimizes the cost function: @@ -33,7 +33,7 @@ y, t, x, uout = lsim(sys,u,t,x0=x0) plot(t,x', lab=["Position" "Velocity"], xlabel="Time [s]") ``` """ -function lqr(A, B, Q, R) +function lqrc(A, B, Q, R) S = care(A, B, Q, R) K = R\B'*S return K @@ -47,28 +47,28 @@ Calculate the optimal Kalman gain See also `LQG` """ -kalman(A, C, R1,R2) = Matrix(lqr(A',C',R1,R2)') +kalman(A, C, R1,R2) = Matrix(lqrc(A',C',R1,R2)') function lqr(sys::AbstractStateSpace, Q, R) if iscontinuous(sys) - return lqr(sys.A, sys.B, Q, R) + return lqrc(sys.A, sys.B, Q, R) else - return dlqr(sys.A, sys.B, Q, R) + return lqrd(sys.A, sys.B, Q, R) end end function kalman(sys::AbstractStateSpace, R1,R2) if iscontinuous(sys) - return Matrix(lqr(sys.A', sys.C', R1,R2)') + return Matrix(lqrc(sys.A', sys.C', R1,R2)') else - return Matrix(dlqr(sys.A', sys.C', R1,R2)') + return Matrix(lqrd(sys.A', sys.C', R1,R2)') end end """ - dlqr(A, B, Q, R) - dlqr(sys, Q, R) + lqrd(A, B, Q, R) + lqrd(sys, Q, R) Calculate the optimal gain matrix `K` for the state-feedback law `u[k] = -K*x[k]` that minimizes the cost function: @@ -89,7 +89,7 @@ C = [1 0] sys = ss(A,B,C,0, Ts) Q = I R = I -L = dlqr(A,B,Q,R) # lqr(sys,Q,R) can also be used +L = lqrd(A,B,Q,R) # lqr(sys,Q,R) can also be used u(x,t) = -L*x # Form control law, t=0:Ts:5 @@ -98,15 +98,15 @@ y, t, x, uout = lsim(sys,u,t,x0=x0) plot(t,x', lab=["Position" "Velocity"], xlabel="Time [s]") ``` """ -function dlqr(A, B, Q, R) +function lqrd(A, B, Q, R) S = dare(A, B, Q, R) K = (B'*S*B + R)\(B'S*A) return K end -function dlqr(sys::AbstractStateSpace, Q, R) +function lqrd(sys::AbstractStateSpace, Q, R) !isdiscrete(sys) && throw(ArgumentError("Input argument sys must be discrete-time system")) - return dlqr(sys.A, sys.B, Q, R) + return lqrd(sys.A, sys.B, Q, R) end """ @@ -116,7 +116,7 @@ end Calculate the optimal Kalman gain for discrete time systems """ -dkalman(A, C, R1,R2) = Matrix(dlqr(A',C',R1,R2)') +dkalman(A, C, R1,R2) = Matrix(lqrd(A',C',R1,R2)') """ place(A, B, p, opt=:c) diff --git a/src/types/lqg.jl b/src/types/lqg.jl index 26cce312a..df377789a 100644 --- a/src/types/lqg.jl +++ b/src/types/lqg.jl @@ -16,7 +16,7 @@ described by `sys`. `qQ` and `qR` can be set to incorporate loop transfer recovery, i.e., ```julia -L = lqr(A, B, Q1+qQ*C'C, Q2) +L = lqrc(A, B, Q1+qQ*C'C, Q2) K = kalman(A, C, R1+qR*B*B', R2) ``` @@ -131,7 +131,7 @@ function _LQG(sys::LTISystem, Q1, Q2, R1, R2, qQ, qR; M = sys.C) n = size(A, 1) m = size(B, 2) p = size(C, 1) - L = lqr(A, B, Q1 + qQ * C'C, Q2) + L = lqrc(A, B, Q1 + qQ * C'C, Q2) K = kalman(A, C, R1 + qR * B * B', R2) # Controller system @@ -158,7 +158,7 @@ function _LQGi(sys::LTISystem, Q1, Q2, R1, R2, qQ, qR; M = sys.C) Ce = [C zeros(p, m)] De = D - L = lqr(A, B, Q1 + qQ * C'C, Q2) + L = lqrc(A, B, Q1 + qQ * C'C, Q2) Le = [L I] K = kalman(Ae, Ce, R1 + qR * Be * Be', R2) diff --git a/test/test_plots.jl b/test/test_plots.jl index 630be4b91..4e3c00159 100644 --- a/test/test_plots.jl +++ b/test/test_plots.jl @@ -22,7 +22,7 @@ function getexamples() stepgen() = stepplot(sys, ts[end], l=(:dash, 4)) impulsegen() = impulseplot(sys, ts[end], l=:blue) - L = lqr(sysss.A, sysss.B, [1 0; 0 1], [1 0; 0 1]) + L = lqrc(sysss.A, sysss.B, [1 0; 0 1], [1 0; 0 1]) lsimgen() = lsimplot(sysss, (x,i)->-L*x, ts; x0=[1;2]) margingen() = marginplot([tf1, tf2], ws) diff --git a/test/test_synthesis.jl b/test/test_synthesis.jl index 95ac316b6..c96e82e17 100644 --- a/test/test_synthesis.jl +++ b/test/test_synthesis.jl @@ -108,28 +108,28 @@ end C = [1 0] Q = I R = I - L = dlqr(A,B,Q,R) + L = lqrd(A,B,Q,R) @test L ≈ [0.5890881713787511 0.7118839434795103] sys = ss(A,B,C,0,Ts) L = lqr(sys, Q, R) @test L ≈ [0.5890881713787511 0.7118839434795103] - L = dlqr(sys, Q, R) + L = lqrd(sys, Q, R) @test L ≈ [0.5890881713787511 0.7118839434795103] B = reshape(B,2,1) # Note B is matrix, B'B is compatible with I - L = dlqr(A,B,Q,R) + L = lqrd(A,B,Q,R) @test L ≈ [0.5890881713787511 0.7118839434795103] Q = eye_(2) R = eye_(1) - L = dlqr(A,B,Q,R) + L = lqrd(A,B,Q,R) @test L ≈ [0.5890881713787511 0.7118839434795103] B = [0;1] # Note B is vector, B'B is scalar and INcompatible with matrix Q = eye_(2) R = eye_(1) - @test_throws MethodError L ≈ dlqr(A,B,Q,R) + @test_throws MethodError L ≈ lqrd(A,B,Q,R) #L ≈ [0.5890881713787511 0.7118839434795103] end From 4af569357e9ce1a2ec9521c2de4215df4aae8fd9 Mon Sep 17 00:00:00 2001 From: Fredrik Bagge Carlson Date: Wed, 20 Oct 2021 13:00:26 +0200 Subject: [PATCH 2/4] care/dare -> arec/ared --- README.md | 2 +- src/ControlSystems.jl | 4 ++-- src/matrix_comps.jl | 10 +++++----- src/synthesis.jl | 4 ++-- test/test_discrete.jl | 2 +- test/test_linalg.jl | 20 ++++++++++---------- 6 files changed, 21 insertions(+), 21 deletions(-) diff --git a/README.md b/README.md index 7beabf6ae..70430658a 100644 --- a/README.md +++ b/README.md @@ -68,7 +68,7 @@ ss, tf, zpk ##### Analysis pole, tzero, norm, hinfnorm, linfnorm, ctrb, obsv, gangoffour, margin, markovparam, damp, dampreport, zpkdata, dcgain, covar, gram, sigma, sisomargin ##### Synthesis -care, dare, lyapd, lqr, lqrc, lqrd, place, leadlink, laglink, leadlinkat, rstd, rstc, dab, balreal, baltrunc +arec, ared, lyapd, lqr, lqrc, lqrd, place, leadlink, laglink, leadlinkat, rstd, rstc, dab, balreal, baltrunc ###### PID design pid, stabregionPID, loopshapingPI, pidplots ##### Time and Frequency response diff --git a/src/ControlSystems.jl b/src/ControlSystems.jl index 979d23309..69af03b5e 100644 --- a/src/ControlSystems.jl +++ b/src/ControlSystems.jl @@ -15,8 +15,8 @@ export LTISystem, isproper, # Linear Algebra balance, - care, - dare, + arec, + ared, lyapd, lqr, lqrd, diff --git a/src/matrix_comps.jl b/src/matrix_comps.jl index 37a7300b9..1ad1eed7c 100644 --- a/src/matrix_comps.jl +++ b/src/matrix_comps.jl @@ -1,4 +1,4 @@ -"""`care(A, B, Q, R)` +"""`arec(A, B, Q, R)` Compute 'X', the solution to the continuous-time algebraic Riccati equation, defined as A'X + XA - (XB)R^-1(B'X) + Q = 0, where R is non-singular. @@ -7,12 +7,12 @@ Algorithm taken from: Laub, "A Schur Method for Solving Algebraic Riccati Equations." http://dspace.mit.edu/bitstream/handle/1721.1/1301/R-0859-05666488.pdf """ -function care(A, B, Q, R) +function arec(A, B, Q, R) G = try B*inv(R)*B' catch y if y isa SingularException - error("R must be non-singular in care.") + error("R must be non-singular in arec.") else throw(y) end @@ -31,7 +31,7 @@ function care(A, B, Q, R) return U21/U11 end -"""`dare(A, B, Q, R)` +"""`ared(A, B, Q, R)` Compute `X`, the solution to the discrete-time algebraic Riccati equation, defined as A'XA - X - (A'XB)(B'XB + R)^-1(B'XA) + Q = 0, where Q>=0 and R>0 @@ -40,7 +40,7 @@ Algorithm taken from: Laub, "A Schur Method for Solving Algebraic Riccati Equations." http://dspace.mit.edu/bitstream/handle/1721.1/1301/R-0859-05666488.pdf """ -function dare(A, B, Q, R) +function ared(A, B, Q, R) if !issemiposdef(Q) error("Q must be positive-semidefinite."); end diff --git a/src/synthesis.jl b/src/synthesis.jl index 4984f388c..dd8936e85 100644 --- a/src/synthesis.jl +++ b/src/synthesis.jl @@ -34,7 +34,7 @@ plot(t,x', lab=["Position" "Velocity"], xlabel="Time [s]") ``` """ function lqrc(A, B, Q, R) - S = care(A, B, Q, R) + S = arec(A, B, Q, R) K = R\B'*S return K end @@ -99,7 +99,7 @@ plot(t,x', lab=["Position" "Velocity"], xlabel="Time [s]") ``` """ function lqrd(A, B, Q, R) - S = dare(A, B, Q, R) + S = ared(A, B, Q, R) K = (B'*S*B + R)\(B'S*A) return K end diff --git a/test/test_discrete.jl b/test/test_discrete.jl index 19b89b7fb..641412600 100644 --- a/test/test_discrete.jl +++ b/test/test_discrete.jl @@ -136,7 +136,7 @@ Bm = conv(B⁺, B⁻) # In this case, keep the entire numerator polynomial of t R,S,T = rstc(B⁺,B⁻,A,Bm,Am,Ao,AR) # Calculate the 2-DOF controller polynomials -Gcl = tf(conv(B,T),zpconv(A,R,B,S)) # Form the closed loop polynomial from reference to output, the closed-loop characteristic polynomial is AR + BS, the function zpconv takes care of the polynomial multiplication and makes sure the coefficient vectores are of equal length +Gcl = tf(conv(B,T),zpconv(A,R,B,S)) # Form the closed loop polynomial from reference to output, the closed-loop characteristic polynomial is AR + BS, the function zpconv takes arec of the polynomial multiplication and makes sure the coefficient vectores are of equal length @test ControlSystems.isstable(Gcl) diff --git a/test/test_linalg.jl b/test/test_linalg.jl index 484511c7e..d4540f528 100644 --- a/test/test_linalg.jl +++ b/test/test_linalg.jl @@ -4,34 +4,34 @@ b = [0 1]' c = [1 -1] r = 3 -@test care(a, b, c'*c, r) ≈ [0.5895174372762604 1.8215747248860816; 1.8215747248860823 8.818839806923107] -@test dare(a, b, c'*c, r) ≈ [240.73344504496302 -131.09928700089387; -131.0992870008943 75.93413176750603] +@test arec(a, b, c'*c, r) ≈ [0.5895174372762604 1.8215747248860816; 1.8215747248860823 8.818839806923107] +@test ared(a, b, c'*c, r) ≈ [240.73344504496302 -131.09928700089387; -131.0992870008943 75.93413176750603] -## Test dare method for non invertible matrices +## Test ared method for non invertible matrices A = [0 1; 0 0]; B0 = [0; 1]; Q = Matrix{Float64}(I, 2, 2); R0 = 1 -X = dare(A, B0, Q, R0); +X = ared(A, B0, Q, R0); # Reshape for matrix expression B = reshape(B0, 2, 1) R = fill(R0, 1, 1) @test norm(A'X*A - X - (A'X*B)*((B'X*B + R)\(B'X*A)) + Q) < 1e-14 -## Test dare for scalars +## Test ared for scalars A = 1.0; B = 1.0; Q = 1.0; R = 1.0; -X0 = dare(A, B, Q, R); +X0 = ared(A, B, Q, R); X = X0[1] @test norm(A'X*A - X - (A'X*B)*((B'X*B + R)\(B'X*A)) + Q) < 1e-14 # Test for complex matrices I2 = Matrix{Float64}(I, 2, 2) -@test dare([1.0 im; im 1.0], I2, I2, I2) ≈ [1 + sqrt(2) 0; 0 1 + sqrt(2)] +@test ared([1.0 im; im 1.0], I2, I2, I2) ≈ [1 + sqrt(2) 0; 0 1 + sqrt(2)] # And complex scalars -@test dare(1.0, 1, 1, 1) ≈ fill((1 + sqrt(5))/2, 1, 1) -@test dare(1.0im, 1, 1, 1) ≈ fill((1 + sqrt(5))/2, 1, 1) -@test dare(1.0, 1im, 1, 1) ≈ fill((1 + sqrt(5))/2, 1, 1) +@test ared(1.0, 1, 1, 1) ≈ fill((1 + sqrt(5))/2, 1, 1) +@test ared(1.0im, 1, 1, 1) ≈ fill((1 + sqrt(5))/2, 1, 1) +@test ared(1.0, 1im, 1, 1) ≈ fill((1 + sqrt(5))/2, 1, 1) ## Test gram, ctrb, obsv a_2 = [-5 -3; 2 -9] From b72d550c23c12065f7db9a1d2c6307cc1baae567 Mon Sep 17 00:00:00 2001 From: Fredrik Bagge Carlson Date: Wed, 20 Oct 2021 13:01:33 +0200 Subject: [PATCH 3/4] export new lqrc --- src/ControlSystems.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/src/ControlSystems.jl b/src/ControlSystems.jl index 69af03b5e..87f4d4064 100644 --- a/src/ControlSystems.jl +++ b/src/ControlSystems.jl @@ -19,6 +19,7 @@ export LTISystem, ared, lyapd, lqr, + lqrc, lqrd, kalman, dkalman, From 1b0b60cd506d42dbd2400a062f74407c9699628a Mon Sep 17 00:00:00 2001 From: Fredrik Bagge Carlson Date: Wed, 20 Oct 2021 13:06:52 +0200 Subject: [PATCH 4/4] add deprecations --- src/ControlSystems.jl | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/ControlSystems.jl b/src/ControlSystems.jl index 87f4d4064..76c5a6f8c 100644 --- a/src/ControlSystems.jl +++ b/src/ControlSystems.jl @@ -171,6 +171,10 @@ include("delay_systems.jl") include("plotting.jl") +@deprecate care arec +@deprecate dare ared +@deprecate dlqr lqrd +@deprecate lqr(A,B,R1,R2) lqrc(A,B,R1,R2) @deprecate num numvec @deprecate den denvec @deprecate norminf hinfnorm