Input/Output Functions

[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:

@variables x(t) y(t) z(t)
@derivatives D'~t

eqs1 = [D(x) ~ σ*(y-x),
       D(z) ~ x*y - β*z]

c1 = IOComponent(Inputs = [y], Outputs = [x, z], equations = eqs1)
c2 = IOComponent(Inputs = [x, z], Outputs = [y], equations = D(y) ~ x*(ρ-z)-y )

sys = IOSystem(c1, c2, [c1.x => c2.y,  c1.z => c2.z, c2.y => c1.y])

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 IOFunctions with signature f(dx, x, p, t, u), and where u takes 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:

e = edge(v_s, v_d)

(v_s and v_d are the variables at the source and destination of edge e)

Aggregate the edges per vertex:

a_v = aggregate(e_s, e_d)

where e_s and e_d are arrays of edges incident on v (they have v as source or destination respectively). Taken together edge and aggregate define a network layer, that has inputs v and outputs a_v.
Define the dynamics at the vertex

vertex!(dv, v, p, t, a_v)

This has output v and input a_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) or maximum(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.

[ChrisRackauckas]

I don't think this is quite the direction to go. #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:

1. Open aliasvars are explicit algebraic inputs, solving @asinghvi17 's problem and the problem here.
2. aliasvars can 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.
3. 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 #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.

[ChrisRackauckas]

This issue is a duplicate of #340 so it would probably be best to continue the discussion there since having two spots could get confusing.

[FHell]

Now that I understand the issue #340 better, I can focus this issue down (might be better for a new issue?).

#340 introduces incompletely specified systems. It would be extremely useful for us to be able to obtain a callable form of these incompletely specified systems with the calling convention f(dx, x, p, t, u). So we need a new DEFunction type that's called IOFunction. This was already mentioned in my initial comment but somewhat tangentially (though it was the core intended output).

Possible usecases for this that I can think of from the top of my head:

• Inputing random perturbations (not SDE) into a dynamical system.
• Designing a control by using a neural network to learn a u(x)
• Causal modelling where we combine IOFunctions (this would allow us to introduce the use of MTK in ND in the short term).

[ChrisRackauckas]

Don't those just need inputs and outputs in the AbstractSystem to work?

[FHell]

Don't those just need inputs and outputs in the AbstractSystem to work?

You can already do this via build_function since it takes any amount of args. build_function(x,p,t,u;...) would make f(dx, x, p, t, u).

I think that's the implementation plan I had in mind, but it's not completely obvious to me from reading the documentation how to put the pieces together.

[ChrisRackauckas]

Given the acausal and structural_simplify stuff, I think (?) we are good without input/output functions.

[FHell]

Agreed. Just as a pointer to people searching this issue. If anyone needs causal modeling with the guarantees (each input only gets one incoming connection) or with function with the calling signature f(dx, x, u, p, t) with u the unspecified inputs, we are developing this as a (thinning) layer on top of MTK at https://wuerfel.io/BlockSystems.jl/dev/