Skip to content
Closed
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion README.md
Original file line number Diff line number Diff line change
Expand Up @@ -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
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
Expand Down
2 changes: 1 addition & 1 deletion docs/src/examples/example.md
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
2 changes: 1 addition & 1 deletion example/dc_motor_lqg_design.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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)

Expand Down
13 changes: 9 additions & 4 deletions src/ControlSystems.jl
Original file line number Diff line number Diff line change
Expand Up @@ -15,11 +15,12 @@ export LTISystem,
isproper,
# Linear Algebra
balance,
care,
dare,
dlyap,
arec,
ared,
lyapd,
lqr,
dlqr,
lqrc,
lqrd,
kalman,
dkalman,
lqg,
Expand Down Expand Up @@ -170,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
Expand Down
18 changes: 9 additions & 9 deletions src/matrix_comps.jl
Original file line number Diff line number Diff line change
@@ -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.
Expand All @@ -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
Expand All @@ -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
Expand All @@ -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
Expand All @@ -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)
Expand All @@ -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
Expand Down Expand Up @@ -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
Expand Down
32 changes: 16 additions & 16 deletions src/synthesis.jl
Original file line number Diff line number Diff line change
@@ -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:
Expand Down Expand Up @@ -33,8 +33,8 @@ 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)
S = care(A, B, Q, R)
function lqrc(A, B, Q, R)
S = arec(A, B, Q, R)
K = R\B'*S
return K
end
Expand All @@ -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:
Expand All @@ -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
Expand All @@ -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)
S = dare(A, B, Q, R)
function lqrd(A, B, Q, R)
S = ared(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

"""
Expand All @@ -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)
Expand Down
6 changes: 3 additions & 3 deletions src/types/lqg.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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)
```

Expand Down Expand Up @@ -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
Expand All @@ -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)

Expand Down
2 changes: 1 addition & 1 deletion test/test_discrete.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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)

Expand Down
20 changes: 10 additions & 10 deletions test/test_linalg.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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]
Expand Down
2 changes: 1 addition & 1 deletion test/test_plots.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
10 changes: 5 additions & 5 deletions test/test_synthesis.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down