Extending a system with a constant parameter to allow parameter to vary

Specific Domains Modelling & Simulations
1

Hello!

Let’s say I have a ModelingToolkit system that has a constant parameter, for example a falling object with a constant mass:

function newton()
    @parameters g = -9.8
    @parameters m = 1.0
    @variables t
    @variables x(t) v(t)
    D = Differential(t)
    eq = [
        D(v) ~ m * g,
        D(x) ~ v
    ]

    @named sys = ODESystem(eq)
end

sys = newton()
@unpack x, v = sys
prob = ODEProblem(sys, [x => 0.0, v => 0.0], (0.0, 10.0), [])

sol = solve(prob, Tsit5())

plot(sol)

I would like to change the system so that the mass can vary over time. (I understand that I could just edit the function above, but let’s assume that the function above is part of a third-party library and I would prefer not to edit it.) For example, I would like m = α t, where m is the mass above, t is time, and α is some constant. I can create this second function in MTK like this:

function varmass()
    @parameters α = 2.0
    @variables t m(t)
    eq = [m ~ α * t]
    @named varmass = NonlinearSystem(eq, [t, m], [α])
end

vm = varmass()

@unpack m, t = vm
prob2 = NonlinearProblem(vm, [m => 1.0, t => 1.0], check_length=false)

sol2 = solve(prob2)

plot(sol2)

To be honest, I’m not sure what the plot that the code above creates means, and it’s not exactly a nonlinear optimization problem, so this is probably not the correct way to do it, but hopefully you get the idea of what I’m trying to do.

The final step would presumably be to combine the two systems:

@named connection = NonlinearSystem([
    sys.m ~ vm.m
], [sys.m, vm.m], [])

connected = compose(sys, vm, connection)

Above I’m just trying to say that the mass of our falling object should no longer be the default parameter value, but instead should vary according to the equation in the second system.

This, however, doesn’t work. The final line gives an error about “Cannot convert an object of type NonlinearSystem to an object of type ODESystem”.

I guess there’s some gap in my understanding of how this all works. Is there a recommended way to do this type of thing? Thank you for any help you can offer!

Chris

2

I think you should be able to just use an ODESystem constructor instead of NonlinearSystem, then use a solver where you can specify the linear solver and it should work out the same

3

Thanks for your response! Doing it with all ODESystems doesn’t solve the problem either:

using ModelingToolkit
using DifferentialEquations
using Plots

function newton()
    @variables t
    @parameters g = -9.8
    @variables m(t) = 1.0
    @variables x(t) v(t)
    D = Differential(t)
    eq = [
        D(v) ~ m * g,
        D(x) ~ v
    ]

    @named sys = ODESystem(eq)
end

function varmass()
    @parameters α = 2.0
    @variables t m(t)
    eq = [m ~ 1 + α * t]
    @named varmass = ODESystem(eq)
end

sys = newton()
vm = varmass()

@named connection = ODESystem([
    sys.m ~ vm.m
])

connected = compose(sys, vm, connection)

states(connected)
equations(connected)

simplified_sys = structural_simplify(connected)

The last line above gives the error: ExtraVariablesSystemException: The system is unbalanced. There are 3 highest order derivative variables and 2 equations.

These are the equations:

julia> equations(connected)
4-element Vector{Equation}:
 Differential(t)(v(t)) ~ g*m(t)
 Differential(t)(x(t)) ~ v(t)
 varmass₊m(t) ~ 1 + t*varmass₊α
 connection₊sys₊m(t) ~ connection₊varmass₊m(t)

So it is correct that there are 3 equations and two variables, but one of the equations should just simplify out.

If I do the same thing but just put all the equations in the same system to start, that is in fact what happens:

function newton2()
    @variables t
    @parameters g = -9.8
    @parameters α = 2.0
    @variables m(t) = 1.0
    @variables x(t) v(t)
    D = Differential(t)
    eq = [
        D(v) ~ m * g,
        D(x) ~ v,
        m ~ 1 + α * t,
    ]

    @named sys = ODESystem(eq)
end

sys2 = newton2()

simplified_sys2 = structural_simplify(sys2)

@unpack x, v, m = simplified_sys2
prob3 = ODEProblem(simplified_sys2, [x=>0, v=>0], (0.0, 10.0), [])
sol3 = solve(prob3, Tsit5())
plot(sol3)
plot!(sol3.t, sol3[m], label=m)

Is there something special I need to do to get the simplification to work across the composed systems? Thanks for your help!

4

This actually looks similar to a problem I’ve been having in MethodOfLines.jl where sometimes there are the same number of equations and states so you’d expect it to work but this error shows up.

Not sure what to do here, maybe @YingboMa has an idea.

One idea, try calling ModelingToolkit.flatten on the composed system before structural_simplify

5

Thanks. ModelingToolkit.flatten doesn’t seem to help. Here is a slightly simplified version of the code above the reproduces the error:

using ModelingToolkit
using DifferentialEquations
using Plots

function newton()
    @variables t
    @variables g(t)
    @variables x(t) v(t)
    D = Differential(t)
    eq = [
        D(v) ~ g,
        D(x) ~ v
    ]

    @named sys = ODESystem(eq)
end

function gravity()
    @variables t
    @variables g(t)
    eq = [g ~ -9.8]
    @named gravity = ODESystem(eq)
end

sys = newton()
grav = gravity()

@named connection = ODESystem([
    sys.g ~ grav.g
])

connected = compose(sys, grav, connection)

states(connected)
equations(connected)

simplified_sys = structural_simplify(connected)
6

Thanks for your help so far! I’ve created an issue over on github to continue the discussion: Unbalanced system error when using "compose" · Issue #2027 · SciML/ModelingToolkit.jl · GitHub