Nonlinear Optimal Control - orthogonal collocation support time points

[dynaTenCo]

I tried to solve a nonlinear optimal control problem in which model is defined as follows:

model = InfiniteModel(Ipopt.Optimizer)
@infinite_parameter(model, t in [0, 36], num_supports = 21, 
                           derivative_method = OrthogonalCollocation(3))

My intention is to discretize the time domain of 0-36 seconds into 5 finite elements of 7.2 s each, with 3 internal nodes in each element.

The optimization problem is solved successfully. However, when I access the support values used to transcribe u(t) along the domain of t, using supports.(u), the result is a uniformly spaced set of points 1.8s apart.
My understanding is that orthogonal collocation selects nodes which are typically not uniformly spaced, depending on the roots of the polynomial approximation using Lobatto quadrature.
Solution values of state variables and control variable seem to be OK.
What is the right way to access the support values for u(t) along the domain of t? Does InfiniteOpt interpolate and return the solution values at equally spaced points of time? Is there a mistake in the model definition (num_supports, collocation nodes etc?)

Full code is given below:

using InfiniteOpt, Ipopt
model = InfiniteModel(Ipopt.Optimizer)

@infinite_parameter(model, t in [0, 36], num_supports = 21, 
                           derivative_method = OrthogonalCollocation(3))

@variable(model, x1, Infinite(t), start = 0)
@variable(model, x2, Infinite(t), start = 0)
@variable(model, u, Infinite(t), start = 0)

@objective(model, Min, integral(γ1*(x1-2.5)^2+γ2*(x2-1)^2+γ3*(u-25)^2, t))

@constraint(model, x1(0) == x10)
@constraint(model, x2(0) == x20)
@constraint(model, u(0) == u0)

@constraint(model, c1, deriv(x1, t) == -k1*x1 - k3*x1^2 + (xF - x1)*u)
@constraint(model, c2, deriv(x2, t) == k1*x1 - k2*x2 - x2*u)

optimize!(model)

u_opt = value.(u)
u_ts = supports.(u)
value.(x1)
value.(x2)

Version info:
Julia 1.6.2
InfiniteOpt v0.5.0
Optim v1.6.0
Ipopt v0.7.0

Thanks in advance.

[pulsipher]

Hi there, great question! By default, we hide the internal supports (and associated variables/constraints) generated by collocation methods. However, these can be recovered using the label = All keyword argument as explained here and here. So in this case, modify the last few lines to:

u_opt = value.(u, label = All)
u_ts = supports.(u, label = All)
value.(x1, label = All)
value.(x2, label = All)

Under the hood, we employ a sophisticated support point management system to ensure, user points, generated points, derivatives, and measures (e.g., integrals) work together seamlessly. This system assigns at least one label to each support value so we can keep track of where they came from and/or how they were generated as explained here. The labels we use are detailed here. By default, only the supports with labels that inherit from PublicLabel are returned. Using All will make it such that all the supports are returned regardless of the labels that they have.

[pulsipher]

The code above isn't quite right for what you intend. When specifying the supports you only specify the ones that will comprise the finite elements (not the internal collocation nodes). Also, when using OrthogonalCollocation the number specified is inclusive of the finite element points since it uses Lobatto quadrature. Hence, you would need to specify OrthogonalCollocation(5) to have 3 internal nodes per finite element (as explained here). So for your example I would do:

@infinite_parameter(model, t in [0, 36], num_supports = 6, derivative_method = OrthogonalCollocation(5))

Note that using num_supports = 6 is equivalent to supports = collect(range(0, 36, length = 6)) in this case which is equivalent to supports = collect(0:7.2:36).

Again note, that currently we do not hold control variables constant over the internal nodes so you might get some unexpected results. However, I think I will have time to address #216 this month.

With regard to the prediction horizon vs control horizon, does the prediction horizon denote the time horizon of the optimal control subproblem solved at each MPC iteration and the control horizon refers to the overall duration of running the MPC (i.e., prediction horizon < control horizon)? Or are you referring to something else?

[dynaTenCo]

Thanks for the info regarding proper definition of supports.

Regarding prediction and control horizons:
In the MPC literature, prediction horizon, p, is the number of future finite time elements the MPC controller must evaluate by model prediction when optimizing the control variables at control interval k.
Control horizon, m, is the number of control variable moves to be optimized at control interval k. The control horizon falls between 1 and the prediction horizon p.
In general, the recommendation seems to be to keep m << p. I think this is mainly just to reduce the cost of computation, since in MPC only the first control move is applied and the rest discarded.
I believe in the case of InfiniteOpt, both prediction horizon and control horizon correspond to the defined time domain. It probably doesn't matter since the solution is pretty fast in Julia. It may become significant for large systems.

[pulsipher]

Thanks for the clarification. The MPC problems I had encountered typically use m = p following your notation. In any case, supporting this directly is probably beyond the scope of InfiniteOpt since it is not a control specific package (hence we don't distinguish between state and control variables). As such, running InfiniteOpt will not be most efficient way to conduct MPC and will likely be insufficient for real-time application with high control frequencies (e.g., autonomous vehicle control); in which case domain specific tools would be better suited (e.g., Gekko). InfiniteOpt is more intended as a convenient high-level sandbox for experimenting with general infinite-dimensional optimization problems.

With that said, one simple workaround would be to constrain the control variables to be constant after a certain time. For example,
with a control variable y(t):

@constraint(model, y(t_p) == y, DomainRestrictions(t => [t_p, t_f]))

where t_p is the final control horizon time and t_f is the final prediction horizon time.

[dynaTenCo]

I am not sure if there is any theoretical benefit in keeping the control horizon less than the prediction horizon. I think it became a practice to reduce the computational complexity. But given the computing power these days, it may not be a concern.
I have used Gekko for the last couple of years and its a fantastic package, with mature features. I was looking for pure Julia alternatives and found the SciML ecosystem (https://sciml.ai/), with DifferentialEquations.jl, ModelingToolkit.jl, DiffEqFlux.jl, etc. It seems to provide a lot more in terms of integrating differential equations with neural networks, which provides powerful modeling capabilities and learning dynamics from data (universal differential equations). However there seems to be some limitations in nonlinear optimal control at present. That's when I came across InfiniteOpt. As you mentioned it's a great sandbox to play around, and has helped greatly improve my understanding of discretization of continuous systems, orthogonal collocation, etc.
It would be great if InfiniteOpt would be composable with SciML ecosystem! Please explore this possibility, so that a model built using ModelingToolkit may be used in InfiniteOpt for optimal control.
Anyway, thanks for your time so far, your comments have been very valuable. I will continue to play around with InfiniteOpt.
Thanks.

[pulsipher]

I am happy to help. I have had a similar experience with the packages you described above. The SciML ecosystem is certainly very rich for dealing with systems of differential equations. It is definitely on my radar to make a bridge in accordance with #44. In fact, I am currently working on an extension package that will convert InfiniteModels over to ModelingToolkit and vice versa. This will facilitate simulation and initialization techniques. It will also allow us to take advantage of the preprocessing techniques (e.g., order reduction) that ModelingToolkit has. The main limitation is that ModelingToolkit has a pretty limited syntax for constrained optimization problems at the moment, so not all the features of an InfiniteModel can readily be transferred over. The other limitation of course is finding the necessary time :)

I appreciate your taking the time to ask questions here. This is helping figure out what we should include on the upcoming FAQ page (#217). This also is inspiring needed feature additions (e.g., #216). So thank you for using InfiniteOpt.jl and for the conversation!