Multi-level ODESystem

[cuihantao]

Hi,

I'm looking to implement a minimal working example in ModelToolkit.jl for a control systems use case. Linear control systems can be constructed with commonly used transfer functions, which are reused to build larger control blocks. Consider a low-pass filter (a lag transfer function) described by

\dot{x} = (K u - x) / T

where u is the input variable, x is the output differential variable, and T and K are parameters.

In my example, I'm looking to chain two systems, each built from two chained lag transfer functions. The input and output of the lag transfer function, as well as the full system, is shown in the image below.

I ran into an error with the following code

using ModelingToolkit

@parameters t T K
@variables u(t) x(t)
@derivatives D'~t

lag = [D(x) ~ (K*u-x)/T]

lag1 = ODESystem(lag, name=:lag1)
lag2 = ODESystem(lag, name=:lag2)
conn1 = [0 ~ lag1.x - lag2.u]
sys1 = ODESystem(conn1,t,[],[],systems=[lag1,lag2])

lag3 = ODESystem(lag, name=:lag3)
lag4 = ODESystem(lag, name=:lag4)
conn2 = [0 ~ lag3.x - lag4.u]
sys2 = ODESystem(conn2,t,[],[],systems=[lag3,lag4])

conn12 = [0 ~ lag2.x - lag3.u,
          0 ~ lag1.u - (t-5)^2]
sys12 = ODESystem(conn12, t,[],[],systems=[sys1,sys2])

# supply parameters

u10 = [lag1.u => 26.0,
       lag1.x => 10.0,
       lag2.u => 10.0,
       lag2.x => 20.0,
      ]
u20 = [lag3.x => 20.0,
       lag3.u => 20.0,
       lag4.x => 10.0,
       lag4.u => 20.0,
       ]

p1  = [lag1.T => 2.0,
      lag1.K => 1.0,
      lag2.T => 5.0,
      lag2.K => 2.0,
      ]
p2 = [lag3.T => 2.0,
      lag3.K => 1.0,
      lag4.T => 5.0,
      lag4.K => 2.0]

tspan = (0.0,10.0)
prob = ODEProblem(sys12,[u10;u20],tspan,[p1;p2])

The error is

ArgumentError: invalid index: nothing of type Nothing

Stacktrace:
 [1] to_index(::Nothing) at ./indices.jl:297
 [2] to_index(::Array{Float64,1}, ::Nothing) at ./indices.jl:274
 [3] to_indices at ./indices.jl:325 [inlined]
 [4] to_indices at ./indices.jl:322 [inlined]
 [5] setindex! at ./abstractarray.jl:1073 [inlined]
 [6] varmap_to_vars(::Array{Pair{Operation,Float64},1}, ::Array{Variable,1}) at /home/hcui7/.julia/packages/ModelingToolkit/TAGTl/src/variables.jl:254
 [7] ODEProblem{true,tType,isinplace,P,F,K,PT} where PT where K where F where P where isinplace where tType(::ODESystem, ::Array{Pair{Operation,Float64},1}, ::Tuple{Float64,Float64}, ::Array{Pair{Operation,Float64},1}; version::Nothing, tgrad::Bool, jac::Bool, checkbounds::Bool, sparse::Bool, simplify::Bool, linenumbers::Bool, parallel::ModelingToolkit.SerialForm, eval_expression::Bool, kwargs::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}) at /home/hcui7/.julia/packages/ModelingToolkit/TAGTl/src/systems/diffeqs/abstractodesystem.jl:260
 [8] ODEProblem{true,tType,isinplace,P,F,K,PT} where PT where K where F where P where isinplace where tType(::ODESystem, ::Array{Pair{Operation,Float64},1}, ::Tuple{Float64,Float64}, ::Array{Pair{Operation,Float64},1}) at /home/hcui7/.julia/packages/ModelingToolkit/TAGTl/src/systems/diffeqs/abstractodesystem.jl:258
 [9] ODEProblem(::ODESystem, ::Array{Pair{Operation,Float64},1}, ::Vararg{Any,N} where N; kwargs::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}) at /home/hcui7/.julia/packages/ModelingToolkit/TAGTl/src/systems/diffeqs/abstractodesystem.jl:231
 [10] ODEProblem(::ODESystem, ::Array{Pair{Operation,Float64},1}, ::Vararg{Any,N} where N) at /home/hcui7/.julia/packages/ModelingToolkit/TAGTl/src/systems/diffeqs/abstractodesystem.jl:231
 [11] top-level scope at In[5]:2

Where do I have errors? Thanks.

[ChrisRackauckas]

I think for you example you really want a causal modeling aspect, so you might want to wait on https://github.com/SciML/ModelingToolkit.jl/pull/400 to be done. @shashi is picking it up so we can pick this back up then. I can find a workaround for you, but it would take a bit of time and I think you'll just be happier when proper inputs/outputs are supported. Probably not the best answer for today but hopefully the best for tomorrow (well, soon)!

[cuihantao]

Chris - Glad to learn the solution is on the way. Not in any hurry, and I'll just wait.

In my application, there could be a lot of instances of sys1, sys2, ..., given in the following. They shares the same ODE structure but with different parameters.

lag1 = ODESystem(lag, name=:lag1)
lag2 = ODESystem(lag, name=:lag2)
conn1 = [0 ~ lag1.x - lag2.u]
sys1 = ODESystem(conn1,t,[],[],systems=[lag1,lag2])

lag3 = ODESystem(lag, name=:lag3)
lag4 = ODESystem(lag, name=:lag4)
conn2 = [0 ~ lag3.x - lag4.u]
sys2 = ODESystem(conn2,t,[],[],systems=[lag3,lag4])

To create variables and assign parameters correctly, I made lag1, lag2, lag3 and lag4 here. My approach here does not seem practical if I have lots of sys.

Is there already a solution for this type of problem? Not sure how exactly it should be called. Maybe like "batch instantiation and connection of ODESystems"?

[ChrisRackauckas]

https://github.com/SciML/ModelingToolkit.jl/issues/341#issuecomment-624760102 is a nice way of doing it, but we should build it around the coming causal system since the way you're doing it (and the way that prototype does it) makes the DAE very large when in reality it should be an ODE, which just makes it numerically harder to solve. So I think what we should do is first finish the causal system and then finish the network system in a way that uses the causality.

Pinging @asinghvi17 so he's aware how this is progressing.

[cuihantao]

Thanks for the explanation. I will watch the development of causal systems. Feel free to close this issue.