Help with initializing a model

[bradcarman]

Take the following model (Note: Hydraulic branch of the ModelingToolkitStandardLibrary is required currently)...

using ModelingToolkitStandardLibrary.Hydraulic.IsothermalCompressible: FixedPressure,
                                                                       Orifice, Volume,
                                                                       HydraulicFluid, HydraulicPort, liquid_density
using ModelingToolkitStandardLibrary.Blocks: Ramp, Constant
using ModelingToolkitStandardLibrary.Mechanical.Translational: Force, MechanicalPort, Mass

using ModelingToolkit
using ModelingToolkit: t_nounits as t, D_nounits as D
using OrdinaryDiffEq

@mtkmodel System begin
    @components begin
        src = FixedPressure(p = 101325)
        res = Orifice()
        act = Volume(direction = -1, area = 0.1)
        # mass = Mass(m=1e-5)
        fluid = HydraulicFluid(density = 1000, bulk_modulus = 2.0e9)
    end
    @equations begin
        connect(src.port, res.port_a)
        connect(res.port_b, act.port)
        # connect(act.flange, mass.flange)

        connect(fluid, src.port)
    end
end

@mtkbuild sys = System()

Note: this model was providing an error (#3138) when including the Mass component, so it's being removed. However this should still work. This model will not have any forces in the actuator, and will therefore have zero pressure, therefore no compressibility, and the flow dm will be dictated by the orifice. This means the system is solving for the equations...

dm = function of the orifice pressure differential (orifice equation)
dm = rho*dx*area (swept volume equation)

I can solve this system and get the correct result.

prob = ODEProblem(sys, [], (0, 1), []) 
sol = solve(prob, Rodas5P())

using Plots
plot(sol; idxs = [sys.act.dx*sys.fluid.ρ*sys.act.area, sys.act.dm], ylims=(6,10))

As can be seen, both the actuator velocity dx and the mass flow dm are constants...

However, the trouble is, how can one specify the initial actuator position x? As shown with the explanation above, the problem is not compressible anymore, and therefore x is eliminated essentially from the problem, I believe this is what causes an issue here. So for example, take...

initsys = ModelingToolkit.generate_initializesystem(sys)
initsys = structural_simplify(initsys) #ERROR: InvalidSystemException: The system is structurally singular!
initprob = NonlinearProblem(complete(initsys), [t=>0, sys.act.x => 1], [])
initsol = solve(initprob; abstol=1e-6)
initsol[sys.act.x] #1.0

I can assemble an initialization problem and set sys.act.x to 1.0 and get a successful solution. However, note this problem cannot be structurally simplified! See #2952. How can one do this with the ODEProblem interface? The system will complain about an unbalanced initialization, which is not my intent, I only want to set an observable variable start value, not necessarily add an initial equation...

prob = ODEProblem(sys, [sys.act.x => 1], (0, 1), []) # Warning: Initialization system is overdetermined. 2 equations for 1 unknowns. 

So, is there another interface for ODEProblem that allows me to set a start value for an observed variable without breaking the initialization system?

[AayushSabharwal]

However, the trouble is, how can one specify the initial actuator position x? As shown with the explanation above, the problem is not compressible anymore, and therefore x is eliminated essentially from the problem, I believe this is what causes an issue here.

act.x is still an unknown in the system? I'm not sure I understand what you mean here.

julia> @mtkbuild sys = System()
Model sys with 2 equations
Unknowns (2):
  act₊x(t)
  act₊dx(t)
Parameters (74):
  src₊port₊p_gas [defaults to fluid₊p_gas]
  res₊port_a₊p_gas [defaults to fluid₊p_gas]
  res₊port_b₊β [defaults to fluid₊β]
  act₊port₊β [defaults to fluid₊β]
  res₊port_b₊μ [defaults to fluid₊μ]
  res₊valve₊port_b₊n [defaults to fluid₊n]
  fluid₊let_gas [defaults to 1]
  res₊valve₊base₊port_b₊ρ [defaults to fluid₊ρ]
  src₊port₊β [defaults to fluid₊β]
  src₊port₊ρ_gas [defaults to fluid₊ρ_gas]
  res₊valve₊port_b₊p_gas [defaults to fluid₊p_gas]

I can assemble an initialization problem and set sys.act.x to 1.0 and get a successful solution. However, note this problem cannot be structurally simplified!

It can be, if you pass fully_determined = false to structural_simplify. Also note that setting sys.act.x to 1.0 in the initialization system is effectively changing the guess of sys.act.x. It doesn't ensure sys.act.x will initialize to 1.0.

only want to set an observable variable start value, not necessarily add an initial equation

This doesn't make sense. How would you expect to set the initial value of an observed variable without adding an initial equation? If the variable is not part of the state realization, it needs to be solved for. E.g. if y ~ 2x was an observed equation, and you provided y => 1 to ODEProblem, it would need to solve for x to initialize the ODEProblem.

So, is there another interface for ODEProblem that allows me to set a start value for an observed variable without breaking the initialization system?

It's not breaking the initialization system. The constraints provided for initialization are more than the number of variables being constrained, and as such the system is overdetermined. This does not necessarily mean initialization will fail.

[bradcarman]

Thank's @AayushSabharwal, I mispoke about sys.act.x. Let me explain my issue in a different way. If I run the following code, I get an initialization that I don't want, I want sys.act.x to start at 1.0...

prob = ModelingToolkit.InitializationProblem(sys, 0.0) 
sol = solve(prob)
unknowns(sys) .=> sol
#=
  act₊x(t) => 0.0
 act₊dx(t) => -0.08028832692241133
=#

So then I try to set sys.act.x to 1.0 and I get the following...

prob = ModelingToolkit.InitializationProblem(sys, 0.0, [sys.act.x => 1.0]) 
sol = solve(prob)
unknowns(sys) .=> sol
#=
  act₊x(t) => 0.08028832692241133
 act₊dx(t) => 0.08028832692241133
=#

How can I get act.x(t) to start at 1.0???

[AayushSabharwal]

unknowns(sys) .=> sol

That's not correct. If you look at sol for the second case, sol.u has only one element so the broadcast unknowns(sys) .=> sol assigns both unknowns the same value. Use unknowns(sys) .=> sol[unknowns(sys)] to get the appropriate value. The initialization problem may not have the same number of unknowns as the original system.

julia> unknowns(sys) .=> sol[unknowns(sys)]
2-element Vector{Pair{SymbolicUtils.BasicSymbolic{Real}, Float64}}:
  act₊x(t) => 1.0
 act₊dx(t) => 0.08028832692241133

[ChrisRackauckas]

sol.u has only one element so the broadcast unknowns(sys) .=> sol assigns both unknowns the same value.

That's kind of nasty 😅 . But yeah it's a consequence specifically of the scalar broadcast rule. Any other case it would error with length differences, it's just if sol.u is a scalar then scalars broadcast to the whole array...

[AayushSabharwal]

In general though, regardless of broadcast rules the operation unknowns(sys) .=> sol is invalid. The Initialization problem may or may not have the same unknowns, or even the same number of unknowns, as the original ODESystem.

[ChrisRackauckas]

Yes agreed, it's just a funny coincidence that the broadcast doesn't throw a nice error in this one instance.

[bradcarman]

Wow, I learned something new, I didn't realize I could broadcast to arrays of different size, I thought that was supposed to throw an error. Anyway...

OK, so I can see the following, if I use a guess, then I have a system that solves for act.x and act.flange.v

prob = ModelingToolkit.InitializationProblem(sys, 0.0; guesses = [sys.act.x => 1.0]) 

julia> prob.f.sys
Model sys with 2 equations
Unknowns (2):
  act₊x(t) [defaults to 0]
  act₊flange₊v(t) [defaults to 0]

sol = solve(prob)

julia> unknowns(sys) .=> sol[unknowns(sys)]
2-element Vector{Pair{SymbolicUtils.BasicSymbolic{Real}, Float64}}:
  act₊x(t) => 1.0
 act₊dx(t) => 0.08028832692241133

But, the system isn't really "solving" for sys.act.x, I can supply any guess value and I get back the same value...

prob = ModelingToolkit.InitializationProblem(sys, 0.0; guesses = [sys.act.x => 12345.0]) 
sol = solve(prob)
sol[sys.act.x] #12345.0

So I believe that act.x should require an initial equation. But if I provide that initial equation, I get a warning...

julia> prob = ModelingToolkit.InitializationProblem(sys, 0.0, [sys.act.x => 1.0])
┌ Warning: Initialization system is overdetermined. 2 equations for 1 unknowns. Initialization will default to using least squares. To suppress this warning pass warn_initialize_determined = false. To make this warning into an error, pass fully_determined = true
└ @ ModelingToolkit C:\Users\bradl\.julia\packages\ModelingToolkit\OYoU3\src\systems\diffeqs\abstractodesystem.jl:1578
NonlinearLeastSquaresProblem with uType Vector{Float64}. In-place: true
u0: 1-element Vector{Float64}:
 0.0

julia> prob.f.sys
Model sys with 2 equations
Unknowns (1):
  act₊dx(t) [defaults to 0]

Solving works and gives the correct result, but this makes me wonder for more complex problems, what is the best practice?

julia> sol = solve(prob)
retcode: Success
u: 1-element Vector{Float64}:
 0.08028832692241133

julia> unknowns(sys) .=> sol[unknowns(sys)]
2-element Vector{Pair{SymbolicUtils.BasicSymbolic{Real}, Float64}}:
  act₊x(t) => 1.0
 act₊dx(t) => 0.08028832692241133

[ChrisRackauckas]

Wow, I learned something new, I didn't realize I could broadcast to arrays of different size, I thought that was supposed to throw an error. Anyway...

It's just the case of a scalar. It's this:

x = rand(4)
x .= 0 # works, turns all values to 0

it's the same as that case of "broadcasting by a scalar works", but you just accidentally hit it since braodcasting a nonlinear system just uses the sol.u which for this case is just a scalar.

It's not intentional that doesn't throw an error, but that's why.

[AayushSabharwal]

Wow, I learned something new, I didn't realize I could broadcast to arrays of different size, I thought that was supposed to throw an error.

You can't, in general. Single-element arrays are an exception:

julia> [1] .+ [1, 2, 3]
3-element Vector{Int64}:
 2
 3
 4

julia> [1, 2] .+ [1, 2, 3]
ERROR: DimensionMismatch: arrays could not be broadcast to a common size: a has axes Base.OneTo(2) and b has axes Base.OneTo(3)

But, the system isn't really "solving" for sys.act.x, I can supply any guess value and I get back the same value...

This happens because the initialization system is singular. What this means in your case is that while the equations contain act.x, if you substitute all the observed equations into equations(sys), the resulting expressions do not contain act.x. Effectively, the residual is independent of act.x. See:

julia> prob = ModelingToolkit.InitializationProblem(sys, 0.0; guesses = [sys.act.x => 1.0])
julia> prob.f.sys
Model sys with 2 equations
Unknowns (2):
  act₊x(t) [defaults to 0]
  act₊flange₊v(t) [defaults to 0]
Parameters (75):
  src₊port₊p_gas [defaults to fluid₊p_gas]
  res₊port_a₊p_gas [defaults to fluid₊p_gas]
  res₊port_b₊β [defaults to fluid₊β]
  act₊port₊β [defaults to fluid₊β]
  t
  res₊port_b₊μ [defaults to fluid₊μ]
  res₊valve₊port_b₊n [defaults to fluid₊n]
  fluid₊let_gas [defaults to 1]
  res₊valve₊base₊port_b₊ρ [defaults to fluid₊ρ]
  src₊port₊β [defaults to fluid₊β]
  src₊port₊ρ_gas [defaults to fluid₊ρ_gas]
⋮

julia> subs = map(observed(prob.f.sys)) do eq
           eq.lhs => eq.rhs
       end |> Dict;

julia> subeqs = map(equations(prob.f.sys)) do eq
           ModelingToolkit.fixpoint_sub(eq, subs)
       end;
julia> subeqs
2-element Vector{Equation}:
 0 ~ (-2res₊valve₊base₊port_a₊ρ*src₊p*ifelse(res₊area₊k > res₊valve₊base₊minimum_area, res₊area₊k, res₊valve₊base₊minimum_area)) / (((0.0001 + (4(res₊valve₊base₊port_a₊ρ^2)*(src₊p^2)) / (ifelse(src₊p > 0, res₊valve₊base₊Cd, res₊valve₊base₊Cd_reverse)^2))^0.25)*ifelse(src₊p > 0, res₊valve₊base₊Cd, res₊valve₊base₊Cd_reverse)) - act₊area*act₊port₊ρ*act₊flange₊v(t)*(1.0001^((1//2)*(-1 + 1 / act₊port₊n)))
 0 ~ 0

Note how the last equation is 0 ~ 0. The first equation also does not contain act.x. This is verified by:

julia> ModelingToolkit.vars(subeqs)
Set{Any} with 10 elements:
  res₊valve₊base₊Cd
  act₊port₊ρ
  res₊valve₊base₊Cd_reverse
  act₊port₊n
  res₊valve₊base₊port_a₊ρ
  act₊area
  src₊p
  res₊valve₊base₊minimum_area
  res₊area₊k
  act₊flange₊v(t)

Hence the residual is independent of act.x, which is why you run into that issue.

So I believe that act.x should require an initial equation. But if I provide that initial equation, I get a warning...

Yes, it will throw a warning but you do need an equation. I believe this is a limitation of the current structural_simplify being unable to inline solve linear systems.

Solving works and gives the correct result, but this makes me wonder for more complex problems, what is the best practice?

Maybe we should have a warning similar to the over- and under- determined warning which says the system is singular and identifies the variables which don't affect the residual. That would help you as a modeler identify problematic initial conditions.

[ChrisRackauckas]

Maybe we should have a warning similar to the over- and under- determined warning which says the system is singular and identifies the variables which don't affect the residual. That would help you as a modeler identify problematic initial conditions.

That would be helpful.

[AayushSabharwal]

Interestingly without fully_determined = false structural_simplify does throw an error complaining its singular, but points to the wrong variable?

julia> structural_simplify(ModelingToolkit.generate_initializesystem(sys; guesses = [sys.act.x => 1.0]))
ERROR: InvalidSystemException: The system is structurally singular! Here are the problematic variables:
 act₊flange₊v(t)

[ChrisRackauckas]

Note I had the PR #2909 which would make such cases throw, but our test cases had a lot of hits for it. Maybe it's something to consider again.