Exercise
Objective: Solve the dynamic optimization benchmark problems 2 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
VIDEO
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 ( )
VIDEO
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
VIDEO
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 ( )
VIDEO
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
VIDEO
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 ( )
VIDEO
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
VIDEO
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 ( )
VIDEO
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
VIDEO
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 ( )
VIDEO
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