Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Error: ArgumentError: xxx is neither an observed nor an unknown variable. #3138

Closed
bradcarman opened this issue Oct 22, 2024 · 33 comments
Closed
Assignees

Comments

@bradcarman
Copy link
Contributor

MWE (note this requires the Hydraulic branch for the ModelingToolkitStandardLibrary): The following code produces the error "ERROR: ArgumentError: mass₊s(t) is neither an observed nor an unknown variable."

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₁)
        connect(res.port₂, act.port)
        connect(act.flange, mass.flange)

        connect(fluid, src.port)
    end
end

@mtkbuild sys = System()

initialization_eqs = [
    sys.act.x ~ 1
    sys.act.dx ~ 0
]

initsys = ModelingToolkit.generate_initializesystem(sys; initialization_eqs)
initsys = structural_simplify(initsys)
initprob = NonlinearProblem(initsys, [t=>0], [])
initsol = solve(initprob; abstol=1e-6)
initsol[sys.mass.s]

The full stack trace is:

Stacktrace:
  [1] build_explicit_observed_function(sys::NonlinearSystem, ts::Num; inputs::Nothing, expression::Bool, eval_expression::Bool, eval_module::Module, output_type::Type, checkbounds::Bool, drop_expr::Function, ps::Vector{…}, return_inplace::Bool, param_only::Bool, op::Type, throw::Bool)
    @ ModelingToolkit C:\Users\bradl\.julia\packages\ModelingToolkit\OYoU3\src\systems\diffeqs\odesystem.jl:501
  [2] #observed#269
    @ C:\Users\bradl\.julia\packages\ModelingToolkit\OYoU3\src\systems\abstractsystem.jl:843 [inlined]
  [3] observed
    @ C:\Users\bradl\.julia\packages\ModelingToolkit\OYoU3\src\systems\abstractsystem.jl:820 [inlined]
  [4] #347
    @ C:\Users\bradl\.julia\packages\ModelingToolkit\OYoU3\src\systems\abstractsystem.jl:1684 [inlined]
  [5] get!(default::ModelingToolkit.var"#347#348"{…}, h::Dict{…}, key::SymbolicUtils.BasicSymbolic{…})
    @ Base .\dict.jl:479
  [6] (::ModelingToolkit.ObservedFunctionCache{NonlinearSystem})(::Num)
    @ ModelingToolkit C:\Users\bradl\.julia\packages\ModelingToolkit\OYoU3\src\systems\abstractsystem.jl:1683
  [7] observed
    @ C:\Users\bradl\.julia\packages\SciMLBase\tEuIM\src\scimlfunctions.jl:4525 [inlined]
  [8] observed
    @ C:\Users\bradl\.julia\packages\SymbolicIndexingInterface\cwAFH\src\index_provider_interface.jl:212 [inlined]
  [9] _getu(sys::SciMLBase.NonlinearSolution{…}, ::SymbolicIndexingInterface.ScalarSymbolic, ::SymbolicIndexingInterface.ScalarSymbolic, sym::Num)
    @ SymbolicIndexingInterface C:\Users\bradl\.julia\packages\SymbolicIndexingInterface\cwAFH\src\state_indexing.jl:155
 [10] getu
    @ C:\Users\bradl\.julia\packages\SymbolicIndexingInterface\cwAFH\src\state_indexing.jl:31 [inlined]
 [11] getindex(A::SciMLBase.NonlinearSolution{…}, sym::Num)
    @ SciMLBase C:\Users\bradl\.julia\packages\SciMLBase\tEuIM\src\solutions\solution_interface.jl:70

The current environment is:

Status `C:\Work\Packages\ModelingToolkitStandardLibrary.jl\Project.toml`
  [d360d2e6] ChainRulesCore v1.25.0
  [2b5f629d] DiffEqBase v6.158.1
  [615f187c] IfElse v0.1.1
⌃ [961ee093] ModelingToolkit v9.45.0
⌃ [0c5d862f] Symbolics v6.13.1
  [37e2e46d] LinearAlgebra
@bradcarman bradcarman added the bug Something isn't working label Oct 22, 2024
@AayushSabharwal
Copy link
Member

I can't reproduce your error. When running the example, I get:

julia> @mtkbuild sys = System()
ERROR: UndefVarError: `Aₒ` not defined in `ModelingToolkitStandardLibrary.Hydraulic.IsothermalCompressible`
Suggestion: check for spelling errors or missing imports.
Stacktrace:
 [1]
   @ ModelingToolkitStandardLibrary.Hydraulic.IsothermalCompressible ~/Julia/SciML/ModelingToolkit.jl/src/systems/model_parsing.jl:908
 [2] __Orifice__
   @ ~/Julia/SciML/ModelingToolkit.jl/src/systems/model_parsing.jl:138 [inlined]
 [3] #_#433
   @ ~/Julia/SciML/ModelingToolkit.jl/src/systems/model_parsing.jl:25 [inlined]
 [4] macro expansion
   @ ~/Julia/SciML/ModelingToolkit.jl/src/systems/abstractsystem.jl:2051 [inlined]
 [5]
   @ Main ~/Julia/SciML/ModelingToolkit.jl/src/systems/model_parsing.jl:138
 [6] __System__
   @ ~/Julia/SciML/ModelingToolkit.jl/src/systems/model_parsing.jl:138 [inlined]
 [7] #_#433
   @ ~/Julia/SciML/ModelingToolkit.jl/src/systems/model_parsing.jl:25 [inlined]
 [8] top-level scope
   @ ~/Julia/SciML/ModelingToolkit.jl/src/systems/abstractsystem.jl:2051
Some type information was truncated. Use `show(err)` to see complete types.

@mgithubk
Copy link

@AayushSabharwal I think that comment is related to the port names in your Orifice model. Can you share the code for your orifice model?

@AayushSabharwal
Copy link
Member

I'm using whatever is on ModelingToolkitStandardLibrary#Hydraulic

@mgithubk
Copy link

I made some changes to it that I haven't pushed yet. Also, I don't have the permissions to push them anymore and need to get them reinstated. I'll let you know when they are.

@mgithubk
Copy link

mgithubk commented Oct 22, 2024

Okay, I created a fork
https://github.com/mgithubk/ModelingToolkitStandardLibrary.jl
It contains updated library components.

There are some changes that need to be made to the model toot though.

#MWE Model Showing Initialization issue with a mass
@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()

initialization_eqs = [
sys.act.x ~ 1
sys.act.dx ~ 0
]

initsys = ModelingToolkit.generate_initializesystem(sys; initialization_eqs)
initsys = structural_simplify(initsys)
initprob = NonlinearProblem(initsys, [t=>0], [])
initsol = solve(initprob; abstol=1e-6)
u0 = unknowns(sys) .=> initsol[unknowns(sys)] # ERROR: ArgumentError: mass₊s(t) is neither an observed nor an unknown variable.

prob = ODEProblem(sys, u0, (0, 1), [])
sol = solve(prob, Rodas5P())
plot(sol; idxs = sys.act.x)

@bradcarman
Copy link
Contributor Author

I updated the Hydraulic branch of ModelingToolkitStandardLibrary to include these changes needed. As @mgithubk the script to produce the error is now...

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

initialization_eqs = [
    sys.act.x ~ 1
    sys.act.dx ~ 0
]

initsys = ModelingToolkit.generate_initializesystem(sys; initialization_eqs)
initsys = structural_simplify(initsys)
initprob = NonlinearProblem(initsys, [t=>0], [])
initsol = solve(initprob; abstol=1e-6)
initsol[sys.mass.s]

@AayushSabharwal
Copy link
Member

Apparently this happens because sys.mass.s doesn't appear in any of the equations in initsys. It's not in the algebraic equations, it doesn't have a default, and isn't in initialization_eqs so it never shows up. As a result, structural_simplify completely ignores it.

@bradcarman
Copy link
Contributor Author

And yet it's an unknown of sys...

julia> unknowns(sys)
6-element Vector{SymbolicUtils.BasicSymbolic{Real}}:
 act₊x(t)
 act₊rho(t)
 mass₊s(t)
 mass₊v(t)
 act₊p(t)
 act₊drho(t)

How is this possible?

@AayushSabharwal
Copy link
Member

I don't understand why it isn't possible 😅 mass.s is an unknown of sys because of the equation D(mass.s) ~ mass.v.

Let me ask this question: how is it meant to solve for the initial value of mass.s, given that the only equation it's involved in is D(mass.s) ~ mass.v and it doesn't have a default? Equivalently, consider this system:

@variables x(t) y(t)
@mtkbuild sys = ODESystem([D(x) ~ y, y ~ 3t], t)

If I now do ODEProblem(sys, [], (0.0, 1.0)) what determines the initial value of x?

@AayushSabharwal
Copy link
Member

AayushSabharwal commented Oct 23, 2024

The initialization system is formed from:

  1. Algebraic equations
  2. Observed equations
  3. Defaults
  4. Initialization equations
  5. variable => value mappings passed to ODEProblem

mass.s appears in none of these. So initialization can't solve for it.

@ChrisRackauckas ChrisRackauckas removed the bug Something isn't working label Oct 23, 2024
@mgithubk
Copy link

@AayushSabharwal Would mass.s be in the defaults mentioned above If we set a guess for s?

@component function Mass(; name, m, g = 0)
pars = @parameters begin
m = m
g = g
end
@nAmed flange = MechanicalPort()

vars = @variables begin
    s(t), _[guess = 0]_
    v(t)
    f(t)
end

eqs = [flange.v ~ v
       flange.f ~ f
       D(s) ~ v
       D(v) ~ f / m + g]

return compose(ODESystem(eqs, t, vars, pars; name = name),
    flange)

end

@AayushSabharwal
Copy link
Member

Guesses are different. Guesses are the u0 used for the NonlinearProblem that is solved during intialization. Defaults are constraints to be respected during initialization.

@mgithubk
Copy link

I don't understand why it isn't possible 😅 mass.s is an unknown of sys because of the equation D(mass.s) ~ mass.v.

Let me ask this question: how is it meant to solve for the initial value of mass.s, given that the only equation it's involved in is D(mass.s) ~ mass.v and it doesn't have a default? Equivalently, consider this system:

@variables x(t) y(t)
@mtkbuild sys = ODESystem([D(x) ~ y, y ~ 3t], t)

If I now do ODEProblem(sys, [], (0.0, 1.0)) what determines the initial value of x?

Would another way of making your point here be that the system has not steady-state value - since there is no limit to x and no force resisting the motion of the mass?

@AayushSabharwal
Copy link
Member

No, steady state values are unrelated.

@variables x(t) y(t)
@mtkbuild sys = ODESystem([D(x) ~ y, y ~ 3t], t; defaults = [x => 2y + 1])

x still doesn't have a steady state, but it has a default which turns into the initialization equation x ~ 2y+1 and provides a method to determine the initial value of x.

@AayushSabharwal
Copy link
Member

AayushSabharwal commented Oct 23, 2024

I think the best question to think about regarding your problem is, for the system in this issue if you were given initial values for all other variables, how would you then calculate the initial value of mass.s?

@mgithubk
Copy link

No, steady state values are unrelated.

@variables x(t) y(t)
@mtkbuild sys = ODESystem([D(x) ~ y, y ~ 3t], t; defaults = [x => 2y + 1])

x still doesn't have a steady state, but it has a default which turns into the initialization equation x ~ 2y+1 and provides a method to determine the initial value of x.

Oh, I think I got it. At t=0 there is a means to determine x.

@mgithubk
Copy link

I think the best question to think about regarding your problem is, for the system in this issue if you were given initial values for all other variables, how would you then calculate the initial value of mass.s?

Well, if you have act.x then you know mass.s

@mgithubk
Copy link

So the model works if you change the initial conditions accordingly:

initialization_eqs = [ sys.act.x ~ 1 sys.mass.s ~ 1 sys.act.dx ~ 0 ]

However, could the initialization check have figured that out for us?

@AayushSabharwal
Copy link
Member

MTK has no reason to assume and no basis to figure out that act.x is initially 1. Is that initial value valid for all possible systems?

@ChrisRackauckas
Copy link
Member

ChrisRackauckas commented Oct 23, 2024

I think the question is, how are you expecting MTK to know that mass.s at time 0 is 1? You're making some assumption that MTK does not know about, and for pedagogical reasons it would be good to know where you believe this information has been stated to the system.

@mgithubk
Copy link

what I mean is that when I initialize, it would be good if I got an error indicated that sys.mass.s is missing.

@ChrisRackauckas
Copy link
Member

It tells you exactly that???

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

initialization_eqs = [
    sys.act.x ~ 1
    sys.act.dx ~ 0
]

ODEProblem(sys, [], (0.0,1.0); initialization_eqs)
ERROR: Initialization incomplete. Not all of the state variables of the
DAE system can be determined by the initialization. Missing
variables:

SymbolicUtils.BasicSymbolic{Real}[mass₊s(t)]

Stacktrace:
  [1] ModelingToolkit.InitializationProblem{…}(sys::ODESystem, t::Float64, u0map::Dict{…}, parammap::SciMLBase.NullParameters; guesses::Dict{…}, check_length::Bool, warn_initialize_determined::Bool, initialization_eqs::Vector{…}, fully_determined::Bool, check_units::Bool, kwargs::@Kwargs{…})
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/CCIoy/src/systems/diffeqs/abstractodesystem.jl:1571
  [2] (ModelingToolkit.InitializationProblem{…})(::ODESystem, ::Float64, ::Vararg{…}; kwargs::@Kwargs{…})
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/CCIoy/src/systems/diffeqs/abstractodesystem.jl:1519
  [3] InitializationProblem
    @ ~/.julia/packages/ModelingToolkit/CCIoy/src/systems/diffeqs/abstractodesystem.jl:1518 [inlined]
  [4] #InitializationProblem#931
    @ ~/.julia/packages/ModelingToolkit/CCIoy/src/systems/diffeqs/abstractodesystem.jl:1507 [inlined]
  [5] InitializationProblem
    @ ~/.julia/packages/ModelingToolkit/CCIoy/src/systems/diffeqs/abstractodesystem.jl:1506 [inlined]
  [6] process_DEProblem(constructor::Type, sys::ODESystem, u0map::Vector{…}, parammap::SciMLBase.NullParameters; implicit_dae::Bool, du0map::Nothing, version::Nothing, tgrad::Bool, jac::Bool, checkbounds::Bool, sparse::Bool, simplify::Bool, linenumbers::Bool, parallel::Symbolics.SerialForm, eval_expression::Bool, eval_module::Module, use_union::Bool, tofloat::Bool, symbolic_u0::Bool, u0_constructor::typeof(identity), guesses::Dict{…}, t::Float64, warn_initialize_determined::Bool, build_initializeprob::Bool, initialization_eqs::Vector{…}, fully_determined::Bool, check_units::Bool, kwargs::@Kwargs{…})
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/CCIoy/src/systems/diffeqs/abstractodesystem.jl:917
  [7] process_DEProblem
    @ ~/.julia/packages/ModelingToolkit/CCIoy/src/systems/diffeqs/abstractodesystem.jl:835 [inlined]
  [8] (ODEProblem{…})(sys::ODESystem, u0map::Vector{…}, tspan::Tuple{…}, parammap::SciMLBase.NullParameters; callback::Nothing, check_length::Bool, warn_initialize_determined::Bool, eval_expression::Bool, eval_module::Module, kwargs::@Kwargs{…})
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/CCIoy/src/systems/diffeqs/abstractodesystem.jl:1107
  [9] ODEProblem
    @ ~/.julia/packages/ModelingToolkit/CCIoy/src/systems/diffeqs/abstractodesystem.jl:1095 [inlined]
 [10] (ODEProblem{true})(::ODESystem, ::Vector{Any}, ::Vararg{Any}; kwargs::@Kwargs{initialization_eqs::Vector{Equation}})
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/CCIoy/src/systems/diffeqs/abstractodesystem.jl:1082
 [11] ODEProblem(::ODESystem, ::Vector{Any}, ::Vararg{Any}; kwargs::@Kwargs{initialization_eqs::Vector{Equation}})
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/CCIoy/src/systems/diffeqs/abstractodesystem.jl:1071
 [12] top-level scope
    @ ~/Desktop/test.jl:273
Some type information was truncated. Use `show(err)` to see complete types.

@mgithubk
Copy link

I get this error when I comment out the initial condition for mass.s:

ERROR: ArgumentError: mass₊s(t) is neither an observed nor an unknown variable. Stacktrace: [1] build_explicit_observed_function(sys::NonlinearSystem, ts::Vector{…}; inputs::Nothing, expression::Bool, eval_expression::Bool, eval_module::Module, output_type::Type, checkbounds::Bool, drop_expr::Function, ps::Vector{…}, return_inplace::Bool, param_only::Bool, op::Type, throw::Bool) @ ModelingToolkit ~/.julia/packages/ModelingToolkit/Jr7hI/src/systems/diffeqs/odesystem.jl:487 [2] observed(sys::NonlinearSystem, sym::Vector{SymbolicUtils.BasicSymbolic{…}}; eval_expression::Bool, eval_module::Module) @ ModelingToolkit ~/.julia/packages/ModelingToolkit/Jr7hI/src/systems/abstractsystem.jl:834 [3] observed @ ~/.julia/packages/ModelingToolkit/Jr7hI/src/systems/abstractsystem.jl:811 [inlined] [4] #343 @ ~/.julia/packages/ModelingToolkit/Jr7hI/src/systems/abstractsystem.jl:1663 [inlined] [5] get!(default::ModelingToolkit.var"#343#344"{…}, h::Dict{…}, key::Vector{…}) @ Base ./dict.jl:479 [6] ObservedFunctionCache @ ~/.julia/packages/ModelingToolkit/Jr7hI/src/systems/abstractsystem.jl:1662 [inlined] [7] observed @ ~/.julia/packages/SciMLBase/tEuIM/src/scimlfunctions.jl:4525 [inlined] [8] observed @ ~/.julia/packages/SymbolicIndexingInterface/cwAFH/src/index_provider_interface.jl:212 [inlined] [9] _getu(sys::SciMLBase.NonlinearSolution{…}, ::SymbolicIndexingInterface.NotSymbolic, elt::SymbolicIndexingInterface.ScalarSymbolic, sym::Vector{…}) @ SymbolicIndexingInterface ~/.julia/packages/SymbolicIndexingInterface/cwAFH/src/state_indexing.jl:260 [10] getu @ ~/.julia/packages/SymbolicIndexingInterface/cwAFH/src/state_indexing.jl:31 [inlined] [11] getindex(A::SciMLBase.NonlinearSolution{…}, sym::Vector{…}) @ SciMLBase ~/.julia/packages/SciMLBase/tEuIM/src/solutions/solution_interface.jl:79 [12] top-level scope @ ~/Documents/GitHub/SimpleActuationMWE/mtksl.jl:123 Some type information was truncated. Use show(err) to see complete types.

@AayushSabharwal
Copy link
Member

Is this with the code in #3138 (comment)?

@AayushSabharwal
Copy link
Member

I think I misunderstood you on the call then. The initializationsystem alone will never tell you what's missing or not. ODEProblem tells you what's missing before constructing the initialization system

@mgithubk
Copy link

It is with this code:

`@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()

initialization_eqs = [
sys.act.x ~ 0.001
sys.act.dx ~ 0
#sys.mass.s ~ 0.001
]

initsys = ModelingToolkit.generate_initializesystem(sys; initialization_eqs)
initsys = structural_simplify(initsys)
initprob = NonlinearProblem(initsys, [t=>0], [])
initsol = solve(initprob; abstol=1e-6)
u0 = unknowns(sys) .=> initsol[unknowns(sys)] # ERROR: ArgumentError: mass₊s(t) is neither an observed nor an unknown variable.

prob = ODEProblem(sys, u0, (0, 1), [])
sol = solve(prob, Rodas5P())
plot(sol; idxs = sys.act.x)`

@AayushSabharwal
Copy link
Member

AayushSabharwal commented Oct 23, 2024

Yeah so you shouldn't be solving the initializationsystem manually. Constructing the ODEProblem and solving it does that for you, and will complain that mass.s is missing an initial condition

@ChrisRackauckas
Copy link
Member

That's because you're not initializing, you're purposefully skipping the standard pipeline to user lower level debugging / development utilities and those obviously cannot error because then there would be no way to look what is being computed. I don't know what to say other than, use the documented standard path instead of the internals if you want the standard errors. ModelingToolkit.generate_initializesystem is a piece of code only documented as a way to recreate the internals and skip the errors in order to understand why the symbolic manipulation is saying what it says (right here https://docs.sciml.ai/ModelingToolkit/stable/tutorials/initialization/#Diving-Deeper:-Constructing-the-Initialization-System).

@ChrisRackauckas
Copy link
Member

Let me bring this to a higher level. The way symbolic-numeric compilation works is:

  1. You generate a simplified symbolic representation of the ODESystem
  2. You generate a symbolic representation of the NonlinearSystem that is solved during initialization
  3. You check that these representations are well-defined
  4. You perform code generation
  5. You run the numerical code

As a debugging procedure, and because it's an open compiler that is built in Julia and thus generally all internals are accessible, we allow people to do (2) in an interactive session, which can be a nice way to ask "why is it not able to solve for x in the initialization?", since you can see the full set of equations, you can see what it simplifies to, and you can run the nonlinear solve independently of the ODE and understand why it's computing the values it does and the residual of those values.

Your question is essentially, "I did (2) so why didn't (3) run?". If the code that creates the symbolic representation errored if the representation is invalid, then it would be impossible to ever investigate why it's invalid. Thus it makes sense to have one function generate the initialization system, and then another function check whether the generated initialization system is fully well-defined and throw errors / warnings. You are skipping "normal MTK" and just calling that first part, so you get what you asked for.

So then the next question is, why are you doing this? What led you to this code? Since the only place in the entire documentation which mentions this is right here https://docs.sciml.ai/ModelingToolkit/stable/tutorials/initialization/#Diving-Deeper:-Constructing-the-Initialization-System, and it says it's a debugging utility to do in order to investigate initialization, I do not understand how with our documentation how you can get to this code without already doing the "normal thing". So I'm really curious, what led you to write or try this code? If it's some tutorial that's outside the docs, we should really clean that up.

@bradcarman
Copy link
Contributor Author

Thanks @ChrisRackauckas and @AayushSabharwal this is clear now. I was using the generate_initializesystem to get better insight to what was happening and try to determine how many initial conditions are needed. The error that mass.s(t) was missing, this is clear now, I was taking this as a different meaning originally. So now I think the workflow is clear, one can do the following to back out how many initial conditions are needed by carefully reading the errors and warnings.

prob = ODEProblem(sys, [], (0,1)) # ERROR: Initialization incomplete. Missing ... SymbolicUtils.BasicSymbolic{Real}[mass₊s(t)]
prob = ODEProblem(sys, [sys.mass.s => 1], (0,1)) # Warning: Initialization system is underdetermined. 2 equations for 4 unknowns.
prob = ODEProblem(sys, [sys.mass.s => 1, sys.act.x => 1], (0,1)) # Warning: Initialization system is underdetermined. 2 equations for 3 unknowns.
prob = ODEProblem(sys, [sys.mass.s => 1, sys.act.x => 1, sys.act.dx => 1], (0,1)) #OK

@ChrisRackauckas
Copy link
Member

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

I wonder how hard in the underdetermined cases here it would be to bubble up some message of what equation is likely missing. The problem is that there's likely an unlimited number of things you could define to make it determined, so it might just not be possible.

@AayushSabharwal
Copy link
Member

If we can get our hands on var_eq_matching and a couple other data structures, we can have a message similar to the error when you don't pass fully_determined = false where it reports unmatched variables. However, that choice is still largely arbitrary.

@ChrisRackauckas
Copy link
Member

Yeah it's just a heuristic so I'm not sure how helpful it ends up being in practice.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

4 participants