Generic types

README.md

@@ -15,9 +15,27 @@ To install, in the Julia REPL:
 ```julia
 Pkg.add("ControlSystems")
 ```
-
 Note that the latest version of this package requires Julia 0.5. Users of Julia 0.4 should use [v0.1.4](https://github.com/JuliaControl/ControlSystems.jl/tree/v0.1.4) instead.
 
+## News
+### 2018-02-01
+- LTISystem types are now more generic and can hold matrices/vectors of arbitrary type. Examples (partly pseudo-code):
+```julia
+ss(1)
+ss(1.)
+ss(1im)
+ss(ForwardDiff.Dual(1.))
+ss(GPUArray([1]))
+ss(SparseMatrix([1])
+```
+Similar for `tf,zpk` etc.
+- Continuous time systems are simulated with continuous time solvers from `OrdinaryDiffEq.jl`
+- Breaking: `lsim(sys, u::Function)` syntax has changed from `u(t,x)` to `u(x,t)` to be consistent with `OrdinaryDiffEq`
+- Breaking: `feedback(P,C)` no longer returns `feedback(P*C)`. The behavior is changed to `feedback(P1, P2) = P1/(1+P1*P2)`.
+- Type `Simulator` provides lower level interface to continuous time simulation.
+- Example [autodiff.jl](https://github.com/JuliaControl/ControlSystems.jl/tree/master/example/autodiff.jl) provides an illustration of how the new generic types can be used for automatic differentiation of a cost function through the continuous-time solver, which allows for optimization of the cost function w.r.t. PID parameters.
+
+
 ## Documentation
 
 All functions have docstrings, which can be viewed from the REPL, using for example `?tf `.

REQUIRE

@@ -2,4 +2,5 @@ julia 0.6.0
 Plots 0.7.4
 Polynomials
 LaTeXStrings
-Requires
+OrdinaryDiffEq
+IterTools

docs/src/lib/constructors.md

@@ -16,6 +16,5 @@ sminreal
 ss
 ss2tf
 tf
-tfg
 zpk
 ```

docs/src/man/creatingtfs.md

@@ -49,26 +49,6 @@ Continuous-time transfer function model
 
 The transfer functions created using this method will be of type `TransferFunction{SisoZpk}`.
 
-## tfg - Generalized Representation
-If you want to work with transfer functions that are not rational functions, it is possible to use the `tfg` representation
-```julia
-tfg(str::String), tfg(str::Expr)
-```
-This function will either convert `str` to an expression or directly accept an `Expr` and create a transfer function.
-### Example:
-```julia
-tfg("1/((s+1)*exp(-sqrt(s)))")
-
-## output
-
-TransferFunction:
-1/((s+1)*exp(-sqrt(s)))
-
-Continuous-time transfer function model
-```
-The transfer functions created using this method will be of type `TransferFunction{SisoGeneralized}`.
-This type will work with some functions like `bodeplot, stepplot` but not others ,like `poles`.
-
 ## Converting between types
 It is sometime useful to convert one representation to another, this is possible using the same functions, for example
 ```julia

example/autodiff.jl

@@ -0,0 +1,114 @@
+using ControlSystems, OrdinaryDiffEq, NLopt, BlackBoxOptim
+p0          = Float64[0.2,0.8,1] # Initial guess
+K(kp,ki,kd) = pid(kp=kp, ki=ki, kd=kd)
+K(p)        = K(p...)
+
+# Define process model
+ζ  = 0.1
+ω  = 1.; ω² = ω^2
+P  = tf(ω²,[1, 2ζ*ω, ω²])*tf(1,[1,1])
+
+Ω  = logspace(-1,2,150)  # Frequency vector to eval constraints
+h  = 0.1 # Sample time for time-domain evaluation
+Tf = 60.  # Time horizon
+t  = 0:h:Tf-h
+
+Ms = 1.4 # Maximum allowed magnitude of sensitivity function
+Mt = 1.4 # Maximum allowed magnitude of complimentary sensitivity function
+
+p  = copy(p0)
+
+function timedomain(p)
+    C     = K(p[1], p[2], p[3])
+    S     = 1/(1+P*C) # Sensitivity fun
+    PS    = ss(P*S)   # TF from load disturbance to output
+    s     = Simulator(PS, (t,x) -> [1]) # Sim. unit step load disturbance
+    ty    = eltype(p) # So that all inputs to solve have same numerical type (ForwardDiff.Dual)
+    x0    = zeros(PS.nx) .|> ty
+    tspan = (ty(0.),ty(Tf))
+    sol   = solve(s, x0, tspan, Tsit5())
+    y     = PS.C*sol(t) # y = C*x
+    C,y
+end
+
+function freqdomain(p)
+    C     = K(p[1], p[2], p[3])
+    S     = 1/(1+P*C) # Sensitivity fun
+    T     = tf(1.) - S# Comp. Sensitivity fun
+    Sw    = vec(bode(S,Ω)[1]) # Freq. domain constraints
+    Tw    = vec(bode(T,Ω)[1]) # Freq. domain constraints
+    Sw,Tw
+end
+
+# Evaluates and plots candidate controller
+function evalsol(p::Vector)
+    C,y = timedomain(p)
+    Sw,Tw = freqdomain(p)
+    plot(t,y', layout=2, show=false)
+    plot!(Ω, [Sw Tw] , lab=["Sw" "Tw"], subplot=2, xscale=:log10, yscale=:log10, show=false)
+    plot!([Ω[1],Ω[end]], [Ms,Ms], c = :black, l=:dash, subplot=2, show=false, lab="Ms")
+    plot!([Ω[1],Ω[end]], [Mt,Mt], c = :purple, l=:dash, subplot=2, lab="Mt")
+    gui()
+    false
+end
+evalsol(res::BlackBoxOptim.OptimizationResults) = evalsol(best_candidate(res))
+
+
+function constraintfun(p)
+    Sw,Tw = freqdomain(p)
+    [maximum(Sw)-Ms; maximum(Tw)-Mt]
+end
+
+
+function costfun(p)
+    C,y = timedomain(p)
+    mean(abs,y) # ~ Integrated absolute error IAE
+end
+
+
+function runopt(p, costfun, constraintfun;
+    f_tol = 1e-5,
+    x_tol = 1e-3,
+    c_tol = 1e-8,
+    f_cfg = ForwardDiff.GradientConfig(costfun, p),
+    g_cfg = ForwardDiff.JacobianConfig(constraintfun, p),
+    lb = zeros(length(p)))
+
+    c1 = constraintfun(p)
+    np = length(p)
+    nc = length(c1)
+
+    function f(p::Vector, grad::Vector)
+        if length(grad) > 0
+            grad .= ForwardDiff.gradient(costfun,p,f_cfg)
+        end
+        costfun(p)
+    end
+
+    function c(result, p::Vector, grad)
+        if length(grad) > 0
+            grad .= ForwardDiff.jacobian(constraintfun,p,g_cfg)'
+        end
+        result .= constraintfun(p)
+    end
+
+    opt = Opt(:LD_SLSQP, np)
+    lower_bounds!(opt, lb)
+    xtol_rel!(opt, x_tol)
+    ftol_rel!(opt, f_tol)
+
+    min_objective!(opt, f)
+    inequality_constraint!(opt, c, c_tol*ones(nc))
+    minf,minx,ret = NLopt.optimize(opt, p)
+end
+
+f_cfg = ForwardDiff.GradientConfig(costfun, p)
+g_cfg = ForwardDiff.JacobianConfig(constraintfun, p)
+@time minf,minx,ret = runopt(1p0, costfun, constraintfun, x_tol=1e-6, c_tol=1e-12, f_cfg=f_cfg, g_cfg=g_cfg)
+evalsol(minx)
+
+# # Optimize costfun using derivative-free method
+# res1 = compare_optimizers(costfun; SearchRange = (0.,2.), NumDimensions = length(p0), MaxTime = 20.0)
+# res2 = bboptimize(costfun, NumDimensions = length(p0), MaxTime = 20.0)
+# evalsol(res2)
+# p = best_candidate(res2)

src/ControlSystems.jl

@@ -1,14 +1,16 @@
 module ControlSystems
 
+const BlasNumber = Union{Base.LinAlg.BlasFloat, Base.LinAlg.BlasInt} # QUESTION: WHy include BlasInt, why not include BlasComplex?
+
 export  LTISystem,
         StateSpace,
         TransferFunction,
         ss,
         tf,
-        tfg,
         zpk,
         ss2tf,
         LQG,
+        primitivetype,
         # Linear Algebra
         balance,
         care,
@@ -57,6 +59,8 @@ export  LTISystem,
         step,
         impulse,
         lsim,
+        solve,
+        Simulator,
         # Frequency Response
         freqresp,
         evalfr,
@@ -67,27 +71,55 @@ export  LTISystem,
         numpoly,
         denpoly
 
-using Plots, LaTeXStrings, Requires
+
+# QUESTION: are these used? LaTeXStrings, Requires, IterTools
+using Polynomials, OrdinaryDiffEq, Plots
 import Base: +, -, *, /, (./), (==), (.+), (.-), (.*), (!=), isapprox, convert, promote_op
 
-include("types/lti.jl")
-include("types/transferfunction.jl")
-include("types/statespace.jl")
-include("types/tf2ss.jl")
-include("types/lqg.jl")
+abstract type AbstractSystem end
+
+include("types/Lti.jl")
+
+abstract type SisoTf{T<:Number} end
+
+# Transfer functions and tranfer function elemements
+include("types/TransferFunction.jl")
+include("types/SisoTfTypes/SisoZpk.jl")
+include("types/SisoTfTypes/SisoRational.jl")
+include("types/SisoTfTypes/promotion.jl")
+include("types/SisoTfTypes/conversion.jl")
+
+include("types/StateSpace.jl")
+
+# Convenience constructors
+include("types/tf.jl")
+include("types/zpk.jl")
+include("types/ss.jl")
 
+include("types/lqg.jl") # QUESTION: is it really motivated to have an LQG type?
+
+include("utilities.jl")
+
+include("types/promotion.jl")
+include("types/conversion.jl")
 include("connections.jl")
-include("discrete.jl")
+
+# Analysis
+include("freqresp.jl")
+include("timeresp.jl")
+
 include("matrix_comps.jl")
 include("simplification.jl")
-include("synthesis.jl")
+
+include("discrete.jl")
 include("analysis.jl")
-include("timeresp.jl")
-include("freqresp.jl")
-include("utilities.jl")
-include("plotting.jl")
+include("synthesis.jl")
+
+include("simulators.jl")
 include("pid_design.jl")
 
+
+
 # The path has to be evaluated upon initial import
 const __CONTROLSYSTEMS_SOURCE_DIR__ = dirname(Base.source_path())
 

src/analysis.jl

@@ -7,24 +7,11 @@ pole(sys::SisoTf) = error("pole is not implemented for type $(typeof(sys))")
 
 @doc """`dcgain(sys)`
 
-Compute the gain of SISO system `sys`.""" ->
-function dcgain(sys::StateSpace, zs::Vector=tzero(sys))
-    !issiso(sys) && error("Gain only defined for siso systems")
-    return ( sys.C*(-sys.A\sys.B) + sys.D)[1]
-end
-
-@doc """`dcgain(sys)`
-
 Compute the dcgain of system `sys`.
 
 equal to G(0) for continuous-time systems and G(1) for discrete-time systems.""" ->
-function dcgain(sys::TransferFunction)
-    !issiso(sys) && error("Gain only defined for siso systems")
-    if sys.Ts > 0
-        return map(s->sum(s.num.a)/sum(s.den.a), sys.matrix)
-    end
-
-    return map(s -> s.num[end]/s.den[end], sys.matrix)
+function dcgain(sys::LTISystem)
+    return sys.Ts > 0 ? evalfr(sys, 1) : evalfr(sys, 0)
 end
 
 @doc """`markovparam(sys, n)`
@@ -78,12 +65,13 @@ function _zpk_kern(sys::TransferFunction, iy::Int, iu::Int)
     return _zpk_kern(s)
 end
 
-function _zpk_kern(s::SisoRational)
-  return roots(s.num), roots(s.den), s.num[1]/s.den[1]
+function _zpk_kern(f::SisoRational)
+    f_zpk = convert(SisoZpk, f)
+    return f_zpk.z, f_zpk.p, f_zpk.k
 end
 
-function _zpk_kern(s::SisoZpk)
-    return s.z, s.p, s.k
+function _zpk_kern(f::SisoZpk)
+    return f.z, f.p, f.k
 end
 
 @doc """`Wn, zeta, ps = damp(sys)`
@@ -142,8 +130,8 @@ end
 #
 # Note that this returns either Vector{Complex64} or Vector{Float64}
 tzero(sys::StateSpace) = tzero(sys.A, sys.B, sys.C, sys.D)
-function tzero(A::Matrix{Float64}, B::Matrix{Float64}, C::Matrix{Float64},
-        D::Matrix{Float64})
+function tzero(A::Matrix{<:BlasNumber}, B::Matrix{<:BlasNumber}, C::Matrix{<:BlasNumber},
+        D::Matrix{<:BlasNumber})
     # Balance the system
     A, B, C = balance_statespace(A, B, C)
 
@@ -176,7 +164,7 @@ end
 Implements REDUCE in the Emami-Naeini & Van Dooren paper. Returns transformed
 A, B, C, D matrices. These are empty if there are no zeros.
 """
-function reduce_sys(A::Matrix{Float64}, B::Matrix{Float64}, C::Matrix{Float64}, D::Matrix{Float64}, meps::Float64)
+function reduce_sys(A::Matrix{<:BlasNumber}, B::Matrix{<:BlasNumber}, C::Matrix{<:BlasNumber}, D::Matrix{<:BlasNumber}, meps::BlasNumber)
     Cbar, Dbar = C, D
     if isempty(A)
         return A, B, C, D
@@ -230,10 +218,10 @@ end
 
 # Determine the number of non-zero rows, with meps as a tolerance. For an
 # upper-triangular matrix, this is a good proxy for determining the row-rank.
-function fastrank(A::Matrix{Float64}, meps::Float64)
+function fastrank(A::AbstractMatrix, meps::Real)
     n, m = size(A, 1, 2)
     if n*m == 0     return 0    end
-    norms = Array{Float64}(n)
+    norms = Vector{eltype(A)}(n)
     for i = 1:n
         norms[i] = norm(A[i, :])
     end
@@ -256,11 +244,11 @@ function margin{S<:Real}(sys::LTISystem, w::AbstractVector{S}; full=false, allMa
     vals = (:wgm, :gm, :wpm, :pm, :fullPhase)
     if allMargins
         for val in vals
-            eval(:($val = Array{Array{Float64,1}}($ny,$nu)))
+            eval(:($val = Array{Array{eltype(sys),1}}($ny,$nu)))
         end
     else
         for val in vals
-            eval(:($val = Array{Float64,2}($ny,$nu)))
+            eval(:($val = Array{eltype(sys),2}($ny,$nu)))
         end
     end
     for j=1:nu
@@ -337,17 +325,17 @@ end
 
 function _allPhaseCrossings(w, phase)
     #Calculate numer of times real axis is crossed on negative side
-    n =  Array{Float64,1}(length(w)) #Nbr of crossed
-    ph = Array{Float64,1}(length(w)) #Residual
+    n =  Vector{eltype(w)}(length(w)) #Nbr of crossed
+    ph = Vector{eltype(w)}(length(w)) #Residual
     for i = eachindex(w) #Found no easier way to do this
         n[i], ph[i] = fldmod(phase[i]+180,360)#+180
     end
     _findCrossings(w, n, ph)
 end
 
 function _findCrossings(w, n, res)
-    wcross = Array{Float64,1}()
-    tcross = Array{Float64,1}()
+    wcross = Vector{eltype(w)}()
+    tcross = Vector{eltype(w)}()
     for i in 1:(length(w)-1)
         if res[i] == 0
             wcross = [wcross; w[i]]