DifferentialEquations, Control Systems, and Linearization

General Usage
,
1

I have a more general question regarding Differential Equations, Control Systems, and Linearization. Thanks to Tamas_Papp, ssfrr, DrPapa, ChrisRackauckas, and longemen3000 for help in a previous thread on ‘ForwardDiff & “complicated” functions?’.

Systems in Control Engineering are often described as

\frac{dx}{dt} = f(x,u,t;p) \\ y = g(x,u,t;p)

[state space form, with obvious generalizations to DAEs, and other types of models. x is normally used for states, u for inputs, and p or \theta for parameters/constants.] Here, u is often considered a time function, but in (state) feedback control, u=u(x,t).

Differential Equations solvers assume that the input function u(x,t) already is specified, and that the differential equation is of form:

\frac{dx}{dt} = h(x,t;p)

which is straightforward to set up by h(x,t;p) = f(x,u(x,t),t;p). This works fine with Julia’s DifferentialEquations package.

The ModelingToolkit.jl package has a tool for finding the Jacobian \frac{\partial h}{\partial x}, and I understand that with some tinkering, it is possible to find \frac{\partial f}{\partial u}. The previous thread also showed how I can find \frac{\partial f}{\partial x}, \frac{\partial f}{\partial u} via the ForwardDiff.jl package.

However, in “real life” control models, u would be implemented as function u(x,t). How do I then find the Jacobian \frac{\partial f}{\partial u}?

More concretely… Suppose I have the “trivial” model:

# Model with input function
function model(dx,x,p,t,u)
   dx[1] = -p[1]*x[1] + p[2]*u(t)
end

To prepare this model for DifferentialEquations.jl, I introduce:

md(dx,x,p,t) = model(dx,x,p,t,u)

and then I define an input function, e.g.:

u(t) = t < 5 ? 0.0 : 1.0

(as an example; a step function).

Of course, I can go through my model and replace u(t) with u, do the linearization with ForwardDiff, etc. And then reinsert u(t) to do the actual simulation after I have found the Jacobians/linearized model (which is often used for controller design). But this seems error prone. And for more complex models, it seems to be a hazzle.

OR: will ForwardDiff.jl or ModelingToolkit.jl handle Jacobians when u in the function call model(dx,x,p,t,u) actually is a function?

2

ModelingToolkit handles functions in its derivative definitions, it just assumes in modelingtoolkitize that all parameters are constants. Automatically detecting when that’s not true is a little difficult. ForwardDiff.jl will have a bit of a difficulty though, since changing u from a function to a number is a difficult translation.

Note that model(dx,x,(p,u),t) (i.e. using destructuring) might be a nice way to write this.

If there’s a lot of specific things that need to be done on this kind of system, which is indeed the case, then a specific ODEControlProblem and DAEControlProblem which takes in a control variable u and specifically uses these tools to correctly build Jacobians would be a good idea. I made DiffEqOptimalControl.jl to do this, but never really got around to it yet. That said, we can definite hash out a few interface ideas here to get something useful implemented. I know a few other people who might want to make use of these.

The Bigger Deal

At a higher level, instead of creating ControlProblem, I think the real problem is that we don’t have a great interface on handling parameters at a very high level. What I mean is, people keep asking about doing things like, mixing cache variables and functions with parameters and then trying to estimate only some parameter vector p. I don’t want to say “it has to be a vector”, though that’s kind of true for our current parameter estimation tools (it’s required for example with the adjoint sensitivity analysis). What this is pointing to is that, what we’re actually missing is a system for specifying parameters.

So here’s the deal. In the documentation we say “parameters”, but what we actually mean is “user data”. What we need is a way to specify what is the “parameters” portion of “user data”. There is a pretty sensible sensible I have in mind. If you give us a Vector, we can assume that the vector is our parameter vector. If you give us any other type, it might make sense to assume the .p portion is the parameters. This would mean that your kind of approach would be like (u=u_func,p = ...) would be an easy way to split the function and we’d know to only estimate/change/etc. p. But if it’s all type-based, then we can do things like

getparameters(userdata::Vector) = userdata
getparameters(userdata) = hasfield(userdata,:p) ? userdata.p : nothing

and then users could define their own getparameters and setparameters! for their type. Then this would make explicit parameter handling with implicit data, like user function, cache variables, etc., much easier to use with this kind of functionality. Then I would pull it out of DiffEq and try to make it a more Julia-wide convention.

This would intersect with a few things I’ve been talking with @DrPapa about in uncertainty quantification, with @jessebett in DiffEqFlux, and @MikeInnes with alternatives to implicit parameterization in Flux. @mauro3 might have some ideas from doing Parameters.jl as well. But in general I think this direction of “pass userdata, but I need information to understand what the parameters p part is!” is in general something that is coming up in a lot of modeling domains, and a solution specific to control problems may be the wrong direction.

6 Likes
3

I have definitely implemented this a few times in my life, but I haven’t come up with a simple scheme. And maybe there is no simple scheme, in that case the options would be to learn and use some complicated DealWithParas.jl package (which someone first has to write of course) or roll my own. Often people end up rolling their own because they can’t be bothered to learn the package. By the time they think it would be good to use the package, they are too entrenched in their hand-rolled solution to do anything about it. At least, that’s how I can see my future self… nonetheless, a solution would be great.

1 Like
4

This can’t really be a “roll your own” kind of thing, because if DiffEq doesn’t “know about it”, then it can’t take the adjoint of your userdata with respect to it.

1 Like
5

Well, I meant handling of parameters in general. But sure, DiffEq will need to have an API for this case, and presumably others.

1 Like
6

I, too, have been dipping my feet into the bog of how to store/represent model user data. The solution that I have ended up using quite a lot was mainly developed to solve a somewhat different issue, but it could still be helpful (at least as a discussion point).

My issue was mainly one of organising large-ish numbers of different parameters sets, estimated for different models and different targets. My solution was to introduce an abstract parameter type and a parameter collection type which mostly utilises overloads to getindex and friends to behave like parameter vectors while still containing data which can be used in other ways (for me this is mainly important for filtering parameter sets out from a parameter set collection).

A main benefit of relying on types with getindex/setindex overloads to access a parameter vector is that they just work with anything that expects a vector of parameters.

I wrote a small package for this which I think has a grand total of two users (me and @Gaussia) and has not been published/released in any way -
DiffEqParameters.jl. It’s very immature and it could do with re-naming but it has been helpful to me.

This may be overkill for the purposes discussed here but I thought it close enough to the issue at hand to warrant mentioning.

5 Likes
7

I believe lens is the best API for this. As I noted before, lens can be used to specify arbitrary nested fields along which you are taking the derivative, for example. This can be done without defining any methods and works with vanilla Julia objects without any use of symbolic math framework (although they are very useful, of course). Because of this flexibility, I use it anywhere I need to tweak model parameters/inputs. See this example in Bifurcations.jl where I specify the parameter of a plain named tuple by a simple lens (param_axis = @lens _.i). SteadyStateFit.jl uses lenses as “input setter” (like the argument u of f(x, u, t; p) in the OP) and the model parameter specifier. I also created a toolkit called Kaleido.jl that can be used to, e.g., add a constraint to a model and re-parametrize a model using TransformVariables.jl.

A simple implementation ControlProblem would be such that it generates an ODE parameter ControlProblemWrapper like below on the fly:

struct ControlProblemWrapper
    f!     # original dynamics
    p      # original parameter
    input  # time-to-input mapping
    lens   # input specifier
end

function f!(du, u, p::ControlProblemWrapper, t)
    q = set(p.p, p.lens, p.input(t))
    p.f!(du, q, t)
end

which can be used as

# Model with input function
function model!(dx,x,p,t)
   dx[1] = -p[1]*x[1] + p[2]*p[3]
end

p1 = ControlProblemWrapper(
    model!,
    (1.0, 2.0, 3.0),  # 3.0 as a place holder
    u(t) = t < 5 ? 0.0 : 1.0,
    (@lens _[3]),
)

Now it is super easy to make other part of the model time-dependent:

p2 = ControlProblemWrapper(
    model!,
    (1.0, 2.0, 3.0),
    u(t) = t < 5 ? 0.0 : 1.0,
    (@lens _[1]),  # make the first parameter time-dependent
)

Note that above lenses are very simple ones. You can use it to specify (@lens _.arbitrary.nested[:complex].fields[1]).

11 Likes
8

Hello, I’m fresh off the boat to using Julia (coming from Python and Stan), and I am looking to do more or less this same type of problem. Has there been any resolution to this discussion?

As a simple case, I have a vector response function, Y = f(u,t,theta;X), where X are system characteristics (like the volume of a tank, say, or inflow concentrations) potentially to be optimized to meet some constraints, and theta are parameters to be learned (like reaction rates, settling velocities) from observed concentration data via calibration (in my case, preferably Bayesian to support uncertainty quantification over the unknown parameters).

I can easily make an ODEProblem work to separate user-specified “inputs” (i.e., X values) and parameters in the derivative function, but this solution throws errors I don’t yet understand how to interpret when passing the same to Turing for parameter estimation. Any tips/pointers would be much appreciated - especially so if there is a settled way of handling this pattern, now. I am working on a simple example to see if the SciML ecosystem is right for my research focus (optimization of wasteload allocations in aquatic systems). So while this is currently just an ODEProblem, it will quickly become a PDE problem, as well.

9

Open a new issue. This is unrelated to above, all of which here are pretty automatic in 2021.

10

Thanks, I did so. For future readers interested in the question I posed here, see this post instead.