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).

[B-LIE]

Try named_ss(sys, ...; initialization_abstol = 1e-8, initialization_reltol = 1e-8) and see if it makes a difference

I tried to use with linearize(...), but got an error message. Reason for trying with linearize is that I it was not 100% clear to how to use the "recipe" with named_ss...

I guess I could replace

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)

with

lsys = named_ss(...)
Gb = balance_statespace(lsys)[1]
Gb = minreal(Gb, 1e-8)

?

[B-LIE]

Btw... here is the steady state gain from turbine opening angle $\alpha_1$ to turbine shaft power:

This mean that there is a point where the "gain" changes sign, as is also shown by the dcgain plot:

[In the dcgain plot, I have not scaled the power to MW, but kept it as W]

---

This is an interesting control challenge, I guess...

• Is it even realistic to handle the change in sign of the gain, or should one operate at angles which avoids the problem?
• The gain (as well as pole/zero) changes dramatically with $\alpha_1$, which probably makes "gain scheduling" an interesting option.

---

A totally different topic... (kind of out-of-place):

• I tend to develop models/code in Jupyter notebooks. This is not an ideal solution (problems with large documents), but I like to be able to mix code with text/math/etc.
• A possible alternative is Pluto notebooks; I have not tested that option. To me, the big advantages seem to be (i) Julia code, (ii) possibility with sliders/interactivity, (iii) handling of versioning, (iv) sharing possibilities. On the flip side: I don't know (a) how big documents one can have in Pluto notebooks, (b) how "problematic" it will be with not being able to re-use variable names, etc. (for me as an "old-fashioned" programmer).
• Any thoughts on this?

[baggepinnen]

Is it even realistic to handle the change in sign of the gain

If you know the operating point $\alpha_1$ with good enough accuracy I guess you can handle this with gain scheduling. Due to the potentially very large relative uncertainty in the gain close to the zero-gain operating point, the controller should probably either not do anything here, or do something to move away from the zero-gain operating point.

Pluto notebooks can handle long (as in a lot of text and code) documents just fine, but they become awkward when computations in some cells take a long time to run. All cells that are affected by a change anywhere are rerun automatically. This is great for reproducibility, and there are mechanisms to mitigate this, but it's still awkward.

how "problematic" it will be with not being able to re-use variable names, etc. (for me as an "old-fashioned" programmer).

this is slightly annoying sometimes, but it also has upsides, and is quite easy to work around by either

• use another variable name
• introduce a local scope with a let block instead of a begin block