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.

[ChrisRackauckas]

So is the solution then provide some standard types like PIDController which are callable functions (u,p,t) to make this easier to scale in a type-stable way?

[baggepinnen]

So is the solution then provide some standard types like PIDController which are callable functions (u,p,t) to make this easier to scale in a type-stable way?

I think this will be hard to do since each and every implementation of a PID controller is slightly different, and most often include things like integrator anti-widup, saturations and filters :/

I have experimented a bit on my own and in my implementations I had controllers always be their own ODESystem and the user had themselves to specify how the controller system interacted with the controlled system. To calculate the control Jacoiban, the user had to manually specify with respect to what the dynamics is differentiated etc. I haven't yet landed on an implementation where this is smooth in all situations though and it's still very much experimental.

[cadojo]

Is there some way to do both function signatures for control inputs? The control function passed as an input to ODEProblem or ODEFunction has to have one input (assumed t) or three inputs (assumed x, p, t). That would let people pass in data interpolations without wrapping them in closures. The downside to this is some added complexity, and possibly confusing users.

[jonniedie]

I’m still loosely of the opinion that putting ctrl(x,p,t) and ctrl(t) under the same interface may prove to be a mistake. The two are generally treated as separate subjects and not meaningfully interchangable.

For example, space missions are designed by solving trajectory optimization problems of the form ctrl(t), but then a model with controls of the form ctrl(x,p,t) is made to follow the optimal trajectory. A lot of the time, it is two completely separate groups within an organization handling both of these problems because the there isn’t much skill overlap between the two. And tooling-wise they are pretty separate problems as well. In one case, you are either building variational equations and solving the BVP or discretizing and solving directly as an optimization problem. In the other case, you are tagging which internal parameters you want to tune your control response to and optimizing over those. Simulink has a really cool feature for this where you just mark the parameters in your model you want to optimize over as “tunable” and it will tune it to meet a set of common controls objectives that you set.

[baggepinnen]

I’m still loosely of the opinion that putting ctrl(x,p,t) and ctrl(t) under the same interface may prove to be a mistake. The two are generally treated as separate subjects and not meaningfully interchangeable.

I agree that they are not meaningfully interchangeable, but it might sometimes be very interesting to co-optimize them. Correct me if I'm wrong (no space experience at all here :P) but in the planning of space missions, the optimal control problem and the feedback control problem has time constants separated by several orders of magnitudes, and treating them as completely different problems makes total sense. Putting this in other words, the feedback controller is so fast that the controlled system can be considered completely disturbance/uncertainty free from the view of the optimal control problem. In a field like robotics this might not at all be the case, where quick movements with flexible structures causes a very tight interaction between the feedback controller that is rejecting disturbances, and the optimal trajectory, and the possibility to co-optimize them would be very valuable. It's not all that common to do so nowadays, at least as far as I know, but if the complexity this entails could be managed by future tools, it would allow us to treat these two as two components of the same joint problem.

I recognize that the ability to co-optimize open and closed-loop control can be separated from the interface to the two problems, which is perhaps what @jonniedie is primarily talking about? As long as the feedback case is represented somehow, and that open and closed-loop controllers can be both simulated and optimized together, I think MTK would open up possibilities that require quite a lot of custom implementation today.

[jonniedie]

the optimal control problem and the feedback control problem has time constants separated by several orders of magnitudes, and treating them as completely different problems makes total sense

Yep, that's exactly it. Space trajectory optimization is typically treats the vehicle as a point-mass because the attitude dynamics are much faster than what matters for the problem. But I guess that also kinda distracts from the point I was making. The real issue is not that they are separate, but that then might need to be used at the same time.

I recognize that the ability to co-optimize open and closed-loop control can be separated from the interface to the two problems, which is perhaps what @jonniedie is primarily talking about?

Yeah, this is mostly what I was getting at. If the interface is the same for both problems, I feel like it would be hard to handle cases where you want to do open-loop optimization with a model that has closed-loop controls inside of it.

[ChrisRackauckas]

• Implement Difference in Symbolics.jl
• Make DiscreteSystem use Difference instead of Differential
• In ODESystem, specialize Difference operators to be non-differential equations and updated using PeriodicCallback
  • Merge when possible
• Add history variables
• Throw parsing error if any terms are in different directions?
• Throw error if delay is forward at ODEProblem construction

@sharanry

[MuradMathematics]

Since this was discussed here I would like to add the copy I posted on Slacks

I don't know if I understood the problem correctly, but wouldn't having the controlled parameters be dictated by a separate subsystem (something along the lines of an OptimizationProblem system) be possible? I might also be biased, but in Control System descriptions the controlled parameters are typically set by a separate system (be it an MPC routine or just a plain LQR, where the Riccati equation would be the "state" of the optimization subsystem and the output would be the optimal K feedback multiplier). In MATLAB you would have such a description as a seperate "block", however I don't know if this would conflict with #1059 . I think I missed the point of the seperation p.p and p.ctrl somewhere while reading, too.

I obviously thought of a seperation between Optimizer and actual Plant, but my ideas might be too narrow for something that should be general.

[ChrisRackauckas]

Yes, that would be possible but it would require full recompilation and codegen every time the controller is changed. This is how it's done in Modelica, and I think that makes the optimal control implementation a bit more difficult. Or even just experimenting with controls.

[ronisbr]

Hi!

I did not fully understand the proposal to handle the problem. I just want say that sometimes the control function can be very complicated. In my case, it has the entire representation of a satellite embedded software (controller and estimator). Currently, I am using DEDataArray to create a workspace that has more than 2,000 parameters. Those parameters are updated inside the control function (DiscreteCallback) to generate some torques that modify the derivatives.

[ChrisRackauckas]

You can do this by connecting a parameteric function / component to an input variable, so essentially the v9 parameter types is all that's needed and a separate interface isn't required.