Silver bullet: An optimal control solver
Silver bullet is a direct transcription finite element method for the numerical solution of one-dimensional optimal control problems. Instead of collocation it uses L2 plenalties. While it is trivial to construct examples for which collocation methods diverge in general, the L2 penalty approach is proven to be convergent for the general one-dimensional optimal control problem with general differential-algebraic equality- and inequality-constraints. While the code is under progress, I present on this page some numerical results.

Origins of the method

In 2016/2017 I had an employment at the Mercedes Formula One team. All Formula One teams use numerical methods for optimal control to minimize lap time. I was working on a replacement for the current numerical optimization methods. Working on these methods during my employment raised my interest for research in new numerical methods for industrial optimal control problems.

Practical methods discretise the optimal control problem in a direct way into a large sparse nonlinear program. This in turn is solved using software for constrained nonlinear optimization. The process of transcription from an infinite-dimensional optimisation problem into a finite-dimensional one involves numerical approximation schemes for the underlying differential-algebraic constraints. During the time of my employment, there was no method known that had a theoretical proof that the numerical solution would converge to the analytic solution of the original problem. This is why after my employment and completion of my master's degree I performed research to finally develop a transcription method that is proven to be convergent.

The method is described here: . Dr Eric Kerrigan from Imperial College London and I wrote a paper that presents a very thorough proof of convergence for this method: . The proof works under very mild assumptions.

The idea of the method is as follows. As common, finite elements are used to parametrise the arcs of the solution. Since collocation can diverge, the path-constraints are treated not with collocation but by integrating their Lebesgue two-norm and adding it as a quadratic penalty to the cost-functional. Inequality constraints are treated with an integral of a logarithmic barrier function that is also added to the cost-functional. The L2 penalty yields symmetric stencils and is apparently stiffly stable. In result of the constraint treatments on the continuous level we obtain a finite-dimensional unconstrained penalty-barrier nonlinear programming problem. This is solved efficiently with an SQP-like numerical optimization method. The arising linear systems have narrow-banded matrices and can be solved fully in parallel, employing a divide-and-conquer technique.


John Betts provides a collection of optimal control problems. Below we solve a small selection of problems described therein, using our method. All problems were solved with isotropic mesh, finite elements of polynomial degree 4 to 10, and all in less than 5 minutes. The computations were performed without parallelism in Matlab 2016b on a Lenovo Yoga from 2015 with Intel i7 4800U processor and 8GB RAM. Since our method does not need to distinguish between states and controls, so do we not either, and call everything "species". The dimensions of the problems are between 1000 and 8000 degrees of freedom for the non-linear program. We plan for a parallel implementation in C++ for a compute cluster. In future we want to solve far bigger problems.

Aly-Chan problem
This problem has four smooth species. u is a fully singular arc that touches its box-constraints [-1,1] only at one time-instance. Numerically this problem is difficult to treat as usually the NLP solver does not know what use to make of the degrees of freedom in u and thus starts ringing wildly. In our case we find that the numerical solution is a smooth approximation to the analytic solution, which is u(t)=-sin(t).

Figure 1: Numerical solution to the Aly-Chan problem.

Brachistochrone problem
The problem consists of four smooth species. The problem is well-scaled. The solution can be rapidly approximated with a few elements of high order. A difficulty with this problem is that the end-time is the variable to minimize. After transforming the system into an interval of fixed length, the constraint equations are nonlinear of type res=x1 * x2 * cos(x3). Our NLP solver requires about 40 iterations, depending on mesh-size and elemental degree.

Figure 2: Numerical solution to the brachistochrone problem.

This is a badly scaled non-linear problem of five smooth species. It stems from a biological model for the computation of a treatment plan for HIV with chemotherapy.

Figure 3: Numerical solution to the chemotherapy problem.

Container crane problem
This is a problem with 8 smooth species. We have no deep understanding of this problem, yet, but were able to verify through mesh convergence that the presented arcs give an accurate solution.

Figure 4: Numerical solution to the container crane problem.

Drug treatment problem
This is a badly scaled problem of four smooth species that arises from a biological model for a drug treatment of HIV. While smooth, three of the four species have non-smooth first derivatives, making it difficult for the finite elements to accurately represent the exact solution.

Figure 5: Numerical solution to the drug treatment problem.

Free flying robot problem
This problem has twelve species with bang-bang controls. We used 80 finite elements of degree 8. Six out of the twelve states are non-smooth. While the arcs of the smooth species are well represented by the numerical solution, we find some non-optimal mesh-related artifacts nearby the discontinuities. This is because the high polynomial degree does not help approximating discontinuities well on a coarse mesh.

Figure 6: Numerical solution of the free flying robot problem.

The free flying robot problem searches for an economic control to reposition a robot in space. The robot has pneumatic propulsions. The amount of used pressure gas is minimized. The animation below visualizes the above control diagram.

Figure 7: Animation of the optimal dynamics for the free flying robot.

An example from raceline optimization

The following is a crude testing model for computing time-optimal trajectories through a track. The objective is to minimize the end-time. There are no-slip path-constraints, initial velocity constraints (start at rest), and start and finish coordinate constraints along the path of the track. The optimal arcs to be found involve raceline (plotted in blue), and arcs of accelerations and velocities.

Figure 8: Raceline along sample track.

Figure 9: Optimal velocities and accelerations along raceline.

The problem statement below consists of three parts: Equations (2)-(5) model the track parametrization. T is final time, s is track parametrization, n is track tangential vector (solved as part of the control problem), lambda is an auxiliary to express the tangent condition of n to the track. gamma is a given curve that returns track center line coordinates with respect to input s of track parametrization. Equations (6)-(8) model the dynamics. Equations (9)-(11) form the actual car model. Equations (12)-(14) form the point constraints. (12)-(13) make sure that the car drives from the start line to the finish line on the track. (14) gives the initial velocity vector. The model is simple and does not consider, e.g., angular orientation of the car. Constraints in (15) involve the track width B and the square of the maximum acceleration u^2_max .

Figure 10: Problem statement for raceline optimization with a simplified car model.

Contact via e-mail: