Status of ControlSystems.jl

control

#1

Hello,
What is the status of the JuliaControl organization and the package ControlSystems.jl?

I want to model state space systems. Is that already implemented in ControlSystems.jl?
I know there is also the package ControlToolbox.jl, that implements it, but it has only one contributor and I am not sure if I should use it.
ControlSystems.jl has also different branches. What is the status of the autodiff branch?

We would like to use less Matlab/Simulink, because we have only one license, but I am not sure what is the best way to proceed with Julia libraries.
We also noticed that while there is a function to implement an LQR controller, the function to implement a Kalman estimator is missing.

Any comments welcome!

Uwe


#2

If you want to model a state space system why wouldn’t you use the excellent StateSpace.jl?


#3

does StateSpace even work on Julia 0.6? Last commit was 2 years ago.


#4

To be more specific: With ControlSystems.jl I can solve the LQR problem (linear quadratic regulator), see:

  lqr(A, B, Q, R)

  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)

  Solve the LQR problem for state-space system sys. Works for both discrete
  and continuous time systems.

  Usage example:

  A = [0 1; 0 0]
  B = [0;1]
  C = [1 0]
  sys = ss(A,B,C,0)
  Q = eye(2)
  R = eye(1)
  L = lqr(sys,Q,R)
  
  u(t,x) = -L*x # Form control law,
  t=0:0.1:5
  x0 = [1,0]
  y, t, x, uout = lsim(sys,u,t,x0)
  plot(t,x, lab=["Position", "Velocity"]', xlabel="Time [s]")

But what I want now is to solve the LQG problem (Linear Quadratic Gaussian).
In Matlab I can use the command lqg for this, but I did not find a way to do this with Julia yet.

Octave does not fully implement this function, but the open source package SciLab does:
https://help.scilab.org/docs/6.0.1/en_US/lqg.html

Any ideas how to implement this in Julia?


#5

Exploit the duality by calling

K = lqr(A',C',Q,R)'

K will be your Kalman gain. Notice the three transposes.


#6

You further have the function kalman and the helper type LQG available, see docs for LGQ below

julia> using ControlSystems

help?> LQG
search: LQG lqg lqgi

  G = LQG(A,B,C,D, Q1, Q2, R1, R2; qQ=0, qR=0, integrator=false)
  G = LQG(sys, args...; kwargs...)

  Return an LQG object that describes the closed control loop around the process sys=ss(A,B,C,D) where the controller is of LQG-type. The
  controller is specified by weight matrices R1,R2 that penalizes state deviations and control signal variance respectively, and
  covariance matrices Q1,Q2 which specify state drift and measurement covariance respectively. This constructor calls lqr and kalman and
  forms the closed-loop system.

  If integrator=true, the resulting controller will have intregral action. This is achieved by adding a model of a constant disturbance on
  the inputs to the system described by A,B,C,D.

  qQ and qR can be set to incorporate loop transfer recovery, i.e.,

  L = lqr(A, B, Q1+qQ*C'C, Q2)
  K = kalman(A, C, R1+qR*B*B', R2)

     Fields
    ≡≡≡≡≡≡≡≡

  When the LQG-object is populated by the lqg-function, the following fields have been made available

    •    L is the feedback matrix, such that A-BL is stable. Note that the length of the state vector (and the width of L) is increased
        by the number of inputs if the option integrator=true.
      
    •    K is the kalman gain such that A-KC is stable
      
    •    sysc is a dynamical system describing the controller u=L*inv(A-BL-KC+KDL)Ky
      

     Functions
    ≡≡≡≡≡≡≡≡≡≡≡

  Several other properties of the object are accessible with the indexing function getindex() and are called with the syntax G[:function].
  The available functions are (some have many alternative names, separated with / )

  -G[:cl] / G[:closedloop] is the closed-loop system, including observer, from reference to output, precompensated to have static gain 1
  (u = −Lx + lᵣr). -G[:S] / G[:Sin] Input sensitivity function -G[:T] / G[:Tin] Input complementary sensitivity function -G[:Sout] Output
  sensitivity function -G[:Tout] Output complementary sensitivity function -G[:CS] The transfer function from measurement noise to control
  signal -G[:DS] The transfer function from input load disturbance to output -G[:lt] / G[:looptransfer] / G[:loopgain]  =  PC -G[:rd] /
  G[:returndifference]  =  I + PC -G[:sr] / G[:stabilityrobustness]  =  I + inv(PC) -G[:sysc] / G[:controller] Returns the controller as a
  StateSpace-system

  It is also possible to access all fileds using the G[:symbol] syntax, the fields are P ,Q1,Q2,R1,R2,qQ,qR,sysc,L,K,integrator

     Example
    ≡≡≡≡≡≡≡≡≡

  qQ = 1
  qR = 1
  Q1 = 10eye(4)
  Q2 = 1eye(2)
  R1 = 1eye(6)
  R2 = 1eye(2)
  
  G = LQG(sys, Q1, Q2, R1, R2, qQ=qQ, qR=qR, integrator=true)
  
  Gcl = G[:cl]
  T = G[:T]
  S = G[:S]
  sigmaplot([S,T],logspace(-3,3,1000))
  stepplot(Gcl)

If you check the implementation of kalman you see that it uses the above described duality

function kalman(sys::StateSpace, R1,R2)
    if iscontinuous(sys)
        return lqr(sys.A', sys.C', R1,R2)'
    else
        return dlqr(sys.A', sys.C', R1,R2)'
    end
end

#7

Not working for me:

julia> using ControlSystems

help?> LQG
search:

Couldn't find LQG
Perhaps you meant lqr, lq, Libc, lcm, log, lu, let, BLAS, all, cld, fld or Val
  No documentation found.

  Binding LQG does not exist.

julia> 

On which branch are you?


#8

Master. kalman exists on latest tagged version.


#9

Ok, I did checkout master. This had huge consequences:

julia> Pkg.checkout("ControlSystems")
INFO: Checking out ControlSystems master...
INFO: Pulling ControlSystems latest master...
INFO: Installing ArrayViews v0.7.0
INFO: Upgrading BenchmarkTools: v0.2.5 => v0.3.1
INFO: Upgrading BinDeps: v0.8.7 => v0.8.8
INFO: Upgrading BinaryProvider: v0.2.3 => v0.3.2
INFO: Upgrading Blosc: v0.3.0 => v0.5.0
INFO: Installing CMakeWrapper v0.1.0
INFO: Upgrading Calculus: v0.3.0 => v0.4.0
INFO: Upgrading Compat: v0.54.0 => v0.68.0
INFO: Upgrading Conda: v0.7.1 => v0.8.1
INFO: Downgrading DataFrames: v0.10.1 => v0.0.0
INFO: Upgrading DataStructures: v0.7.4 => v0.8.3
INFO: Upgrading DocStringExtensions: v0.4.3 => v0.4.4
INFO: Installing FixedSizeArrays v0.2.5
INFO: Upgrading Formatting: v0.3.0 => v0.3.2
INFO: Upgrading HDF5: v0.8.8 => v0.9.2
INFO: Upgrading JSON: v0.17.1 => v0.17.2
INFO: Upgrading MacroTools: v0.4.0 => v0.4.1
INFO: Upgrading MathProgBase: v0.7.0 => v0.7.1
INFO: Installing Missings v0.2.9
INFO: Downgrading NullableArrays: v0.1.2 => v0.1.0
INFO: Installing Options v0.2.6
INFO: Downgrading Plots: v0.17.2 => v0.10.3
INFO: Upgrading Polynomials: v0.3.1 => v0.3.2
INFO: Upgrading PyCall: v1.15.0 => v1.16.1
INFO: Upgrading Roots: v0.4.3 => v0.6.0
INFO: Upgrading Showoff: v0.1.1 => v0.2.0
INFO: Upgrading SpecialFunctions: v0.3.8 => v0.4.0
INFO: Upgrading StaticArrays: v0.7.0 => v0.7.1
INFO: Downgrading StatsBase: v0.19.5 => v0.6.10
INFO: Installing VersionParsing v1.1.1
INFO: Removing Contour v0.4.0
INFO: Removing DataArrays v0.6.2
INFO: Removing GR v0.31.0
INFO: Removing GZip v0.3.0
INFO: Removing SortingAlgorithms v0.2.1

I think tagging a new version would be very useful.

So you suggest to use master and not the branch autodiff?

Thanks for you good work!

Uwe


#10

Well, master is not working for me:

julia> using ControlSystems

WARNING: deprecated syntax "typealias BBox Measures.Absolute2DBox" at /home/ufechner/.julia/v0.6/Plots/src/Plots.jl:109.
Use "const BBox = Measures.Absolute2DBox" instead.

WARNING: deprecated syntax "typealias AVec AbstractVector" at /home/ufechner/.julia/v0.6/Plots/src/types.jl:6.
Use "const AVec = AbstractVector" instead.

WARNING: deprecated syntax "typealias AMat AbstractMatrix" at /home/ufechner/.julia/v0.6/Plots/src/types.jl:7.
Use "const AMat = AbstractMatrix" instead.

WARNING: deprecated syntax "typealias KW Dict{Symbol,Any}" at /home/ufechner/.julia/v0.6/Plots/src/types.jl:8.
Use "const KW = Dict{Symbol,Any}" instead.

WARNING: deprecated syntax "abstract AbstractBackend" at /home/ufechner/.julia/v0.6/Plots/src/types.jl:12.
Use "abstract type AbstractBackend end" instead.

WARNING: deprecated syntax "abstract AbstractPlot{T<:AbstractBackend}" at /home/ufechner/.julia/v0.6/Plots/src/types.jl:13.
Use "abstract type AbstractPlot{T<:AbstractBackend} end" instead.

WARNING: deprecated syntax "abstract AbstractLayout" at /home/ufechner/.julia/v0.6/Plots/src/types.jl:14.
Use "abstract type AbstractLayout end" instead.

WARNING: deprecated syntax "typealias SubplotMap Dict{Any,Subplot}" at /home/ufechner/.julia/v0.6/Plots/src/types.jl:66.
Use "const SubplotMap = Dict{Any,Subplot}" instead.

WARNING: deprecated syntax "typealias P2 FixedSizeArrays.Vec{2,Float64}" at /home/ufechner/.julia/v0.6/Plots/src/components.jl:4.
Use "const P2 = FixedSizeArrays.Vec{2,Float64}" instead.

WARNING: deprecated syntax "typealias P3 FixedSizeArrays.Vec{3,Float64}" at /home/ufechner/.julia/v0.6/Plots/src/components.jl:5.
Use "const P3 = FixedSizeArrays.Vec{3,Float64}" instead.

WARNING: deprecated syntax "abstract AbstractSurface" at /home/ufechner/.julia/v0.6/Plots/src/components.jl:511.
Use "abstract type AbstractSurface end" instead.
ERROR: 
WARNING: deprecated syntax "abstract Functor{N}".
Use "abstract type Functor{N} end" instead.

WARNING: deprecated syntax "abstract FixedArray{T,NDim,SIZE}".
Use "abstract type FixedArray{T,NDim,SIZE} end" instead.

WARNING: deprecated syntax "abstract MutableFixedArray{T,NDim,SIZE}<:FixedArray{T,NDim,SIZE}".
Use "abstract type MutableFixedArray{T,NDim,SIZE}<:FixedArray{T,NDim,SIZE} end" instead.

WARNING: deprecated syntax "typealias MutableFixedVector{T,CARDINALITY} MutableFixedArray{T,1,Tuple{CARDINALITY}}".
Use "MutableFixedVector{T,CARDINALITY} = MutableFixedArray{T,1,Tuple{CARDINALITY}}" instead.

WARNING: deprecated syntax "typealias MutableFixedMatrix{T,M,N} MutableFixedArray{T,2,Tuple{M,N}}".
Use "MutableFixedMatrix{T,M,N} = MutableFixedArray{T,2,Tuple{M,N}}" instead.

WARNING: deprecated syntax "typealias FixedVector{CARDINALITY,T} FixedArray{T,1,Tuple{CARDINALITY}}".
Use "FixedVector{CARDINALITY,T} = FixedArray{T,1,Tuple{CARDINALITY}}" instead.

WARNING: deprecated syntax "typealias FixedMatrix{Row,Column,T} FixedArray{T,2,Tuple{Row,Column}}".
Use "FixedMatrix{Row,Column,T} = FixedArray{T,2,Tuple{Row,Column}}" instead.

WARNING: deprecated syntax "abstract FixedVectorNoTuple{CARDINALITY,T}<:FixedVector{CARDINALITY,T}".
Use "abstract type FixedVectorNoTuple{CARDINALITY,T}<:FixedVector{CARDINALITY,T} end" instead.
LoadError: LoadError: LoadError: UndefVarError: scale not defined
Stacktrace:
 [1] include_from_node1(::String) at ./loading.jl:576
 [2] include(::String) at ./sysimg.jl:14
 [3] include_from_node1(::String) at ./loading.jl:576
 [4] eval(::Module, ::Any) at ./boot.jl:235
 [5] _require(::Symbol) at ./loading.jl:490
 [6] require(::Symbol) at ./loading.jl:405
 [7] include_from_node1(::String) at ./loading.jl:576
 [8] eval(::Module, ::Any) at ./boot.jl:235
 [9] _require(::Symbol) at ./loading.jl:490
 [10] require(::Symbol) at ./loading.jl:405
while loading /home/ufechner/.julia/v0.6/Plots/src/components.jl, in expression starting on line 187
while loading /home/ufechner/.julia/v0.6/Plots/src/Plots.jl, in expression starting on line 115
while loading /home/ufechner/.julia/v0.6/ControlSystems/src/ControlSystems.jl, in expression starting on line 70

julia> 

Any idea?


#11

autodiff is a development branch, so it’s not expected to be stable. In fact, we’re currently rearranging a lot in that particular branch.

I’m not sure about the changes to the packages upon checkout of master, and a new release will likely not be tagged until we have sorted out https://github.com/JuliaControl/ControlSystems.jl/pull/118

Is latest tagged version not working for you? In either case, finding the kalman gain should be straight forward using any version.


#12

Well, there is no function lqg in the stable branch. Furthermore the stable branch downgrades a lot of my other packages, probably because it requires an old version of plots.jl. Just a bit annoying. And master doesn’t work. :disappointed_relieved:


#13

Then I guess the best bet is to use master and resolve any issues locally until a new stable version is tagged. Work on ControlSystems.jl has rather low priority for most of the contributers since we make limited use of it ourselves. Any contributions fixing issues is of course welcome.


#14

Ok, I managed to install master. Next problem: The example in the help of LQG is not working. What is missing?

using ControlSystems
  qQ = 1
  qR = 1
  Q1 = 10eye(4)
  Q2 = 1eye(2)
  R1 = 1eye(6)
  R2 = 1eye(2)
  
  G = LQG(sys, Q1, Q2, R1, R2, qQ=qQ, qR=qR, integrator=true)
  
  Gcl = G[:cl]
  T = G[:T]
  S = G[:S]
  sigmaplot([S,T],logspace(-3,3,1000))
  stepplot(Gcl)

julia> G = LQG(sys, Q1, Q2, R1, R2, qQ=qQ, qR=qR, integrator=true)
ERROR: UndefVarError: sys not defined

I see, that sys is not defined, but how can I define it?

Uwe


#15

sys seems to be the system you want to design a controller for. The help ?LQG states “Return an LQG object that describes the closed control loop around the process sys=ss(A,B,C,D)”

The dimensions of Q1,Q2,R1,R2 would naturally depend on the dimensions of the StateSpace system sys.