Integral Objective (Luus)

Objective: Set up and solve the Luus optimal control benchmark problem. Create a program to optimize and display the results. Estimated Time: 30 minutes

The Luus optimal control problem has an integral objective function. Integrals are a natural expression of many systems where the accumulation of a quantity is maximized or minimized.

$$\min_u \frac{1}{2} \int_0^2 x_1^2(t) \, dt$$ $$\mathrm{subject \; to}$$ $$\frac{dx_1}{dt}=u$$ $$x_1(0) = 1$$ $$-1 \le u(t) \le 1$$

Integral expressions are reformulated in differential and algebraic form by defining a new variable.

$$x_2 = \frac{1}{2} \int_0^2 x_1^2(t) \, dt$$

The new variable and integral are differentiated and included as an additional equation. The problem then becomes a minimization of the new variable `x_2` at the final time.

$$\min_u x_2\left(t_f\right)$$ $$\mathrm{subject \; to}$$ $$\frac{dx_1}{dt}=u$$ $$\frac{dx_2}{dt} = \frac{1}{2} x_1^2(t)$$ $$x_1(0) = 1$$ $$x_2(0) = 0$$ $$t_f = 2$$ $$-1 \le u(t) \le 1$$

The initial condition of the integral starts at zero and becomes the integral in the time range of 0 to 2. The value that is minimized is at the final point in the time horizon of the optimal control problem.

Solutions

Luus Problem in MATLAB and Python and Video

import numpy as np
import matplotlib.pyplot as plt
from gekko import GEKKO

m = GEKKO()

nt = 101
m.time = np.linspace(0,2,nt)

# Variables
x1 = m.Var(value=1)
x2 = m.Var(value=0)
u = m.Var(value=0,lb=-1,ub=1)

p = np.zeros(nt)
p[-1] = 1.0
final = m.Param(value=p)

# Equations
m.Equation(x1.dt()==u)
m.Equation(x2.dt()==0.5*x1**2)

# Objective Function
m.Obj(x2*final)

m.options.IMODE = 6
m.solve()
# m.solve(remote=False) # for local solve

plt.figure(1)
plt.plot(m.time,x1.value,'k-',lw=2,label=r'$x_1$')
plt.plot(m.time,x2.value,'b-',lw=2,label=r'$x_2$')
plt.plot(m.time,u.value,'r--',lw=2,label=r'$u$')
plt.legend(loc='best')
plt.xlabel('Time')
plt.ylabel('Value')
plt.show()

[$[Get Code]]

# this solution uses the m.integral function and produces the same answer
import numpy as np
import matplotlib.pyplot as plt
from gekko import GEKKO

m = GEKKO()

nt = 101
m.time = np.linspace(0,2,nt)

# Variables
x1 = m.Var(value=1)
u = m.Var(value=0,lb=-1,ub=1)

p = np.zeros(nt)
p[-1] = 1.0
final = m.Param(value=p)

# Equations
m.Equation(x1.dt()==u)
x2 = m.Intermediate(m.integral(0.5*x1**2))

# Objective Function
m.Obj(m.integral(0.5*x1**2)*final)

m.options.IMODE = 6
m.solve()
# m.solve(remote=False) # for local solve

plt.figure(1)
plt.plot(m.time,x1.value,'k-',lw=2,label=r'$x_1$')
plt.plot(m.time,x2.value,'b-',lw=2,label=r'$x_2$')
plt.plot(m.time,u.value,'r--',lw=2,label=r'$u$')
plt.legend(loc='best')
plt.xlabel('Time')
plt.ylabel('Value')
plt.show()

[$[Get Code]]

References

Beal, L., GEKKO for Dynamic Optimization, 2018, URL: https://gekko.readthedocs.io/en/latest/
Hedengren, J. D. and Asgharzadeh Shishavan, R., Powell, K.M., and Edgar, T.F., Nonlinear Modeling, Estimation and Predictive Control in APMonitor, Computers and Chemical Engineering, Volume 70, pg. 133–148, 2014. Article

Dynamic Optimization