## Optimal Control Benchmark Problems

#### Exercise

Objective: Solve the dynamic optimization benchmark problems2 and more dynamic optimization benchmark problems. Complete the 9 exercises as shown in the Jupyter Notebook link below. For each problem, create a program to optimize and display the results. Estimated Time (each): 10-30 minutes

#### Example 1a

• Nonlinear, unconstrained, minimize final state

$$\min_{u(t)} x_2 \left( t_f \right)$$ $$\mathrm{subject \; to}$$ $$\frac{dx_1}{dt}=u$$ $$\frac{dx_2}{dt}=x_1^2 + u^2$$ $$x(0) = [1 \; 0]^T$$ $$t_f=1$$

#### Example 1b

• Nonlinear, unconstrained, minimize final state with terminal constraint

$$\min_{u(t)} x_2 \left( t_f \right)$$ $$\mathrm{subject \; to}$$ $$\frac{dx_1}{dt}=u$$ $$\frac{dx_2}{dt}=x_1^2 + u^2$$ $$x(0) = [1 \; 0]^T$$ $$x_1 \left( t_f \right)=1$$ $$t_f=1$$

#### Solutions to Benchmarks 1a and 1b

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

m = GEKKO()

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

# Variables
x1 = m.Var(value=1)
x2 = m.Var(value=0)
u = m.Var(value=-0.75)

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

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

# Objective Function
m.Obj(x2*final)

m.options.IMODE = 6
m.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()

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

m = GEKKO()

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

# Variables
x1 = m.Var(value=1)
x2 = m.Var(value=0)
u = m.Var(value=-0.48)

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

# Equations
m.Equation(x1.dt()==u)
m.Equation(x2.dt()==x1**2 + u**2)
m.Equation(final*(x1-1)==0)

# Objective Function
m.Obj(x2*final)

m.options.IMODE = 6
m.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()

#### Example 2

• Nonlinear, constrained, minimize final state

$$\min_{u(t)} x_4 \left( t_f \right)$$ $$\mathrm{subject \; to}$$ $$\frac{dx_1}{dt}=x_2$$ $$\frac{dx_2}{dt}=-x_3 \, u + 16 \, t - 8$$ $$\frac{dx_3}{dt}=u$$ $$\frac{dx_4}{dt}=x_1^2+x_2^2+0.005 \left(x_2 + 16 \, t -8 -0.1x_3\,u^2\right)^2$$ $$x(0) = [0 \; -1 \; -\sqrt{5} \; 0]^T$$ $$-4 \le u \le 10$$ $$t_f=1$$

#### Solution to Benchmark 2

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

m = GEKKO()

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

# Parameters
u = m.MV(value=9,lb=-4,ub=10)
u.STATUS = 1
u.DCOST = 0

# Variables
t = m.Var(value=0)
x1 = m.Var(value=0)
x2 = m.Var(value=-1)
x3 = m.Var(value=-np.sqrt(5))
x4 = m.Var(value=0)

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

# Equations
m.Equation(t.dt()==1)
m.Equation(x1.dt()==x2)
m.Equation(x2.dt()==-x3*u+16*t-8)
m.Equation(x3.dt()==u)
m.Equation(x4.dt()==x1**2+x2**2 \
+0.005*(x2+16*t-8-0.1*x3*(u**2))**2)

# Objective Function
m.Obj(x4*final)

m.options.IMODE = 6
m.options.NODES = 4
m.options.MV_TYPE = 1
m.options.SOLVER = 3
m.solve()

print(m.path)

print('Objective = min x4(tf): ' + str(x4[-1]))

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

#### Example 3

• Tubular reactor with parallel reaction

$$\max_{u(t)} x_2 \left( t_f \right)$$ $$\mathrm{subject \; to}$$ $$\frac{dx_1}{dt}=-\left(u+0.5u^2\right) x_1$$ $$\frac{dx_2}{dt}=u \, x_1$$ $$x(0) = [1 \; 0]^T$$ $$0 \le u \le 5$$ $$t_f=1$$

#### Solution to Benchmark 3

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

m = GEKKO()

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

# Parameters
u = m.MV(value=1,ub=5,lb=0)
u.STATUS = 1

# Variables
x1 = m.Var(value=1)
x2 = m.Var(value=0)

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

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

# Objective Function
m.Obj(-x2*final)

m.options.IMODE = 6
m.solve()

print('Objective: ' + str(x2[-1]))

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()

#### Example 4

• Batch reactor with consecutive reactions A->B->C

$$\max_{T(t)} x_2 \left( t_f \right)$$ $$\mathrm{subject \; to}$$ $$\frac{dx_1}{dt}=-k_1 \, x_1^2$$ $$\frac{dx_2}{dt}=k_1 \, x_1^2 - k_2 \, x_2$$ $$k_1 = 4000 \, \exp{\left(-\frac{2500}{T}\right)}$$ $$k_2 = 6.2e5 \, \exp{\left(-\frac{5000}{T}\right)}$$ $$x(0) = [1 \; 0]^T$$ $$298 \le T \le 398$$ $$t_f=1$$

#### Solution to Benchmark 4

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

m = GEKKO()

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

# Parameters
T = m.MV(value=362,ub=398,lb=298)
T.STATUS = 1
T.DCOST = 0

# Variables
x1 = m.Var(value=1)
x2 = m.Var(value=0)

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

# Intermediates
k1 = m.Intermediate(4000*m.exp(-2500/T))
k2 = m.Intermediate(6.2e5*m.exp(-5000/T))

# Equations
m.Equation(x1.dt()==-k1*x1**2)
m.Equation(x2.dt()==k1*x1**2 - k2*x2)

# Objective Function
m.Obj(-x2*final)

m.options.IMODE = 6
m.solve()

print('Objective: ' + str(x2[-1]))

plt.figure(1)

plt.subplot(2,1,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.ylabel('Value')
plt.legend(loc='best')

plt.subplot(2,1,2)
plt.plot(m.time,T.value,'r--',lw=2,label=r'$T$')
plt.legend(loc='best')
plt.xlabel('Time')
plt.ylabel('Value')

plt.show()

#### Example 5

• Catalytic reactor with A<->B->C

$$\max_{u(t)} \left(1 - x_1 \left( t_f \right) - x_2 \left( t_f \right) \right)$$ $$\mathrm{subject \; to}$$ $$\frac{dx_1}{dt}=u \left(10 \, x_2 - x_1 \right)$$ $$\frac{dx_2}{dt}=-u \left(10 \, x_2 - x_1 \right)-\left(1-u\right) x_2$$ $$x(0) = [1 \; 0]^T$$ $$0 \le u \le 1$$ $$t_f=12$$

#### Solution to Benchmark 5

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

m = GEKKO()

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

# Parameters
u = m.MV(value=1,ub=1,lb=0)
u.STATUS = 1
u.DCOST = 0

# Variables
x1 = m.Var(value=1)
x2 = m.Var(value=0)

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

# Equations
m.Equation(x1.dt()==u*(10*x2-x1))
m.Equation(x2.dt()==-u*(10*x2-x1)-(1-u)*x2)

# Objective Function
m.Obj(-final*(1-x1-x2))

m.options.IMODE = 6
m.solve()

print('Objective: ' + str(1-x1[-1]-x2[-1]))

plt.figure(1)

plt.subplot(2,1,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.ylabel('Value')
plt.legend(loc='best')

plt.subplot(2,1,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()

#### References

1. 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
2. M. Čižniar, M. Fikar, M.A. Latifi: A MATLAB Package for Dynamic Optimisation of Processes, 7th International Scientific – Technical Conference – Process Control 2006, June 13 – 16, 2006, Kouty nad Desnou, Czech Republic. Article