Restructuring ControlSystems.jl and its types

[mfalt]

There has been a long going discussion on how to restructure the toolbox to allow for better extendability and performance. The discussion has so far mainly been between @mfalt @baggepinnen @aytekinar and @neveritt. It is very valuable to me that we are able to agree on a solution that can help the package become better and to avoid any unnecessary forks of the project (I think this wold be really bad for the community).

I therefore tried to summarize the proposals on the wiki HERE. Anyone is welcome to edit and add new proposals and Pros/Cons on each of the points. I have tried to be as impartial as possible when stating the proposals.

The more views the better, I think anyone who is currently using or planning to use this toolbox should feel free to add their view! I am therefore also tagging some users that have previously shown interest in the package: @ufechner7, @jcrist, @denglerchr, @phelipe, @jleny, I hope you don't mind.

The main difference between our views is whether or not Proposal 2 is a good option. The main point of this proposal is to have both a SISO and a MIMO type visible to the user. The SISO type is currently only an internal representation of a function that is not directly available to the user.

I will try to wait until a few people have had the chance to weigh in before I state my opinions.

[olof3]

Some thoughts:

1. It would be good to have compatibility with a decent DSP package, perhaps DSP. For example in order to include designed filters into control loops. This is slightly annoying to do in MatLab. I guess this could be achieved by a
   a) Joint core package with basic Signals and Systems functionality
   b) Conversion of DSP types to ControlSystems types

2. Transfer functions should be able to have complex coefficients (and operate on complex signals). This is not uncommon in DSP applications (also supported by DSP) but it also comes up in a handful of control applications.

3. It would be helpful with support for systems with internal time-delays, it is a bit frustrating to be forced to apply pade when considering feedback for systems with time-delays. There are two articles by Gahinet describing how this can be done. Would also be nice to be able to do exp(-s).

4. An alternative to having separate functions for discrete and continuous time would be to have a type argument that is either p (differential operator) or q (shift operator). Don't know how this stacks up performance-wise, but the type trees would be half the size.

I don't feel the possible benefits of separate SISO and MIMO motivates the added complications..

[mfalt]

Thanks for the input @olof3, I updated the Wiki to reflect that complex numbers should not be a problem to allow in this structure. I also updated and corrected some things in the text and trees.

I think your point 3. is nice to have in the future, but it would have to be a new type I think, something I hope won't be too hard to implement.

The problem with 4. which is similar to what we have today, is that there would be extra information (sample time) in the discrete case that is not needed in the continuous. The nice thing about splitting the types is that different methods can be invoked automatically based on the Discrete/Continuous property when needed, and if not, then the functions are implemented for their ancestors.

[jleny]

I second the idea of trying to have types providing a common interface with signal processing functions, improving on Matlab in that aspect, and you might be able to attract more contributors that way. Not sure how feasible this is philosophically though (you might have to change the name of the package again 😄). In that case, exposing SISO could be useful because a lot of DSP seems to only care about that (from what I see is available in Matlab for example). It might also be useful in your plans to add sysid functions.

On the other hand just skimming through the proposal it seems to be a lot of work that has to be balanced against more visible short term benefits of say adding new functionalities, and it's important I think to maintain a type structure that is not too complicated and does not discourage occasional contributors, and of course a perhaps higher-level interface that is not too far from Matlab. This is just my opinion, I haven't spent enough time studying the proposal to say that it goes against these principles. In any case, documentation will be very important if the type tree becomes more complicated.

[baggepinnen]

Regarding conversion between DSP transfer functions and tfs from CS
I do not think it is restrictive to not expose SISO-functions to the user. When converting from DSP, their SISO would fit in a 1x1 CS::MIMO

Conversion from CS::MIMO would not be supported unless it was the special case 1x1, in which case it would be converted to a regular DSP::SISO, If DSP has support for MIMO, CS::MIMO would be converted to the corresponing DSP::MIMO.

[mfalt]

It would be nice if we could avoid an explicit split of continuous and discrete systems in the type tree, especially since the split would occur for all concrete types e.g in TF and SS. One way to do avoid this would be to let the type of the sample time in a system be a type parameter, i.e LTISystem{S<:Union{Real,Void}}, so that for example TransferFunction could be defined

type TransferFunction{T<:Union{Real,Void}, S<:SifoTf}
    Ts::T
    mat::Array{S,2}
end

It seems that there would be no performance loss, and that the types are the same size when T is Void as when the Ts field is omitted completely. And dispatch could be done by defining for example

ContT = Void
DiscT = Real

and then

bode{T<:DiscT, S<:SisoTf}(sys::TransferFunction{T,S}) = "some code for disc systems"
bode{T<:ContT, S<:SisoTf}(sys::TransferFunction{T,S}) = "some code for cont systems"

when we want different behaviour for continuous and discrete. And if we force Ts to be for example Float64, then we could do

typealias TransferFunctionCont{S} TransferFunction{Void,S}
typealias TransferFunctionDisc{S} TransferFunction{Float64,S}

which would enable for example

bode(sys::TransferFunctionDisc) = "some code for disc systems"
bode(sys::TransferFunctionCont) = "some code for cont systems"

All of this could be invisible to the user, except for when running for example

> typeof(sys)
   TransferFunction{Float64, SisoRational{Float64}}
   #or
   TransferFunction{Void, SisoRational{Float64}}

which is not completely clear to the user. The output of just

> sys

would however still be the human readable form.

[ilayn]

Hi there,

I don't have a dog in the race but I'm really happy to see that you guys took this initiative and made a really nice tool. I'm doing the same in Python for some time now under the name of harold and I'm pretty much interested in making algorithms as open as possible so that inter-translation between different languages are much easier.

What I have been experiencing on Python side is that inheritance from a common object is not necessarily a better structured way. For example, freqresp of SS and TF objects are almost mutually exclusive. Similarly, zeros, and few other properties are not even related any more. Similarly, I think the discrete/continuous distinction is something matlab poisoned us over the years. I am not saying that the algorithms don't change depending on the type but the model types only differ by Ts. I created two classes like this and they have class properties

G.SamplingSet

which returns either Z or R strings. There are some more tricks around this but the gist is this. Similarly filter/DSP representation is simply representing the arrays backwards. But the core functionalities stay the same (SOS representations might require some care though). Otherwise I guess the following mess is around the corner

Anyways, super congrats for the toolbox !

P.S. : I've just updated the SciPy Riccati solvers (PR is being reviewed) and I'll do the same to Lyapunov solvers so in case you are interested, you can update them for better accuracy. The tests are included and CAREX/DAREX examples are used. I guess that is something for the core library. Riccatis and Lyapunovs are in harold but I'll just yank them out to SciPy since they are needed for non-control people too. Since it is Python, I guess it would be useful to you guys maybe.

[ma-laforge]

RE: Polynomials (Wiki):

We plan to replace all uses of polynomials to be of the types Polynomials.jl
instead of the current internal type. It should be possible to create systems
 using polynomials, for example sys = tf(Poly([1,2]), Poly([1,2,1])).

I would personally prefer to see Polynomials.jl provide some sort of Ratio type (name TBD) that would define a ratio of polynomials. I think it would be better located there. It would be nice if the module would also overload either the /(::Poly,::Poly), or //(::Poly,::Poly) operator to make it more readable for one to build these ratio structures.

If that were done, one would instead create sys = tf(Poly([1,2]) // Poly([1,2,1])). In this case, TransferFunction could simply be a wrapper object to assist with function dispatch / type inheritance (so that TransferFunction <: AbstractSystem, or something similar).

[ma-laforge]

RE: Proposal 1 - Discrete/Continuous time

Any --> System -->  LTISystem --> TransferFunction{S<:SisoFunction} --> DTf{S}
                              \                                     \-> CTf{S}

Background: I have not put in tons of hours thinking about the practical nature of identifying the >>system itself<< as being discrete vs continuous. My gut feeling is that continuous vs discrete is an orthogonal entity (proccess type/simulation type/type of time/type of execution/... not sure what this TimeType should be called). Thus, the TimeType would be used at a higher level - like when you actually excite the system with an input.

EDITED: (previous argument likely not valid). I don't really have a practical example of why TimeType should not be associated with the system. It just seems a bit strange to me - even though there are likely some practical advantages - especially in Julia's type/dispatch environment.

In any case: if the control group decides to go this route, it might be better to tweak this type hierarchy a bit. I propose something like (feel free to improve type names, etc):

abstract SystemTimeType
immutable ContinuousTime <: SystemTimeType; end
immutable DiscreteTime <: SystemTimeType
    Ts::Float64
end

abstract System{T<:SystemTimeType}
abstract LTISystem{T<:SystemTimeType} <: System{T}
type TransferFunction{T<:SystemTimeType, S<:SisoFunction} <: LTISystem{T}; ...

This topology is a bit more like @mfalt 's suggestion above (using T<:Union{Real,Void}), but extends SystemTypeType up to LTISystem & System. This allows generic functions to be dispatched on ::LTISystem{DiscreteTime}, for example - whether the system is described using a ratio of polynomials, the state-space formalism, or a block-level description. I also find it is a little bit cleaner than making Ts of Void type for continuous time (not that it is a bad choice).

[ma-laforge]

RE: Remove array syntax for creating MIMO systems

Specifically: proposing to change how MIMO systems are created

[will] break[] the syntax that users are used to from example MATLAB
where [1sys 2sys; 3sys 4sys] creates a MIMO system

I am not an expert at dealing with MIMO systems with computational packages.

However, if you wish to provide compatibility with Matlab: have you considered creating a seperate module (ex: MatlabControl.jl)?

This way, you can build your system in a way that is more appropriate for the Julia language and its users, without worrying too much about scaring away Matlab users.

With a little ingenuity, there is probably a way to provide a Matlab-compatible wrapper module for those who don't know if they want to invest time to learn how to use Julia's core modules.

[ma-laforge]

@mfalt:

There has been a long going discussion on how to restructure the toolbox to allow
for better extendability and performance.
[...] The more views the better

@olof3:

It would be good to have compatibility with a decent DSP package, perhaps DSP.
 [...] I guess this could be achieved by a
a) Joint core package with basic Signals and Systems functionality

If the JuliaControl group actually wants to go this route, it might be good to consider looking at the DSP topic found in Julia's new discourse website:
https://discourse.julialang.org/c/domain/dsp

My guess is that we could propose that StefanKarpinski merge contol systems with DSP, as a single topic:

• DSP / Control?
• Signals & systems?

To be quite honest, I am not a big expert on the taxonomy of these subjects. I am aware there is alot of overlap between the two subjects, but I don' quite know what the "parent" field is called (i.e. where we can put in the common functionality).

NOTE: I personally believe the reason why we are used to seeing these two subjects divided in existing numerical packages is mostly a marketing issue: you want to be able to sell a DSP license to those working in digital signal processing, and a Control System license to the control people.... You might even get to sell both to many others that don't want to stick to one "domain".

[mfalt]

I think after a lot of discussion, much of it off-line, we have started to converge to a structure of the types that we like, similar to Proposal 1

Below I will outline one such design that I think could work well. These are the criteria I think we should meet:

• Dispatch on Continuous/ Discrete, TransferFunction/StateSpace
• Able to handle systems with complex coefficients
• Able to handle multiple internal representations of transfer functions, for example "zero pole gain" and "rational function"
• Simplicity for the user while maintaining computational efficiency
• Handle sparse systems

My proposal consist of one internal type tree for functions that represent the frequency responses "SisoFunction", and one for the systems that the user handles. I like the proposal by @ma-laforge for taking care of the discrete/continuous using something like

abstract SystemTimeType
immutable ContinuousTime <: SystemTimeType; end
immutable DiscreteTime <: SystemTimeType
    Ts::Float64
end

For the following tree RC<:Number encodes if the system is represented using Real or Complex numbers and of which accuracy. T<:SystemTimeType as described above, and for FransferFunction we use the SisoFunction to encode the frequency response:

Any --> System -->  LTISystem{RC<:Number, T<:SystemTimeType} --> TransferFunction{RC, T, S<:SisoFunction{RC}}
                                                            \ 
                                                             \-> StateSpace{RC, T}

Any --> SisoFunction{RC<:Number} --> SisoTf{RC} -->  SisoRational{RC}
                                                \->  SisoZpk{RC}

An implementation of Transferfunction and StateSpace could then look like

type TransferFunction{RC<:Number, T<:SystemTimeType, S<:SisoFunction{RC}} <: LTISystem
    matrix::Matrix{S}
    time::T
    nu::Int64
    ny::Int64
end

type StateSpace{RC<:Number, T<:SystemTimeType} <: LTISystem
    A::Matrix{RC}
    B::Matrix{RC}
    C::Matrix{RC}
    D::Matrix{RC}
    time::T
    nx::Int
    nu::Int
    ny::Int
end

It then remains to handle sparse systems. One proposal is to have AbstractTransferfunction (which I think we want anyway( and AbstractStateSpace and then have separate subtypes where "Matrix" above is replaced with SparseMatrixCSC{S,Int64} and SparseMatrixCSC{RC,Int64}. An alternative, which is possible now that we have diagonal dispatch is to do

type TransferFunction{RC<:Number, T<:SystemTimeType, S<:SisoFunction, M<:AbstractMatrix{S}} <: LTISystem
    matrix::M
    time::T
    nu::Int64
    ny::Int64
end

We would obviously need to consider how these are presented to the user, but it allows for great flexibility in working with structured systems. We could define aliases for commonly used types such as

const DenseStateSpace{RC,T} = StateSpace{RC,T,M,M,M,M} where {RC,T,M<:Matrix{RC}}
const SparseStateSpace{RC,T} = StateSpace{RC,T,M,M,M,M} where {RC,T,M<:SparseMatrixCSC{RC,Int64}}

or whatever the correct syntax is at the moment.
And for the user the syntax and output could still be simple, for example

> julia Ts = 2.0
> A = sprandn(50,50,0.1); B = ...;
> ss(A,B,C,D, Ts)
50x50 SparseStateSpace{Float64,DiscreteTime}:
Ts= 2.0
A = ...
B = ...

Anyone is more than welcome to comment on this design and help out with the implementation. There is a bit of work to be done, but I think we could greatly improve the structure and usability of the code. Another question that needs answering: Do we still care about carrying around names of inputs and outputs? It seems to create a bit of clutter in the code, and I personally never use it.
Comments @baggepinnen, @olof3 ?

[olof3]

After some more off-line discussions we seem to have further converged to something like this.

The thought behind it is that it should be straight-forward, while being flexible enough for most things.

Feedback and discussion are more than welcome.

[ma-laforge]

@olof3:

1. Without trying to build a solution for a concrete control problem at the moment, it is very difficult for me to give useful information about the ideal type hierarchy.

2. My guess is that @baggepinnen might have the same problem: he might not have an immediate need for a more complex type system.

3. That said, I am trying to follow the discussion/contribute if I can, and am hoping to see something come out of this.

Comment on "Toolbox 2"

1. I think the best solution for "dealing with continuous vs discrete time" will be more evident once you start developing an actual solution to a practical problem. The good thing I can see at the moment is that in all 3 suggestions, the "time type" is always a parameter of the system (👍). That means Julia will be able to dispatch directly on time.

2. It might be best if you implemented this toolbox as a separate module (probably making use of ControlSystems.jl as a base) so you can get the freedom to hammer out a good solution.

3. Personally, I think it would be nice to see some sort of a pointer to this "unofficial" (under development/with unstable API) toolbox from the JuilaControl group. Maybe there could be a JuliaControl/ExperimentalAPI repository that could provide links to such projects. This way, users could more easily find your alternative solution than they would searching on the main Juila package page.

[baggepinnen]

I'm mostly on board with the discussion carried out between Olof and Mattias (haven't caught up on the latest developments though ^^). I should mention that I share office with Mattias and we both share corridor with Olof, so a lot of discussions are going on offline. Since we are all at the same department, your view is all the more appreciated @ma-laforge, hopefully you do not share the same inevitable biases we do in Lund.

[baggepinnen]

To answer

Questions:
Is there any need to have discrete systems with sampling times of another type than Float64? I can't really come up with anything.. Symbolic computations of discrete system frequency response, is the only thing that I can really come up with.

Yes there might be. A rational sample time might allow for higher accuracy computations for some problems? In any case, simplicity is important as well, so it might be okay if this exotic use case require an extra type implementation.

On another point. I can visualize many benefits of having abstract types on many levels in the tree. The tree becomes more complex, which is a drawback, but the implementation of a new special-case type becomes much easier since many methods can specialize on a nearby abstract type, thus eliminating the need to reimplement more methods than necessary. Consider, e.g., the propsoed type

struct StateSpace{CoD, T, SoM} <: RationalLtiSystem{CoD, T, SoM}
    A::Matrix{T}
    B::Matrix{T}
    C::Matrix{T}
    D::Matrix{T}
    Ts::Float64
end

and say I want to implement a special-case StateSpace type where the fields are SMatrixes. My special type can be <: RationalLtiSystem, but I may then not make use of all methods that would work for arbitrary StateSpace-like types with A,B,C,D matrices that are <: AbstractMatrix

[ma-laforge]

@baggepinnen: I should mention that I share office with Mattias and we both share corridor with Olof, so a lot of discussions are going on offline.

Thanks for the clarification... I just wanted to make sure Mattias & Olof did not feel like there was no value in improving the API.

I myself have often noticed that the more basic API (presently implemented in ControlSystems.jl) could benefit from Julia's advanced type system. I find this "classical" API is not all that natural when it comes to expressing your solution... so I have a tendancy to pull my mind away from the problem at hand in order to focus on "making the software components work".

[ma-laforge]

@baggepinnen : Yes there might be. A rational sample time might allow for higher accuracy computations for some problems?

Good point. I myself have often found myself wanting to express "sampling frequency" (either rad or Hz) instead of sampling time.... but if you want to tackle this issue, I think that solution 2 would be best:

You could have 3 concete "Discrete time subtypes" (2 of them parameterized):

abstract type AbstractDiscreteTime <: SystemTimeType
struct DiscreteSamplingTime <: AbstractDiscreteTime
    Ts::Float64
end
struct DiscreteSamplingFreq{UNIT} <: AbstractDiscreteTime
    Fs::Float64
end

#Type aliases:
const DiscreteSamplingFreqRad = DiscreteSamplingFreq{:rad_s}
const DiscreteSamplingFreqHz = DiscreteSamplingFreq{:Hz}

...Of course I find these type names a bit verbose... It would be nice to find clear, unambiguous alternatives.

RE: Precision:

I am still uncertain whether or not using fractional time or frequency to store sampling time/frequency is actually necessary in order to get good numerical accuracy.

However, I have noticed that it can be better to describe a sampled system by specifying the time window and number of samples (instead of sampling time/frequency). Not having this level of expressivity can cause headaches when developping a solution.

[ma-laforge]

On a related note:

I often find it awkward to express poles in terms of radians/s or Hz (I would prefer being able to use either without an explicit 2pi multiplication). As a response, I have developped the lightweight Pole structure for my own experimental "SignalProcessing.jl" library:

abstract type AbastractFrequency{Symbol} end #Id
const Hertz = AbastractFrequency{:Hz}
const RadiansPerSecond = AbastractFrequency{:rad}

immutable Pole{T<:Number,AbastractFrequency}
    v::T
end

#Aliases (Still not sure if these are useful):
const PoleHz{T<:Number} = Pole{T,Hertz}
const PoleRad{T<:Number} = Pole{T,RadiansPerSecond}

Library Link: https://github.com/ma-laforge/SignalProcessing.jl

I can now express pole locations more naturally. For example, if I have a time vector, tCT, I can create a step response with the following line of Julia code:

u = step(tCT, Pole(3/pulsewidth,:rad), tdel=tmax/4)

(See implementation here)

[mfalt]

This issue is mostly done, closing in favor of smaller remaining issues such as #28