## 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

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-',LineWidth=2,label=r'$x_1$')
plt.plot(m.time,x2.value,'b-',LineWidth=2,label=r'$x_2$')
plt.plot(m.time,u.value,'r--',LineWidth=2,label=r'$u$')
plt.legend(loc='best')
plt.xlabel('Time')
plt.ylabel('Value')
plt.show()

# 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-',LineWidth=2,label=r'$x_1$')
plt.plot(m.time,x2.value,'b-',LineWidth=2,label=r'$x_2$')
plt.plot(m.time,u.value,'r--',LineWidth=2,label=r'$u$')
plt.legend(loc='best')
plt.xlabel('Time')
plt.ylabel('Value')
plt.show()

#### References

1. Beal, L., GEKKO for Dynamic Optimization, 2018, URL: https://gekko.readthedocs.io/en/latest/
2. 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