New descriptor system type

[andreasvarga]

I am struggling to define a new extension of the existing descriptor state space structure, to cover new functionality with structured or sparse system matrices. Recall that the basic definition used in DescriptorSystems is

struct DescriptorStateSpace{T, ET <: ETYPE{T}, MT<:AbstractMatrix{T}} <: AbstractDescriptorStateSpace
    A::MT
    E::ET
    B::MT
    C::MT
    D::MT
    Ts::Float64
    function DescriptorStateSpace{T}(A::MT, E::ETYPE{T},
                                     B::MT, C::MT, D::MT,  Ts::Real) where {T, MT<:AbstractMatrix{T}}
        dss_validation(A, E, B, C, D, Ts)
        new{T, typeof(E), MT}(A, E, B, C, D, Float64(Ts))
    end
end

which restricts the matrices A, B, C, and D to have the same type (they can be, for example, all sparse matrices), while E may have its own type (e.g., sparse, Diagonal, UpperTriangular, ... etc.).

I started to use the following definition of an extended data structure, which ensures full flexibility regarding the types of system matrices

struct DescriptorStateSpaceExt{T} <: AbstractDescriptorStateSpace
    A::AbstractMatrix
    E::Union{AbstractMatrix,UniformScaling}
    B::AbstractMatrix
    C::AbstractMatrix
    D::AbstractMatrix
    Ts::Float64
    function DescriptorStateSpaceExt{T}(A::AbstractMatrix{T}, E::Union{AbstractMatrix{T},UniformScaling{Bool}},
                                     B::AbstractMatrix{T}, C::AbstractMatrix{T}, D::AbstractMatrix{T},  Ts::Real) where {T}
        dss_validation(A, E, B, C, D, Ts)
        new(A, E, B, C, D, Float64(Ts))
    end
end

Using this data type I extended several functions such as gbalmr, gh2norm, ghanorm, evalfr, c2d, 'inv' (WIP), to handle large scale matrix models with the A and E matrices sparse or banded. These functions parallel and extend the functionality of the MATLAB's Control System Toolbox for sparse models. To support all operations with descriptor systems involving this data type involves however a substantial effort to cover all possible cases. This is still a WIP and very error prone.

To handle sparse models, I was practically forced to include SparseArrays among the dependencies of DescriptorSystems and to use explicitly the inquiry function issparse to detect sparse models. This implies that type testing occurs at runtime and there is no support during compilation for sparse models. A typical example is the function iszero to detect if a descriptor system has a null transfer function matrix. Here, the rank computations involve calling the function rank with different tolerance parameters for dense and sparse models.

A remedy would be to define another data type for sparse descriptor system state space objects to cover exclussively the functions involving sparse models. A variant with more flexibility would be to define this object with only A and E enforced to be sparse, while B, C and D just AbstractMatrices. This would prevent conversions which are almost always necessary in computations involving these matrices.

I would appreciate any opinion regarding the above aspects.

@baggepinnen What do you think what is the best balance between flexibility and efficiency?

[baggepinnen]

In ControlSystems we have the following (abridged) type to represent systems where each matrix may be of a unique type

struct HeteroStateSpace{TE <: TimeEvolution, AT<:AbstractMatrix,BT<:AbstractVecOrMat,CT<:AbstractMatrix,DT<:AbstractMatrix} <: AbstractStateSpace{TE}
    A::AT
    B::BT
    C::CT
    D::DT
end

Note how all fields A,B,C,D have their own type parameter. This makes the struct fully concretely typed, i.e., what type each matrix has is part of the type parameters. An example

HeteroStateSpace{Continuous, SMatrix{5, 5, Float64, 25}, SMatrix{5, 1, Float64, 5}, SMatrix{4, 5, Float64, 20}, SMatrix{4, 1, Float64, 4}}

note how all matrices are of different types, since StaticArrays.jl encode the size of the matrix in the type of the matrix.

Using such a system type, you could statically (at compile time) check whether or not a system object has sparse matrices using something like this

function is_sparse_system(sys::HeteroStateSpace{<:Any, AT,BT,CT,DT}) where {AT,BT,CT,DT}
    AT <: SparseMatrixCSC
end

This can resolve at compile time since the return value may be computed based on the type alone.

Another more explicit approach is to define the methods

is_sparse_system(sys::HeteroStateSpace) = false # default case
is_sparse_system(sys::HeteroStateSpace{<:Any, <: SparseMatrixCSC}) = true

where the first one is default, and the second one is invoked only for systems where the A matrix has type SparseMatrixCSC. This latter approach is also useful to dispatch to different implementation depending on what type the matrices have, e.g.,

gbalmr(sys::HeteroStateSpace) = ... # default dense implementation
gbalmr(sys::HeteroStateSpace{<:Any, <: SparseMatrixCSC}) = ... # Customized sparse implementation

[andreasvarga]

Thanks for the inputs. I decided to define a general descriptor state-space object DescriptorStateSpaceExt which is similar to the HeteroStateSpace object. However, for handling sparse models, I defined a dedicated SparseDescriptorStateSpace object, for which I implemented most of the usual operations (without concatennations). I hope this is a more efficient way to handle these modes, although I have no strong arguments for this. The DescriptorStateSpaceExt object is actually superfluous for my intended functionality, the only motivation for this structure was to allow large scale models with A and E banded matrices, as provided by the BandedMatrices package. Unfortunately, this package is only seldom updated and many bug reports are simply ignored. So I will ponder to remove this object, to avoid troubles with BandedMatrices.

[baggepinnen]

The DescriptorStateSpaceExt object is actually superfluous for my intended functionality, the only motivation for this structure was to allow large scale models with A and E banded matrices, as provided by the BandedMatrices package

Does this system type have the algorithms implemented in Julia? I'm asking, since this package, as well as MatrixPencils.jl, explicitly calls out to BLAS in several places, effectively making it impossible to use modern Julia features such as automatic differentiation and computations with BigFloats etc. I've noticed this since switching to using MatrixPencils.jl for zero computations in ControlSystems, we can no longer operate on systems with BigFloat coefficients.

[andreasvarga]

There are no special algorithms implemented for this system type (yet). For the SparseDescriptorStateSpace type, I implemented several special functions, to overcome calling the MatrixPencils functions. But some functions there work for general AbstractMatices, as for example the balancing routines.

Making MatrixPencils.jl generic, is a delicate issue, because of the use of tolerances at each rank decision. Removing completely the LAPACK calls is not the best option, since, in my opinion, most computations are still peformed using Float64 data, and therefore the code is best suited for this type of data. When making MatrixEquations generic, I kept the codes for BlasFloat and implemented new ones in Julia for AbstractMatrices. I think a similar approach would be necessary for MatrixPencils too. But this is a huge effort!

[andreasvarga]

I registered a new version of MatrixPencils.jl, where I implemented several new generic functions involved in the computation of zeros (and other similar functionality). With the new functions, computing the zeros with BigFloat data should work (I hope you can confirm this). If you use other functions which you would like to be extended in similar way please let me know.

[baggepinnen]

Wow thanks a lot, I'll be sure to try it out in the upcoming days!

[baggepinnen]

I've played around a bit with the new generic number support and it generally works well. The one missing thing I've still found is lsminreal which calls sklf_right! / sklf_left!

[andreasvarga]

I registered a new version of MatrixPencils.jl, where all basic structural reductions are now generic. I hope you can use the new version to compute minimal realizations.

[baggepinnen]

Thanks, it works very well, this is an amazing capability!

[andreasvarga]

I mostly finished making MatrixPencil.jl generic. The only part which is still not generic is formed by functions using the generalized Schur decomposition, for which a generic software is still not available in the GenericSchur.jl package (but possibly will be available in the future). To make the spectral decomposition generic I needed to include MatrixEquations.jl as a dependency. I think this is not a big issue since the two packages are usually used together in control related software. I am myself delighted of the new capabilities provided by the DescriptorSystems.jl. Thanks for the motivating impulses!