more balanced state-space realization when converting from tf

[olof3]

This gives very poorly conditioned state-space realization of Padé approximations of small delays. Which in turns gives poorly balanced Padé approximations of DelayLtiSystems with small delays.

I think that we have discussed this poor balancing before.

julia> ss(pade(1e-3, 3))
StateSpace{Continuous, Float64}
A = 
  0.0      1.0                       0.0
  0.0      0.0                       1.0
 -1.2e11  -5.999999999999999e7  -12000.0
B = 
 0.0
 0.0
 1.0
C = 
 2.4e11  0.0  24000.0
D = 
 -1.0

Continuous-time state-space model

[baggepinnen]

It is poorly conditioned, but it has exactly the same information as the transfer function. How does balancing influence future computations? Any operation performed on the system will lose accuracy, but perhaps it's useful if subsequent computations are improved by balancing.

[olof3]

For a lot of computations, related to, e.g., eigenvalues/matrix exponentials there is a preliminary balancing step, so in many cases the poor conditioning isn't noticeable.

I think the main thing is that one wants the states to have roughly the same nominal magnitude, for balanced truncation, etc. But this is also for convenient plotting and when one wants to inspect the matrices.

The following is along the lines of what I think should be done. The implementation of gebal could probably be specialized to the case of companion forms.

sys1 = ss(pade(1e-7, 4))

ilo, ihi, scale = LAPACK.gebal!('B', copy(sys1.A))
S = Diagonal(scale)

sys2 = ss(S\sys1.A*S, S\sys1.B, sys1.C*S, sys1.D)

Compare for example

x1 = step(sys1)[3]
x2 = step(sys2)[3]

In the case below, the poorly balanced version gives a LAPACK error

s = tf("s")
P = 1/s
gram(feedback(P, sys1), :o)
gram(feedback(P, sys2), :o)

[baggepinnen]

How does the gebal! scaling relate to what we do in balance_transform? Should we incorporate gebal into balance_transform, replace balance_transform or have them both?

[olof3]

How does the gebal! scaling relate to what we do in balance_transform? Should we incorporate gebal into balance_transform, replace balance_transform or have them both?

I wasn't aware of balance_transform. After looking into it, there is eventually a call to gebal! there as well, but it looks at the B and C matrices as well, which seems a little bit strange to me, but perhaps there is a good reason for it.

I have a feeling that there might be a way to exploit the companion form structure to get a very cheap implementation of the balancing, but gebal! should be good enough for a start.

[baggepinnen]

Fixed, balancing is now on by default when converting tf -> ss

julia> ss(pade(1e-3, 3))
StateSpace{Continuous, Float64}
A = 
     0.0               4096.0                  0.0
     0.0                  0.0               8192.0
 -3576.2786865234375  -7324.218749999999  -12000.0
B = 
   0.0
   0.0
 128.0
C = 
 55.87935447692871  0.0  187.5
D = 
 -1.0

Continuous-time state-space model