SciML / ModelingToolkit.jl

An acausal modeling framework for automatically parallelized scientific machine learning (SciML) in Julia. A computer algebra system for integrated symbolics for physics-informed machine learning and automated transformations of differential equations
https://mtk.sciml.ai/dev/
Other
1.41k stars 204 forks source link

`remake` fails with `ERROR: ArgumentError: ... is neither an observed nor an unknown variable.` #2757

Closed bradcarman closed 2 months ago

bradcarman commented 4 months ago

Using remake on the below model causes an error about "neither an observed nor an unknown variable"

using ModelingToolkit
using ModelingToolkit: t_nounits as t, D_nounits as D
using ModelingToolkitStandardLibrary.Mechanical.Translational
using ModelingToolkitStandardLibrary.Blocks
using DifferentialEquations

function MassSpringDamper(;name, mass, gravity, damping, stiffness, initial_position)

    pars = @parameters begin
        mass=mass
        gravity=gravity
        stiffness=stiffness
    end

    vars = []

    systems = @named begin
        d = Damper(; d=damping)
        m = Mass(;m=mass, g=gravity, s=initial_position)
        s = Spring(;k=stiffness)
        port_m = MechanicalPort()
        port_sd = MechanicalPort()
    end

    eqs = [
        connect(m.flange, s.flange_a, d.flange_a, port_m)
        connect(s.flange_b, d.flange_b, port_sd)
    ]

    return ODESystem(eqs, t, vars, pars; systems, name)
end

function System(; name)

    pars = @parameters begin
        gravity = 0

        wheel_mass = 25 #kg
        wheel_stiffness = 1e2
        wheel_damping = 1e4

        car_mass = 1000 #kg
        suspension_stiffness = 1e4
        suspension_damping = 10

        human_and_seat_mass = 100
        seat_stiffness = 1
        seat_damping = 1

        car_velocity = 1

        Kp=1
        Kd=1
        Ki=1
    end

    vars = @variables begin
        x(t)=0
        err(t)=0
        derr(t)=0
        dderr(t)=0
        active_force(t)=0
        dactive_force(t)=0
    end

    systems = @named begin
        wheel = MassSpringDamper(; mass=4*wheel_mass, gravity, damping=wheel_damping, stiffness=wheel_stiffness, initial_position=0.5)
        car_and_suspension = MassSpringDamper(; mass=car_mass, gravity, damping=suspension_damping, stiffness=suspension_stiffness, initial_position=1)
        seat = MassSpringDamper(; mass=4*human_and_seat_mass, gravity, damping=seat_damping, stiffness=seat_stiffness, initial_position=1.5)
        road = Position()
        force = Force()
    end

    eqs = [
        D(x) ~ car_velocity
        road.s.u ~ 0 #YF(x)

        connect(road.flange, wheel.port_sd)
        connect(wheel.port_m, car_and_suspension.port_sd)
        connect(car_and_suspension.port_m, seat.port_sd)
        connect(seat.port_m, force.flange)

        # controller
        err ~ seat.m.s - 1.5
        D(err) ~ derr
        D(derr) ~ dderr
        D(active_force) ~ dactive_force
        dactive_force ~ Kp*derr + Kp*Ki*err + Kp*Kd*dderr

        force.f.u ~ active_force
    ]

    return ODESystem(eqs, t, vars, pars; systems, name)
end

@mtkbuild sys = System()
prob = ODEProblem(sys, [], (0, 1), [sys.Kp=>1, sys.Ki=>1, sys.Kd=>1, sys.seat_stiffness=>1])
prob′ = remake(prob, p=[sys.Kp=>1])  # Error 

This gives the error..

ERROR: ArgumentError: car_and_suspension₊suspension_damping is neither an observed nor an unknown variable.
Stacktrace:
  [1] build_explicit_observed_function(sys::ODESystem, ts::SymbolicUtils.BasicSymbolic{…}; inputs::Nothing, expression::Bool, output_type::Type, checkbounds::Bool, drop_expr::typeof(RuntimeGeneratedFunctions.drop_expr), ps::Vector{…}, op::Type, throw::Bool)
    @ ModelingToolkit C:\Users\bradl\.julia\packages\ModelingToolkit\cAhZr\src\systems\diffeqs\odesystem.jl:445
  [2] build_explicit_observed_function
    @ C:\Users\bradl\.julia\packages\ModelingToolkit\cAhZr\src\systems\diffeqs\odesystem.jl:376 [inlined]
  [3] #708
    @ C:\Users\bradl\.julia\packages\ModelingToolkit\cAhZr\src\systems\diffeqs\abstractodesystem.jl:456 [inlined]
  [4] get!(default::ModelingToolkit.var"#708#721"{…}, h::Dict{…}, key::SymbolicUtils.BasicSymbolic{…})
    @ Base .\dict.jl:479
  [5] #863#generated_observed
    @ C:\Users\bradl\.julia\packages\ModelingToolkit\cAhZr\src\systems\diffeqs\abstractodesystem.jl:455 [inlined]
  [6] observed
    @ C:\Users\bradl\.julia\packages\SciMLBase\JUp1I\src\scimlfunctions.jl:4385 [inlined]
  [7] observed
    @ C:\Users\bradl\.julia\packages\SymbolicIndexingInterface\A1VUA\src\index_provider_interface.jl:100 [inlined]
  [8] _getu(sys::ODEProblem{…}, ::SymbolicIndexingInterface.ScalarSymbolic, ::SymbolicIndexingInterface.NotSymbolic, sym::SymbolicUtils.BasicSymbolic{…})
    @ SymbolicIndexingInterface C:\Users\bradl\.julia\packages\SymbolicIndexingInterface\A1VUA\src\state_indexing.jl:108
  [9] getu(sys::ODEProblem{…}, sym::SymbolicUtils.BasicSymbolic{…})
    @ SymbolicIndexingInterface C:\Users\bradl\.julia\packages\SymbolicIndexingInterface\A1VUA\src\state_indexing.jl:32
 [10] (::SciMLBase.var"#675#679"{ODEProblem{…}, SymbolicIndexingInterface.ProblemState{…}})(::Pair{Any, Any})
    @ SciMLBase .\none:0
 [11] iterate
    @ .\generator.jl:47 [inlined]
 [12] Dict{Any, Any}(kv::Base.Generator{Dict{…}, SciMLBase.var"#675#679"{…}})
    @ Base .\dict.jl:83
 [13] anydict(d::Base.Generator{Dict{…}, SciMLBase.var"#675#679"{…}})
    @ SciMLBase C:\Users\bradl\.julia\packages\SciMLBase\JUp1I\src\remake.jl:385
 [14] _updated_u0_p_symmap(prob::ODEProblem{…}, u0::Vector{…}, ::Val{…}, p::Dict{…}, ::Val{…})
    @ SciMLBase C:\Users\bradl\.julia\packages\SciMLBase\JUp1I\src\remake.jl:498
 [15] _updated_u0_p_internal(prob::ODEProblem{…}, ::Missing, p::Vector{…}; interpret_symbolicmap::Bool, use_defaults::Bool)
    @ SciMLBase C:\Users\bradl\.julia\packages\SciMLBase\JUp1I\src\remake.jl:399
 [16] _updated_u0_p_internal
    @ C:\Users\bradl\.julia\packages\SciMLBase\JUp1I\src\remake.jl:388 [inlined]
 [17] #updated_u0_p#688
    @ C:\Users\bradl\.julia\packages\SciMLBase\JUp1I\src\remake.jl:548 [inlined]
 [18] updated_u0_p
    @ C:\Users\bradl\.julia\packages\SciMLBase\JUp1I\src\remake.jl:529 [inlined]
 [19] remake(prob::ODEProblem{…}; f::Missing, u0::Missing, tspan::Missing, p::Vector{…}, kwargs::Missing, interpret_symbolicmap::Bool, use_defaults::Bool, _kwargs::@Kwargs{})
    @ SciMLBase C:\Users\bradl\.julia\packages\SciMLBase\JUp1I\src\remake.jl:95
 [20] top-level scope
    @ c:\Work\Packages\ActiveSuspension\mwe.jl:105
Some type information was truncated. Use `show(err)` to see complete types.

Note: the error occurs regardless of the parameter that is set in remake

Environment:

pkg> st
  [0c46a032] DifferentialEquations v7.13.0
  [961ee093] ModelingToolkit v9.15.0
  [16a59e39] ModelingToolkitStandardLibrary v2.7.1
AayushSabharwal commented 4 months ago

The problem here is that in MassSpringDamper, there is no @parameters damping = damping like there is for the other parameters. This is because a variable can't have a "double ParentScope" metadata.

When @named wheel = MassSpringDamper(; mass=4*wheel_mass, gravity, damping=wheel_damping, stiffness=wheel_stiffness, initial_position=0.5) is done, wheel_damping gets its SymScope metadata set to ParentScope, indicating that it is a variable from the parent of the system it is being passed to (and thus shouldn't be namespaced). If it is aliased to a new parameter inside MassSpringDamper via @parameters damping = damping, then damping doesn't have the SymScope metadata, so when it is passed to Damper, it gets assigned ParentScope. Without the aliasing,ParentScope can't be applied twice to wheel_damping which is why Damper thinks the variable is from MassSpringDamper and not System.

The long explanation is only to document the error and find a solution. It is definitely not something I'd expect a user to be concerned with, and we should either throw an error in such a scenario or handle it properly. I think DelayParentScope should be used here.

baggepinnen commented 3 months ago

Here's another example that throws the same error, but I'm not sure if the same explanation as to why applies

using ModelingToolkit
using ModelingToolkitStandardLibrary: Blocks
gravity = 9.81
gain_u1 = 0.89 / 1.4
d0 = 70
d1 = 17
n0 = 55
@named motor_dynamics = Blocks.FirstOrder(T = 0.001)
x0 = [0.85, 1, π/12, π/2]
@parameters t
@variables  u(t)[1:2]=0
u = collect(u)
@variables y(t)=0.85 v(t)=1 ϕ(t)=π/12 ω(t)=π/2
D = Differential(t)
eqs = [
    D(y) ~ v,
    D(v) ~ -gravity + gain_u1 * cos(ϕ)*(motor_dynamics.output.u + gravity/gain_u1),
    D(ϕ) ~ ω,
    D(ω) ~ -d0 * ϕ - d1 * ω + n0 * u[2],
    motor_dynamics.input.u ~ u[1]
]

@named model = ODESystem(eqs, t; systems=[motor_dynamics])

linearize(model, collect(u), [y, v, ϕ, ω])
julia> linearize(model, collect(u), [y, v, ϕ, ω])
ERROR: ArgumentError: v(t) is neither an observed nor an unknown variable.
Stacktrace:
  [1] build_explicit_observed_function(sys::NonlinearSystem, ts::Vector{…}; inputs::Nothing, expression::Bool, output_type::Type, checkbounds::Bool, drop_expr::typeof(RuntimeGeneratedFunctions.drop_expr), ps::Vector{…}, return_inplace::Bool, op::Type, throw::Bool)
    @ ModelingToolkit ~/.julia/dev/ModelingToolkit/src/systems/diffeqs/odesystem.jl:451
  [2] build_explicit_observed_function(sys::NonlinearSystem, ts::Vector{SymbolicUtils.BasicSymbolic{Real}})
    @ ModelingToolkit ~/.julia/dev/ModelingToolkit/src/systems/diffeqs/odesystem.jl:378
  [3] observed(sys::NonlinearSystem, sym::Vector{SymbolicUtils.BasicSymbolic{Real}})
    @ ModelingToolkit ~/.julia/dev/ModelingToolkit/src/systems/abstractsystem.jl:534
  [4] _getu(sys::NonlinearSystem, ::SymbolicIndexingInterface.NotSymbolic, ::SymbolicIndexingInterface.ScalarSymbolic, sym::Vector{…})
    @ SymbolicIndexingInterface ~/.julia/packages/SymbolicIndexingInterface/A1VUA/src/state_indexing.jl:158
  [5] getu(sys::NonlinearSystem, sym::Vector{SymbolicUtils.BasicSymbolic{Real}})
    @ SymbolicIndexingInterface ~/.julia/packages/SymbolicIndexingInterface/A1VUA/src/state_indexing.jl:32
  [6] linearization_function(sys::ODESystem, inputs::Vector{…}, outputs::Vector{…}; simplify::Bool, initialize::Bool, op::Dict{…}, p::SciMLBase.NullParameters, zero_dummy_der::Bool, initialization_solver_alg::NonlinearSolve.GeneralizedFirstOrderAlgorithm{…}, kwargs::@Kwargs{})
    @ ModelingToolkit ~/.julia/dev/ModelingToolkit/src/systems/abstractsystem.jl:1965
  [7] linearization_function
    @ ~/.julia/dev/ModelingToolkit/src/systems/abstractsystem.jl:1867 [inlined]
  [8] linearize(sys::ODESystem, inputs::Vector{…}, outputs::Vector{…}; op::Dict{…}, t::Float64, allow_input_derivatives::Bool, zero_dummy_der::Bool, kwargs::@Kwargs{})
    @ ModelingToolkit ~/.julia/dev/ModelingToolkit/src/systems/abstractsystem.jl:2325
  [9] linearize(sys::ODESystem, inputs::Vector{Num}, outputs::Vector{Num})
    @ ModelingToolkit ~/.julia/dev/ModelingToolkit/src/systems/abstractsystem.jl:2321
 [10] top-level scope
    @ REPL[3]:1
Some type information was truncated. Use `show(err)` to see complete types.
baggepinnen commented 2 months ago

My example in the post above no longer errors

baggepinnen commented 2 months ago

Example in OP now errors with

julia> prob = ODEProblem(sys, [], (0, 1), [sys.Kp=>1, sys.Ki=>1, sys.Kd=>1, sys.seat_stiffness=>1])
ERROR: Could not evaluate value of parameter car_and_suspension₊d₊d. Missing values for variables in expression car_and_suspension₊suspension_damping.
Stacktrace:
ChrisRackauckas commented 2 months ago

And I think that error is a good statement for what @AayushSabharwal said, so I think this is good to go.