Difference between revisions of "Lotka Volterra fishing problem"

From mintOC
Jump to: navigation, search
m
(Source Code)
 
(44 intermediate revisions by 8 users not shown)
Line 4: Line 4:
 
|nw        = 1
 
|nw        = 1
 
|nre      = 3
 
|nre      = 3
}}
+
}}<!-- Do not insert line break here or Dimensions Box moves up in the layout...
  
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 '''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 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.
  
Line 16: Line 16:
 
== Mathematical formulation ==
 
== Mathematical formulation ==
  
For <math>t \in [t_0, t_f]</math> almost everywhere the mixed-integer optimal control problem is given by
+
The mixed-integer optimal control problem is given by
  
 +
<p>
 
<math>
 
<math>
\begin{array}{llcl}
+
\begin{array}{llclr}
 
  \displaystyle \min_{x, w} & x_2(t_f)  \\[1.5ex]
 
  \displaystyle \min_{x, w} & x_2(t_f)  \\[1.5ex]
  \mbox{s.t.} & \dot{x}_0(t) & = & x_0(t) - x_0(t) x_1(t) - \; c_0 x_0(t) \; w(t), \\
+
  \mbox{s.t.}  
  & \dot{x}_1(t) & = & - x_1(t) + x_0(t) x_1(t) - \; c_1 x_1(t) \; w(t),  \\
+
& \dot{x}_0 & = & x_0 - x_0 x_1 - \; c_0 x_0 \; w, \\
  & \dot{x}_2(t) & = & (x_0(t) - 1)^2 + (x_1(t) - 1)^2,  \\[1.5ex]
+
  & \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, \\
 
  & x(0) &=& (0.5, 0.7, 0)^T, \\
 
  & w(t) &\in&  \{0, 1\}.
 
  & w(t) &\in&  \{0, 1\}.
 
\end{array}  
 
\end{array}  
 
</math>
 
</math>
 +
</p>
  
 
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>.
 
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>.
Line 46: Line 49:
 
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.
 
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<bibref>Sager2008</bibref> 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>.
+
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">
 
<gallery caption="Reference solution plots" widths="180px" heights="140px" perrow="2">
Line 58: Line 61:
 
== Source Code ==
 
== Source Code ==
  
=== C ===
+
Model descriptions are available in
  
The differential equations in C code:
+
* [[:Category:ACADO | ACADO code]] at [[Lotka Volterra fishing problem (ACADO)]]
<source lang="cpp">
+
* [[:Category:AMPL | AMPL code]] at [[Lotka Volterra fishing problem (AMPL)]]
/* steady state with u == 0 */
+
* [[:Category:APMonitor | APMonitor code]] at [[Lotka Volterra fishing problem (APMonitor)]]
double ref0 = 1, ref1 = 1;
+
* [[:Category:Bocop | Bocop code]] at [[Lotka Volterra fishing problem (Bocop)]]
 
+
* [[:Category:Casadi | Casadi code]] at [[Lotka Volterra fishing problem (Casadi)]]
/* Biomass of prey */
+
* [[:Category:Gekko | GEKKO Python code]] at [[Lotka Volterra fishing problem (GEKKO)]]
rhs[0] =  xd[0] - xd[0]*xd[1] - p[0]*u[0]*xd[0];
+
* [[:Category:JModelica | JModelica code]] at [[Lotka Volterra fishing problem (JModelica)]]
/* Biomass of predator */
+
* [[:Category:Julia/JuMP | JuMP code]] at [[Lotka Volterra fishing problem (JuMP)]]
rhs[1] = - xd[1] + xd[0]*xd[1] - p[1]*u[0]*xd[1];
+
* [[:Category:Muscod | Muscod code]] at [[Lotka Volterra fishing problem (Muscod)]]
/* Deviation from reference trajectory */
+
* [[:Category:switch | switch code]] at [[Lotka Volterra fishing problem (switch)]]
rhs[2] = (xd[0]-ref0)*(xd[0]-ref0) + (xd[1]-ref1)*(xd[1]-ref1);
+
* [[:Category:TomDyn/PROPT | PROPT code]] at [[Lotka Volterra fishing problem (TomDyn/PROPT)]]
</source>
+
* [[:Category:Julia/JuMP | Julia code]] at [[Lotka Volterra fishing problem (Julia Neural Network solve)]]
 
+
=== AMPL ===
+
 
+
The model in AMPL code for a fixed control discretization grid with a collocation method. We need a model file lotka_ampl.mod,
+
<source lang="AMPL">
+
# ----------------------------------------------------------------
+
# Lotka Volterra fishing problem with collocation (explicit Euler)
+
# (c) Sebastian Sager
+
# ----------------------------------------------------------------
+
 
+
param T    > 0;
+
param nt  > 0;
+
param nu  > 0;
+
param ntperu > 0;
+
param c1  > 0;
+
param c2  > 0;
+
param ref1 > 0;
+
param ref2 > 0;
+
param dt := T / (nt-1);
+
set I:= 0..nt;
+
set U:= 0..nu-1;
+
 
+
param uidx {I};
+
 
+
var x {I, 1..2} >= 0;
+
var w {U} >= 0, <= 1 binary;
+
 
+
minimize Deviation:
+
  0.5 * dt * ( (x[0,1] - ref1)*(x[0,1] - ref1) + (x[0,2] - ref2)*(x[0,2] - ref2) )
+
+ 0.5 * dt * ( (x[nt,1] - ref1)*(x[nt,1] - ref1) + (x[nt,2] - ref2)*(x[nt,2] - ref2) )
+
+ dt * sum {i in I diff {0,nt} } ( (x[i,1] - ref1)*(x[i,1] - ref1)
+
                                + (x[i,2] - ref2)*(x[i,2] - ref2) ) ;
+
 
+
subj to ODE_DISC_1 {i in I diff {0}}:
+
  x[i,1] = x[i-1,1] + dt * ( x[i-1,1] - x[i-1,1]*x[i-1,2] - x[i-1,1]*c1*w[uidx[i-1]] );
+
+
subj to ODE_DISC_2 {i in I diff {0}}:
+
  x[i,2] = x[i-1,2] + dt * ( - x[i-1,2] + x[i-1,1]*x[i-1,2] - x[i-1,2]*c2*w[uidx[i-1]] );
+
</source>
+
a data file lotka_ampl.dat,
+
<source lang="AMPL">
+
# ------------------------------------
+
# Data: Lotka Volterra fishing problem
+
# ------------------------------------
+
 
+
# Algorithmic parameters
+
param ntperu := 1;
+
param nu := 100;
+
param nt := 100;
+
 
+
# Problem parameters
+
param T := 12.0;
+
param c1 := 0.4;
+
param c2 := 0.2;
+
param ref1 := 1.0;
+
param ref2 := 1.0;
+
 
+
# Initial values differential states
+
let x[0,1] := 0.5;
+
let x[0,2] := 0.7;
+
fix x[0,1];
+
fix x[0,2];
+
 
+
# Initial values control
+
let {i in U} w[i] := 0.0;
+
 
+
param mysum;
+
 
+
# Set indices of controls corresponding to time points
+
for {i in 0..nu-1} {
+
for {j in 0..ntperu-1} {
+
let uidx[i*ntperu+j] := i;
+
}
+
}
+
let uidx[nt] := nu-1;
+
</source>
+
and a running script lotka_ampl.run,
+
<source lang="AMPL">
+
# ----------------------------------------------------
+
# Solve Lotka Volterra fishing problem via collocation
+
# ----------------------------------------------------
+
 
+
model ampl_lotka.mod;
+
data ampl_lotka.dat;
+
 
+
option solver bonmin;
+
 
+
solve;
+
</source>
+
  
 
== Variants ==
 
== Variants ==
Line 166: Line 80:
 
There are several alternative formulations and variants of the above problem, in particular
 
There are several alternative formulations and variants of the above problem, in particular
  
* a prescribed time grid for the control function <bibref>Sager2006</bibref>,
+
* 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 <bibref>Sager2005</bibref>,
+
* 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>,
+
* 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>, see [[Lotka Volterra multi-arcs problem]]
* different fishing control functions for the two species,
+
* different fishing control functions for the two species, see [[Lotka Volterra Multimode fishing problem]]
 +
* different fishing control functions that fish an absolute value from the two species, see [[Lotka Volterra absolute fishing problem]]
 +
* a terminal constrained formulation, where a violation is penalized via slack variables, see [[Lotka Volterra (terminal constraint violation)]]
 
* different parameters and start values.
 
* different parameters and start values.
  
 
== Miscellaneous and Further Reading ==
 
== Miscellaneous and Further Reading ==
The Lotka Volterra fishing problem was introduced by Sebastian Sager in a proceedings paper <bibref>Sager2006</bibref> and revisited in his PhD thesis <bibref>Sager2005</bibref>. These are also the references to look for more details.
+
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.
  
<!--Testing Graphviz
 
<graphviz border='frame' format='svg'>
 
digraph G {Hello->World!}
 
</graphviz>
 
-->
 
 
== References ==
 
== References ==
<bibreferences/>
+
<biblist />
  
 
<!--List of all categories this page is part of. List characterization of solution behavior, model properties, ore presence of implementation details (e.g., AMPL for AMPL model) here -->
 
<!--List of all categories this page is part of. List characterization of solution behavior, model properties, ore presence of implementation details (e.g., AMPL for AMPL model) here -->
Line 188: Line 99:
 
[[Category:Tracking objective]]
 
[[Category:Tracking objective]]
 
[[Category:Chattering]]
 
[[Category:Chattering]]
[[Category:AMPL]]
+
[[Category:Sensitivity-seeking arcs]]
 
[[Category:Population dynamics]]
 
[[Category:Population dynamics]]
 +
 +
 +
<!--Testing Graphviz
 +
<graphviz border='frame' format='svg'>
 +
digraph G {Hello->World!}
 +
</graphviz>
 +
-->

Latest revision as of 16:16, 20 September 2021

Lotka Volterra fishing problem
State dimension: 1
Differential states: 3
Discrete control functions: 1
Interior point equalities: 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 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 ODE model. The interior point equality conditions fix the initial values of the differential states.

The optimal integer control functions shows 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


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

Here the differential states (x_0, x_1) 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 \min \; x_2(t_f). The decision, whether the fishing fleet is actually fishing at time t is denoted by w(t).

Parameters

These fixed values are used within the model.


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

Reference Solutions

If the problem is relaxed, i.e., we demand that w(t) be in the continuous interval [0, 1] instead of the binary choice \{0,1\}, the optimal solution can be determined by means of 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 x_2(t_f) = 1.34408. As follows from MIOC theory [Sager2011d]Author: S. Sager
How published: University of Heidelberg
Month: August
Note: Habilitation
Title: On the Integration of Optimization Approaches for Mixed-Integer Nonlinear Optimal Control
Url: http://mathopt.de/PUBLICATIONS/Sager2011d.pdf
Year: 2011
Link to Google Scholar
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 x_2(t_f) = 1.38276, and one with 56 switches and x_2(t_f) = 1.34416.


Source Code

Model descriptions are available in

Variants

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

  • a prescribed time grid for the control function [Sager2006]Address: Heidelberg
    Author: S. Sager; H.G. Bock; M. Diehl; G. Reinelt; J.P. Schl\"oder
    Booktitle: Recent Advances in Optimization
    Editor: A. Seeger
    Note: ISBN 978-3-5402-8257-0
    Pages: 269--289
    Publisher: Springer
    Series: Lectures Notes in Economics and Mathematical Systems
    Title: Numerical methods for optimal control with binary control functions applied to a Lotka-Volterra type fishing problem
    Volume: 563
    Year: 2009
    Link to Google Scholar
    , see also Lotka Volterra fishing problem (AMPL),
  • a time-optimal formulation to get into a steady-state [Sager2005]Address: Tönning, Lübeck, Marburg
    Author: S. Sager
    Editor: ISBN 3-89959-416-9
    Publisher: Der andere Verlag
    Title: Numerical methods for mixed--integer optimal control problems
    Url: http://mathopt.de/PUBLICATIONS/Sager2005.pdf
    Year: 2005
    Link to Google Scholar
    ,
  • the usage of a different target steady-state, as the one corresponding to  w(t) = 1 which is (1 + c_1, 1 - c_0), see Lotka Volterra multi-arcs problem
  • different fishing control functions for the two species, see Lotka Volterra Multimode fishing problem
  • different fishing control functions that fish an absolute value from the two species, see Lotka Volterra absolute fishing problem
  • a terminal constrained formulation, where a violation is penalized via slack variables, see Lotka Volterra (terminal constraint violation)
  • different parameters and start values.

Miscellaneous and Further Reading

The Lotka Volterra fishing problem was introduced by Sebastian Sager in a proceedings paper [Sager2006]Address: Heidelberg
Author: S. Sager; H.G. Bock; M. Diehl; G. Reinelt; J.P. Schl\"oder
Booktitle: Recent Advances in Optimization
Editor: A. Seeger
Note: ISBN 978-3-5402-8257-0
Pages: 269--289
Publisher: Springer
Series: Lectures Notes in Economics and Mathematical Systems
Title: Numerical methods for optimal control with binary control functions applied to a Lotka-Volterra type fishing problem
Volume: 563
Year: 2009
Link to Google Scholar
and revisited in his PhD thesis [Sager2005]Address: Tönning, Lübeck, Marburg
Author: S. Sager
Editor: ISBN 3-89959-416-9
Publisher: Der andere Verlag
Title: Numerical methods for mixed--integer optimal control problems
Url: http://mathopt.de/PUBLICATIONS/Sager2005.pdf
Year: 2005
Link to Google Scholar
. These are also the references to look for more details.

References

There were no citations found in the article.