named_ss & linearization -- strange result

[B-LIE]

I've got a slightly complicated system (a hydro power system from turbine opening angle $\alpha_1$ to grid frequency $f_\mathrm{e}$). I'm trying to linearize the system at various opening angles. At the nominal angle, I run a 5 % step increase in $\alpha_1$, and get an increase in $f_\mathrm{e}$:

This should indicate a positive steady state gain in the system.

First, I use linearize -- the MTK linearization function. I get a linear model (A,B,C,D) with input derivative, so I create a linear model (A_1,B_1,C_1,D_1) from the time derivative of the input to the output:

mats, sysL_2 = linearize(hp_ns_LoL, [hp_ns_LoL.u_v], [hp_ns_LoL.f_e]; op=op_pt, allow_input_derivatives=true)
sys_2 = ss(mats...)
# 
nx = size(sys_2.A,1)
A_1 = [sys_2.A sys_2.B[:,1]; zeros(1,nx) 0]
B_1 = [sys_2.B[:,2]; 1]
C_1 = [sys_2.C sys_2.D[:,1]]
D_1 = sys_2.D[:,2]
sys_1 = ss(A_1,B_1,C_1,D_1)

Finally, I use that the input is L(D(u)) = s*u, i.e., I multiply the transfer function above with "s" to get u as input:

ss(zpk(sys_1)*tf("s"))

which gives the model in zpk form as:

This system has a dcgain of +225.

---

Next, I try to use named_ss from ControlSystemsMTK -- which is where the weird result comes in:

named_ss(hp_ns_LoL, [hp_ns_LoL.u_v], [hp_ns_LoL.f_e]; op=op_pt, allow_input_derivatives=true)

The result in zpk form is:

[The model from MTK had similar false zeros/poles, but these disappeared when I multiplied with tf("s")...]

Observe here:

• This system has close to a pure integrator
• This system appears to have the wrong sign

If I "guess" that this system really has derivative of u as input and do the same trick as above, I get:

which has a dcgain of ca. -230 -- which has the wrong sign.

Questions:

1. Is there something I do incorrectly? (OK -- I didn't include every step)
2. It is weird that it seems like the result of named_ss has the derivative of the input as input (my system clearly does not have a true integrator; a steady state is reached)
3. The sign is wrong, as far as I can see.

What are the explanations? Is this a bug, or do I do something wrong?

Here is the tuple of matrices that linearize returns:

(A = [0.0 2.778983834717109e8 1.4122312296634873e6 0.0; 0.0 0.0 0.0 0.037848975765016724; 0.0 24.837541148074962 0.12622006230897712 0.0; -0.0 -4.620724819774693 -0.023481719514324866 -0.6841991610512456], B = [-5.042589978197361e8 0.0; -0.0 0.0; -45.068824982602656 -0.0; 8.384511049369085 54.98555939873381], C = [0.0 0.0 0.954929658551372 0.0], D = [0.0 0.0])

I'm on MTK v9.64.1, ControlSystems v1.11.2, ControlSystemsMTK v2.3.1.

[B-LIE]

Even stranger... with the following matrices, I get the same transfer function with the two methods ("my" elimination of input derivative -- you suggested this method to me a while ago; the named_ss method):

(A = [-0.0075449237853825925 1.6716817118020731e-6 0.0; 1864.7356343162514 -0.4131578457122937 0.0; 0.011864343330426718 -2.6287085638214332e-6 0.0], B = [0.0 0.0; 0.0 52566.418015009294; 0.0 0.3284546792274811], C = [1.4683007399899215e8 0.0 0.0], D = [-9.157636303058283e7 0.0])

[baggepinnen]

With the matrices from the first post, these two methods both produce the same transfer function

lsys = (A = [0.0 2.778983834717109e8 1.4122312296634873e6 0.0; 0.0 0.0 0.0 0.037848975765016724; 0.0 24.837541148074962 0.12622006230897712 0.0; -0.0 -4.620724819774693 -0.023481719514324866 -0.6841991610512456], B = [-5.042589978197361e8 0.0; -0.0 0.0; -45.068824982602656 -0.0; 8.384511049369085 54.98555939873381], C = [0.0 0.0 0.954929658551372 0.0], D = [0.0 0.0])

G = ControlSystemsMTK.causal_simplification(lsys, [1=>2])
G2 = (ss(lsys.A, lsys.B[:,1], lsys.C, lsys.D[:,1])) + (ss(lsys.A, lsys.B[:,2], lsys.C, lsys.D[:,2]))*tf('s')
ControlSystemsBase.bodeplot([G, G2], size=(800,1200))

causal_simplification is the function that is called inside named_ss to handle the input derivatives

[baggepinnen]

same with your last set of matrices

[B-LIE]

causal_simplification is the function that is called inside named_ss to handle the input derivatives

I found that function, but didn't understand how to specify the second input argument.

[baggepinnen]

The second argument is a map from input indices to the indices of their derivatives, in this case 1 is the input and 2 is its derivative

[baggepinnen]

Look at the code I posted that calls the function

[B-LIE]

Yea, I just saw it.

[B-LIE]

OK -- if I run dcgain() on the two systems, I get:

([224.94741728183658;;], [-5.495512471427165e18;;])

where the first one (the correct one, I think) is from my removal of the input derivative, while the second one is what I get from named_ss.

Maybe I was just lucky? The model from named_ss appears to give the correct Bode plot in the frequency you display (which is perhaps the one of interest), but there is a "false" almost-integration.

[baggepinnen]

The numerics of this system appear to be extremely sensitive, minreal with the default tolerance does not produce the correct result, I have to lower the tolerance quite a bit. dcgain evaluates the frequency response in 0, but you can try dcgain(sys, 1e-6) to evaluate in a slightly positive frequency to avoid issues with the "almost integrators"

[baggepinnen]

By the way, try your best to avoid converting the statespace system to a transfer function, transfer functions are terrible from a numerics perspective. What do you get if you call dcgain(sys, 1e-6) where sys is the statespace object returned from named_ss?

[B-LIE]

What do you get if you call dcgain(sys, 1e-6) where sys is the statespace object returned from named_ss

If I just do dcgain(sys), then I get

1×1 Matrix{Float64}:
 -5.495512471427165e18

If I do dcgain(sys, 1e-6), I get

1×1 Matrix{Float64}:
 224.5610584029963

[baggepinnen]

With this PR

• enable tzeros for BigFloat ControlSystems.jl#976

the zero computation performs a numerical balancing of the statespace system before calling the zero routine in MatrixPencils, this improves the numerics for this system. It also makes possible to compute zeros in BigFloat precision by

using GenericSchur # required
tzeros(big(1.0)sys)

multiplying by big(1.0) gives the system extended precision. The issue with your particular system is that the poles and zeros at the origin are not always cancelling, and different ways of simplifying the system thus produces different result near the origin:

If you look at the statespace matrices and see large coefficients, like the ~1e8 in the B matrix here

julia> ss(lsys...)
StateSpace{Continuous, Float64}
A = 
  0.0   2.778983834717109e8   1.4122312296634873e6   0.0
  0.0   0.0                   0.0                    0.037848975765016724
  0.0  24.837541148074962     0.12622006230897712    0.0
 -0.0  -4.620724819774693    -0.023481719514324866  -0.6841991610512456
B = 
  -5.042589978197361e8   0.0
  -0.0                   0.0
 -45.068824982602656    -0.0
   8.384511049369085    54.98555939873381
C = 
 0.0  0.0  0.954929658551372  0.0
D = 
 0.0  0.0

Continuous-time state-space model

it's often good for numerics to perform balancing

julia> balance_statespace(ss(lsys...))[1]
StateSpace{Continuous, Float64}
A = 
 0.0  1060.0982035511433          5.387234610227536        0.0
 0.0     0.0                      0.0                    620.117618934034
 0.0    24.837541148074962        0.12622006230897712      0.0
 0.0    -0.00028202666136320145  -1.4332104195754923e-6   -0.6841991610512456
B = 
 -60.112356879679695      0.0
   0.0                    0.0
  -1.408400780706333      0.0
   1.5992185686815423e-5  0.00010487663154360545
C = 
 0.0  0.0  30.557749073643905  0.0
D = 
 0.0  0.0

Continuous-time state-space model

this balancing is now done automatically in tzeros (also in bodeplot), but this does not help dcgain in this case.

---

This PR

• add option to simplify causality in a different way #73

adds an option to simplify the causality using the sys[:,uinds] + sys[:,duinds]*tf('s') method rather than the descriptor-system based method that is currently used.

---

With these two PRs merged, your particular system can be robustly handled by

lsys = (A = [0.0 2.778983834717109e8 1.4122312296634873e6 0.0; 0.0 0.0 0.0 0.037848975765016724; 0.0 24.837541148074962 0.12622006230897712 0.0; -0.0 -4.620724819774693 -0.023481719514324866 -0.6841991610512456], B = [-5.042589978197361e8 0.0; -0.0 0.0; -45.068824982602656 -0.0; 8.384511049369085 54.98555939873381], C = [0.0 0.0 0.954929658551372 0.0], D = [0.0 0.0])
G = ControlSystemsMTK.causal_simplification(lsys, [1=>2])
Gb = balance_statespace(G)[1]
Gb = minreal(Gb, 1e-8)

julia> dcgain(Gb)
1×1 Matrix{Float64}:
 224.94742230326355

[baggepinnen]

Another thing you might want to try is to tighten the tolerance used for root finding when linearizing, this might be the source of the numerical issues. Try named_ss(sys, ...; initialization_abstol = 1e-8, initialization_reltol = 1e-8) and see if it makes a difference

[B-LIE]

Thanks -- and sorry for slow response. Today is my exercise day (in the morning), and after lunch, the cafes I've been to have had rotten wifi and virtually zero cell phone connection (so no tethering).