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.
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.
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!
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
I came to think about this thread again, the relation between the initial state x_0 and the initial output y_0 and it’s derivatives in a free response (absent control input) is given by the observability matrix \mathcal{O}. It is sometimes useful to use this property for initialization, e.g.
y0 = [1,0,0,0]
x0 = obsv(P) \ y0
1 Like