use timeevol types to dispatch matrix equations

README.md

@@ -50,19 +50,19 @@ A documentation website is available at [http://juliacontrol.github.io/ControlSy
 
 Some of the available commands are:
 ##### Constructing systems
-ss, tf, zpk
+`ss, tf, zpk`
 ##### Analysis
-poles, tzeros, norm, hinfnorm, linfnorm, ctrb, obsv, gangoffour, margin, markovparam, damp, dampreport, zpkdata, dcgain, covar, gram, sigma, sisomargin
+`poles, tzeros, 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
+`are, lyap, lqr, place, leadlink, laglink, leadlinkat, rstd, rstc, dab, balreal, baltrunc`
 ###### PID design
-pid, stabregionPID, loopshapingPI, pidplots
+`pid, stabregionPID, loopshapingPI, pidplots`
 ##### Time and Frequency response
-step, impulse, lsim, freqresp, evalfr, bode, nyquist
+`step, impulse, lsim, freqresp, evalfr, bode, nyquist`
 ##### Plotting
-lsimplot, stepplot, impulseplot, bodeplot, nyquistplot, sigmaplot, marginplot, gangoffourplot, pidplots, pzmap, nicholsplot, pidplots, rlocus, leadlinkcurve
+`lsimplot, stepplot, impulseplot, bodeplot, nyquistplot, sigmaplot, marginplot, gangoffourplot, pidplots, pzmap, nicholsplot, pidplots, rlocus, leadlinkcurve`
 ##### Other
-minreal, sminreal, c2d
+`minreal, sminreal, c2d`
 ## Usage
 
 This toolbox works similar to that of other major computer-aided control

docs/src/examples/example.md

@@ -21,7 +21,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       = lqr(Discrete,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

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 = lqr(p60, 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)
 

src/ControlSystems.jl

@@ -15,13 +15,9 @@ export  LTISystem,
         # Linear Algebra
         balance,
         balance_statespace,
-        care,
-        dare,
-        dlyap,
+        are,
         lqr,
-        dlqr,
         kalman,
-        dkalman,
         lqg,
         lqgi,
         covar,

src/matrix_comps.jl

@@ -1,42 +1,52 @@
-"""`care(A, B, Q, R)`
+"""
+    are(::Continuous, 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.
 
 Uses `MatrixEquations.arec`.
 """
-function care(A, B, Q, R)
+function are(::ContinuousType, A::AbstractMatrix, B, Q, R)
     arec(A, B, R, Q)[1]
 end
 
-care(A::Number, B::Number, Q::Number, R::Number) = care(fill(A,1,1),fill(B,1,1),fill(Q,1,1),fill(R,1,1))
-
 """
-    dare(A, B, Q, R; kwargs...)
+    are(::Discrete, A, B, Q, R; kwargs...)
 
 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
 
 Uses `MatrixEquations.ared`. For keyword arguments, see the docstring of `ControlSystems.MatrixEquations.ared`
 """
-function dare(A, B, Q, R; kwargs...)
+function are(::DiscreteType, A::AbstractMatrix, B, Q, R; kwargs...)
     ared(A, B, R, Q; kwargs...)[1]
 end
 
-dare(A::Number, B::Number, Q::Number, R::Number) = dare(fill(A,1,1),fill(B,1,1),fill(Q,1,1),fill(R,1,1))
+are(t::TimeEvolType, A::Number, B::Number, Q::Number, R::Number) = are(t, fill(A,1,1),fill(B,1,1),fill(Q,1,1),fill(R,1,1))
+
+@deprecate care(args...; kwargs...) are(Continuous, args...; kwargs...)
+@deprecate dare(args...; kwargs...) are(Discrete, args...; kwargs...)
 
 """
-    dlyap(A, Q; kwargs...)
+    lyap(A, Q; kwargs...)
 
 Compute the solution `X` to the discrete Lyapunov equation
 `AXA' - X + Q = 0`.
 
 Uses `MatrixEquations.lyapd`. For keyword arguments, see the docstring of `ControlSystems.MatrixEquations.lyapd`
 """
-function dlyap(A::AbstractMatrix, Q; kwargs...)
+function lyap(::DiscreteType, A::AbstractMatrix, Q; kwargs...)
     lyapd(A, Q; kwargs...)
 end
 
+LinearAlgebra.lyap(::ContinuousType, args...; kwargs...) = lyap(args...; kwargs...)
+LinearAlgebra.lyap(::DiscreteType, args...; kwargs...) = dlyap(args...; kwargs...)
+
+plyap(::ContinuousType, args...; kwargs...) = MatrixEquations.plyapc(args...; kwargs...)
+plyap(::DiscreteType, args...; kwargs...) = MatrixEquations.plyapd(args...; kwargs...)
+
+@deprecate dlyap(args...; kwargs...) lyap(Discrete, args...; kwargs...)
+
 
 """
     U = grampd(sys, opt; kwargs...)
@@ -53,11 +63,10 @@ function grampd(sys::AbstractStateSpace, opt::Symbol; kwargs...)
     if !isstable(sys)
         error("gram only valid for stable A")
     end
-    func = iscontinuous(sys) ? MatrixEquations.plyapc : MatrixEquations.plyapd
     if opt === :c
-        func(sys.A, sys.B; kwargs...)
+        plyap(sys.timeevol, sys.A, sys.B; kwargs...)
     elseif opt === :o
-        func(sys.A', sys.C'; kwargs...)
+        plyap(sys.timeevol, sys.A', sys.C'; kwargs...)
     else
         error("opt must be either :c for controllability grammian, or :o for
                 observability grammian")
@@ -147,10 +156,9 @@ function covar(sys::AbstractStateSpace, W)
     if !isstable(sys)
         return fill(Inf,(size(C,1),size(C,1)))
     end
-    func = iscontinuous(sys) ? plyapc : plyapd
     Wc = cholesky(W).L
     Q1 = sys.nx == 0 ? B*Wc : try
-        func(A, B*Wc)
+        plyap(sys.timeevol, A, B*Wc)
     catch e
         @error("No solution to the Lyapunov equation was found in covar.")
         rethrow(e)
@@ -178,7 +186,7 @@ covar(sys::TransferFunction, W) = covar(ss(sys), W)
 # Note: the H∞ norm computation is probably not as accurate as with SLICOT,
 # but this seems to be still reasonably ok as a first step
 """
-`..  norm(sys, p=2; tol=1e-6)`
+    norm(sys, p=2; tol=1e-6)
 
 `norm(sys)` or `norm(sys,2)` computes the H2 norm of the LTI system `sys`.
 
@@ -206,7 +214,7 @@ LinearAlgebra.norm(sys::TransferFunction, p::Real=2; tol=1e-6) = norm(ss(sys), p
 
 
 """
-`   (Ninf, ω_peak) = hinfnorm(sys; tol=1e-6)`
+    Ninf, ω_peak = hinfnorm(sys; tol=1e-6)
 
 Compute the H∞ norm `Ninf` of the LTI system `sys`, together with a frequency
 `ω_peak` at which the gain Ninf is achieved.
@@ -234,7 +242,7 @@ hinfnorm(sys::AbstractStateSpace{<:Discrete}; tol=1e-6) = _infnorm_two_steps_dt(
 hinfnorm(sys::TransferFunction; tol=1e-6) = hinfnorm(ss(sys); tol=tol)
 
 """
-`   (Ninf, ω_peak) = linfnorm(sys; tol=1e-6)`
+    Ninf, ω_peak = linfnorm(sys; tol=1e-6)
 
 Compute the L∞ norm `Ninf` of the LTI system `sys`, together with a frequency
 `ω_peak` at which the gain `Ninf` is achieved.

src/synthesis.jl

@@ -1,101 +1,51 @@
 """
-    lqr(A, B, Q, R, args...; kwargs...)
+    lqr(sys, Q, R)
+    lqr(Continuous, A, B, Q, R, args...; kwargs...)
+    lqr(Discrete, A, B, Q, R, args...; kwargs...)
 
 Calculate the optimal gain matrix `K` for the state-feedback law `u = -K*x` that
 minimizes the cost function:
 
-J = integral(x'Qx + u'Ru, 0, inf).
-
-For the continuous time model `dx = Ax + Bu`.
-
-`lqr(sys, Q, R)`
+J = integral(x'Qx + u'Ru, 0, inf) for the continuous-time model `dx = Ax + Bu`.
+J = sum(x'Qx + u'Ru, 0, inf) for the discrete-time model `x[k+1] = Ax[k] + Bu[k]`.
 
 Solve the LQR problem for state-space system `sys`. Works for both discrete
 and continuous time systems.
 
-The `args...; kwargs...` are sent to the Riccati solver, allowing specification of cross-covariance etc. See `?MatrixEquations.arec` for more help.
+The `args...; kwargs...` are sent to the Riccati solver, allowing specification of cross-covariance etc. See `?MatrixEquations.arec / ared` for more help.
 
-See also `LQG`
-Usage example:
+# Examples
+Continuous time
 ```julia
 using LinearAlgebra # For identity matrix I
 using Plots
 A = [0 1; 0 0]
-B = [0;1]
+B = [0; 1]
 C = [1 0]
 sys = ss(A,B,C,0)
 Q = I
 R = I
-L = lqr(sys,Q,R)
+L = lqr(sys,Q,R) # lqr(Continuous,A,B,Q,R) can also be used
 
 u(x,t) = -L*x # Form control law,
 t=0:0.1:5
 x0 = [1,0]
 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, args...; kwargs...)
-    S, _, K = arec(A, B, R, Q, args...; kwargs...)
-    return K
-end
-
-"""
-    kalman(A, C, R1, R2)
-    kalman(sys, R1, R2)
-
-Calculate the optimal Kalman gain
-
-The `args...; kwargs...` are sent to the Riccati solver, allowing specification of cross-covariance etc. See `?MatrixEquations.arec/ared` for more help.
-
-See also `LQG`
-"""
-kalman(A, C, R1,R2, args...; kwargs...) = Matrix(lqr(A',C',R1,R2, args...; kwargs...)')
-
-function lqr(sys::AbstractStateSpace, Q, R, args...; kwargs...)
-    if iscontinuous(sys)
-        return lqr(sys.A, sys.B, Q, R, args...; kwargs...)
-    else
-        return dlqr(sys.A, sys.B, Q, R)
-    end
-end
-
-function kalman(sys::AbstractStateSpace, R1, R2, args...; kwargs...)
-    if iscontinuous(sys)
-        return Matrix(lqr(sys.A', sys.C', R1,R2, args...; kwargs...)')
-    else
-        return Matrix(dlqr(sys.A', sys.C', R1,R2, args...; kwargs...)')
-    end
-end
 
-
-"""
-    dlqr(A, B, Q, R, args...; kwargs...)
-    dlqr(sys, Q, R, args...; kwargs...)
-
-Calculate the optimal gain matrix `K` for the state-feedback law `u[k] = -K*x[k]` that
-minimizes the cost function:
-
-J = sum(x'Qx + u'Ru, 0, inf).
-
-For the discrte time model `x[k+1] = Ax[k] + Bu[k]`.
-
-See also `lqg`
-
-The `args...; kwargs...` are sent to the Riccati solver, allowing specification of cross-covariance etc. See `?MatrixEquations.ared` for more help.
-
-Usage example:
+Discrete time
 ```julia
 using LinearAlgebra # For identity matrix I
 using Plots
 Ts = 0.1
 A = [1 Ts; 0 1]
 B = [0;1]
 C = [1 0]
-sys = ss(A,B,C,0, Ts)
+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 = lqr(Discrete, A,B,Q,R) # lqr(sys,Q,R) can also be used
 
 u(x,t) = -L*x # Form control law,
 t=0:Ts:5
@@ -104,25 +54,40 @@ 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, args...; kwargs...)
-    S, _, K = ared(A, B, R, Q, args...; kwargs...)
+function lqr(::ContinuousType, A, B, Q, R, args...; kwargs...)
+    S, _, K = arec(A, B, R, Q, args...; kwargs...)
     return K
 end
 
-function dlqr(sys::AbstractStateSpace, Q, R, args...; kwargs...)
-    !isdiscrete(sys) && throw(ArgumentError("Input argument sys must be discrete-time system"))
-    return dlqr(sys.A, sys.B, Q, R, args...; kwargs...)
+function lqr(::DiscreteType, A, B, Q, R, args...; kwargs...)
+    S, _, K = ared(A, B, R, Q, args...; kwargs...)
+    return K
 end
 
+@deprecate lqr(A::AbstractMatrix, args...; kwargs...)  lqr(Continuous, A, args...; kwargs...)
+@deprecate dlqr(args...; kwargs...)  lqr(Discrete, args...; kwargs...)
+
+
+
 """
-    dkalman(A, C, R1, R2)
-    dkalman(sys, R1, R2)
+    kalman(Continuous, A, C, R1, R2)
+    kalman(Discrete, A, C, R1, R2)
+    kalman(sys, R1, R2)
 
-Calculate the optimal Kalman gain for discrete time systems
+Calculate the optimal Kalman gain
 
-The `args...; kwargs...` are sent to the Riccati solver, allowing specification of cross-covariance etc. See `?MatrixEquations.ared` for more help.
+The `args...; kwargs...` are sent to the Riccati solver, allowing specification of cross-covariance etc. See `?MatrixEquations.arec/ared` for more help.
 """
-dkalman(A, C, R1,R2, args...; kwargs...) = Matrix(dlqr(A',C',R1,R2, args...; kwargs...)')
+kalman(te, A, C, R1,R2, args...; kwargs...) = Matrix(lqr(te, A',C',R1,R2, args...; kwargs...)')
+
+function lqr(sys::AbstractStateSpace, Q, R, args...; kwargs...)
+    return lqr(sys.timeevol, sys.A, sys.B, Q, R, args...; kwargs...)
+end
+
+function kalman(sys::AbstractStateSpace, R1, R2, args...; kwargs...)
+    return Matrix(lqr(sys.timeevol, sys.A', sys.C', R1,R2, args...; kwargs...)')
+end
+
 
 """
     place(A, B, p, opt=:c)

src/types/Lti.jl

@@ -80,21 +80,12 @@ Base.propertynames(sys::LTISystem, private::Bool=false) =
 common_timeevol(systems::LTISystem...) = common_timeevol(timeevol(sys) for sys in systems)
 
 
-"""`isstable(sys)`
+"""
+    isstable(sys)
 
 Returns `true` if `sys` is stable, else returns `false`."""
-function isstable(sys::LTISystem)
-    if iscontinuous(sys)
-        if any(real.(poles(sys)).>=0)
-            return false
-        end
-    else
-        if any(abs.(poles(sys)).>=1)
-            return false
-        end
-    end
-    return true
-end
+isstable(sys::LTISystem{Continuous}) = all(real.(poles(sys)) .< 0)
+isstable(sys::LTISystem{<:Discrete}) = all(abs.(poles(sys)) .< 1)
 
 # Fallback since LTISystem not AbstractArray
 Base.size(sys::LTISystem, i::Integer) = size(sys)[i]

src/types/TimeEvolution.jl

@@ -61,3 +61,7 @@ common_timeevol(x::Continuous, ys::Continuous...) = Continuous()
 isapprox(x::TimeEvolution, y::TimeEvolution, args...; kwargs...) = false
 isapprox(x::Discrete, y::Discrete, args...; kwargs...) = isapprox(x.Ts, y.Ts, args...; kwargs...)
 isapprox(::Continuous, ::Continuous, args...; kwargs...) = true
+
+const TimeEvolType = Union{<:TimeEvolution, Type{<:TimeEvolution}}
+const DiscreteType = Union{Discrete, Type{Discrete}}
+const ContinuousType = Union{Continuous, Type{Continuous}}