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

**URL:** <https://discourse.julialang.org/t/ann-descriptorsystems-jl-a-package-for-manipulation-of-generalized-state-space-descriptor-system-representations/62566>\
**Category:** Package Announcements\
**Tags:** announcement, control\
**Created:** [June 8, 2021, 9:20am UTC](https://discourse.julialang.org/t/ann-descriptorsystems-jl-a-package-for-manipulation-of-generalized-state-space-descriptor-system-representations/62566 "2021-06-08T09:20:45Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![andreasvarga](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/andreasvarga/32/11634_2.png) [@andreasvarga](https://discourse.julialang.org/u/andreasvarga)\
**Post date:** [June 8, 2021, 9:20am UTC](https://discourse.julialang.org/t/ann-descriptorsystems-jl-a-package-for-manipulation-of-generalized-state-space-descriptor-system-representations/62566/1 "2021-06-08T09:20:45Z")

</div>

I’m pleased to announce the release of v1.0 of the package [DescriptorSystems.jl](https://github.com/andreasvarga/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](https://github.com/andreasvarga/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](https://github.com/andreasvarga/DescriptorSystems.jl) package is based on the tools for matrix pencil manipulations provided in the [MatrixPencils.jl](https://github.com/andreasvarga/MatrixPencils.jl) package, solvers for Lyapunov, Sylvester and Riccati matrix equations available in the [MatrixEquations.jl](https://github.com/andreasvarga/MatrixEquations.jl) package, and polynomial and rational function manipulation tools provided by the [Polynomials.jl](https://github.com/JuliaMath/Polynomials.jl) package.

The functionality of [DescriptorSystems.jl](https://github.com/andreasvarga/DescriptorSystems.jl) and [MatrixEquations.jl](https://github.com/andreasvarga/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](https://github.com/andreasvarga/DescriptorSystemTools) collection of tools for MATLAB, based on Fortran implementations available in the _Systems and Control Library_ [SLICOT](https://github.com/SLICOT/SLICOT-Reference/). 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](https://doi.org/10.1007/978-1-4471-5102-9_100054-1) (see [here](https://arxiv.org/abs/1902.00009) for an extended and updated version).

**Acknowledgements**

I would like to thank [John Verzani](https://github.com/jverzani) for the helpful work in implementing the _rational function_ object within the [Polynomials.jl](https://github.com/JuliaMath/Polynomials.jl) package.

---

<div class="post-metadata">

**Author:** ![ufechner7](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ufechner7/32/51363_2.png) [@ufechner7](https://discourse.julialang.org/u/ufechner7)\
**Post date:** [June 8, 2021, 11:22am UTC](https://discourse.julialang.org/t/ann-descriptorsystems-jl-a-package-for-manipulation-of-generalized-state-space-descriptor-system-representations/62566/2 "2021-06-08T11:22:23Z")

</div>

How is the relationship to [GitHub - JuliaControl/ControlSystems.jl: A Control Systems Toolbox for Julia](https://github.com/JuliaControl/ControlSystems.jl) ?  
It also supports state space systems: [Creating Systems · ControlSystems.jl](http://juliacontrol.github.io/ControlSystems.jl/latest/man/creating_systems/#Creating-State-Space-Systems)

---

<div class="post-metadata">

**Author:** ![andreasvarga](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/andreasvarga/32/11634_2.png) [@andreasvarga](https://discourse.julialang.org/u/andreasvarga)\
**Post date:** [June 8, 2021, 11:52am UTC](https://discourse.julialang.org/t/ann-descriptorsystems-jl-a-package-for-manipulation-of-generalized-state-space-descriptor-system-representations/62566/3 "2021-06-08T11:52:03Z")

</div>

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).

---

<div class="post-metadata">

**Author:** ![ufechner7](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ufechner7/32/51363_2.png) [@ufechner7](https://discourse.julialang.org/u/ufechner7)\
**Post date:** [June 8, 2021, 11:57am UTC](https://discourse.julialang.org/t/ann-descriptorsystems-jl-a-package-for-manipulation-of-generalized-state-space-descriptor-system-representations/62566/4 "2021-06-08T11:57:05Z")

</div>

In which aspects is your approach more general?

---

<div class="post-metadata">

**Author:** ![andreasvarga](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/andreasvarga/32/11634_2.png) [@andreasvarga](https://discourse.julialang.org/u/andreasvarga)\
**Post date:** [June 8, 2021, 12:16pm UTC](https://discourse.julialang.org/t/ann-descriptorsystems-jl-a-package-for-manipulation-of-generalized-state-space-descriptor-system-representations/62566/5 "2021-06-08T12:16:17Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![imcinerney](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/imcinerney/32/9882_2.png) [@imcinerney](https://discourse.julialang.org/u/imcinerney)\
**Post date:** [June 8, 2021, 12:21pm UTC](https://discourse.julialang.org/t/ann-descriptorsystems-jl-a-package-for-manipulation-of-generalized-state-space-descriptor-system-representations/62566/6 "2021-06-08T12:21:41Z")

</div>

> [@ufechner7](#):
>
> How is the relationship to [GitHub - JuliaControl/ControlSystems.jl: A Control Systems Toolbox for Julia](https://github.com/JuliaControl/ControlSystems.jl) ?

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

---

<div class="post-metadata">

**Author:** ![imcinerney](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/imcinerney/32/9882_2.png) [@imcinerney](https://discourse.julialang.org/u/imcinerney)\
**Post date:** [June 8, 2021, 12:38pm UTC](https://discourse.julialang.org/t/ann-descriptorsystems-jl-a-package-for-manipulation-of-generalized-state-space-descriptor-system-representations/62566/7 "2021-06-08T12:38:09Z")

</div>

@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

---

<div class="post-metadata">

**Author:** ![mfalt](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mfalt/32/12467_2.png) [@mfalt](https://discourse.julialang.org/u/mfalt)\
**Post date:** [June 8, 2021, 1:26pm UTC](https://discourse.julialang.org/t/ann-descriptorsystems-jl-a-package-for-manipulation-of-generalized-state-space-descriptor-system-representations/62566/8 "2021-06-08T13:26:03Z")

</div>

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](https://github.com/olof3/ControlMatrixEquations.jl)  
I am not sure if the plan is to integrate it with MatrixEquations.jl

---

<div class="post-metadata">

**Author:** ![jlperla](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jlperla/32/34332_2.png) [@jlperla](https://discourse.julialang.org/u/jlperla)\
**Post date:** [June 8, 2021, 4:08pm UTC](https://discourse.julialang.org/t/ann-descriptorsystems-jl-a-package-for-manipulation-of-generalized-state-space-descriptor-system-representations/62566/9 "2021-06-08T16:08:25Z")

</div>

This all seems great. Is [GitHub - olof3/ControlMatrixEquations.jl: Solvers for Sylvester, Lyapunov, and Riccati Equations](https://github.com/olof3/ControlMatrixEquations.jl) 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?

---

<div class="post-metadata">

**Author:** ![olof3](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/olof3/32/23875_2.png) [@olof3](https://discourse.julialang.org/u/olof3)\
**Post date:** [June 8, 2021, 8:02pm UTC](https://discourse.julialang.org/t/ann-descriptorsystems-jl-a-package-for-manipulation-of-generalized-state-space-descriptor-system-representations/62566/10 "2021-06-08T20:02:04Z")

</div>

> [@imcinerney](#):
>
> 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).

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.

---

<div class="post-metadata">

**Author:** ![olof3](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/olof3/32/23875_2.png) [@olof3](https://discourse.julialang.org/u/olof3)\
**Post date:** [June 8, 2021, 8:29pm UTC](https://discourse.julialang.org/t/ann-descriptorsystems-jl-a-package-for-manipulation-of-generalized-state-space-descriptor-system-representations/62566/11 "2021-06-08T20:29:21Z")

</div>

> [@andreasvarga](#):
>
> I’m pleased to announce the release of v1.0 of the package [DescriptorSystems.jl](https://github.com/andreasvarga/DescriptorSystems.jl) for manipulation of descriptor systems.

Great work!

---

<div class="post-metadata">

**Author:** ![olof3](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/olof3/32/23875_2.png) [@olof3](https://discourse.julialang.org/u/olof3)\
**Post date:** [June 8, 2021, 8:58pm UTC](https://discourse.julialang.org/t/ann-descriptorsystems-jl-a-package-for-manipulation-of-generalized-state-space-descriptor-system-representations/62566/12 "2021-06-08T20:58:46Z")

</div>

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](https://github.com/olof3/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](https://en.wikipedia.org/wiki/Bartels%E2%80%93Stewart_algorithm) 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

---

<div class="post-metadata">

**Author:** ![imcinerney](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/imcinerney/32/9882_2.png) [@imcinerney](https://discourse.julialang.org/u/imcinerney)\
**Post date:** [June 9, 2021, 12:57am UTC](https://discourse.julialang.org/t/ann-descriptorsystems-jl-a-package-for-manipulation-of-generalized-state-space-descriptor-system-representations/62566/13 "2021-06-09T00:57:15Z")

</div>

> [@olof3](#):
>
> 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?

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):

```julia
>> 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:

```julia
  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:

```julia
>> 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).

---

<div class="post-metadata">

**Author:** ![olof3](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/olof3/32/23875_2.png) [@olof3](https://discourse.julialang.org/u/olof3)\
**Post date:** [June 9, 2021, 8:32am UTC](https://discourse.julialang.org/t/ann-descriptorsystems-jl-a-package-for-manipulation-of-generalized-state-space-descriptor-system-representations/62566/14 "2021-06-09T08:32:13Z")

</div>

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’

---

<div class="post-metadata">

**Author:** ![BLI](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bli/32/37206_2.png) [@BLI](https://discourse.julialang.org/u/BLI)\
**Post date:** [June 9, 2021, 8:49am UTC](https://discourse.julialang.org/t/ann-descriptorsystems-jl-a-package-for-manipulation-of-generalized-state-space-descriptor-system-representations/62566/15 "2021-06-09T08:49:38Z")

</div>

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?

---

<div class="post-metadata">

**Author:** ![andreasvarga](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/andreasvarga/32/11634_2.png) [@andreasvarga](https://discourse.julialang.org/u/andreasvarga)\
**Post date:** [June 9, 2021, 8:51am UTC](https://discourse.julialang.org/t/ann-descriptorsystems-jl-a-package-for-manipulation-of-generalized-state-space-descriptor-system-representations/62566/16 "2021-06-09T08:51:11Z")

</div>

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
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

```

---

<div class="post-metadata">

**Author:** ![andreasvarga](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/andreasvarga/32/11634_2.png) [@andreasvarga](https://discourse.julialang.org/u/andreasvarga)\
**Post date:** [June 9, 2021, 9:00am UTC](https://discourse.julialang.org/t/ann-descriptorsystems-jl-a-package-for-manipulation-of-generalized-state-space-descriptor-system-representations/62566/17 "2021-06-09T09:00:34Z")

</div>

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`) .

---

<div class="post-metadata">

**Author:** ![andreasvarga](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/andreasvarga/32/11634_2.png) [@andreasvarga](https://discourse.julialang.org/u/andreasvarga)\
**Post date:** [June 15, 2021, 2:28pm UTC](https://discourse.julialang.org/t/ann-descriptorsystems-jl-a-package-for-manipulation-of-generalized-state-space-descriptor-system-representations/62566/18 "2021-06-15T14:28:38Z")

</div>

I thank you for your comments on [MatrixEquations.jl](https://github.com/andreasvarga/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](https://github.com/SLICOT/SLICOT-Reference) Fortran library (now freely available, under a BSD license).

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

```julia
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:

```julia

>> 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.

```

---

<div class="post-metadata">

**Author:** ![BambOoxX](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bambooxx/32/22179_2.png) [@BambOoxX](https://discourse.julialang.org/u/BambOoxX)\
**Post date:** [June 16, 2021, 7:39am UTC](https://discourse.julialang.org/t/ann-descriptorsystems-jl-a-package-for-manipulation-of-generalized-state-space-descriptor-system-representations/62566/19 "2021-06-16T07:39:41Z")

</div>

@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 ?

---

<div class="post-metadata">

**Author:** ![andreasvarga](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/andreasvarga/32/11634_2.png) [@andreasvarga](https://discourse.julialang.org/u/andreasvarga)\
**Post date:** [June 16, 2021, 10:09am UTC](https://discourse.julialang.org/t/ann-descriptorsystems-jl-a-package-for-manipulation-of-generalized-state-space-descriptor-system-representations/62566/20 "2021-06-16T10:09:39Z")

</div>

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

[Next page](https://discourse.julialang.org/t/ann-descriptorsystems-jl-a-package-for-manipulation-of-generalized-state-space-descriptor-system-representations/62566.md?page=2)