OptControl.jl: An Interpreter for Optimal Control Problem

Jingyi Yang¹, Yuebao Yang², Mingtao Li¹¹1Corresponding author
Xi’an Jiaotong University

Abstract

OptControl.jl²²2https://ai4energy.github.io/OptControl.jl/dev/(OptControl) implements that modeling optimal control problems with symbolic algebra system based on Julia language, and generates the corresponding numerical optimization codes to solve them with packages from Julia. OptControl does not define a data type, but generates a solution script by handling Julia strings and runs the script automatically. It also provides an interface to save script files. Meanwhile, OptControl supports component-based modeling, which makes it easy to build the optimal control problem of complex systems. All of OptControl’s dependency packages come from ecosystem of Julia.

1 Motivation

The essence of optimal control problem is an optimization problem. More accurately, it is a functional extreme value problem. From the perspective of implementation, the analytical solution of optimal control can be obtained only in certain cases, such as linear systems. However, practical problems usually are nonlinear systems or some much complex systems, and the analytical solutions are generally difficult to obtain. Therefore, the numerical solution that can calculate the result is a sharp weapon. Although the numerical solution will have some deviations in the results, the biased results can also give people some inspiration to understand the real world.

Generally, the numerical solution of optimal control problem can be transformed into numerical optimization problem. Numerical optimization problems can be solved by using JuMP.jl(JuMP). JuMP is a modeling language for mathematical optimization embedded in Julia[1].

Using JuMP to solve an optimal control problem can be divided into five steps

•

Step 1: Abstract problems from real world,
•

Step 2: Obtain state space model,
•

Step 3: Build numerical optimization model,
•

Step 4: Generate jump model,
•

Step 5: Solve the model.

In fact, JuMP finish Step 4, and the job of Step 5 is finished by a specific solver and JuMP provides an interface to Step 5. Therefore, it can also be considered that JuMP covers Step 4 and Step 5. Then the first three steps are left to the user, as shown in the Figure 1.

Refer to caption — Figure 1: Steps to solve an optimal control problem.

When dealing with the problem, the first step involves an abstract from real world to mathematical expression, which can only be completed by the advanced human brain. Can steps 2 and 3 be automatically finished? This is exactly what OptControl wants to accomplish. OptControl focuses on automated implementation.

OptControl focuses on automatic implementation, that is, how to automatically build the optimal problem, how to automatically build the jump optimization model and how to call solver to solve automatically. It is not involved that how to solve an optimal problem, how to build a symbolic system and so on . OptControl aims to integrate resources (various software packages in Julia ecosystem) and completing the above five steps automatically.

So OptControl can be regarded as an interpreter to implement the transformation from state space model to optimal control problem.It has three features:

1.

Accept state space model in form of symbolic system built from Symbolics.jl or ModelingToolkit.jl,
2.

Automatically generate the script solving optimization problem expressed in JuMP model and run it automatically,
3.

Provide an interface to save script files for users to modify freely

2 Framework of OptControl.jl

OptControl’s ability is gradually improved.

2.1 Build optimization model

The first thing to complete is the solution from Step 3 to Step 5, as shown in Figure 2. The function of this solution is:

•

generateJuMPcodes——dealing with linear systems,
•

generateNLJuMPcodes——dealing with nonlinear systems.

generateJuMPcodes and generateNLJuMPcodes accept the equation of state expressed in symbolic form with Symbolics. Symbolics.jl(Symbolics) is a symbolic algebra system with high performance and can be extended in user language[2].

The symbolic state space model can be transformed into Julia function which is used to obtain the discrete form of the model. And next taking the state in the discrete model as the variable of jump system. Finally, the JuMP model can be built and problem can be solved by solver.

2.2 Obtain state space model from ODESystem

Further, we want to automatically implement Steps 2 to 5 as shown in Figure 3. The solutions from step 2 to step 5 need to use the acausal component-based model in ModelingToolkit.jl(ModelingToolkit)[3]. ODESystem is a Julia data type in ModelingToolkit for modeling ordinary differential equation system(ODESystem).

The function of this solution is:

•

generateMTKcodes—— accept the ODESystem

The differential equation governing the system in odesystem is actually the state space form of the optimal control problem. The difference between them is that some of variables in the odesystem system are state-variables in the optimal control problem, while others are control-variables. In other words, the state space equation in the optimal control is essentially the same as the differential equation in ODESystem, but the difference is that some variables have special meanings in optimal control problem.

generateMTKcodes uses generate_function in ModelingToolkit to generate Julia function. Then the next step is the same as that in 2.1. Therefore, the solution ideas in 2.1 and 2.2 are the same, because they all generate Julia functions. The generated functions in 2.1 are from Symbolics while the functions in 2.2 are from ODESystem in ModelingToolkit, which is the main difference.

3 Math work in OptControl

3.1 Meaning of Simulation and Control

The state space equations in the optimal control problem are:

\boldsymbol{\dot{x}}(t)=\boldsymbol{A}\boldsymbol{x}(t)+\boldsymbol{B}\boldsymbol{u}(t)=f[\boldsymbol{x}(t),\boldsymbol{u}(t),t]

Where $\boldsymbol{x(t)}$ is the state-variable vector of the system, $\boldsymbol{u(t)}$ is the control-variable vector of the system. They are all functions of the independent variable $t$ , which means that they change with time. In fact, in the control problem, the coefficient matrix $\boldsymbol{A},\boldsymbol{B}$ can also change with time, expressed as $\boldsymbol{A(t)},\boldsymbol{B(t)}$ .

Omitting the $t$ , the above equation can be simplified as:

\dot{\boldsymbol{x}}=\boldsymbol{A}\boldsymbol{x}+\boldsymbol{B}\boldsymbol{u}=f(\boldsymbol{x},\boldsymbol{u},t)

If we consider mathematics without considering physical meaning of $\boldsymbol{x},\boldsymbol{u}$ , it is an ordinary differential equation problem with respect to time. If $\boldsymbol{u}$ is not artificially determined, it evolves itself in the system, which will be a problem of solving differential equations, as Equation 3 shows.

\dot{\boldsymbol{x}}=f(\boldsymbol{x},\boldsymbol{p},t)

Where $\boldsymbol{p}$ is the parameter vector of the system.

Solving the differential equation means dynamic simulation of the system in the real world. Therefore, the essence of control problem and dynamic simulation problem is the same. The difference is whether some variables in the problem can be changed by human intervention. It can also be argued that the dynamic simulation problem is that we want to see how the system evolves, and the control problem is that we want the system to evolve according to our expectations. Because we have expectations, we need human intervention in the system, which is the meaning of $\boldsymbol{u}$ . $\boldsymbol{u}$ is a mathematical form of intervention for system. If we set $\boldsymbol{u}$ without changing it, the problem is still a dynamic simulation problem, because the $\boldsymbol{u}$ does not interfere so that human expectations are not reflected in the system.

The same essence is why OptControl can accept the system from ModelingToolkit. ModelingToolkit was originally a tool for constructing dynamic simulation problems and ODESystem describes the simulation model of a dynamic system, which does not have an interface that can manually intervene in the system - the control variable $\boldsymbol{u}$ . We can build an ODESystem and observe how the system changes, without being able to control its evolution direction from beginning to end (in fact, occasional intervention can be achieved through the Callback function, but it is far from the "control"). But we can turn a simulation problem into a control problem with a few changes. Just add meanings to some variables in ODESystem, and the job is done, as the Equation 4 shows.

\dot{\boldsymbol{x}}=f(\boldsymbol{x},\boldsymbol{p},t)\Rightarrow\dot{\boldsymbol{x}}=f(\boldsymbol{x},\boldsymbol{u},t)

3.2 Mathematical Form of Optimization

In the previous section we explored the meaning of $\boldsymbol{u}$ . Then the remaining question is how to describe the optimization with math[4]. The whole optimization problem can be divided into two parts, the optimal in the control process and the optimal in the control final state.Equation 5 represents an objective for the optimal[5].

\min\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \Phi(\boldsymbol{x}(t_{f}),t_{f})+\int_{t_{0}}^{t_{f}}L[\boldsymbol{x}(t),\boldsymbol{u}(t),t]dt

Where $\Phi(\boldsymbol{x}(t_{f}),t_{f})$ expresses an expectation of the final state. $\int_{t_{0}}^{t_{f}}L[\boldsymbol{x}(t),\boldsymbol{u}(t),t]dt$ expresses an expectation in the control process. The combination of Equation 2 and Equation 5 becomes the general form of the governing equation for the optimal control problem.

\begin{matrix}\min\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \Phi(\boldsymbol{x}(t_{f}),t_{f})+\int_{t_{0}}^{t_{f}}L[\boldsymbol{x}(t),\boldsymbol{u}(t),t]dt\\ s.t.\hskip 28.45274pt\dot{\boldsymbol{x}}=f[\boldsymbol{x}(t),\boldsymbol{u}(t),t]\end{matrix}

3.3 Numerical optimization model

Equation 6 is in continuous form and needs to be discretized if numerical optimization methods are used. The discrete method uses Euler’s method[6], then there are:

\begin{matrix}\min\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \Phi(\boldsymbol{x}(t_{f}),t_{f})+\sum_{i=1}^{n}L(\boldsymbol{x}_{i},\boldsymbol{u}_{i},t_{i})\\ s.t.\hskip 11.38092pt\boldsymbol{x}_{i+1}=\boldsymbol{x}_{i}+f(\boldsymbol{x}_{i},\boldsymbol{u}_{i},t_{i})*dt\end{matrix}

If the backward Euler’s rule is used:

\begin{matrix}\min\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \Phi(\boldsymbol{x}(t_{f}),t_{f})+\sum_{i=1}^{n}L(\boldsymbol{x}_{i},\boldsymbol{u}_{i},t_{i})\\ s.t.\hskip 11.38092pt\boldsymbol{x}_{i+1}=\boldsymbol{x}_{i}+f(\boldsymbol{x}_{i+1},\boldsymbol{u}_{i+1},t_{i+1})*dt\end{matrix}

In addition, there are many discrete methods, such as trapezoidal method, Adams method[7] and so on. Once Equation 6 or Equation 7 is obtained, the next step is to use JuMP to build the corresponding JuMP model, which is the form the solver can solve. These steps are done automatically by OptControl.

4 Case Study

4.1 Case1: Optimal Control of Linear System

To solve the following linear optimal control problem:

\begin{matrix}min\int_{0}^{2}u^{2}dt\\ s.t.\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \dot{\boldsymbol{x}}=\begin{bmatrix}0&1\\ 0&0\end{bmatrix}\boldsymbol{x}+\begin{bmatrix}0\\ 1\end{bmatrix}u\\ \boldsymbol{x}(0)=\begin{bmatrix}1\\ 1\end{bmatrix},\boldsymbol{x}(2)=\begin{bmatrix}0\\ 0\end{bmatrix}\end{matrix}

The steps to solve this problem using OptControl are:

1.

Use ModelingToolkit or Symbolics to model the system,
2.

Set states such as initial state and final state,
3.

Call generateJuMPcodes to solve.

⬇

using OptControl, Statistics, ModelingToolkit

@variables t u x[1:2]

f = [0 1; 0 0] * x + [0, 1] * u

L = 0.5 * u^2

t0 = [1.0, 1.0]

tf = [0.0, 0.0]

tspan = (0.0, 2.0)

N = 100

sol = generateJuMPcodes(L, f, x, u, tspan, t0, tf; N=N)

In this optimal problem, the analytical solution of $x_{1}$ is

x_{1}(t)=0.5*t^{3}-1.75*t^{2}+t+1

Comparing the analytical solution and the numerical solution, it can be found the mean square error is 2.696E-6.

Even if the result under this error cannot be used, it is of great reference significance and can give people inspiration. Even if this result is not adopted, it can give people inspiration to solve problem.

4.2 Case2: Optimal Control of Nonlinear system

To solve the following nonlinear optimal control problem

\begin{array}[]{c}min\int_{0}^{2}u^{2}dt\\ s.t.\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \dot{\boldsymbol{x}}=\begin{bmatrix}exp&cos\\ sin&1\end{bmatrix}\boldsymbol{x}+\begin{bmatrix}0\\ 1\end{bmatrix}u\\ \boldsymbol{x}(0)=\begin{bmatrix}1\\ 1\end{bmatrix},\boldsymbol{x}(1)=\begin{bmatrix}0\\ 0\end{bmatrix}\end{array}

Use ModelingToolkit or Symbolics to define symbolic variables, giving initial and final states. Call generateNLJuMPcodes to get the result.

⬇

using OptControl, ModelingToolkit, Test

@variables t u x[1:2]

f = [exp(x[1]) + cos(x[2]), sin(x[1]) + x[2]] + [1, 0] * u

L = u^2

t0 = [1.0, 1.0]

tf = [0.0, 0.0]

tspan = (0.0, 2.0)

N = 100

sol = generateNLJuMPcodes(L, f, x, u, tspan, t0, tf; N=N)

4.3 Case3: Optimal Control of RC Circuit System

This is a simple circuit system. The supply voltage is 1 volt, the resistance is 1 ohm, and the capacitance is 1 farad, as Figure 4 shows.

The optimal control problem we constructed is how to change the voltage to make the capacitor voltage change from 1V to 3V within 1s, while satisfying the goal of the voltage as low as possible in the whole process. This has a physical meaning and can illustrate how optimal control can be applied in component-based systems.

The components of the code are taken from the documentation of ModelingToolkit³³3https://mtk.sciml.ai/dev/tutorials/spring_mass/. If it is a simulation problem, it should be solved by calling DifferentialEquations.jl[8] after the ODESystem is simplified. Now it is an optimal control problem, so we set the optimization objective, define the relevant parameters, and call generateMTKcodes to solve.

⬇

using OptControl, ModelingToolkit, Test

# Components ......

# Define and Simplify System

R = 1.0

C = 1.0

V = 1.0

@named resistor = Resistor(R=R)

@named capacitor = Capacitor(C=C)

@named source = ConstantVoltage(V=V)

@named ground = Ground()

rc_eqs = [

connect(source.p, resistor.p)

connect(resistor.n, capacitor.p)

connect(capacitor.n, source.n)

connect(capacitor.n, ground.g)

]

@named _rc_model = ODESystem(rc_eqs, t)

@named rc_model = compose(_rc_model,

[resistor, capacitor, source, ground])

sys = structural_simplify(rc_model)

# Build Optimal Control Problem and Solve

L = 0.5 * (source.V^2)

t0 = [1.0]

tf = [3.0]

tspan = (0.0, 1.0)

N = 100

sol = OptControl.generateMTKcodes(L, sys, states(sys), [source.V], tspan, t0, tf;N=N)

5 Conclusion

OptControl implements the automatic construction from the state space equations to the optimal control problem and the automatic construction from the ordinary differential equation system of the ModelingToolkit to the system optimal control problem. The core of the problem lies in the choice of design control variables $\boldsymbol{u}$ . When discretizing, OptControl provides an interface to select the discretization methods. In future work, more discrete methods will be added.

Another important feature of OptControl is that it does not solve the problem directly, but generates a solution script and runs it. This means that OptControl is like a conductor, it decomposes the problem, and then calls the package in the Julia ecosystem to solve the problem. OptControl provides an interface to get the script file, which means that when the script does not meet your needs, you can modify it directly. You can add whatever code you need to the generated script. If you are not familiar with JuMP’s modeling language, then you can learn some advanced usage of JuMP through the generated scripts. If you also want to choose some different solvers, modify the script.

In the future, OptControl may provide more interfaces. But it will not change the role as a conductor. OptControl will always focus on automatically generating solutions to optimal control problems, rather than developing a modeling language like ModelingToolkit and JuMP, nor developing solving algorithms like the solvers called by JuMP. The original intention of OptControl is to break through barriers, integrate tools, and solve optimal control problems conveniently and quickly.

References

[1] Iain Dunning, Joey Huchette, and Miles Lubin. Jump: A modeling language for mathematical optimization. SIAM Review, 59(2):295--320, 2017.
[2] Shashi Gowda, Yingbo Ma, Alessandro Cheli, Maja Gwozdz, Viral B Shah, Alan Edelman, and Christopher Rackauckas. High-performance symbolic-numerics via multiple dispatch. arXiv preprint arXiv:2105.03949, 2021.
[3] Yingbo Ma, Shashi Gowda, Ranjan Anantharaman, Chris Laughman, Viral Shah, and Chris Rackauckas. Modelingtoolkit: A composable graph transformation system for equation-based modeling. arXiv preprint arXiv:2103.05244, 2021.
[4] Nocedal Jorge and J Wright Stephen. Numerical optimization, 2006.
[5] Jack Macki and Aaron Strauss. Introduction to optimal control theory. Springer Science & Business Media, 2012.
[6] Richard L Burden, J Douglas Faires, and Annette M Burden. Numerical analysis. Cengage learning, 2015.
[7] Kai Diethelm, Neville J Ford, and Alan D Freed. Detailed error analysis for a fractional adams method. Numerical algorithms, 36(1):31--52, 2004.
[8] Christopher Rackauckas and Qing Nie. Differentialequations.jl--a performant and feature-rich ecosystem for solving differential equations in julia. Journal of Open Research Software, 5(1), 2017.