Parameter representation for controls

[ChrisRackauckas]

The underlying problem is that we need to represent the generated f(du,u,p,t) in a way where p can hold the necessary parameters while also being a valid object for the optimal control. Without that, #1059 will not be able to generate an ODEProblem. But what kind of ODEProblem should it be generating? This issue is to specify the interface for the forward simulation and describe why that needs to be the case.

Proposed code

In a modified version of https://diffeqflux.sciml.ai/dev/examples/optimal_control/ where dx[1] also has a parameter, i.e.

function dxdt_(dx,x,p,t)
    x1, x2 = x
    dx[1] = p.p[1] * x[2]
    dx[2] = p.ctrl[1](t)^3
end

i.e. the parameters are held in some kind of struct:

ParameterController <: AbstractParameterData
  p::P
  ctrl::C
end

In this form, we need that each ctrl is a ctrl(t) which internally holds its parameterization. Then we would just build it like:

ODEProblem(sys,[x1 => ...], [p1 => ...],ctrls = [ctrl1 => DataInterpolations.ConstantInterpolation(data,data_t)])

Sounds good? Okay now the hard part.

Parameter Representation for Optimal Control

Now we want to optimize p. That can mean multiple things: optimize p.p, optimize p.ctrl, or optimize p in full. How do we make this happen? Well we can have p itself act as a vector. If trainctrl == false, then p[i] = p.p[i]. If trainp == false then p[i] == ...?. So let's talk about that part.

We need to make the parameter interface recursive. ctrl should be something that follows the interface that it acts like a parameter vector, so that ctrl[i] is something that acts like a parameter vector, so that ctrl[i][j] is something that acts like a parameter vector. If that's the case, we just map i to continue over each next parameter vector.

DataInterpolations.jl was setup in 2018 with this in mind (gotta appreciate that forward thinking 🤔 😉 ) so if you did ConstantInterpolation{true}(u,t), the u and t would be trainable, and ConstantInterpolation(u,t) = ConstantInterpolation{false}(u,t) would just make the values u trainable and the node points t of the controller would be fixed. To make this work, we would just need to have length(ConstantInterpolation(u,t)) work, and then paramlength(p) = paramlength(p.p) + paramlength(p.ctrl). Now if you have multiple controls, p.ctrl can be a tuple of controls (or vector, depending on the type stability), and then paramlength(x::Tuple) = sum(paramlength,x). Then length(p::ParameterController ) = paramlength(p).

Given this indexing, pvec = Vector(p) would "just work" via Julia's generic dispatches. However, we would still need a way to change a vector back into p, which we should implement by overloading ArrayInterface.restructure(p,pvec) (though it might just work by default?). Once this is the case, then all of GalacticOptim, adjoints, DiffEqFlux, etc. should be able to train the values of the parameters and the controllers in this form.

Pieces that are made to be constant would just be things that are kept out of the vectorization of the parameter struct.

Where to put it

Since this is a generally nice interface, it might make sense to have this be a package that we point people towards which defines a nice structure to parameters. This is also the general solution to a long-standing issue SciML/DiffEqFlux.jl#178 and SciML/SciMLSensitivity.jl#80, it might be good to call it SciMLParameters.jl or something and make it a full interface so that it's documented and people know how to make their parameters "fit" with all of the packages in the ecosystem. Having this would work make p = TupleParams(p1,p2) work in DiffEqFlux, if it's using this paramlength where on tuple it knows to just extend. ModelingToolkit.jl would then just use this new "SciML common parameter interface"

[cadojo]

One design feature I'd propose is the ability to specify control laws that take states, non-control parameters, and time as inputs. SciML is machine-learning focused as the name suggests, so perhaps the use case I'm about to describe isn't a top use-case for primary SciML users. That said, I know (as a user) I'm primarily interested in the SciML ecosystem as an aerospace controls engineer. The spacecraft separation dynamics use-case for DifferentialEquations.jl is very close to the kinds of use-cases I'm interested in.

Since a new parameter design is up for discussion here, I thought I'd describe one desired use case as an example, and use it to motivate a tweak to the Proposed Code section above.

In short – I think it would be nice to have the control functions take x, p.p, t as inputs.

Flight Control Use-case

As a graduate student, I worked with public models for a sub-scale model aircraft that NASA uses for flight control research. Researchers published a paper which described the model aircraft dynamics, and presented nonlinear flight control analysis. The general equations of motion are shown here. I'm now realizing that LaTeX math doesn't show here in GitHub-flavored Markdown!

All symbols in the equations of motion can be divided into the following groups. Note that I'm using the terms "variables" and "parameters" with their colloquial definitions, not with ModelingToolkit syntax.

• state variables (airspeed, angle of attack, pitch rate, pitch angle)
• constant parameters (vehicle mass, etc.)
• non-constant parameters (aerodynamic coefficients, which are interpolated from pre-logged flight-data as described by the paper linked above)
• control parameters (control surface deflection, thrust magnitude (0% to 100%)

Link with ModelingToolkit

Today, the groups above map to the following ModelingToolkit components.

• state variables map to ModelingToolkit states
• constant parameters map to ModelingToolkit parameters
• non-constant parameters map to ModelingToolkit forcing functions, as specified with @register
• control parameters will map to the subset of ModelingToolkit parameters specified with the controls keyword if Control parameter specification for AbstractODESystem, DiscreteSystem #1059 merges

Addition to Proposed Code

The @register forcing-function functionality seems to only allow for time-varying parameters (or forcing functions) in your equations of motion. My understanding of aerodynamic coefficient data (I've never actually seen any tables, I'm just getting this from reading papers) is the non-constant parameters (the aerodynamic coefficients) are affected by time, and/or the state variables.

Similarly, it would be helpful if the parameter field in the ODEFunction had the following functionality.

function dxdt_(dx,x,p,t)
    x1, x2 = x
    dx[1] = p.p[1] * x[2]
    dx[2] = p.ctrl[1](x,p.p,t)^3
end

I've never personally designed subsystems with ModelingToolkit, so it's possible that the functionality I'm describing is already possible by nesting a NonlinearSystem in an ODESystem.

[JTaets]

I 100% agree with @cadojo on making the controller state and time dependent.
I am a nobody, but would like to add a remark. Ignore this if it is irrelevant.

Usually p.ctrl[1](x,p.p,t) contains zero order holds (=discontinuities), Because practically a computer can only calculate a new control at multiples of some delta_t. (i.e. it is a hybrid system)

Currently, if you want this to be differentiable with DiffEqFlux (afaik), you need to augment the state with the control vector, and have dx[...] = 0. Then the control can be updated by a PeriodicCallback.

I just want to make sure that this is taken into account if possible. See #894.

[ChrisRackauckas]

Usually p.ctrl1 contains zero order holds (=discontinuities), Because practically a computer can only calculate a new control at multiples of some delta_t. (i.e. it is a hybrid system)

That's a control that should be in a differencing term, very different.

I've never personally designed subsystems with ModelingToolkit, so it's possible that the functionality I'm describing is already possible by nesting a NonlinearSystem in an ODESystem.

You can. I wonder if having this in this form too is good, or a bad form of redundancy.

[cadojo]

You can. I wonder if having this in this form too is good, or a bad form of redundancy.

I think a benefit for providing an abstract control parameter form is you don't have to re-define your system every time you try a new controller design.

If I export a model from a package, say PolynomialGTM.jl, that model would be able to specify that "these two parameters are controls, and can be used to affect the dynamics of the system".

Without modifying the model, users could specify ctrl = [(x,p,t) -> (<whatever your control law is>)] as an argument to ODEProblem or ODEFunction.

If we don't add some functionality to specify state-varying and time-varying control laws, the user would have to construct a new model of the system for every new control law, right?

Also, leaving controls as a symbolic parameter in the dynamics allows for functionality like calculate_control_jacobian, as introduced and discussed in #1059.

[ChrisRackauckas]

If we don't add some functionality to specify state-varying and time-varying control laws, the user would have to construct a new model of the system for every new control law, right?

Good point. Should it be a function of p though, or should the control just hold its p values?

Also, leaving controls as a symbolic parameter in the dynamics allows for functionality like calculate_control_jacobian, as introduced and discussed in #1059.

Another good point.

[cadojo]

Edit: I'm realizing I didn't answer the question "should the control just hold the p values". I would argue that parameters can change as a result of discrete callbacks, so the control should take p as an argument so that the control law reflects all callbacks. But I'm sure there's another way to do it as well.

Should it be a function of p though, or should the control just hold its p values?

Now that we're talking this out, maybe controls should be entirely separate from parameters. Folks often talk about controls as "control parameters" because they typically don't have their own dynamics (when they do, they can be appended to the state vector). I wanted to avoid a breaking change, so I designed the new controls specification such that all controls have to be a subset of the parameters. This initially felt natural to me, because of the "control parameters" term that folks use.

Maybe that's not the right approach though. When you write a control law as a function, that law f is often a function of states, parameters, and time: f(x, p, t). Of course, we don't want the controls to recursively depend on themselves.

It would be great for this change to break as few things as possible, of course. Maybe controls for a system still default to an empty array, but instead of controls being a subset of parameters, they're defined separately by the user, and they're implemented separately by ModelingToolkit with a design similar to the new ParameterController <: AbstractParameterData type designed above.

For a spring-mass-damper, today a user would write...

@parameters t fₑ d k
@variables x(t) ẋ(t)
δ = Differential(t)

eqs = [
    δ(x) ~ ẋ,
    δ(ẋ)~ - d*ẋ - k*x + fₑ
]

model = ODESystem(eqs, t, [x, ẋ], [fₑ, d, k]; name = :HarmonicOscillator)

If we want fₑ to be our control, maybe we do the following...

@parameters t d k
@variables x(t) ẋ(t)
@controls fₑ
δ = Differential(t)

eqs = [
    δ(x) ~ ẋ,
    δ(ẋ)~ - d*ẋ - k*x + fₑ
]

# signature could look like ODESystem(eqs, iv = iv(eqs), states = states(eqs), params = params(eqs), ctrls = ctrls(eqs)
model = ODESystem(eqs, t, [x, ẋ], [d, k], [fₑ]; name = :HarmonicOscillator)

For ODEProblem and ODEFunction, users pass in an array of functions which all map (x, p, t) to a scalar value, and the array of functions is the same size as controls(sys).

The generated ODEFunction could look like...

function dxdt_(dx,x,p,t)
    # p.p == (d, k)
    # p.c == (fₑ) where fₑ = (x,p,t) -> (some scalar output)
    dx[1] = x[2]
    dx[2] = -p.p[1] * x[2] - p.p[1] * x[1] + p.c(x,p.p,t)
end

[ChrisRackauckas]

Now that we're talking this out, maybe controls should be entirely separate from parameters. Folks often talk about controls as "control parameters" because they typically don't have their own dynamics (when they do, they can be appended to the state vector). I wanted to avoid a breaking change, so I designed the new controls specification such that all controls have to be a subset of the parameters. This initially felt natural to me, because of the "control parameters" term that folks use.

And my proposal already does exactly that division between control parameters and other parameters. Those parameters are just stored inside of the control functions, and then the adjoint procedure would have a way to include/exclude them. I guess the only remaining question is then, since they have internal parameters, whether we want them to be functions p.c[i](x,p.p,t) or p.c[i](t). I can see the advantage of the former.

[jonniedie]

So one thing I think we should make explicit here is what we are all talking about when we say "controls". For most engineers, the term used by itself mostly means closed-loop (feedback) control. Most of the controls-related work in ModelingToolkit has been focused on open-loop optimal control that is solved offline. Both are important to have an interface for in ModelingToolkit, but they're mostly separate things and I'm not sure giving them the same interface is what we actually want.

For aerospace (and a lot of other industries), you'll need both--sometimes within the same model. That can go both ways: an offline trajectory optimization that uses a model with internal feedback control or a closed-loop feedback system that solves for open-loop optimal trajectories online (albeit at a coarser time step).

@cadojo it looks like what you're looking for here is a more structured way to do closed-loop control. I usually handle this in ModelingToolkit with alias variables (for example, fₑ=k_p*x for a proportional controller, then fₑ is just a variable that will get alias-eliminated by structural_simplify). But this has a few problems that you and @JTaets mentioned:

For one, it isn't easy to calculate the control Jacobian and linear control methods require you to handle the control Jacobian separately from the state Jacobian.

The other issue is discrete-time variables. Since your controller will eventually need to run on digital hardware, control laws are almost always written as difference equations of the discretely sampled continuous state variables. I usually handle this right now by making continuous state variables with zero dynamics and then update them in a callback. This gets a little clunky, though. It would be really nice to have a discrete independent variable that you could use to write difference equations. For example:

@parameters t
@parameters d k
@variables x(t) ẋ(t) fₑ(t)
δ = Differential(t)
n = DiscreteSampleTime(t; dt=0.01)

eqs = [
    δ(x) ~ ẋ,
    δ(ẋ) ~ - d*ẋ - k*x + fₑ
    fₑ(n+1) ~ -x(n) + x(n-1) + fₑ(n)   # not a meaningful control law, just an example
]

[ChrisRackauckas]

The other issue is discrete-time variables. Since your controller will eventually need to run on digital hardware, control laws are almost always written as difference equations of the discretely sampled continuous state variables. I usually handle this right now by making continuous state variables with zero dynamics and then update them in a callback. This gets a little clunky, though. It would be really nice to have a discrete independent variable that you could use to write difference equations. For example:

That's part of #894. You're close. My plan is to instead have that like:

@parameters t
@parameters d k
@variables x(t) ẋ(t) fₑ(t)
δ = Differential(t)
D = Difference(t; dt=0.01)

eqs = [
    δ(x) ~ ẋ,
    δ(ẋ) ~ - d*ẋ - k*x + fₑ
    D(fₑ) ~ -x + fₑ   # not a meaningful control law, just an example
]

And then we want to add delay handling, like:

eqs = [
    δ(x) ~ ẋ,
    δ(ẋ) ~ - d*ẋ - k*x(t-1) + fₑ # delay differential equation
    D(fₑ) ~ -x + x(t-dt) + fₑ   # delay term in the update equation
]

But that handling of the n thing is kind of hard, which is why it's not completed yet. But that is one of the high priority things (@YingboMa thoughts on the delays here?)

[jonniedie]

Won't the delay in the difference equation be sensitive to floating point in the case where you are getting the delayed state of a discrete variable (because you're asking it to retrieve a value right at the discontinuity)? Or are delayed discrete variables intended to be handled with higher-order difference equations?

[ChrisRackauckas]

Won't the delay in the difference equation be sensitive to floating point in the case where you are getting the delayed state of a discrete variable (because you're asking it to retrieve a value right at the discontinuity)? Or are delayed discrete variables intended to be handled with higher-order difference equations?

Symbolically that should be fine. I think this would be a bit easier to parse than the n form, since a Difference operator would work very similarly to the Differential parsing.

I guess the remaining question is just whether to have the controls be ctrl(t) or ctrl(u,p,t). The only issue with the latter is that people might think that generally using different functions there is a good idea, when it is not. With ctrl(t) you might think "oh, I should be passing interpolators like ConstantInterpolation(u,t), which then makes it all the same type. If someone starts to say u1 => param1 * t which we lower to (u,p,t)->p[5] * t and then put it all into a tuple, we can quickly end up in compile time trouble if there's hundreds of new functions every time. But maybe that is just best handled through documentation and tutorials leading people towards the right libraries of control laws.

[baggepinnen]

I guess the remaining question is just whether to have the controls be ctrl(t) or ctrl(u,p,t).

I can offer an additional control-oriented perspective:
In the field of control, optimal control in the sense of ctrl(t), i.e., open-loop control, is very much a niche. I would venture a guess that only 1 out of 1000 applications of something an engineer would call automatic control is actually referring to ctrl(t), as opposed to the far more common situation of the control signal being a function of measurements of the state (or some function of the state in general). Do applications of open-loop optimal control exist? Yes of course, but they are vastly outnumbered by simple feedback controllers.

Even in the case one actually wants to compute an open-loop optimal control trajectory, there is most likely a feedback controller in the picture making sure the system is actually following the computed trajectory. In robotics, where trajectory optimization is fairly common (perhaps not so much in industry, but in academia at least), the interaction with the servo/joint controllers can not be neglected, and an optimized trajectory will always be executed by a feedback-controlled system in the end. Simulating the interaction between the feedback controller and the feedforward system/trajectory optimizer is very important to make any use of optimal control at all.