ControlSystemsMTK: Easy Way to Initialize Output?

Specific Domains Modelling & Simulations
1

Is there a convenient way to the initialize the unknowns in a ODESystem created by ControlSystemsMTK based on a desired initial output?

I needed a simple time delay model for something, and I figured using pade from ControlSystems and converting it using ControlSystemsMTK would be fine. And it works, but setting the initial states for the converted system isn’t effortless.

using ControlSystemsBase, ControlSystemsMTK, Plots, ModelingToolkit, OrdinaryDiffEq
delay_tf = pade(0.8, 2, 5)
@named delay = ODESystem(delay_tf)

@Independent_variables t
connections = [delay.input.u ~ 1.0 + 0.2 * sin(t)]
sys = compose(ODESystem(connections, t, name=:bob), delay)
sys = structural_simplify(sys)

u0 = [delay.x[5] => 20.0] # Close enough for now
prob = ODEProblem(sys, u0, (0.0, 10.0))
sol = solve(prob, Tsit5())

display(plot(sol.t, sol[delay.output.u]))
sol[delay.output.u][1] == 20.0 # is false

I tried a few simple attempts at leaning on MTK’s initialization methods but didn’t find anything quickly.

My last thought was maybe the system was in a canonical form, I could figure out how to initialize the states from that, but it looks like the equations in the created model don’t follow a canonical form.

2

Specifying the output instead of one of the state variables should work?

u0 = [delay.output.u => 20.0]

If you want to obtain the delay on a canonical form, try

julia> ss(delay_tf, balance=false)
StateSpace{Continuous, Float64}
A = 
 -18.75  -187.49999999999997  -1171.8749999999998  -4394.531249999999  -7690.429687499997
   1.0      0.0                   0.0                  0.0                 0.0
   0.0      1.0                   0.0                  0.0                 0.0
   0.0      0.0                   1.0                  0.0                 0.0
   0.0      0.0                   0.0                  1.0                 0.0
B = 
 1.0
 0.0
 0.0
 0.0
 0.0
C = 
 0.0  0.0  117.18749999999997  -1757.8124999999993  7690.429687499997
D = 
 0.0

Continuous-time state-space model

but I recommend against going this route, the solver will not like this realization and will likely take longer time to perform the simulation if you avoid balancing.

3

I did try [delay.output.u => 20.0] and few variations thereof but they usually came back with an overdetermined warning.

using ControlSystemsBase, ControlSystemsMTK, ModelingToolkit, OrdinaryDiffEq
delay_tf = pade(0.8, 2, 5)
@named delay = ODESystem(delay_tf)

@independent_variables t
connections = [delay.input.u ~ 1.0 + 0.2 * sin(t)]
sys = compose(ODESystem(connections, t, name=:bob), delay)
sys = structural_simplify(sys)

u0 = [delay.output.u => 20.0]
prob = ODEProblem(sys, u0, (0.0, 10.0))
sol = solve(prob, Tsit5())

using Plots
display(plot(sol.t, sol[delay.output.u]))
sol[delay.output.u][1] == 20.0

julia> include("mwe.jl")
┌ Warning: Initialization system is overdetermined. 1 equations for 0 unknowns. Initialization will default to using least squares. `SCCNonlinearProblem` can only be used for initialization of fully determined systems and hence will not be used here. To suppress this warning pass warn_initialize_determined = false. To make this warning into an error, pass fully_determined = true
└ @ ModelingToolkit ~/.julia/packages/ModelingToolkit/KpwHT/src/systems/diffeqs/abstractodesystem.jl:1517
false

julia> sol
retcode: InitialFailure
Interpolation: specialized 4th order "free" interpolation
t: 1-element Vector{Float64}:
 0.0
u: 1-element Vector{Vector{Float64}}:
 [0.0, 0.0, 0.0, 0.0, 0.0]

I’ll check out the unbalanced version. Thanks!

4

you can unset some variables by

u0 = [delay.output.u => 20.0; collect(delay.x[:] .=> nothing)]

this should prevent it from being over-determined

This works

u0 = [delay.output.u => 20.0; delay.x[5] => nothing]#; delay.x[4] => nothing; delay.x[3] => nothing]
prob = ODEProblem(sys, u0, (0.0, 10.0))
sol = solve(prob, Tsit5())

using Plots
display(plot(sol.t, sol[delay.output.u]))
sol[delay.output.u][1] ≈ 20.0
# true