Problems with minreal for StateSpace systems

[olof3]

1. Calling minreal on an unobservable/uncontrollable statespace realization throws an error. These would be the main cases where minreal is meaningful.

minreal(ss([-1 0; 0 -2], [1; 0], [1 0], 0))

Some additional test cases that currently fail:
```julia
minreal(ss([-1 0; 0 -2], [1; 0], [1 1], 0))
minreal(ss([-1 0; 0 -2], [1; 1], [1 0], 0))

# Discrete time
minreal(ss([0.5 0; 0 0.2], [1; 0], [1 0], 0, 1))
minreal(ss([0.5 0; 0 0.2], [1; 0], [1 1], 0, 1))
minreal(ss([0.5 0; 0 0.2], [1; 1], [1 0], 0, 1))

2. Calling minreal on a discrete-time (minimal) statespace realization returns a continuous time statespace realization.

   minreal(ss(0.5, 1, 1, 0, 1))

[mfalt]

2. Should definitely be fixed as soon as possible, and shouldn't be too hard.

I am not sure what algorithm we should use to handle 1.

[olof3]

Here is a survey on minimal realizations

De Schutter, B. (2000). "Minimal state-space realization in linear system theory: an overview". Journal of computational and applied mathematics

There are a few different options to choose from. Perhaps the one in the reference by van Dooren is the way to go. An alternative that I've made some quick experiments with, and which seems promising, is to rely on the balanced realization from

Tombs, M. S., & Postlethwaite, I. (1987). "Truncated balanced realization of a stable non-minimal state-space system." International Journal of Control

If it works out I think that it could be good enough for now.

I'll try to get a tentative PR together soon.

[baggepinnen]

@olof3 Did you ever get along to trying something out here? I'm finding that minreal fails to be useful too often :(

[olof3]

I don't quite remember, here is as far as I got https://github.com/JuliaControl/ControlSystems.jl/compare/master...olof3:balreal_work

I think I came to the conclusion that balreal using the approach in "Truncated balanced realization of a stable non-minimal state-space system" by Tombs & Postlethwaite could also be used for minreal, but there is quite a zoo of different methods. Not sure what is the best option.

For best results of that approach one should use Hammarling's approach to find the Cholesky factors of the solutions to the Lyapunov equations. In the branch I do the Cholesky factorization explicitly. If I remember correctly the code in the branch does not work satisfactorily but I think it was somewhat better than what we currently have.

If you have a broken example I could see if it works and at least add it to the tests.

[baggepinnen]

The following is a double-mass model with a spring and damper in between two masses

A = [
      0 0 1 0
      0 0 0 1
      -2 2 -1 1
      3 -3 1 -1
]
B = [0, 0, 1.0, 0]
C = [1 -1 0 0]
G = ss(A, B, C, 0)

since I've chosen the C matrix as the difference between the two position states, it should have zero gain for zero frequency (the two positions converge to the same position). Our toolbox does not like this system at all! OrdinaryDiffEq.jl can not simulate it properly either, close to stationarity after an impulse, there will be some numerical issues

impulseplot(G, 50)

[olof3]

Hmm, interesting. I would guess that cancellation of marginally stable poles is a well-studied topic, although I don't know much about it. A quickfix for this particular instance would be ss(minreal(zpk(sys))), but I guess you are after a more general solution?

One possible approach could be the one in "Balanced realization and model reduction for unstable systems." by Zhou, Kemin, Gregory Salomon, and Eva Wu It seems to rely on separating the model in its stable and unstable parts so we would need some functionality for that as well.

[olof3]

Another ad-hoc fix would be

sys = ss(A - 1e-7I, B, C, 0)
sys_red = balreal(minreal(sys, 1e-5))

[ma-laforge]

Bump. I might be experiencing the same issue.

My system

A=[
 1.0  0.0   0.0         0.0  -0.00067556;
 1.0  1.0  -0.00790937  0.0  -0.00712761;
 1.0  1.0   0.992091    0.0  -0.0592255;
 0.0  0.0   1.0         1.0  -0.247138;
 0.0  0.0   1.0         1.0   0.194932;
]
B=[
 0.00067556  -0.00067556;
 0.0         -0.00712761;
 0.0         -0.0592255;
 0.0         -0.244344;
 0.0         -0.802274;
]
C=[0.0  0.0  0.0  0.0  1.0]
D=[0.0  1.0]

I think it is valid (the rest of my code works), but I'm not 100% sure.

minreal()

Anyhow, when I try to take minreal():

tol=1e-3
sys = ss(A, B, C, D, 1)
mr_sys = minreal(sys, tol)

I get:

The number of columns of B (2) and D (1) are not equal

(Failing in state_space_validation)

Sorry, I'm not an expert in state-space representations/transformations, so I don't have much more to add on this.

[olof3]

Thanks for finding that bug @ma-laforge . It seems unrelated to this issue, so I opened another one.

[baggepinnen]

All examples in this issue works with https://andreasvarga.github.io/MatrixPencils.jl/dev/lstools.html#MatrixPencils.lsminreal we could just use that