Merge pull request #118 from JuliaControl/autodiff

Generic types, major rewriting, bugfixes and several tests, finally!

.travis.yml

@@ -7,7 +7,7 @@ script:
     - if [[ -a .git/shallow ]]; then git fetch --unshallow; fi
     - julia -e 'Pkg.clone(pwd()); Pkg.build("ControlSystems")'
     #- julia -e 'ENV["PYTHON"] = ""; Pkg.clone("PyPlot"); Pkg.build("PyPlot")'
-    - julia -e 'Pkg.clone("https://github.com/JuliaControl/ControlExamplePlots.jl.git");'
+    #- julia -e 'Pkg.clone("https://github.com/JuliaControl/ControlExamplePlots.jl.git");'
     - julia -e 'Pkg.test("ControlSystems"; coverage=true)'
 after_success:
     - julia -e 'cd(Pkg.dir("ControlSystems")); Pkg.add("Coverage"); using Coverage; Coveralls.submit(Coveralls.process_folder())'

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

@@ -1,5 +1,6 @@
-julia 0.6.0
+julia 0.6.0 0.7.0
 Plots 0.7.4
-Polynomials
+Polynomials 0.3.0
 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)

example/symbolic_computations.jl

@@ -0,0 +1,86 @@
+using Revise
+using ControlSystems
+using SymPy
+
+# Some basic demonstations of working with symbolic LTI systems
+# Functionality is rather limited, and for complicated expressions the
+# printing is awful. Anyway, let's get started...
+
+
+# Need to modify some functions, to work with the Sym type
+# Some of these changes are on their way into the respective packages
+Base.complex(::Type{SymPy.Sym}) = Sym
+Base.eigvals(a::Matrix{Sym}) = Vector{Sym}(collect(keys(call_matrix_meth(a, :eigenvals))))
+
+function Polynomials.roots(p::Polynomials.Poly{SymPy.Sym})
+    x = Sym("x")
+    return SymPy.solve(p(x), x)
+end
+
+function SymPy.simplify(G::TransferFunction{ControlSystems.SisoZpk{Sym,Sym}})
+    z, p, k = zpkdata(G)
+
+    return zpk(map(x->simplify.(x), z), map(x->simplify.(x), p), map(x->simplify.(x), k))
+end
+
+function ControlSystems.impulse(sys::StateSpace{Sym})
+    t = symbols("t", real=true)
+    return simplify.(sys.C * exp(sys.A*t) * sys.B)
+end
+
+
+# Define a symbolic parameter
+a = symbols("a", real=true)
+
+# Define a statespace and a trasnfer function
+sys = ss([1 a; a 1], [0; 1], [1 0], 0)
+
+s = tf("s")
+G = (s/(5a)) / ((s/a)^2 + s/(5a) + 1)
+
+
+# Simple conversions
+@edit zpk(sys)
+tf(sys)
+
+# The coefficient looks a bit complicated, but simplifying gives..
+z, p, k = zpkdata(tf(sys))
+simplify.(k[1])
+
+zpkdata(G)
+ss(G)
+
+
+# Controllability/observability matrices
+Mc = ctrb(sys)
+Mo = obsv(sys)
+
+
+## Compute the L∞ norm (or actually  L∞ norm squared) for two symbolic systems
+w = symbols("w", real=true)
+
+sys_to_consider = sys # G
+
+sys_fr = simplify(evalfr(sys_to_consider, im*w)[1])
+sys_fr_mag = simplify(abs(sys_fr)^2)
+
+n, _ = fraction( diff(sys_fr_mag, w) )
+roots = SymPy.solve(SymPy.Poly(n), w)
+
+real_roots = roots[SymPy.imag(roots) .== 0]
+
+maximum([subs(sys_fr_mag, w => r) for r in real_roots])
+
+
+# Compute the impulse resonse of some systems (on statespace form)
+impulse(sys)[1]
+simplify(impulse(ss(G))[1])
+
+rosenbrock = [1/(s+1) 1/(s+a); 1/(s+1) 1/(s+1)]
+ss(rosenbrock)
+impulse(ss(rosenbrock))
+
+
+# Analytic impulse response
+sys2 = ss([-2.5 0;1 1.5],[1;3],[1 2],Sym(2.5))
+impulse(sys2)[1]

src/ControlSystems.jl

@@ -5,10 +5,10 @@ export  LTISystem,
         TransferFunction,
         ss,
         tf,
-        tfg,
         zpk,
         ss2tf,
         LQG,
+        primitivetype,
         # Linear Algebra
         balance,
         care,
@@ -57,6 +57,8 @@ export  LTISystem,
         step,
         impulse,
         lsim,
+        solve,
+        Simulator,
         # Frequency Response
         freqresp,
         evalfr,
@@ -67,27 +69,59 @@ export  LTISystem,
         numpoly,
         denpoly
 
-using Plots, LaTeXStrings, Requires
+
+# QUESTION: are these used? LaTeXStrings, Requires, IterTools
+using Polynomials, OrdinaryDiffEq, Plots, LaTeXStrings
 import Base: +, -, *, /, (./), (==), (.+), (.-), (.*), (!=), isapprox, convert, promote_op
+import Base.LinAlg: BlasFloat
+
+abstract type AbstractSystem end
+
+include("types/Lti.jl")
+
+include("types/SisoTf.jl")
+
+# Transfer functions and tranfer function elemements
+include("types/TransferFunction.jl")
+include("types/SisoTfTypes/polyprint.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/lti.jl")
-include("types/transferfunction.jl")
-include("types/statespace.jl")
-include("types/tf2ss.jl")
-include("types/lqg.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")
 
+include("plotting.jl")
+
+
+
 # The path has to be evaluated upon initial import
 const __CONTROLSYSTEMS_SOURCE_DIR__ = dirname(Base.source_path())