Optimal Control Benchmark Problems

Exercise

Objective: Solve the dynamic optimization benchmark problems² 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

DynOpt Benchmarks on GitHub or on Google Colab.

DynOpt Benchmark Solutions on GitHub or on Google Colab.

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

Dynamic Optimization Benchmark 1a and 1b in MATLAB and Python

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

[$[Get Code]]

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

[$[Get Code]]

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.0005 \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

Dynamic Optimization Benchmark 2 in MATLAB and Python

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.0005*(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()

[$[Get Code]]

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

Dynamic Optimization Benchmark 3 in MATLAB and Python

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

[$[Get Code]]

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

Dynamic Optimization Benchmark 4 in MATLAB and Python

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

[$[Get Code]]

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

Dynamic Optimization Benchmark 5 in MATLAB and Python

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

[$[Get Code]]

References

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

Dynamic Optimization