Problem with Minimal Realization in `ControlSystems.jl` package...

[B-LIE]

Today, minreal works on Windows 11/ Julia v1.9. THANKS!

The result from minreal is, however, "wrong" -- or perhaps extremely sensitive to tolerances...

Here is my system eigenvalues -- I insert copyable matrices for the system below:

From physics, there should be 4 states (independent differential equations) and 2 algebraic equations.

• If I apply the minreal function to the system with default tolerances, the system is reduced to a 3 state ODE.
• If I specify the relative tolerance rtol, the system minimal realization has 3 states all the way to rtol=1e-16.
• If I specify relative tolerance at rtol=1e-17, the minimal realization gives 1 state
• If I specify relative tolerance rtol at 1e-18 and smaller, the minimal realization gives 6 states

A = [-6.537773175952662 0.0 0.0 0.0 -9.892378564622923e-9 0.0; 0.0 -6.537773175952662 0.0 0.0 0.0 -9.892378564622923e-9; 2.0163803998106024e8 2.0163803998106024e8 -0.006223894167415392 -1.551620418759878e8 0.002358202548321148 0.002358202548321148; 0.0 0.0 5.063545034365582e-9 -0.4479539754649166 0.0 0.0; -2.824060629317756e8 2.0198389074625736e8 -0.006234569427701143 -1.5542817673286995e8 -0.7305736722226711 0.0023622473513548576; 2.0198389074625736e8 -2.824060629317756e8 -0.006234569427701143 -1.5542817673286995e8 0.0023622473513548576 -0.7305736722226711]

B = [0.004019511633336128; 0.004019511633336128; 0.0; 0.0; 297809.51426114445; 297809.51426114445]

C = [0.0 0.0 0.0 1.0 0.0 0.0]

D = [0.0]

linsys = ss(A,B,C,D)

[B-LIE]

Hm. Fredrik Bagge Carlson made me aware of that there probably is a pole-zero cancellation, so that the minimal realization probably should have 3 states.

I now notice that:

• the original system has 2 real poles and 2 complex poles, and and 2 stray poles that should be removed
• the zero cancels one of the real poles, so the reduced system should contain one real pole and two complex conjugate poles.
• However, using minreal leads to 3 real poles.

[baggepinnen]

Here's the original discussion https://github.com/JuliaControl/ControlSystems.jl/issues/835

[baggepinnen]

For context, ControlSystems.minreal is a thin wrapper that just calls MatrixPencils.lsminreal. The system is poorly balanced, calling balance_statespace before computing the minimal realization computes the right thing

using ControlSystemsBase, Plots
A = [-6.537773175952662 0.0 0.0 0.0 -9.892378564622923e-9 0.0; 0.0 -6.537773175952662 0.0 0.0 0.0 -9.892378564622923e-9; 2.0163803998106024e8 2.0163803998106024e8 -0.006223894167415392 -1.551620418759878e8 0.002358202548321148 0.002358202548321148; 0.0 0.0 5.063545034365582e-9 -0.4479539754649166 0.0 0.0; -2.824060629317756e8 2.0198389074625736e8 -0.006234569427701143 -1.5542817673286995e8 -0.7305736722226711 0.0023622473513548576; 2.0198389074625736e8 -2.824060629317756e8 -0.006234569427701143 -1.5542817673286995e8 0.0023622473513548576 -0.7305736722226711]

B = [0.004019511633336128; 0.004019511633336128; 0.0; 0.0; 297809.51426114445; 297809.51426114445]

C = [0.0 0.0 0.0 1.0 0.0 0.0]

D = [0.0]

linsys = ss(A,B,C,D)
bsys, T = balance_statespace(linsys)
rsys = minreal(bsys)

but calling minreal on the unbalanced model

rsys = minreal(linsys)

results in an inaccurate result.

Full reproducing code follows. Note how the Bode curve of the "original" unbalanced model is inaccurate as well

using ControlSystemsBase, Plots
A = [-6.537773175952662 0.0 0.0 0.0 -9.892378564622923e-9 0.0; 0.0 -6.537773175952662 0.0 0.0 0.0 -9.892378564622923e-9; 2.0163803998106024e8 2.0163803998106024e8 -0.006223894167415392 -1.551620418759878e8 0.002358202548321148 0.002358202548321148; 0.0 0.0 5.063545034365582e-9 -0.4479539754649166 0.0 0.0; -2.824060629317756e8 2.0198389074625736e8 -0.006234569427701143 -1.5542817673286995e8 -0.7305736722226711 0.0023622473513548576; 2.0198389074625736e8 -2.824060629317756e8 -0.006234569427701143 -1.5542817673286995e8 0.0023622473513548576 -0.7305736722226711]

B = [0.004019511633336128; 0.004019511633336128; 0.0; 0.0; 297809.51426114445; 297809.51426114445]

C = [0.0 0.0 0.0 1.0 0.0 0.0]

D = [0.0]

linsys = ss(A,B,C,D)
rsys = minreal(linsys)

bsys, T = balance_statespace(linsys)
brsys = minreal(bsys)

w = exp10.(LinRange(-2, 2, 200))

bodeplot(linsys, w, l=4, lab="Original")
bodeplot!(rsys, w, l=2, lab="Minimal of original")
bodeplot!(bsys, w, l=4, lab="Balanced")
bodeplot!(brsys, w, l=2, lab="Minimal of balanced")

[B-LIE]

With doing balance_statespace before minreal, the minimal realization seems to give exactly correct "physical" eigenvalues.

I find it puzzling that there is a difference in the bodeplot of the original system linsys and the reduced system brsys:

• there is a marked difference in steady state gain
• there is a marked difference in gain overshoot
• there is a marked difference in phase

I would have guessed that these should be identical more or less.

• Perhaps it is necessary to tweak some tolerances?

[baggepinnen]

I think that the problem is simply explained by poor conditioning of the problem matrices, the coefficients are huge

julia> norm(linsys.A, Inf), norm(linsys.B, Inf), norm(linsys.C, Inf)
(2.824060629317756e8, 297809.51426114445, 1.0)

whereas after balancing

julia> norm(bsys.A, Inf), norm(bsys.B, Inf), norm(bsys.C, Inf)
(6.537773175952662, 0.032156093066689026, 0.0625)

[andreasvarga]

Ais an extremely badly conditioned matrix for eigenvalue computations. The following example illustrates the significant difference between the eigenvalues computed with the function eigvals (which uses scaling by default) and the eigenvalues provided via the function schur (which works always without scaling).

using LinearAlgebra
A = [-6.537773175952662 0.0 0.0 0.0 -9.892378564622923e-9 0.0; 
0.0 -6.537773175952662 0.0 0.0 0.0 -9.892378564622923e-9; 
2.0163803998106024e8 2.0163803998106024e8 -0.006223894167415392 -1.551620418759878e8 0.002358202548321148 0.002358202548321148; 
0.0 0.0 5.063545034365582e-9 -0.4479539754649166 0.0 0.0; 
-2.824060629317756e8 2.0198389074625736e8 -0.006234569427701143 -1.5542817673286995e8 -0.7305736722226711 0.0023622473513548576; 
2.0198389074625736e8 -2.824060629317756e8 -0.006234569427701143 -1.5542817673286995e8 0.0023622473513548576 -0.7305736722226711];
F = schur(A);
[F.values eigvals(A,scale=false) eigvals(A) ]

julia> [F.values eigvals(A,scale=false) eigvals(A) ]
6×3 Matrix{ComplexF64}:
    -7.39155+0.0im          -7.39155+0.0im          -7.27071+0.0im
    -6.75547+0.0im          -6.75547+0.0im          -6.68319+0.0im
   -0.233155+0.726709im      -0.3739+0.0im         -0.518487-0.924514im
   -0.233155-0.726709im    -0.233155-0.726709im    -0.518487+0.924514im
     -0.3739+0.0im         -0.233155+0.726709im  1.79694e-17+0.0im
 -0.00364769+0.0im       -0.00364769+0.0im       4.63885e-16+0.0im

Many computations in MatrixPencils are based on the Schur form and therefore may potentially face accuracy problems for badly scaled data, although the results are computed with backward stable algorithms.

julia> A ≈ F.vectors * F.Schur * F.vectors'
true

So, the preliminary scaling of model data may be an escape for badly scaled models.

[andreasvarga]

A fourth order model can be also computed using an additive spectral decomposition applied to the already scaled model (here is an example using DescriptorSystems).

using DescriptorSystems
 sys = dss(bsys.A,bsys.B,bsys.C,bsys.D)
sys1, sys2 = gsdec(sys, job="stable",smarg = -1.e-7);

We have sys = sys1 + sys2, where sys1 has order four, while sys2 has order two and is practically zero:

julia> iszero(sys2)
true

[andreasvarga]

I implemented some functions addressing the scaling issues for both standard and descriptor systems. With these functions, the above problems can be appropriately addressed. Functions supporting scaling of system models, will be soon available in the DescriptorSystems package too. So, I will close this issue being solved.