mintOC - User contributions [en]

Time optimal car problem (GEKKO)

2022-08-28T13:23:10Z

JohnHedengren: Update with latest Gekko v1.0.5

This page contains a solution of the [[Time optimal car problem]] in [https://gekko.readthedocs.io/en/latest/ GEKKO] Python format. The GEKKO package is available with '''pip install gekko'''. The Python code uses orthogonal collocation and a simultaneous optimization method. The end point constraints are imposed as hard constraints.

<source lang="Python">
from gekko import GEKKO
import matplotlib.pyplot as plt
import numpy as np

# set up the gekko model
m = GEKKO()

# set up the time (minimize the time with time scaling)
m.time = np.linspace(0, 1, 100)

# set up the variables
Z1 = m.Var(value=0, ub=330, lb=0)
Z2 = m.Var(value=0, ub=33, lb=0)
m.fix_final(Z2, 0)
m.fix_final(Z1, 300)

# set up the value we modify over the horizon
tf = m.FV(value=500, lb=0.1)
tf.STATUS = 1

# set up the MV
u = m.MV(integer=True, lb=-2, ub=1)
u.STATUS = 1

# set up the equations
m.Equation(Z1.dt() / tf == Z2)
m.Equation(Z2.dt() / tf == u)

# set the objective
m.Obj(tf)

# set up the options
m.options.IMODE = 6
m.options.SOLVER = 1

# solve
m.solve(disp=False)

# print the time
print("Total time taken: " + str(tf.NEWVAL))

# plot the results
plt.figure()
plt.subplot(211)
plt.plot(np.linspace(0,1,100)*tf.NEWVAL, Z1, label=r'$Z_1$')
plt.plot(np.linspace(0,1,100)*tf.NEWVAL, Z2, label=r'$Z_2$')
plt.ylabel('Z')
plt.legend()
plt.subplot(212)
plt.plot(np.linspace(0,1,100)*tf.NEWVAL, u, label=r'$u$')
plt.ylabel('u')
plt.xlabel('Time')
plt.legend()
plt.show()
</source>

== Results with APOPT (MINLP) ==

An MINLP solution is calculated with APOPT with an objective function value of <math>30.3077</math>. Thanks to [https://www.linkedin.com/in/tjasperson/ Tanner Jasperson] for providing a solution.

[[File:Time_optimal_car_GEKKO.png]]

[[Category:Gekko]]

File:Batch reactor GEKKO.png

2019-04-11T13:49:25Z

JohnHedengren:

Batch reactor (GEKKO)

2019-04-11T13:49:07Z

JohnHedengren: Add gekko solution for batch reactor.

This page contains a solution of the [[Batch reactor]] problem in [https://gekko.readthedocs.io/en/latest/ GEKKO] Python format. The model in Python code for a fixed control discretization grid uses orthogonal collocation and a simultaneous optimization method. The GEKKO package is available with '''pip install gekko'''.

<source lang="Python">
import numpy as np
import matplotlib.pyplot as plt
from gekko import GEKKO

m = GEKKO()

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

# Variables
x1 = m.Var(value=1)
x2 = m.Var(value=0)
T = m.MV(value=398,lb=298,ub=398)
T.STATUS = 1
T.DCOST = 1e-6
k1 = m.Intermediate(4000*m.exp(-2500/T))
k2 = m.Intermediate(6.2e5*m.exp(-5000/T))

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

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

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

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

print('Optimal x2(tf): ' + str(x2.value[-1]))

plt.figure(1)
plt.subplot(2,1,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[-1],x2.value[-1],'o',\
color='orange',MarkerSize=5,\
label=r'$x_2\left(t_f\right)$')
plt.legend()
plt.ylabel('Mole Fraction')
plt.subplot(2,1,2)
plt.plot(m.time,T.value,'r--',\
LineWidth=2,label='Temperature')
plt.ylabel(r'$T/;(^oC)$')
plt.legend(loc='best')
plt.xlabel('Time')
plt.ylabel('Value')
plt.show()
</source>

A solution is calculated with an objective function value of x2(tf) = 0.6108.

[[File:Batch_reactor_GEKKO.png]]

[[Category:Gekko]]

Batch reactor

2019-04-11T13:42:22Z

JohnHedengren: <math> </math> rendering appears to be broken

{{Dimensions
|nd = 1
|nx = 2
|nu = 1
|nc = 2
|nre = 2
}}

This batch reactor problem describes the consecutive reaction of some substance A via substance B into a desired product C.

The system is interacted with via the control function <math> T(t) </math> which stands for the temperature.
The goal is to produce as much of substance B (which can then be converted into product C) as possible within the time limit.

== Mathematical formulation ==

The optimal control problem is given by


<math>
\begin{array}{llcl}
\displaystyle \max_{x, u} & x_2(t_f) \\[1.5ex]
\mbox{s.t.} & \dot{x}_1 & = & -k_1 x_1^2.\\
& \dot{x}_2 & = & k_1 x_1^2 - k_2 x_2,\\
& k_1 & = & 4000 \; e^{(-2500/T(t))}, \\
& k_2 & = & 620000 \; e^{(-5000/T(t))}, \\[1.5ex]
& x(0) &=& (1, 0)^T, \\
& T(t) &\in& [298, 398].
\end{array}
</math>


<math> x_1(t) </math> and <math> x_2(t) </math> represent the concentrations of A and B at timepoint <math> t </math> respectively. The control function <math> T(t) </math> represents the temperature.

== Parameters ==

The starting time and end time are given by <math> [t_0, t_f] = [0, 1] </math>.

== Reference Solutions ==

This solution was computed using JuMP with a collocation method and 300 discretization points. The differential equations were solved using the explicit Euler Method. The source code can be found at [[Batch reactor (JuMP)]].

The optimal objective value of the problem is x2(tf) = 0.611715.

<gallery caption="Reference solution plots" widths="180px" heights="140px" perrow="2">
File:batch_reactor_states_plot.png| Optimal states <math> x_1(t)</math> and <math>x_2(t)</math>.
File:batch_reactor_control_plot.png| Optimal control <math>u(t)</math>.
</gallery>

== Source Code ==

Model descriptions are available in
* [[:Category:Gekko | GEKKO Python code]] at [[Batch reactor (GEKKO)]]
* [[:Category: Julia/JuMP | JuMP code]] at [[Batch reactor (JuMP)]]

== References ==
The problem can be found in [http://www.tomopt.com/docs/TOMLAB_PROPT.pdf the Tomlab PROPT guide] or in [http://www.kirp.chtf.stuba.sk/moodle/mod/resource/view.php?id=5464 the Dynopt guide].


[[Category:MIOCP]]
[[Category:ODE model]]
[[Category:Chemical engineering]]

Batch reactor

2019-04-11T13:41:03Z

JohnHedengren: Add link to GEKKO solution

{{Dimensions
|nd = 1
|nx = 2
|nu = 1
|nc = 2
|nre = 2
}}

This batch reactor problem describes the consecutive reaction of some substance A via substance B into a desired product C.

The system is interacted with via the control function <math> T(t) </math> which stands for the temperature.
The goal is to produce as much of substance B (which can then be converted into product C) as possible within the time limit.

== Mathematical formulation ==

The optimal control problem is given by


<math>
\begin{array}{llcl}
\displaystyle \max_{x, u} & x_2(t_f) \\[1.5ex]
\mbox{s.t.} & \dot{x}_1 & = & -k_1 x_1^2.\\
& \dot{x}_2 & = & k_1 x_1^2 - k_2 x_2,\\
& k_1 & = & 4000 \; e^{(-2500/T(t))}, \\
& k_2 & = & 620000 \; e^{(-5000/T(t))}, \\[1.5ex]
& x(0) &=& (1, 0)^T, \\
& T(t) &\in& [298, 398].
\end{array}
</math>


<math> x_1(t) </math> and <math> x_2(t) </math> represent the concentrations of A and B at timepoint <math> t </math> respectively. The control function <math> T(t) </math> represents the temperature.

== Parameters ==

The starting time and end time are given by <math> [t_0, t_f] = [0, 1] </math>.

== Reference Solutions ==

This solution was computed using JuMP with a collocation method and 300 discretization points. The differential equations were solved using the explicit Euler Method. The source code can be found at [[Batch reactor (JuMP)]].

The optimal objective value of the problem is <math> x_2(t_f) = 0.611715 </math>.

<gallery caption="Reference solution plots" widths="180px" heights="140px" perrow="2">
File:batch_reactor_states_plot.png| Optimal states <math> x_1(t)</math> and <math>x_2(t)</math>.
File:batch_reactor_control_plot.png| Optimal control <math>u(t)</math>.
</gallery>

== Source Code ==

Model descriptions are available in
* [[:Category:Gekko | GEKKO Python code]] at [[Batch reactor (GEKKO)]]
* [[:Category: Julia/JuMP | JuMP code]] at [[Batch reactor (JuMP)]]

== References ==
The problem can be found in [http://www.tomopt.com/docs/TOMLAB_PROPT.pdf the Tomlab PROPT guide] or in [http://www.kirp.chtf.stuba.sk/moodle/mod/resource/view.php?id=5464 the Dynopt guide].


[[Category:MIOCP]]
[[Category:ODE model]]
[[Category:Chemical engineering]]

Oil shale pyrolysis (GEKKO)

2019-03-15T02:06:33Z

JohnHedengren:

This page contains a solution of the [[Oil Shale Pyrolysis]] problem in [https://gekko.readthedocs.io/en/latest/ GEKKO] Python format. The model in Python code for a fixed control discretization grid uses orthogonal collocation and a simultaneous optimization method. The GEKKO package is available with '''pip install gekko'''.

<source lang="Python">
from gekko import GEKKO
import numpy as np
import matplotlib.pyplot as plt

m=GEKKO()

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

# Constants
R = m.Const(value=1.9858775e-3) # Gas constant [kcal/mol*K]

# Final Time = FV
tf = m.FV(value=1,lb=0.1,ub=20)
tf.STATUS = 1

# MV
T = m.MV(value=425, ub=475, lb=425) # degC
T.STATUS = 1
T.DCOST = 0

# Variables
x1 = m.Var(value = 1) # kerogen
x2 = m.Var(value = 0) # pyrolytic bitumen
x3 = m.Var(value = 0) # oil and gas
x4 = m.Var(value = 0) # carbon residue

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

# Parameters
# Frequency factor
a1 = m.Param(value=np.exp(8.86))
a2 = m.Param(value=np.exp(24.25))
a3 = m.Param(value=np.exp(23.67))
a4 = m.Param(value=np.exp(18.75))
a5 = m.Param(value=np.exp(20.7))
# Activation envergy [Kcal/g-mole]
b1 = m.Param(value=20.3)
b2 = m.Param(value=37.4)
b3 = m.Param(value=33.8)
b4 = m.Param(value=28.2)
b5 = m.Param(value=31.0)

# Intermediates
k1 = m.Intermediate(a1 * m.exp(-b1/(R*(T+273.15))))
k2 = m.Intermediate(a2 * m.exp(-b2/(R*(T+273.15))))
k3 = m.Intermediate(a3 * m.exp(-b3/(R*(T+273.15))))
k4 = m.Intermediate(a4 * m.exp(-b4/(R*(T+273.15))))
k5 = m.Intermediate(a5 * m.exp(-b5/(R*(T+273.15))))

# Equations
m.Equation(x1.dt()/tf == -k1*x1 - (k3 + k4 + k5) * x1*x2)
m.Equation(x2.dt()/tf == k1*x1 - k2*x2 + k3*x1*x2)
m.Equation(x3.dt()/tf == k2*x2 + k4*x1*x2)
m.Equation(x4.dt()/tf == k5*x1*x2)

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

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

plt.figure(1)

tm = m.time * tf.value[0]

plt.subplot(3,1,1)
plt.plot(tm,x1.value,'k:',LineWidth=2,label=r'$x_1$=Kerogen')
plt.plot(tm,x2.value,'b-',LineWidth=2,label=r'$x_2$=Pyrolytic Bitumen')
plt.ylabel('Fraction')
plt.legend(loc='best')

plt.subplot(3,1,2)
plt.plot(tm,x3.value,'g--',LineWidth=2,label=r'$x_3$=Oil & Gas')
plt.plot(tm,x4.value,'r-.',LineWidth=2,label=r'$x_4$=Carbon Residue')
plt.ylabel('Fraction')
plt.legend(loc='best')

plt.subplot(3,1,3)
plt.plot(tm,T.value,'r-',LineWidth=2,label='Temperature')
plt.legend(loc='best')
plt.xlabel('Time')
plt.ylabel(r'Temp ($^o$C)')

plt.show()
</source>

A solution is calculated with an objective function value of <math>x_2(t_f) = 0.35062</math>.

[[File:Oil_Shale_Pyrolysis_GEKKO.png]]

Thanks to [https://www.linkedin.com/in/junho-park-90b86967/ Junho Park] for providing the solution.

[[Category:Gekko]]

Oil shale pyrolysis (GEKKO)

2019-03-15T02:04:03Z

JohnHedengren:

2019-03-14T17:12:04Z

JohnHedengren:

Double Tank (GEKKO)

2019-03-14T17:11:36Z

JohnHedengren: Create results page with Python GEKKO

This page contains a solution of the [[Double Tank]] in [https://gekko.readthedocs.io/en/latest/ GEKKO] Python format. The GEKKO package is available with '''pip install gekko'''. The Python code uses orthogonal collocation and a simultaneous optimization method. The integral is converted to a differential equation through differentiation and definition of a new variable x3.

<source lang="Python">
import numpy as np
import matplotlib.pyplot as plt
from gekko import GEKKO

m = GEKKO() # create GEKKO model

# Add 0.01 as first step
# 0,0.01,0.1,0.2,0.3,...9.9,10.0)
m.time = np.insert(np.linspace(0,10,201),1,0.01)

# change solver options
m.solver_options = ['minlp_gap_tol 0.001',\
'minlp_maximum_iterations 10000',\
'minlp_max_iter_with_int_sol 100',\
'minlp_branch_method 1',\
'minlp_integer_tol 0.001',\
'minlp_integer_leaves 0',\
'minlp_maximum_iterations 200']

k1 = 2
k2 = 3

last = m.Param(np.zeros(202))
last.value[-1] = 1

sigma=m.MV(value=1,lb=1,ub=2,integer=True)
x1 = m.Var(value=2)
x2 = m.Var(value=2)
x3 = m.Var(value=0)
sigma.STATUS = 1

m.Obj(last*x3)

m.Equations([x1.dt() == sigma - m.sqrt(x1),\
x2.dt() == m.sqrt(x1) - m.sqrt(x2),\
x3.dt() == k1*(x2-k2)**2])

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

plt.figure(1)
plt.step(m.time,sigma.value,'r-',label=r'$\sigma$ (1/2)')
plt.plot(m.time,x1.value,'k-',label=r'$x_1$')
plt.plot(m.time,x2.value,'g-',label=r'$x_2$')
plt.xlabel('Time')
plt.ylabel('Variables')
plt.legend(loc='best')
plt.show()
</source>

== Results with APOPT (MINLP) ==

An MINLP solution is calculated with APOPT with an objective function value of <math>4.767757</math>.

[[File:Double_Tank_GEKKO.png]]

[[Category:Gekko]]

Double Tank

Lotka Volterra fishing problem (GEKKO)

2019-03-13T19:17:01Z

JohnHedengren:

This page contains a solution of the MIOCP [[Lotka Volterra fishing problem]] in [https://gekko.readthedocs.io/en/latest/ GEKKO] Python format. The model in Python code for a fixed control discretization grid using orthogonal collocation and a simultaneous optimization method. The GEKKO package is available with '''pip install gekko'''.

<source lang="Python">
import numpy as np
import matplotlib.pyplot as plt
from gekko import GEKKO

m = GEKKO() # create GEKKO model

# Add 0.01 as first step
# 0,0.01,0.1,0.2,0.3,...11.9,12.0)
m.time = np.insert(np.linspace(0,12,121),1,0.01)

# change solver options
m.solver_options = ['minlp_gap_tol 0.001',\
'minlp_maximum_iterations 10000',\
'minlp_max_iter_with_int_sol 100',\
'minlp_branch_method 1',\
'minlp_integer_tol 0.001',\
'minlp_integer_leaves 0',\
'minlp_maximum_iterations 200']

c0 = 0.4
c1 = 0.2

last = m.Param(np.zeros(122))
last.value[-1] = 1

x0 = m.Var(value=0.5,lb=0)
x1 = m.Var(value=0.7,lb=0)
x2 = m.Var(value=0.0,lb=0)
w = m.MV(value=0,lb=0,ub=1,integer=True)
w.STATUS = 1

m.Obj(last*x2)

m.Equations([x0.dt() == x0 - x0*x1 - c0*x0*w,\
x1.dt() == - x1 + x0*x1 - c1*x1*w,\
x2.dt() == (x0-1)**2 + (x1-1)**2])

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

plt.figure(1)
plt.step(m.time,w.value,'r-',label='w (0/1)')
plt.plot(m.time,x0.value,'b-',label=r'$x_0$')
plt.plot(m.time,x1.value,'k-',label=r'$x_1$')
plt.plot(m.time,x2.value,'g-',label=r'$x_2$')
plt.xlabel('Time')
plt.ylabel('Variables')
plt.legend(loc='best')
plt.show()
</source>

An MINLP solution is calculated with [https://apopt.com APOPT] with an objective function value of <math>x_2(t_f) = 1.349497</math>.

[[File:Volterra_fishing_GEKKO.png]]

[[Category:Gekko]]

2017-11-20T22:13:04Z

JohnHedengren:

Lotka Volterra fishing problem (APMonitor)

2017-11-20T22:12:26Z

JohnHedengren: Results section

This page contains a solution of the MIOCP [[Lotka Volterra fishing problem]] in [http://www.apmonitor.com APMonitor] Python format. A MATLAB version is also available from the [http://apmonitor.com/do/index.php/Main/DiscreteVariables Dynamic Optimization Course] as [http://apmonitor.com/do/uploads/Main/lotka_volterra_fishing.zip Example 3 (lotka_volterra_fishing.zip)].

=== APMonitor ===

The model in Python code for a fixed control discretization grid using orthogonal collocation and a simultaneous optimization method. The APMonitor package is available with '''pip install APMonitor''' or from the [https://github.com/APMonitor/apm_python/blob/master/apm.py APMonitor Python Github repository].

<source lang="Python">
import numpy as np
import matplotlib.pyplot as plt

# retrieve apm.py from
# https://raw.githubusercontent.com/APMonitor/apm_python/master/apm.py
# or
# http://apmonitor.com/wiki/index.php/Main/PythonApp
# from apm import *

# pip install with 'pip install APMonitor'
from APMonitor.apm import *

# local APMonitor servers are available for Windows or Linux
# http://apmonitor.com/wiki/index.php/Main/APMonitorServer
# with clients in Python, MATLAB, and Julia

# write model
model = '''
! apopt MINLP solver options (see apopt.com)
File apopt.opt
minlp_maximum_iterations 1000 ! minlp iterations
minlp_max_iter_with_int_sol 50 ! minlp iterations if integer solution is found
minlp_as_nlp 0 ! treat minlp as nlp
nlp_maximum_iterations 200 ! nlp sub-problem max iterations
minlp_branch_method 1 ! 1 = depth first, 2 = breadth first
minlp_gap_tol 0.001 ! covergence tolerance
minlp_integer_tol 0.001 ! maximum deviation from whole number to be considered an integer
minlp_integer_leaves 0 ! create soft (1) integer leaves or hard (2) integer leaves with branching
End File

Constants
c0 = 0.4
c1 = 0.2

Parameters
last

Variables
x0 = 0.5 , >= 0
x1 = 0.7 , >= 0
x2 = 0.0 , >= 0
int_w = 0 , >= 0 , <= 1

Intermediates
w = int_w

Equations
minimize last * x2

$x0 = x0 - x0*x1 - c0*x0*w
$x1 = - x1 + x0*x1 - c1*x1*w
$x2 = (x0-1)^2 + (x1-1)^2
'''
fid = open('lotka_volterra.apm','w')
fid.write(model)
fid.close()

# write data file
time = np.linspace(0,12,121)
time = np.insert(time, 1, 0.01)
last = np.zeros(122)
last[-1] = 1.0
data = np.vstack((time,last))
np.savetxt('data.csv',data.T,delimiter=',',header='time,last',comments='')

# specify server and application name
s = 'http://byu.apmonitor.com'
#s = 'http://127.0.0.1/' # for local APMonitor server
a = 'lotka'

apm(s,a,'clear all')
apm_load(s,a,'lotka_volterra.apm')
csv_load(s,a,'data.csv')

apm_option(s,a,'nlc.imode',6) # Nonlinear control / dynamic optimization
apm_option(s,a,'nlc.nodes',3)

apm_info(s,a,'MV','int_w') # M or MV = Manipulated variable - independent variable over time horizon
apm_option(s,a,'int_w.status',1) # Status: 1=ON, 0=OFF
apm_option(s,a,'int_w.mv_type',0) # MV Type = Zero Order Hold

apm_option(s,a,'nlc.solver',1) # 1 = APOPT

# solve
output = apm(s,a,'solve')
print(output)

# retrieve solution
y = apm_sol(s,a)

plt.figure(1)
plt.step(y['time'],y['int_w'],'r-',label='w (0/1)')
plt.plot(y['time'],y['x0'],'b-',label=r'$x_0$')
plt.plot(y['time'],y['x1'],'k-',label=r'$x_1$')
plt.plot(y['time'],y['x2'],'g-',label=r'$x_2$')
plt.xlabel('Time')
plt.ylabel('Variables')
plt.legend(loc='best')
plt.show()
</source>

== Results with APOPT (MINLP) ==

An MINLP solution is calculated with [https://apopt.com|APOPT] with an objective function value of <math>x_2(t_f) = 1.36</math>. APOPT requires 52 NLP solutions to find an integer solution (111.0 seconds of processing time).

[[File:Volterra_fishing_APMonitor.png]]

<pre>
----------------------------------------------------------------
APMonitor, Version 0.7.9
APMonitor Optimization Suite
----------------------------------------------------------------

--------- APM Model Size ------------
Each time step contains
Objects : 0
Constants : 2
Variables : 5
Intermediates: 1
Connections : 0
Equations : 5
Residuals : 4

Number of state variables: 2178
Number of total equations: - 2057
Number of slack variables: - 0
---------------------------------------
Degrees of freedom : 121

----------------------------------------------
Dynamic Control with APOPT Solver
----------------------------------------------
Iter: 1 I: 0 Tm: 12.99 NLPi: 93 Dpth: 0 Lvs: 2 Obj: 1.34E+00 Gap: NaN
Iter: 2 I: 0 Tm: 2.82 NLPi: 18 Dpth: 1 Lvs: 3 Obj: 1.34E+00 Gap: NaN
Iter: 3 I: 0 Tm: 3.39 NLPi: 30 Dpth: 2 Lvs: 4 Obj: 1.34E+00 Gap: NaN
Iter: 4 I: 0 Tm: 10.37 NLPi: 131 Dpth: 3 Lvs: 5 Obj: 1.34E+00 Gap: NaN
Iter: 5 I: 0 Tm: 1.81 NLPi: 10 Dpth: 4 Lvs: 6 Obj: 1.34E+00 Gap: NaN
Iter: 6 I: 0 Tm: 7.93 NLPi: 102 Dpth: 5 Lvs: 7 Obj: 1.34E+00 Gap: NaN
Iter: 7 I: 0 Tm: 3.22 NLPi: 28 Dpth: 6 Lvs: 8 Obj: 1.34E+00 Gap: NaN
Iter: 8 I: 0 Tm: 1.85 NLPi: 12 Dpth: 7 Lvs: 9 Obj: 1.34E+00 Gap: NaN
Iter: 9 I: 0 Tm: 10.24 NLPi: 140 Dpth: 8 Lvs: 10 Obj: 1.34E+00 Gap: NaN
Iter: 10 I: 0 Tm: 2.04 NLPi: 9 Dpth: 9 Lvs: 11 Obj: 1.35E+00 Gap: NaN
Iter: 11 I: 0 Tm: 6.50 NLPi: 86 Dpth: 10 Lvs: 12 Obj: 1.35E+00 Gap: NaN
Iter: 12 I: 0 Tm: 1.94 NLPi: 17 Dpth: 11 Lvs: 13 Obj: 1.35E+00 Gap: NaN
Iter: 13 I: 0 Tm: 2.17 NLPi: 18 Dpth: 12 Lvs: 14 Obj: 1.35E+00 Gap: NaN
Iter: 14 I: 0 Tm: 3.56 NLPi: 44 Dpth: 13 Lvs: 15 Obj: 1.35E+00 Gap: NaN
Iter: 15 I: 0 Tm: 2.34 NLPi: 23 Dpth: 14 Lvs: 16 Obj: 1.35E+00 Gap: NaN
Iter: 16 I: 0 Tm: 2.46 NLPi: 26 Dpth: 15 Lvs: 17 Obj: 1.35E+00 Gap: NaN
Iter: 17 I: 0 Tm: 2.72 NLPi: 26 Dpth: 16 Lvs: 18 Obj: 1.35E+00 Gap: NaN
Iter: 18 I: 0 Tm: 3.99 NLPi: 60 Dpth: 17 Lvs: 19 Obj: 1.35E+00 Gap: NaN
Iter: 19 I: 0 Tm: 1.91 NLPi: 10 Dpth: 18 Lvs: 20 Obj: 1.35E+00 Gap: NaN
Iter: 20 I: 0 Tm: 1.23 NLPi: 7 Dpth: 19 Lvs: 21 Obj: 1.35E+00 Gap: NaN
Iter: 21 I: 0 Tm: 1.37 NLPi: 8 Dpth: 20 Lvs: 22 Obj: 1.35E+00 Gap: NaN
Iter: 22 I: 0 Tm: 1.37 NLPi: 12 Dpth: 21 Lvs: 23 Obj: 1.35E+00 Gap: NaN
Iter: 23 I: 0 Tm: 1.34 NLPi: 11 Dpth: 22 Lvs: 24 Obj: 1.36E+00 Gap: NaN
Iter: 24 I: 0 Tm: 1.30 NLPi: 8 Dpth: 23 Lvs: 25 Obj: 1.36E+00 Gap: NaN
Iter: 25 I: 0 Tm: 1.33 NLPi: 14 Dpth: 24 Lvs: 26 Obj: 1.36E+00 Gap: NaN
Iter: 26 I: 0 Tm: 1.13 NLPi: 7 Dpth: 25 Lvs: 27 Obj: 1.36E+00 Gap: NaN
Iter: 27 I: 0 Tm: 0.97 NLPi: 7 Dpth: 26 Lvs: 28 Obj: 1.36E+00 Gap: NaN
Iter: 28 I: 0 Tm: 0.99 NLPi: 6 Dpth: 27 Lvs: 29 Obj: 1.36E+00 Gap: NaN
Iter: 29 I: 0 Tm: 0.93 NLPi: 6 Dpth: 28 Lvs: 30 Obj: 1.36E+00 Gap: NaN
Iter: 30 I: 0 Tm: 0.66 NLPi: 5 Dpth: 29 Lvs: 31 Obj: 1.36E+00 Gap: NaN
Iter: 31 I: 0 Tm: 0.73 NLPi: 5 Dpth: 30 Lvs: 32 Obj: 1.36E+00 Gap: NaN
Iter: 32 I: 0 Tm: 0.66 NLPi: 5 Dpth: 31 Lvs: 33 Obj: 1.36E+00 Gap: NaN
Iter: 33 I: 0 Tm: 0.70 NLPi: 5 Dpth: 32 Lvs: 34 Obj: 1.36E+00 Gap: NaN
Iter: 34 I: 0 Tm: 0.67 NLPi: 5 Dpth: 33 Lvs: 35 Obj: 1.36E+00 Gap: NaN
Iter: 35 I: 0 Tm: 0.82 NLPi: 9 Dpth: 34 Lvs: 36 Obj: 1.36E+00 Gap: NaN
Iter: 36 I: 0 Tm: 0.77 NLPi: 8 Dpth: 35 Lvs: 37 Obj: 1.36E+00 Gap: NaN
Iter: 37 I: 0 Tm: 0.75 NLPi: 8 Dpth: 36 Lvs: 38 Obj: 1.36E+00 Gap: NaN
Iter: 38 I: 0 Tm: 0.69 NLPi: 6 Dpth: 37 Lvs: 39 Obj: 1.36E+00 Gap: NaN
Iter: 39 I: 0 Tm: 0.71 NLPi: 6 Dpth: 38 Lvs: 40 Obj: 1.36E+00 Gap: NaN
Iter: 40 I: 0 Tm: 0.81 NLPi: 9 Dpth: 39 Lvs: 41 Obj: 1.36E+00 Gap: NaN
Iter: 41 I: 0 Tm: 0.70 NLPi: 8 Dpth: 40 Lvs: 42 Obj: 1.36E+00 Gap: NaN
Iter: 42 I: 0 Tm: 0.69 NLPi: 6 Dpth: 41 Lvs: 43 Obj: 1.36E+00 Gap: NaN
Iter: 43 I: 0 Tm: 0.67 NLPi: 5 Dpth: 42 Lvs: 44 Obj: 1.36E+00 Gap: NaN
Iter: 44 I: 0 Tm: 0.60 NLPi: 4 Dpth: 43 Lvs: 45 Obj: 1.36E+00 Gap: NaN
Iter: 45 I: 0 Tm: 0.71 NLPi: 6 Dpth: 44 Lvs: 46 Obj: 1.36E+00 Gap: NaN
Iter: 46 I: 0 Tm: 0.59 NLPi: 4 Dpth: 45 Lvs: 47 Obj: 1.36E+00 Gap: NaN
Iter: 47 I: 0 Tm: 0.69 NLPi: 6 Dpth: 46 Lvs: 48 Obj: 1.36E+00 Gap: NaN
Iter: 48 I: 0 Tm: 0.65 NLPi: 5 Dpth: 47 Lvs: 49 Obj: 1.36E+00 Gap: NaN
--Integer Solution: 1.36E+00 Lowest Leaf: 1.34E+00 Gap: 1.35E-02
Iter: 49 I: 0 Tm: 0.57 NLPi: 3 Dpth: 48 Lvs: 48 Obj: 1.36E+00 Gap: 1.35E-02
Iter: 50 I: 0 Tm: 1.29 NLPi: 8 Dpth: 48 Lvs: 47 Obj: 1.36E+00 Gap: 1.35E-02
Warning: best integer solution returned after maximum MINLP iterations
Adjust minlp_max_iter_with_int_sol 50 in apopt.opt to change limit
Successful solution

---------------------------------------------------
Solver : APOPT (v1.0)
Solution time : 111.364099999999 sec
Objective : 1.36258198934523
Successful solution
---------------------------------------------------
</pre>

[[Category:APMonitor]]

Lotka Volterra fishing problem (APMonitor)

2017-11-20T22:03:56Z

JohnHedengren: APMontor solution of Volterra fishing problem

This page contains a solution of the MIOCP [[Lotka Volterra fishing problem]] in [http://www.apmonitor.com APMonitor] Python format. A MATLAB version is also available from the [http://apmonitor.com/do/index.php/Main/DiscreteVariables Dynamic Optimization Course] as [http://apmonitor.com/do/uploads/Main/lotka_volterra_fishing.zip Example 3 (lotka_volterra_fishing.zip)].

=== APMonitor ===

The model in Python code for a fixed control discretization grid using orthogonal collocation and a simultaneous optimization method. The APMonitor package is available with '''pip install APMonitor''' or from the [https://github.com/APMonitor/apm_python/blob/master/apm.py APMonitor Python Github repository].

<source lang="Python">
import numpy as np
import matplotlib.pyplot as plt

# retrieve apm.py from
# https://raw.githubusercontent.com/APMonitor/apm_python/master/apm.py
# or
# http://apmonitor.com/wiki/index.php/Main/PythonApp
# from apm import *

# pip install with 'pip install APMonitor'
from APMonitor.apm import *

# local APMonitor servers are available for Windows or Linux
# http://apmonitor.com/wiki/index.php/Main/APMonitorServer
# with clients in Python, MATLAB, and Julia

# write model
model = '''
! apopt MINLP solver options (see apopt.com)
File apopt.opt
minlp_maximum_iterations 1000 ! minlp iterations
minlp_max_iter_with_int_sol 50 ! minlp iterations if integer solution is found
minlp_as_nlp 0 ! treat minlp as nlp
nlp_maximum_iterations 200 ! nlp sub-problem max iterations
minlp_branch_method 1 ! 1 = depth first, 2 = breadth first
minlp_gap_tol 0.001 ! covergence tolerance
minlp_integer_tol 0.001 ! maximum deviation from whole number to be considered an integer
minlp_integer_leaves 0 ! create soft (1) integer leaves or hard (2) integer leaves with branching
End File

Constants
c0 = 0.4
c1 = 0.2

Parameters
last

Variables
x0 = 0.5 , >= 0
x1 = 0.7 , >= 0
x2 = 0.0 , >= 0
int_w = 0 , >= 0 , <= 1

Intermediates
w = int_w

Equations
minimize last * x2

$x0 = x0 - x0*x1 - c0*x0*w
$x1 = - x1 + x0*x1 - c1*x1*w
$x2 = (x0-1)^2 + (x1-1)^2
'''
fid = open('lotka_volterra.apm','w')
fid.write(model)
fid.close()

# write data file
time = np.linspace(0,12,121)
time = np.insert(time, 1, 0.01)
last = np.zeros(122)
last[-1] = 1.0
data = np.vstack((time,last))
np.savetxt('data.csv',data.T,delimiter=',',header='time,last',comments='')

# specify server and application name
s = 'http://byu.apmonitor.com'
#s = 'http://127.0.0.1/' # for local APMonitor server
a = 'lotka'

apm(s,a,'clear all')
apm_load(s,a,'lotka_volterra.apm')
csv_load(s,a,'data.csv')

apm_option(s,a,'nlc.imode',6) # Nonlinear control / dynamic optimization
apm_option(s,a,'nlc.nodes',3)

apm_info(s,a,'MV','int_w') # M or MV = Manipulated variable - independent variable over time horizon
apm_option(s,a,'int_w.status',1) # Status: 1=ON, 0=OFF
apm_option(s,a,'int_w.mv_type',0) # MV Type = Zero Order Hold

apm_option(s,a,'nlc.solver',1) # 1 = APOPT

# solve
output = apm(s,a,'solve')
print(output)

# retrieve solution
y = apm_sol(s,a)

plt.figure(1)
plt.step(y['time'],y['int_w'],'r-',label='w (0/1)')
plt.plot(y['time'],y['x0'],'b-',label=r'$x_0$')
plt.plot(y['time'],y['x1'],'k-',label=r'$x_1$')
plt.plot(y['time'],y['x2'],'g-',label=r'$x_2$')
plt.xlabel('Time')
plt.ylabel('Variables')
plt.legend(loc='best')
plt.show()
</source>

== Results ==

=== IPOPT ===
The relaxed solution calculated by IPOPT (Ipopt 3.12, Linux x86_64, default settings, 4 GHz quadcore, Linux 4.2.0-23-generic) has an objective function value of <math>x_2(t_f) = 1.34428</math>. IPOPT requires 22 iterations (6.062 seconds of processing time). The following is a Gnuplot compatible tabular version of the solution:

<pre>
# t, x_0, x_1, u
0, 0.5, 0.7, 5.18436e-09
0.12, 0.519594, 0.659989, 5.42459e-09
0.24, 0.542429, 0.623852, 5.79496e-09
0.36, 0.568601, 0.591429, 6.30563e-09
0.48, 0.59823, 0.562571, 6.97415e-09
0.6, 0.631456, 0.537143, 7.827e-09
0.72, 0.66843, 0.515028, 8.90236e-09
0.84, 0.709312, 0.496134, 1.02545e-08
0.96, 0.75426, 0.480401, 1.19611e-08
1.08, 0.803421, 0.4678, 1.41346e-08
1.2, 0.856919, 0.458344, 1.69417e-08
1.32, 0.914842, 0.452091, 2.06377e-08
1.44, 0.97722, 0.449154, 2.56285e-08
1.56, 1.044, 0.449708, 3.2592e-08
1.68, 1.11503, 0.454003, 4.27322e-08
1.8, 1.18998, 0.462372, 5.83658e-08
1.92, 1.26837, 0.475251, 8.44574e-08
2.04, 1.34942, 0.493189, 1.33425e-07
2.16, 1.43209, 0.516862, 2.44977e-07
2.28, 1.51492, 0.547087, 6.24464e-07
2.4, 1.59605, 0.584816, 0.715281
2.52, 1.61776, 0.61833, 0.999999
2.64, 1.61116, 0.6499, 0.999999
2.76, 1.59844, 0.68229, 1
2.88, 1.57963, 0.714939, 1
3, 1.55497, 0.747196, 1
3.12, 1.52487, 0.778343, 0.999999
3.24, 1.48993, 0.807625, 0.999999
3.36, 1.4509, 0.83429, 0.999999
3.48, 1.40865, 0.857635, 0.999998
3.6, 1.36412, 0.877047, 0.999997
3.72, 1.31827, 0.892039, 0.999991
3.84, 1.27203, 0.90228, 0.999646
3.96, 1.22628, 0.907616, 0.536112
4.08, 1.20755, 0.919653, 0.581667
4.2, 1.18502, 0.928525, 0.503883
4.32, 1.16606, 0.936899, 0.464153
4.44, 1.14851, 0.944178, 0.418428
4.56, 1.13278, 0.950658, 0.378617
4.68, 1.11859, 0.956378, 0.341493
4.8, 1.10582, 0.961433, 0.3073
4.92, 1.09438, 0.965904, 0.276462
5.04, 1.08411, 0.969849, 0.248219
5.16, 1.07492, 0.973333, 0.222481
5.28, 1.0667, 0.976412, 0.19935
5.4, 1.05936, 0.97913, 0.178262
5.52, 1.0528, 0.981533, 0.159402
5.64, 1.04696, 0.983656, 0.142349
5.76, 1.04175, 0.985532, 0.127041
5.88, 1.0371, 0.987191, 0.1133
6, 1.03297, 0.988657, 0.100985
6.12, 1.02929, 0.989954, 0.0899604
6.24, 1.02601, 0.991101, 0.0800981
6.36, 1.0231, 0.992117, 0.0712934
6.48, 1.02051, 0.993015, 0.0634308
6.6, 1.01821, 0.993811, 0.0564127
6.72, 1.01617, 0.994515, 0.0501573
6.84, 1.01435, 0.995139, 0.0445824
6.96, 1.01274, 0.995691, 0.0396177
7.08, 1.01131, 0.99618, 0.0351979
7.2, 1.01003, 0.996613, 0.0312714
7.32, 1.0089, 0.996997, 0.027765
7.44, 1.0079, 0.997337, 0.0246626
7.56, 1.00701, 0.997639, 0.0218975
7.68, 1.00622, 0.997905, 0.0194423
7.8, 1.00552, 0.998142, 0.0172603
7.92, 1.0049, 0.998351, 0.0153224
8.04, 1.00435, 0.998537, 0.0136015
8.16, 1.00386, 0.998701, 0.0120735
8.28, 1.00342, 0.998847, 0.0107171
8.4, 1.00303, 0.998976, 0.00951328
8.52, 1.00269, 0.99909, 0.00844491
8.64, 1.00239, 0.999191, 0.00749702
8.76, 1.00212, 0.999281, 0.00665598
8.88, 1.00188, 0.99936, 0.00590975
9, 1.00167, 0.99943, 0.00524776
9.12, 1.00148, 0.999492, 0.00466052
9.24, 1.00131, 0.999547, 0.00413914
9.36, 1.00117, 0.999596, 0.00367206
9.48, 1.00103, 0.999637, 0.00326159
9.6, 1.00092, 0.999675, 0.00290302
9.72, 1.00082, 0.999707, 0.00258504
9.84, 1.00072, 0.999736, 0.00230262
9.96, 1.00064, 0.999762, 0.00205186
10.08, 1.00057, 0.999785, 0.00182951
10.2, 1.00051, 0.999805, 0.00163302
10.32, 1.00045, 0.999822, 0.00146039
10.44, 1.0004, 0.999838, 0.00131011
10.56, 1.00036, 0.999851, 0.00118097
10.68, 1.00032, 0.999863, 0.00107188
10.8, 1.00028, 0.999872, 0.000981659
10.92, 1.00025, 0.99988, 0.00090944
11.04, 1.00022, 0.999886, 0.000854893
11.16, 1.00019, 0.99989, 0.000818284
11.28, 1.00017, 0.999892, 0.000801181
11.4, 1.00014, 0.999891, 0.000807844
11.52, 1.00012, 0.999886, 0.000848734
11.64, 1.00009, 0.999878, 0.000951928
11.76, 1.00006, 0.999864, 0.00121852
11.88, 1.00002, 0.999838, 0.00267036
12, 0.999915, 0.999769, 0
</pre>

[[Category:APMonitor]]

Category:APMonitor

2017-11-20T21:43:46Z

JohnHedengren: /* Availibility */ minor clarification

'''Advanced process monitor (APMonitor)''', is a modeling language for mixed integer and differential algebraic equations. It is a free web-service or locally installed service for solving this class of problems. APMonitor is suited for large-scale problems and allows solutions of linear programming, integer programming, nonlinear programming, nonlinear mixed integer programming, dynamic simulation, moving horizon estimation, and nonlinear model predictive control. APMonitor does not solve the problems directly, but calls nonlinear programming solvers such as APOPT, BPOPT, IPOPT, MINOS, and SNOPT.

==Programming Language Integration==

Julia, MATLAB, and Python are mathematical programming languages that have APMonitor integration through web-service APIs. The interfaces are built-in optimization toolboxes or modules to both load and process solutions of optimization problems. APMonitor is an object-oriented modeling language and optimization suite that relies on programming languages to load, run, and retrieve solutions. APMonitor models and data are compiled at run-time and translated into objects that are solved by an optimization engine such as APOPT (MINLP) or IPOPT (NLP). The optimization engine is not specified by APMonitor, allowing several different optimization engines to be switched out. The simulation or optimization mode is also configurable to reconfigure the model for dynamic simulation, nonlinear model predictive control, moving horizon estimation or general problems in mathematical optimization.

==High Index DAEs==

The highest order of a derivative that is necessary to return a DAE to ODE form is called the ''differentiation index''. A standard way for dealing with high-index DAEs is to differentiate the equations to put them in index-1 DAE or ODE form with Pantelides algorithm. However, this approach can cause a number of undesirable numerical issues such as instability. While the syntax is similar to other modeling languages, APMonitor solves DAEs of any index without rearrangement or differentiation.

== Availibility ==
The APMonitor Python package can be installed with '''pip install APMonitor''' or from the following Python code.

<source lang="python">
# Install APMonitor
import pip
pip.main(['install','APMonitor'])
</source>

Similar interfaces are available for MATLAB and Julia with minor differences in APMonitor syntax between the languages. APMonitor is freely available as a web service or a locally installed web server for Linux or Windows.

== Supported Problem Classes ==

* NLP, LP, QP, MILP, MINLP
* Optimal Control Problems with ODEs or DAEs
* Minimum Time Problems
* Optimal Design
* Model Calibration
* Parameter Estimation

== Algorithms ==

Dynamic optimization uses direct collocation solved with large-scale nonlinear programming solvers. APMonitor provides exact first and second derivatives of continuous functions to the solvers through automatic differentiation and in sparse matrix form.

== References ==

* [https://apmonitor.com APMonitor Optimization Suite]
* [https://apmonitor.com/do Dynamic Optimization Course]
* [https://en.wikipedia.org/wiki/APMonitor APMonitor (Wikipedia)]


[[Category: Implementation]]

Category:APMonitor

2017-11-20T21:40:26Z

JohnHedengren: Creation of APMonitor page

'''Advanced process monitor (APMonitor)''', is a modeling language for mixed integer and differential algebraic equations. It is a free web-service or locally installed service for solving this class of problems. APMonitor is suited for large-scale problems and allows solutions of linear programming, integer programming, nonlinear programming, nonlinear mixed integer programming, dynamic simulation, moving horizon estimation, and nonlinear model predictive control. APMonitor does not solve the problems directly, but calls nonlinear programming solvers such as APOPT, BPOPT, IPOPT, MINOS, and SNOPT.

==Programming Language Integration==

Julia, MATLAB, and Python are mathematical programming languages that have APMonitor integration through web-service APIs. The interfaces are built-in optimization toolboxes or modules to both load and process solutions of optimization problems. APMonitor is an object-oriented modeling language and optimization suite that relies on programming languages to load, run, and retrieve solutions. APMonitor models and data are compiled at run-time and translated into objects that are solved by an optimization engine such as APOPT (MINLP) or IPOPT (NLP). The optimization engine is not specified by APMonitor, allowing several different optimization engines to be switched out. The simulation or optimization mode is also configurable to reconfigure the model for dynamic simulation, nonlinear model predictive control, moving horizon estimation or general problems in mathematical optimization.

==High Index DAEs==

The highest order of a derivative that is necessary to return a DAE to ODE form is called the ''differentiation index''. A standard way for dealing with high-index DAEs is to differentiate the equations to put them in index-1 DAE or ODE form with Pantelides algorithm. However, this approach can cause a number of undesirable numerical issues such as instability. While the syntax is similar to other modeling languages, APMonitor solves DAEs of any index without rearrangement or differentiation.

== Availibility ==
The APMonitor Python package can be installed with '''pip install APMonitor''' or from the following Python code.

<source lang="python">
# Install APMonitor
import pip
pip.main(['install','APMonitor'])
</source>

Similar interfaces are available for MATLAB and Julia with minor differences in syntax between the languages. APMonitor is freely available as a web service or a locally installed web server for Linux or Windows.

== Supported Problem Classes ==

* NLP, LP, QP, MILP, MINLP
* Optimal Control Problems with ODEs or DAEs
* Minimum Time Problems
* Optimal Design
* Model Calibration
* Parameter Estimation

== Algorithms ==

Dynamic optimization uses direct collocation solved with large-scale nonlinear programming solvers. APMonitor provides exact first and second derivatives of continuous functions to the solvers through automatic differentiation and in sparse matrix form.

== References ==

* [https://apmonitor.com APMonitor Optimization Suite]
* [https://apmonitor.com/do Dynamic Optimization Course]
* [https://en.wikipedia.org/wiki/APMonitor APMonitor (Wikipedia)]


[[Category: Implementation]]

Lotka Volterra fishing problem

2017-11-20T20:39:13Z

JohnHedengren: /* Source Code */ Added APMonitor, alphebetized list

{{Dimensions
|nd = 1
|nx = 3
|nw = 1
|nre = 3
}}The '''Lotka Volterra fishing problem''' looks for an optimal fishing strategy to be performed on a fixed time horizon to bring the biomasses of both predator as prey fish to a prescribed steady state. The problem was set up as a small-scale benchmark problem.
The well known [http://en.wikipedia.org/wiki/Lotka_volterra Lotka Volterra equations] for a predator-prey system have been augmented by an additional linear term, relating to fishing by man. The control can be regarded both in a relaxed, as in a discrete manner, corresponding to a part of the fleet, or the full fishing fleet.

The mathematical equations form a small-scale [[:Category:ODE model|ODE model]]. The interior point equality conditions fix the initial values of the differential states.

The optimal integer control functions shows [[:Category:Chattering|chattering]] behavior, making the Lotka Volterra fishing problem an ideal candidate for benchmarking of algorithms.

== Mathematical formulation ==

The mixed-integer optimal control problem is given by


<math>
\begin{array}{llclr}
\displaystyle \min_{x, w} & x_2(t_f) \\[1.5ex]
\mbox{s.t.}
& \dot{x}_0 & = & x_0 - x_0 x_1 - \; c_0 x_0 \; w, \\
& \dot{x}_1 & = & - x_1 + x_0 x_1 - \; c_1 x_1 \; w, \\
& \dot{x}_2 & = & (x_0 - 1)^2 + (x_1 - 1)^2, \\[1.5ex]
& x(0) &=& (0.5, 0.7, 0)^T, \\
& w(t) &\in& \{0, 1\}.
\end{array}
</math>


Here the differential states <math>(x_0, x_1)</math> describe the biomasses of prey and predator, respectively. The third differential state is used here to transform the objective, an integrated deviation, into the Mayer formulation <math>\min \; x_2(t_f)</math>. The decision, whether the fishing fleet is actually fishing at time <math>t</math> is denoted by <math>w(t)</math>.

== Parameters ==

These fixed values are used within the model.

<math>
\begin{array}{rcl}
[t_0, t_f] &=& [0, 12],\\
(c_0, c_1) &=& (0.4, 0.2).
\end{array}
</math>

== Reference Solutions ==

If the problem is relaxed, i.e., we demand that <math>w(t)</math> be in the continuous interval <math>[0, 1]</math> instead of the binary choice <math>\{0,1\}</math>, the optimal solution can be determined by means of [http://en.wikipedia.org/wiki/Pontryagin%27s_minimum_principle Pontryagins maximum principle]. The optimal solution contains a singular arc, as can be seen in the plot of the optimal control. The two differential states and corresponding adjoint variables in the indirect approach are also displayed. A different approach to solving the relaxed problem is by using a direct method such as collocation or Bock's direct multiple shooting method. Optimal solutions for different control discretizations are also plotted in the leftmost figure.

The optimal objective value of this relaxed problem is <math>x_2(t_f) = 1.34408</math>. As follows from MIOC theory <bib id="Sager2011d" /> this is the best lower bound on the optimal value of the original problem with the integer restriction on the control function. In other words, this objective value can be approximated arbitrarily close, if the control only switches often enough between 0 and 1. As no optimal solution exists, two suboptimal ones are shown, one with only two switches and an objective function value of <math>x_2(t_f) = 1.38276</math>, and one with 56 switches and <math>x_2(t_f) = 1.34416</math>.

<gallery caption="Reference solution plots" widths="180px" heights="140px" perrow="2">
Image:lotkaRelaxedControls.png| Optimal relaxed control determined by an indirect approach and by a direct approach on different control discretization grids.
Image:lotkaindirektStates.png| Differential states and corresponding adjoint variables in the indirect approach.
Image:lotka2Switches.png| Control and differential states with only two switches.
Image:lotka56Switches.png| Control and differential states with 56 switches.
</gallery>

== Source Code ==

Model descriptions are available in

* [[:Category:ACADO | ACADO code]] at [[Lotka Volterra fishing problem (ACADO)]]
* [[:Category:AMPL | AMPL code]] at [[Lotka Volterra fishing problem (AMPL)]]
* [[:Category:APMonitor | APMonitor code]] at [[Lotka Volterra fishing problem (APMonitor)]]
* [[:Category:Bocop | Bocop code]] at [[Lotka Volterra fishing problem (Bocop)]]
* [[:Category:Casadi | Casadi code]] at [[Lotka Volterra fishing problem (Casadi)]]
* [[:Category:JModelica | JModelica code]] at [[Lotka Volterra fishing problem (JModelica)]]
* [[:Category:Julia/JuMP | JuMP code]] at [[Lotka Volterra fishing problem (JuMP)]]
* [[:Category:Muscod | Muscod code]] at [[Lotka Volterra fishing problem (Muscod)]]
* [[:Category:switch | switch code]] at [[Lotka Volterra fishing problem (switch)]]
* [[:Category:TomDyn/PROPT | PROPT code]] at [[Lotka Volterra fishing problem (TomDyn/PROPT)]]

== Variants ==

There are several alternative formulations and variants of the above problem, in particular

* a prescribed time grid for the control function <bib id="Sager2006" />, see also [[Lotka Volterra fishing problem (AMPL)]],
* a time-optimal formulation to get into a steady-state <bib id="Sager2005" />,
* the usage of a different target steady-state, as the one corresponding to <math> w(t) = 1</math> which is <math>(1 + c_1, 1 - c_0)</math>,
* different fishing control functions for the two species,
* different parameters and start values.

== Miscellaneous and Further Reading ==
The Lotka Volterra fishing problem was introduced by Sebastian Sager in a proceedings paper <bib id="Sager2006" /> and revisited in his PhD thesis <bib id="Sager2005" />. These are also the references to look for more details.

== References ==
<biblist />


[[Category:MIOCP]]
[[Category:ODE model]]
[[Category:Tracking objective]]
[[Category:Chattering]]
[[Category:Sensitivity-seeking arcs]]
[[Category:Population dynamics]]