ChrisRackauckas
commented
4 years ago I don't think this is quite the direction to go. https://github.com/SciML/ModelingToolkit.jl/issues/340 is better because it makes it at the AbstractSystem level, i.e. we don't need an IOSystem that is different from the other systems if we do this correctly. Instead what we really need is just a way to add causal models to the standard systems.
@variables x(t) y(t) z(t) J(t)
@derivatives D'~t
eqs = [D(x) ~ σ*(y-J),
D(z) ~ x*y - β*z]
vars = [x,z]
aliasvars = [y,J]
aliaseqs = [J ~ x - 0.05sin(x)]
c1 = ODESystem(eqs,t,vars,[],aliasvars = aliasvars, aliaseqs = aliaseqs , name = :c1)
c2 = ODESystem([D(y) ~ J*(ρ-z)-y],t,[y],aliasvars = [J,z], name= :c2)
connectioneqs = [c1.y ~ c2.y, c2.J ~ c1.J, c2.z ~ c1.z]
connected = ODESystem(connectioneqs,t,[],[],name = :connected)Notice how this solves a few birds in one stone:
- Open
aliasvarsare explicit algebraic inputs, solving @asinghvi17 's problem and the problem here. aliasvarscan get an explicit expression, which makes it an explicit algebraic variable defined by other pieces, which is what you call an output. So this simplifies the lowering process by having these together.- At the higher level, the connections just add new connection equations, just like the acausal system. We can lower this to just be an algebraic equation without a differential part if we notice that it's an
aliasvar, so again no extra machinery required and this'll work on all equation types. A well-formed system that is numerically computable then just has, at it's highest level, all aliasvars must have an aliaseq associated with it. These are then just dropped in above the DE calculation and computed in the let expansion.
We can then put a ControlSystem above this setup that makes it nicer to write such control models with a u part etc. but it would be good to have the compilation process all target ODESystem so that way we have a common spot to put all of our optimization and transformation logic. If I'm not mistaken, this handles all of your concerns and is what is proposed in https://github.com/SciML/ModelingToolkit.jl/issues/340.
The only thing that isn't handled is this open "reduce by sum" thing. Now that you all have seen the syntax, maybe you all might have a suggestion.
FHell
It would be very useful to have Input/Output dynamical systems based on ModelingToolkit, we intend to start developing these and use this issue as a place to discuss and gather input on them, and discuss the development. I hope this is appropriate.
The fundamental idea is a slight simplification of the notions of control theory. Every component has input variable, output variables and equations
Inputs = variables whose dynamics are not determined by the component equations Outputs = variables whose dynamics are determined by the component equations
This is additional information that can not be extracted easily from the equations. In some cases it can be merely convention, which variable we think of as determining the other.
Now we can combine components by identifying outputs and inputs. I imagine a syntax like this, but I haven't prototyped it, so I have no idea whether this is actually a convenient interface in practice:
Trying to build an IOSystem should error when there are open Inputs that do not occur as Outputs (that would mean there are variables with undetermined dynamics). I believe this is straightforward to implement in ModelingToolkit, I simply haven't had the time to digest the internals of MTK enough to just implement it.
An advantage of this is that it should be possible to build complex components from simpler ones iteratively, that you can connect several inputs to an output in a straightforward manner, and that we don't turn an ODE composition into a DAE that we then need to simplify again.
Finally, for us it would be necessary to generate not just complete DEFunctions, but also incomplete
IOFunctionswith signaturef(dx, x, p, t, u), and whereutakes the unspecified inputs.The main use case we have in mind is to use ModelingToolkit in Network Dynamics. We have iterated on the design of ND.jl a number of times now.
The interface we are working on right now (mostly to allow with dealing with different layers, and to make the fusion with MTK more straightforward) is composed of three functions, vertex dynamics, aggregator, edge dynamics. Conceptually they interact like this:
Calculate the state of each edge:
(
v_sandv_dare the variables at the source and destination of edge e)Aggregate the edges per vertex:
where
e_sande_dare arrays of edges incident onv(they have v as source or destination respectively). Taken togetheredgeandaggregatedefine a network layer, that has inputsvand outputsa_v. Define the dynamics at the vertexThis has output
vand inputa_v, thus together with the network layer it defines a complete system.NB: The current interface of Network Dynamics has aggregate and vertex rolled into one function
vertex!(dv, v, e_s, e_d, p, t,). Thus while the above sketched version is already implemented, we want to allow for other types of layers as well. The goal is to multiple different parallel layers, which turned out to be important for our use case. One other use case is to simply pass through input/output variables along directed edges.Aggregate function will typically be something like
sum(e_s) - sum(e_d)ormaximum(e_s). This is difficult to express in MTK when I last looked, and this type of terms are contained in all cases we have looked at, limiting the applicability of MTK for us.The new interface allows a clean use of MTK. We can define vertex and edge as the IOFunctions (or variants thereof) I outlined above, and the aggregate function tells us how edges combine at nodes.
Now in this we don't build one large ODE System and compile it using MTK, but build smaller ones and combine them in a network dynamics object. If there are any suggestions for going the other way, essentially building things like aggregation among varying number of variables into MTK equations, that would be a very interesting alternative design that we would also be happy to invest in. We already saw that this is in some sense possible, you can successfully run modelingtoolkitize on a networkdynamics.jl right hand side function and get an MTK expression for the whole network.