[ANN] DescriptorSystems.jl: A package for manipulation of generalized state-space (descriptor) system representations

Package Announcements
,
1

I’m pleased to announce the release of v1.0 of the package DescriptorSystems.jl for manipulation of descriptor systems.

Descriptor systems are the most general models employed to describe linear time-invariant systems using first order linear differential-algebraic equations (DAEs) of the form

    Edx(t)/dt = Ax(t) + Bu(t) ,

    y(t)      = Cx(t) + Du(t) ,

where x(t), u(t) and y(t) are the system state vector, input vector and output vector, respectively, the matrices A, E, B, C and D can have real or complex elements, with A and E square and E possibly singular, and t is the continuous time variable. Discrete descriptor systems, described by difference-algebraic equations, can be modelled in a similar way using x(t+1) instead of dx(t)/dt and t being the discrete time variable. Standard linear time-invariant systems correspond to E = I.

Continuous-time descriptor models frequently arise in the modeling of interconnected systems involving linear DAEs and are also common in modeling constrained mechanical systems (e.g., contact problems) or electrical circuits. Discrete-time descriptor representations are commonly used to model economic processes. Also, several system representations can be converted to descriptor system models, as for example, input-output system descriptions based on rational and polynomial matrices, second and higher order linear DAEs, or discrete linear time-periodic system models.

The DescriptorSystems.jl collection allows the numerically reliable operation on and manipulation of rational and polynomial matrices via their matrix pencil based linearizations as descriptor system realizations. The implementation of the DescriptorSystems.jl package is based on the tools for matrix pencil manipulations provided in the MatrixPencils.jl package, solvers for Lyapunov, Sylvester and Riccati matrix equations available in the MatrixEquations.jl package, and polynomial and rational function manipulation tools provided by the Polynomials.jl package.

The functionality of DescriptorSystems.jl and MatrixEquations.jl packages covers and extends in many aspects the functionality of the computational tools available in the Control System Toolbox of MATLAB and is similar to that of the DSTOOLS collection of tools for MATLAB, based on Fortran implementations available in the Systems and Control Library SLICOT. A short account on the role of the descriptor system representation as basis for reliable numerical computations for computer aided control system analysis and design can be found in the Encyclopedia of Systems and Control (2019) article (see here for an extended and updated version).

Acknowledgements

I would like to thank John Verzani for the helpful work in implementing the rational function object within the Polynomials.jl package.

16 Likes
2

How is the relationship to GitHub - JuliaControl/ControlSystems.jl: A Control Systems Toolbox for Julia ?
It also supports state space systems: Creating Systems · ControlSystems.jl

1 Like
3

I focussed mainly on implementing numerically reliable computational methods for control systems analysis and design, covering both standard systems (as those employed in ControlSystems.jl, but without allowing for time delays) as well as the more general descriptor systems models (similar to that defined in the MATLAB Control System Toolbox).

2 Likes
4

In which aspects is your approach more general?

5

The underlying model is more general, because the E matrix can be also singular. So, systems which are combinations of differential and algebraic equations can be handled.

2 Likes
6

The state space systems in ControlSystems.jl can’t support a non-identity E matrix, so they can’t model descriptor systems.

2 Likes
7

@andreasvarga excellent work! I had actually been wanting descriptor systems a few months ago when I was working on an analysis method and needed to do the conjugate transpose of a state space system and find its H$_infty$ norm (and that of its inverse). Currently, I have that implemented as a brute-force search of the unit circle, but with this package I can hopefully switch that to be an explicit computation.

One thing I am wondering is if some collaboration with the ControlSystems.jl package would be good to define operations between the types (and conversions) - e.g. implement the complex adjoint of a StateSpace type and have it return a descriptor system when needed. Maybe @baggepinnen can chime in on this since he is one of the main devs working on ControlSystems.jl

8

I definitely think we should try to combine the functionality in some way. Either by providing simple conversions, a common interface, or even incorporating the functionality into ControlSystems.jl. We have been looking at MatrixEquations.jl, but it didn’t feel like it suited our needs perfectly. I don’t remember why, @olof3 can probably explain better. He has been working on this smaller package that suits our needs GitHub - olof3/ControlMatrixEquations.jl: Solvers for Sylvester, Lyapunov, and Riccati Equations
I am not sure if the plan is to integrate it with MatrixEquations.jl

9

This all seems great. Is GitHub - olof3/ControlMatrixEquations.jl: Solvers for Sylvester, Lyapunov, and Riccati Equations and related packages intended to provide solutions for reasonably large sylvester and ricatti equations (in the thousands), or is it intended to be used for fairly small and dense ones?

10

Why do you need the descriptor systems for representing the conjugate transpose of a state-space system? Have you checked out the hinfnorm function in ControlSystems.jl, or do you want to compute the H-inf norm of a descriptor system?

We should definitely add the adjoint of a StateSpace in ControlSystems.jl, I’ve been thinking about it but never got around to do it.

11

Great work!

12

Yes, looking at the code in MatrixEquations.jl (thousands of lines of what seems like copy-pasted Fortran code) it did not really seem like a sustainable dependency.

I had a go in ControlMatrixEquations.jl to see if I could make something more readable, and have working algorithms for Sylvester/Lyapunov and Riccati Equations, which is what we primarily need for ControlSystems.jl. There are details that are yet to be done, but otherwise the solvers are in working order. It relies on the standard algorithms, Bartels–Stewart’s for the linear equations and Arnold, W. F., & Laub, A. J. (1984) for Riccati equations. As far as I can tell, these are the same algorithms that are implemented in MatrixEquations.jl.

There is a lot of code-duplication in MatrixEquations, which makes it difficult to maintain or optimize the code in a meaningful way. Here is an example in terms of performance (sizes 500x500)

julia> @btime MatrixEquations.lyapd($A, $Q)
449.960 ms (734641 allocations: 600.41 MiB)

julia> @btime ControlMatrixEquations.lyapd($A, $Q)
178.803 ms (26 allocations: 13.51 MiB)

In terms of numerics, the behavior seems similar for the linear solvers, as would be expected.

There are a few more algorithms in MatrixEquations.jl, e.g., systems of linear equations, and Hammarlings method for Lyapunov equations (lyapchol in Matlab). The latter one is certainly useful in many cases.

Some other ambitions for ControlMatrixEquations was to have better support generic Number types and to make it easy to switch between different algorithms for matrix equations (if the matrix equations community catches on and adds packages with solvers for sparse problems or with optimized cache usage, etc).

I really liked the symmetrical naming in MatrixEquations, with lyapc, lyapd, etc. but I thought that it would be even nicer with even more symmetrical naming. I also think that there are some other design decisions that perhaps could be discussed. @andreasvarga, please feel free to reach out in a PM if you would be interested in collaborating on this. Cheers

1 Like
13

The descriptor system tends to be the easier form to represent the conjugate transpose system in because it only copies the matrices into new places instead of performing computations on them. Take this example (forgive the Matlab notation, but it also has native descriptor system support that we can compare against):

>> A = [1, 0, 1;
0, 1, 1;
0, 0, 1]
A =
     1     0     1
     0     1     1
     0     0     1
>> B = [0; 0; 1]
B =
     0
     0
     1
>> C = eye(3)
C =
     1     0     0
     0     1     0
     0     0     1
>> D = [0]
D =
     0
>> sys = ss(A, B, C, D, 0.1)

If you do the conjugate transpose in Matlab (sys'), then you get the system:

  A = 
       x1  x2  x3  x4
   x1   1   0   0   0
   x2   0   1   0   0
   x3   0   0   1   0
   x4   0   0   0   1
 
  B = 
       u1  u2  u3
   x1  -1   0   0
   x2   0  -1   0
   x3   0   0  -1
   x4   0   0   0
 
  C = 
       x1  x2  x3  x4
   y1   0   0   0   1
 
  D = 
       u1  u2  u3
   y1   0   0   0
 
  E = 
       x1  x2  x3  x4
   x1   1   0   0   0
   x2   0   1   0   0
   x3   1   1   1   0
   x4   0   0   1   0
 
Sample time: 0.1 seconds
Discrete-time state-space model.

So you can see it built a descriptor system to represent the conjugate transpose. You could also try to form it as an explicit system (by removing the E matrix), but the resulting system isn’t as nice:

>> ss( sys', 'explicit' )

  A = 
               x1          x2          x3
   x1      0.9697  -2.971e-18      0.1714
   x2   1.388e-17           1           0
   x3   -0.005357  -1.165e-16        1.03
 
  B = 
               u1          u2          u3
   x1   -6.13e-17  -7.356e-17     -0.4924
   x2      -2.828       2.828           0
   x3       2.872       2.872    -0.08704
 
  C = 
               x1          x2          x3
   y1       1.908  -6.283e-17      0.6963
 
  D = 
       u1  u2  u3
   y1   1   1  -1
 
Sample time: 0.1 seconds
Discrete-time state-space model.

In my application I specifically need the Hinf norm of a system (and that of its inverse) that has the form sys'*Q*sys + R, which is best expressed as a descriptor system as I showed above (since it has that conjugate transpose on the system).

Edit: I guess I should have been clearer - I am using the conjugate transpose of the system (so I basically need the adjoint system).

14

Thanks a lot for sharing this, I hadn’t come across adjoints for discrete-time systems, where this is applicable. In that case it certainly seems preferable with a descriptor system, just as you say, (or an equation moving backwards in time).

If you are only interested in the hinfnorm of a discrete-time StateSpace system, I still think hinfnorm in ControlSystems.jl might be what you are looking for. It implements
P. Bongers, O. Bosgra, M. Steinbuch, 1991, ‘L∞-norm calculation for generalized state space systems in continuous and discrete time’

15

Very nice! – I have been looking forward to this package.

One question: does the package support reducing the DAE to an ODE with the same I/O properties?

16

This example gives me the opportunity to make some remarks on this issue. Going from a descriptor representation to a standard representation (as used in ControlSystems.jl) involves the inversion of the A matrix. Therefore, it is wise to be avoided, unless A is a well-conditioned matrix. But, if A is singular, as for example in case of systems with time-delays, then this conversion is not even possible, and working with the descriptor representation is the only option. This is also true, when inverting systems with singular feedthrough matrix D (e.g., in case of strictly proper systems).

For illustration I am giving a very simple example and show how things are handled in DescriptorSystems.jl.

julia> using DescriptorSystems
julia> sys = dss(0.,1,1,0,Ts=1)  # this constructs the (standard) state-space realization of 1/z with A = 0 (i.e., singular)
DescriptorStateSpace{Float64}

State matrix A:
1×1 Matrix{Float64}:
 0.0

Input matrix B:
1×1 Matrix{Float64}:
 1.0

Output matrix C:
1×1 Matrix{Float64}:
 1.0

Feedthrough matrix D:
1×1 Matrix{Float64}:
 0.0

Sample time: 1.0 second(s).
Discrete-time state-space model.

julia> sysc = sys'     # this builds the conjugate (descriptor) system, with singular E
DescriptorStateSpace{Float64}

State matrix A:
2×2 Matrix{Float64}:
 1.0  0.0
 0.0  1.0

Descriptor matrix E:
2×2 Matrix{Float64}:
 0.0  1.0
 0.0  0.0

Input matrix B:
2×1 Matrix{Float64}:
  0.0
 -1.0

Output matrix C:
1×2 Matrix{Float64}:
 1.0  0.0

Feedthrough matrix D:
1×1 Matrix{Float64}:
 0.0

Sample time: 1.0 second(s).
Discrete-time state-space model.

julia> dss2rm(sysc)   # its transfer function is z (as expected) 
1×1 Matrix{RationalTransferFunction{Float64, :z, Polynomials.Polynomial{Float64, :z}, 1.0}}
From input 1 to output 1
1.0*z

Sample time: 1.0 second(s).
Discrete-time rational transfer function matrix model.

julia> gpole(sysc)   # sysc has a pole at infinity, thus it is unstable (here MATLAB fails to provide the right information)
1-element Vector{Float64}:
 Inf

julia> glinfnorm(sysc)[1]   # the L-infinity norm can be computed and is one (both 1/z and z are all-pass)
1.0
julia> iszero(sysc-inv(sys))    # incidentally, the conjugate is also the inverse
true
1 Like
17

Yes. There are two functions for this: gss2ss which is what you asked for, and dss2ss, which also determines a consistent initial condition and the state mapping matrices. The later, is mainly intended for time response computations, when the state history of the original representation is also requested (see the function timeresp) .

2 Likes
18

I thank you for your comments on MatrixEquations.jl, which motivated me to try to further enhance the performance of Lyapunov solvers.
I registered today a new version of MatrixEquations, where the updated solvers for Lyapunov-type matrix equations exhibit a substantial increase of execution speed (up to two times),
while the memory allocation behaviour has been drastically improved (up to 50 times reduction of the allocated memory). Further enhancements are certainly possible (e.g., by storing
intermediary results in working arrays), but the expected increase in speed (probably less than 10%) does not justify, in my opinion, the use of additional storage.
Old Fortran programmers, like me, would never waste unnecessarily memory for insignificant increase in speed!
Also, the usage of solvers for small order equations based on static arrays led in my experiments to a marginal increase of speed (and involves certainly extra storage allocations).
An alternative to achieve the limiting speed would be to integrate the solvers available in the SLICOT Fortran library (now freely available, under a BSD license).

Here are some performance figures for the improved solvers, run on 500x500 examples:

julia> n = 500; a = rand(n,n); e = rand(n,n); q = rand(n,n); q = q'*q;

julia> @btime ControlMatrixEquations.lyapd(a,q);
  225.434 ms (26 allocations: 13.51 MiB)

julia> @btime MatrixEquations.lyapd(a,q);
  241.155 ms (27 allocations: 11.60 MiB)

julia> @btime ControlMatrixEquations.lyapd(a,q,e);
  810.124 ms (40 allocations: 24.98 MiB)

julia> @btime MatrixEquations.lyapd(a,e,q);
  820.314 ms (37 allocations: 19.26 MiB)

I would like to add that the performance of the Matlab solver dlyap, which is based on the Fortran codes from SLICOT, is probably the limit of achievable speed
(a nice challenge for further speedup of your code).
Here are, for comparison purposes, the time measures computed with Matlab 2020b:


>> a = rand(500); e = rand(500); q = rand(500); q = q'*q;
>> tic; x = dlyap(a,q); toc
Elapsed time is 0.115986 seconds.
>> tic; x = dlyap(a,e,q); toc
Elapsed time is 0.215429 seconds.
5 Likes
19

@andreasvarga, I read a bit the contents of your package, very interesting ! Nice work.
It appears your provide tools to build univariate rational transfer functions. Do you also support building multivariate ones ? Also, are there some tools included in the package the identify the numerator and denominator coefficients of the transfer function from numerical or experimental data ?

20

No, there are no tools for multivariate rational functions and also there are no tools for model identification.