Changes in MTK initialization??

[B-LIE]

Question: What has changed with MTK initialization?

I'm running a model to steady-state, picking out the steady state value, and try to use that in linearization. My code worked yesterday, but this morning I updated Julia packages (it's been a couple of weeks since last update), and now my code does not work.

So - either a bug has been introduced in MTK, or the method to initialize models has changed.

Here is what is going on:

1. Run model to steady state:

tspan_ss = (0,1e3)   # "long" time to reach steady state
u_u_v(t) = u_v_const(t)   # u_u_v(t) is a register_symbolic function that provides input in the model
#
prob_ns_oL_ss = ODEProblem(hp_ns_oL, [], tspan_ss)
sol = solve(prob_ns_oL_ss,Rodas5())
x_ss = sol[end]   # Picks out the final value, interpreted as the steady state

2. Creating model for linearization:

@named hp_ns_LoL = HydroPower_ns(; case = :u_i)   # case is a structural_parameter, 
                                                                                     # used to remove input equation.
hp_ns_LoL = complete(hp_ns_LoL)

3. Specifying operational point for linearization:

op_x = Dict(unknowns(hp_ns_oL) .=> x_ss)  # specifying unknown values
op_u = Dict(hp_ns_LoL.u_v=>alpha_1_n/alpha_1_max) # "u" is the input
#
op_pt = merge(op_x, op_u)  # Operating point

4. Carrying out linearization -- this worked yesterday, but fails now:

mats, _ = 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...)

This leads to the following error message:

Warning: Initialization system is overdetermined. 3 equations for 0 unknowns. Initialization will default to using least squares. `SCCNo...

NOTE to "clarify" my model...

• In the model @variables section, I give initial values to sufficient number of "states" to ensure that the initial unknowns can be computed. And normal simulation works fine.
• I have had some problems using the remake function lately. If I do remake(model_name, u0 = x_ss), it seems like I get an error message related to overdetermined initialization problem. Is there a conflict with the default initial values I have given to the model in the @variables section?
• From the documentation of remake, there appears to have been made some changes??? There is also references to setp and setu, without examples of use (are these functions? keywords? or what?), as well as a replace (function?).
• Has the above possible changes also led to setting the operational point by op = op_pt malfunctioning?

[ChrisRackauckas]

Please give runnable code.

[B-LIE]

OK -- since I'm not able at short notice to create a "simple" case, I'll include the model that I'm actually working on, which is a hydropower system with a Francis turbine and a rotating aggregate (consisting of the rotational mass of the turbine and the generator).

I'll try move the needed code from my Jupyter notebook to a *.jl file and attach it. OK -- I was not allowed to upload a *.jl file, so I uploaded it as a *.txt file. The attached code runs in VSCode. The problem appears in the next-to-last line.

HydroPower_test.txt

[AayushSabharwal]

Initialization in linearization was very broken before. The biggest one was that if you passed x => 1.5 to linearization_function's operating point, and x happened to be an observed variable in the system (something a user has almost no control over) then x would be hardlocked to 1.5 for all calls to linearize with that linearization function. It was also implemented differently to all other initialization, which led to more issues. Now it hooks into the same infrastructure as everything else.

This also means that the operating point needs to be consistent with the constraints of the system. In your example, the initialization system has an equation

0 ~ Vd_p(t) - Vd_i(t)

which comes from the identical equation

        Vd_i ~ Vd_p

Using the additional equations

        Vd_d ~ Vd_p
        Vd_d ~ A_d*v_d

and the default Vd_i => Vd_n which is also a constraint for initialization, the equation turns into

Vd_n ~ A_d * v_d

in your system specification. Given your operating point values, the LHS is 24.3 and the RHS is 24.50493628701447 which is a significant violation of the constraints. Thus, the linearize call fails.

I plan to write a documentation page on this, but here's a little helper for debugging initialization issues with linearization. In the meantime, it's recommended to first go through https://docs.sciml.ai/ModelingToolkit/dev/tutorials/initialization/ for an overview of how this works. Note that for linearization, initialization_eqs provided to the system are not constraints.

Next, create the linearization function as

linfun, ssys = linearize(hp_ns_LoL, [hp_ns_LoL.u_v], [hp_ns_LoL.f_e]; op=op_pt, allow_input_derivatives=true)

linfun is a documented struct. The initialization problem is iprob = linfun.prob.f.initialization_data.initializeprob, and its system can be obtained using iprob.f.sys. It's worth inspecting this system to see what the constraints are, which equations are making it overdetermined, etc. Initialization can be run using SciMLBase.get_initial_values to see which constraints are violated. We plan on adding a higher-level user API for this to make debugging easier.

[AayushSabharwal]

With some digging, the linearization succeeds if u_v is the only entry in op

[B-LIE]

With some digging, the linearization succeeds if u_v is the only entry in op

Yes, but then I don't change the initial "state", do I? This will not allow me to simulate the system to steady state, and then change the initial state so that I can start a new simulation from steady state?

The only thing this does is that I do the linearization in the operating point consisting of the default variable value (given by Vd_i) and the chosen input u_v -- which is not steady state??

[B-LIE]

I plan to write a documentation page on this

This would be very useful.

[baggepinnen]

@B-LIE unrelated to this issue, but are you sure you're not affected by #2967 in your model? You have if statements with more than one equation in each branch

        if system == :w_aggr
            Dt(K_a) ~ Wd_ts - Wd_af - Wd_e
            K_a ~ J_a*w^2/2
            Wd_af ~ k_ba*w*abs(w)/2   # k_ba*w^2/2
        elseif system == :no_aggr
            w ~ w_n
        end 

[AayushSabharwal]

You can change the state, but it should either be consistent with the defaults, or the defaults that shouldn't be respected should be removed. If var => val is a default, passing var => nothing as part of the operating point will remove this default. It's also possible your initial conditions satisfy the constraints, but not to a high enough tolerance. Simulating the model further, might help with this, or increasing the tolerance used for initialization by passing initializealg = SciMLBase.OverrideInit(; abstol = ..., reltol = ...) to linearize.

[B-LIE]

are you sure you're not affected by #2967 in your model?

Hm... I thought I could have any number of equations in each branch... Can I have only one?

If so, can I use the alternative if statement:

system == :w_aggr ? Dt(K_a) ~ Wd_ts - Wd_af - Wd_e : nothing

or something? It is shorter...

[B-LIE]

You can change the state, but it should either be consistent with the defaults,

OK. So if I change op_x to the following,

op_x = merge( Dict(unknowns(hp_ns_oL) .=> x_ss), Dict(hp_ns_oL.Vd_i => nothing) )

then it should work?

---

OK - I tested this, and now it works, but I still get a warning:

Warning: Initialization system is overdetermined. 2 equations for 0 unknowns. Initialization will default to using least squares.

[B-LIE]

Next, create the linearization function as

linfun, ssys = linearize(hp_ns_LoL, [hp_ns_LoL.u_v], [hp_ns_LoL.f_e]; op=op_pt, allow_input_derivatives=true)

is a documented struct. The initialization problem is iprob = linfun.prob.f.initialization_data.initializeprob, and its system can be obtained using iprob.f.sys.

This doesn't work - linfun is a tuple of matrices. Are you referring to some future version of MTK that is in the development stage?

[AayushSabharwal]

My bad, the function should be linearization_function not linearize

[B-LIE]

OK... if I add that hp_ns_oL.Vd_i should be initialized to nothing, I get :

julia> op_pt = merge(op_x, op_u)
Dict{Num, Union{Nothing, Float64}} with 6 entries:
  w(t)      => 55.9432
  Vd_i(t)   => nothing
  Vd_dˍt(t) => 8.44046e-15
  u_v(t)    => 0.233333
  K_a(t)    => 3.12964e8
  v_d(t)    => 0.927487
julia> linfun, ssys = linearization_function(hp_ns_LoL, [hp_ns_LoL.u_v], [hp_ns_LoL.f_e]; op=op_pt, allow_input_derivatives=true)
Warning: Initialization system is overdetermined. 2 equations for 0 unknown
julia> iprob = linfun.prob.f.initialization_data.initializeprob
julia> iprob.f.sys

with response:

If I add that hp_ns_oL.w should be initialized to nothing, linearization_function responds with:

 Warning: Initialization system is overdetermined. 2 equations for 1 unknowns...

& a lengthy stacktrace. The initialization system is now found as:

If I instead add hp_ns_oL.K_a to nothing, linearization_function responds with
Warning: Initialization system is overdetermined. 1 equations for 0 unknowns.

The initialization system is now:

My guess why it is better to set K_a to nothing... if I keep w, then I just need to square w to find K_a. But if I keep K_a, I need to take the square root of w.

• does this make sense?

I have no idea what is the final missing overspecification. Suggestions are appreciated!

[AayushSabharwal]

I'll demonstrate how I went about debugging this problem. We know that this operating point works:

op_x = merge( Dict(unknowns(hp_ns_oL) .=> x_ss), Dict(hp_ns_oL.Vd_i => nothing) )

But is overdetermined, with 2 extra equations. So there are two variables that don't need initial values, which we've given initial values for. If we call

mats, ssys = linearize(hp_ns_LoL, [hp_ns_LoL.u_v], [hp_ns_LoL.f_e]; op=op_pt, allow_input_derivatives=true);

ssys is hp_ns_LoL with the inputs converted to parameters and structurally simplified. If we didn't have a functional operating point, linearization_function returns the same system as its second return value. This is the system that is initialized. A good candidate for extra initial conditions are the algebraic variables. Looking at equations(ssys) and unknowns(ssys):

julia> equations(ssys)
4-element Vector{Equation}:
 Differential(t)(K_a(t)) ~ Wd_ts(t) - Wd_e(t) - Wd_af(t)
 Differential(t)(v_d(t)) ~ Vd_dˍt(t) / A_d
 0 ~ -K_a(t) + (1//2)*J_a*(w(t)^2)
 0 ~ (m_d*Vd_dˍt(t)) / A_d - M_dˍt(t)

julia> unknowns(ssys)
4-element Vector{SymbolicUtils.BasicSymbolic{Real}}:
 K_a(t)
 v_d(t)
 w(t)
 Vd_dˍt(t)

This tells us w and Vd_d_t are algebraic variables. We have provided initial conditions for them, which is likely what makes the system overdetermined. To fix this, we do:

op_x = merge( Dict(unknowns(hp_ns_oL)[1:2] .=> x_ss[1:2]), Dict(hp_ns_oL.Vd_i => nothing, hp_ns_oL.w => nothing) )

Since w and Vd_d_t are also the algebraic variables of hp_ns_oL and we don't want initial values for them. Note hp_ns_oL.w => nothing, which overrides the default of w in hp_ns_LoL. However, to solve for these two algebraic variables the system we're linearizing needs to provide guesses for them, which hp_ns_LoL does not. Ideally linearize should accept a guesses keyword argument to provide extra guesses, but it doesn't. I'll add that. In the interim, we can do:

for x in unknowns(hp_ns_oL)[3:4]
    guesses(hp_ns_LoL)[x] = 1.0
end

These guesses are the guesses for the nonlinear solver which solves for w and Vd_d_t. I chose 1.0 at random, and it seems to work fine. In other cases, the guess may affect the outcome (similar to how the guess affects the result of a newton solve).

Now that we're all set up:

mats, _ = linearize(hp_ns_LoL, [hp_ns_LoL.u_v], [hp_ns_LoL.f_e]; op=op_pt, allow_input_derivatives=true)

this doesn't throw a warning and successfully solves.