Edit on GitHub

Linear Optimization Solution

Linear Optimization solution uses JuMP.jl to solve the Problems. Build the problem and call generateJuMPcodes, work done!

OptControl.generateJuMPcodes — Function
generateJuMPcodes(L, F, state, u, tspan, t0, tf) -> Any
generateJuMPcodes(L, F, state, u, tspan, t0, tf, Φ) -> Any
generateJuMPcodes(
    L,
    F,
    state,
    u,
    tspan,
    t0,
    tf,
    Φ,
    tf_constraint;
    state_ub,
    state_lb,
    u_ub,
    u_lb,
    N,
    discretization,
    model_name,
    writeFilePath,
    tolerance
) -> Any

generateJuMPcodes will generate the solution code that using JuMP.jl to solve problem

A general form of optimal control problem:

\[min (\Phi(\boldsymbol{x}(t_f),t_f)+\int_{t_0}^{t_f} L[\boldsymbol{x}(t),\boldsymbol{u}(t),t]dt) \\ s.t. \hspace{0.4cm} \dot{\boldsymbol{x}} =f[\boldsymbol{x}(t),\boldsymbol{u}(t),t]\]

args:

  • L: above L

\[L[\boldsymbol{x}(t),\boldsymbol{u}(t),t]\]

  • F: above f

\[f[\boldsymbol{x}(t),\boldsymbol{u}(t),t]\]

  • state: states(variable) in F
  • u: control variable u in F
  • tspan: time field
  • t0: initial value at tspan[1], length of t0 must be equal to length of state
  • tf: final(end) value at tspan[2], length of tf must be equal to length of state
  • Φ: above Φ, default: nothing

\[\Phi(\boldsymbol{x}(t_f),t_f)\]

  • tf_constraint: special requirements for end time, default: nothing. Example:

\[x_{1}+x_{2}=0\]

  • state_ub: state's upper limit, length of state_ub must be equal to length of state, default: nothing
  • state_lb: state's lower limit, length of state_lb must be equal to length of state, default: nothing
  • u_ub: u's upper limit, length of u_ub must be equal to length of u, default: nothing
  • u_lb: u's lower limit, length of u_lb must be equal to length of u, default: nothing
  • N: Number of discrete, default: 1000

\[dt = (endTime - startTime) / N\]

  • discretization: discretization methods, default: "trapezoidal"
  • model_name: name of JuMP model, default: "model"
  • writeFilePath: path of generated code, default: nothing
  • tolerance: acceptable_tol of Ipopt, default: 1.0e-6

Example:

using OptControl
@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)
source

Let's see some examples. If you are not familar with ModelingToolkit.jl, you can see the Symbolics.jl document (which ModelingToolkit.jl based on) and try some tests about symbolic computation

using OptControl,ModelingToolkit
@variables x[1:2] y[1:2]
print(scalarize(rand(1:10,2,2)*x+rand(1:10,2,2)*y))
Num[6x[1] + 6x[2] + 9y[1] + y[2], 3y[1] + 5x[2] + 8x[1] + 10y[2]]

PS: scalarize is from Symbolics.jl

If you need, pass a name to writeFilePath and do some changes in script.

Example 1: Fixed value

To solve:

\[min \int_{0}^{2} u^2dt \newline s.t. ~~~~~ \dot{\boldsymbol{x}} =\begin{bmatrix}0&1 \newline 0&0\end{bmatrix}\boldsymbol{x}+ \begin{bmatrix}0 \newline 1 \end{bmatrix}u \newline \boldsymbol{x}(0) = \begin{bmatrix} 1 \newline 1 \end{bmatrix}, \boldsymbol{x}(2)=\begin{bmatrix} 0 \newline 0 \end{bmatrix}\]

Just define variables and build functions. Call generateJuMPcodes and get the results.

The analytical solution of $x_1$ is

\[x_1(t) = 0.5*t^3-1.75*t^2+t+1\]

and we can campare the difference between them by using Mean Square Error(MSE).

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)
xs = collect(range(tspan[1], tspan[2], length=N))
an = @.(0.5 * xs^3 - 1.75 * xs^2 + xs + 1)
mean((an - sol[1][:, 1]).^2)
2.696022442477049e-6

Example 2: Free End

To solve:

\[min \int_{0}^{2} u^2dt \newline s.t. ~~~~~ \dot{\boldsymbol{x}} =\begin{bmatrix}0&1 \newline 0&0\end{bmatrix}\boldsymbol{x}+ \begin{bmatrix}0 \newline 1 \end{bmatrix}u \newline \boldsymbol{x}(0) = \begin{bmatrix} 1 \newline 1 \end{bmatrix}, \boldsymbol{x}(1)=\begin{bmatrix} 0 \newline free \end{bmatrix}\]

If variable is free, use nothing.

The analytical solution of $x_1$ is

\[x_1(t) = t^3-3.0*t^2+t+1\]

and get MSE.

using OptControl, Statistics, ModelingToolkit, 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, nothing]
tspan = (0.0, 1.0)
N = 100
sol = generateJuMPcodes(L, f, x, u, tspan, t0, tf; N=N)
xs = collect(range(tspan[1], tspan[2], length=N))
an = @.(xs^3 - 3.0 * xs^2 + xs + 1)
mean((an - sol[1][:, 1]).^2)
4.6689960309698826e-7

Example 3: End Constraint

To solve:

\[min \int_{0}^{2} u^2dt \newline s.t. ~~~~~ \dot{\boldsymbol{x}} =\begin{bmatrix}0&1 \newline 0&0\end{bmatrix}\boldsymbol{x}+ \begin{bmatrix}0 \newline 1 \end{bmatrix}u \newline \boldsymbol{x}(0) = \begin{bmatrix} 1 \newline 1 \end{bmatrix}, \boldsymbol{x}(1)=\begin{bmatrix} free \newline free \end{bmatrix} \\ x_1(1)+x_2(1)=0\]

Passing special end constraint to parameter tf_constraint.

The analytical solution of $x_1$ is

\[x_1(t) = -1/14*t^2*(t - 6)\]

and get MSE.

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 = [0.0, 0.0]
tf = [nothing, nothing]
tspan = (0.0, 1.0)
tf_con = x[1] + x[2] - 1.0
N = 100
sol = generateJuMPcodes(L, f, x, u, tspan, t0, tf, nothing, tf_con; N=N)
xs = collect(range(tspan[1], tspan[2], length=N))
an = @.(-1 / 14 * xs^2 * (xs - 6))
mean((an - sol[1][:, 1]).^2)
2.3177456677511903e-6

PS:x[1] + x[2] - 1.0 above is $x_1+x_2-1.0$. There nothing to do with independent variable $t$, because tf_constraint means constraint at the end time.

Example 4: Multiple End Constraint

To solve:

\[min \int_{0}^{2} u^2dt \newline s.t. ~~~~~ \dot{\boldsymbol{x}} =\begin{bmatrix}0&1 \newline 0&0\end{bmatrix}\boldsymbol{x}+ \begin{bmatrix}0 \newline 1 \end{bmatrix}u \newline \boldsymbol{x}(0) = \begin{bmatrix} 1 \newline 1 \end{bmatrix}, \boldsymbol{x}(2)=\begin{bmatrix} free \newline free \end{bmatrix} \\ x_1(2)+x_2(2)=0\\ x_1(2)-x_2(2)=0 \]

Actually, solution of those two constraint equations are $x_1=0,x_2=0$Essentially, it's the same as Example 1.

The analytical solution of $x_1$ is

\[x_1(t) = 0.5*t^3-1.75*t^2+t+1\]

and get MSE.

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 = [nothing, nothing]
tspan = (0.0, 2.0)
tf_con = [x[1] + x[2], x[1] - x[2]]
N = 100
sol = generateJuMPcodes(L, f, x, u, tspan, t0, tf, nothing, tf_con; N=N)
xs = collect(range(tspan[1], tspan[2], length=N))
an = @.(0.5 * xs^3 - 1.75 * xs^2 + xs + 1)
mean((an - sol[1][:, 1]).^2)
2.696022442477024e-6

Example 5: Multiple x and u

To solve:

\[min \int_{0}^{2} u^2dt \newline s.t. ~~~~~ \dot{\boldsymbol{x}} =\begin{bmatrix}0&1 \newline 0&0\end{bmatrix}\boldsymbol{x}+ \begin{bmatrix}1&0 \newline 0&1 \end{bmatrix}\boldsymbol{u} \newline \boldsymbol{x}(0) = \begin{bmatrix} 1 \newline 1 \end{bmatrix}, \boldsymbol{x}(2)=\begin{bmatrix} 0.0 \newline free \end{bmatrix} \]

The analytical solution of $u_1,u_2$ is

\[\begin{matrix} u_1=-\frac{9}{14}\\u_2=\frac{9}{14}*t-\frac{9}{7} \end{matrix} \]

and get MSE.

using OptControl, Statistics, ModelingToolkit
@variables t u[1:2] x[1:2]
f = [0 1; 0 0] * x + [1 0; 0 1] * u
L = 0.5 * (u[1]^2 + u[2]^2)
t0 = [1.0, 1.0]
tf = [0.0, nothing]
tspan = (0.0, 2.0)
N = 100
sol = generateJuMPcodes(L, f, x, u, tspan, t0, tf; N=N)
xs = collect(range(tspan[1], tspan[2], length=N))
an_u2 = @.(9 / 14 * xs - 9 / 7)
res1 = mean((-9 / 14 .- sol[2][:, 1]).^2)
res2 = mean((an_u2 - sol[2][:, 2]).^2)
(res1,res2)
(0.002166218961534756, 0.004150451971114561)

Example 6: Add variable limit

To solve:

\[min ~~x_2(1) \newline s.t. ~~~~~ \dot{\boldsymbol{x}} =\begin{bmatrix}-1&0 \newline 1&0\end{bmatrix}\boldsymbol{x}+ \begin{bmatrix}1 \newline 0 \end{bmatrix}u \newline \boldsymbol{x}(0) = \begin{bmatrix} 1 \newline 1 \end{bmatrix}, \boldsymbol{x}(1)=\begin{bmatrix} free \newline free \end{bmatrix} \\ -1.0 \leqslant\boldsymbol{u}\leqslant1.0\]

The analytical solution of $x_1$ is

\[x_1(t) = 2*e^{-t} - 1\]

and get MSE.

using OptControl, Statistics, ModelingToolkit
@variables t u x[1:2]
f = [-1 0; 1 0] * x + [1, 0] * u
L = 0
t0 = [1.0, 0.0]
tf = [nothing, nothing]
tspan = (0.0, 1.0)
N = 100
Φ = x[2]
uub = [1.0]
ulb = [-1.0]
sol = generateJuMPcodes(L, f, x, u, tspan, t0, tf, Φ, nothing;N=N, u_ub=uub, u_lb=ulb)
xs = collect(range(tspan[1], tspan[2], length=N))
an = @.(2 * exp(-xs) - 1)
mean((an - sol[1][:, 1]).^2)
3.244295042930682e-5

In this example, give L a value 0. Because L has nothing to do with optimization objective.