minreal does not give a minimal realisation

[hurak]

When forming a matrix of transfer functions as in the following code

using ControlSystems
G = tf(1, [2, 3])
W1 = tf(4, [5, 6])
W2 = tf(7, [8, 9])
W3 = tf(10, [11, 12])
P = [W1 -W1*G; 0 W2; 0 W3*G; 1 -G]

the corresponding state-space realization is of order 7. This is expected, since the matrix or transfer functions contains 7 first-order transfer functions. Indeed,

julia> P = ss(P);

julia> size(P.A,1)
7

Some of the transfer functions defining P are repeated, hence we want to get a minimal realization

julia> Pmin = minreal(P);

julia> size(Pmin.A,1)
5

This I do not expect, because the number of distinct first-order systems is 4, which is also the order of a minimal realization that I would expect.

By the way, this is also the order of a minimal realization returned by Matlab:

>> G = tf(1, [2, 3]);
W1 = tf(4, [5, 6]);
W2 = tf(7, [8, 9]);
W3 = tf(10, [11, 12]);
P = [W1 -W1*G; 0 W2; 0 W3*G; 1 -G];
P = ss(P);
size(P.A,1)
Pmin = minreal(P);
size(Pmin.A,1)
ans =
     7
3 states removed.
ans =
     4

I will perhaps dig deeper, but at this point there is at least a difference between ControlSystems.jl and Matlab' Control Systems Toolbox.

[baggepinnen]

The mineral function has a number of numerical tolerance arguments you can experiment with, the full documentation of the low level algorithm is here
https://andreasvarga.github.io/MatrixPencils.jl/dev/lstools.html#MatrixPencils.lsminreal

For example:

julia> minreal(Ps, 1e-12).nx
4

[baggepinnen]

This I do not expect, because the number of distinct first-order systems is 4

This reasoning does not always work, the following system has two identical transfer functions but the order 2 realization when converting to ss is minimal

G = tf(1, [1, 1]);
P = [G 0; 0 G]

[baggepinnen]

It's a bit strange that the default tolerance of minreal didn't reduce fully, the Hankel singular values are quite clearly indicating order 4,

julia> hsvd(Ps)
7-element Vector{Float64}:
 0.48310407996435456
 0.31581425998117507
 0.06428150268585676
 0.0004283956439530704
 1.8425536932581322e-17
 1.6516497592076183e-18
 8.002211358937973e-26

balreal also picks up on this

julia> balreal(Ps)[1].nx
4

@andreasvarga Here's a statepace system where lsminreal with default tolerance gives a surprising result

    A = [-1.2 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 -2.7 -1.8 0.0 0.0 0.0 0.0; 0.0 1.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 -1.125 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 -2.590909090909091 -1.6363636363636365 0.0; 0.0 0.0 0.0 0.0 1.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 -1.5]
    B = [1.0 0.0; 0.0 1.0; 0.0 0.0; 0.0 1.0; 0.0 1.0; 0.0 0.0; 0.0 1.0]
    C = [0.8 0.0 -0.4 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.875 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.45454545454545453 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 -0.5]
    lsminreal(A, B, C)

julia> MatrixPencils.lsminreal(A,B,C)[1]
5×5 Matrix{Float64}:
 -2.56995      -0.0605025     1.58008       0.13812       1.40523e-15
 -0.0605025    -1.55492       0.100473     -0.0481364    -3.26223e-16
 -0.994        -0.0871747    -0.0253933     0.0946747     4.18236e-16
 -0.108018      0.311851      0.0908398    -1.14064       2.38888e-16
 -2.84643e-16  -1.70876e-16   2.26594e-16  -6.83558e-18  -1.125

are we missing something?

[andreasvarga]

It is nothing special with this case. The resulting system is just not minimal. This is because the default tolerances used within lsminreal are:
(atol, rtol) = (0.0, 1.7763568394002505e-15)
and probably too hard for this simple example.

If you use either atol = 1.e-13 or rtol = 1.e-13 you get a minimal realization.

As a general remark: I was cautious to set more "generous" default tolerances, to avoid wrong truncations. But you are free to set the default rtol larger when calling lsminreal in ControlSystems.

[baggepinnen]

According to MW their default tolerance is

julia> sqrt(eps())
1.4901161193847656e-8

which would explain the difference @hurak sees. I'll stick with the default tolerance we have to not cause any unnecessary disruption.

Thanks!