Commit 67720393 authored by YABO Agustin's avatar YABO Agustin
Browse files

Examples fixed

parent 3e55c5b3
......@@ -11,9 +11,9 @@ dependencies:
- matplotlib
- numba
- numpy
- bocop=3.1.0=py37h39e3cac_8
- nodejs
- bocop=3.1.1
- nutopy
- nodejs
- ipywidgets
- jupyterlab
- swig
......
This diff is collapsed.
......@@ -19,6 +19,7 @@ constant.2 0.5
constant.3 0.07
# Time discretisation
ode.discretization midpoint_implicit
time.steps 2000
# Bounds for constraints (initial conditions for our case)
......@@ -44,4 +45,13 @@ state.0.init 0.1
state.1.init 0.1
state.2.init 0.1
state.3.init 0.1
control.0.init 0.4679139313552532
\ No newline at end of file
control.0.init 0.4679139313552532
# Ipopt
ipoptIntOption.print_level 5
ipoptIntOption.max_iter 1000
ipoptStrOption.mu_strategy adaptive
ipoptNumOption.tol 1e-12
# Misc
ad.retape 0
......@@ -19,6 +19,7 @@ constant.2 0.3
constant.3 0.07
# Time discretisation
ode.discretization midpoint_implicit
time.steps 2000
# Bounds for constraints (initial conditions for our case)
......@@ -45,3 +46,13 @@ state.1.init 0.1
state.2.init 0.1
state.3.init 0.1
control.0.init 0.5524
# Ipopt
ipoptIntOption.print_level 5
ipoptIntOption.max_iter 1000
ipoptStrOption.mu_strategy adaptive
ipoptNumOption.tol 1e-12
# Misc
ad.retape 0
......@@ -19,6 +19,7 @@ constant.2 0.5
constant.3 0.07
# Time discretisation
ode.discretization midpoint_implicit
time.steps 2000
# Bounds for constraints (initial conditions for our case)
......@@ -44,4 +45,13 @@ state.0.init 0.1
state.1.init 0.1
state.2.init 0.1
state.3.init 0.1
control.0.init 0.5
\ No newline at end of file
control.0.init 0.5
# Ipopt
ipoptIntOption.print_level 5
ipoptIntOption.max_iter 1000
ipoptStrOption.mu_strategy adaptive
ipoptNumOption.tol 1e-12
# Misc
ad.retape 0
......@@ -20,6 +20,7 @@ constant.3 0.07
# Time discretisation
ode.discretization midpoint_implicit
time.steps 2000
# Bounds for constraints (initial conditions for our case)
......@@ -37,7 +38,6 @@ state.0.lowerbound 0
state.1.lowerbound 0
state.2.lowerbound 0
state.3.lowerbound 0
control.0.lowerbound 0
control.0.upperbound 1
......@@ -46,5 +46,13 @@ state.0.init 0.1
state.1.init 0.1
state.2.init 0.1
state.3.init 0.1
control.0.init 0.5
# Ipopt
ipoptIntOption.print_level 5
ipoptIntOption.max_iter 1000
ipoptStrOption.mu_strategy adaptive
ipoptNumOption.tol 1e-12
control.0.init 0.5
\ No newline at end of file
# Misc
ad.retape 0
......@@ -19,6 +19,7 @@ constant.2 0.3
constant.3 0.07
# Time discretisation
ode.discretization midpoint_implicit
time.steps 2000
# Bounds for constraints (initial conditions for our case)
......@@ -44,4 +45,13 @@ state.0.init 0.1
state.1.init 0.1
state.2.init 0.1
state.3.init 0.1
control.0.init 0.5
\ No newline at end of file
control.0.init 0.5
# Ipopt
ipoptIntOption.print_level 5
ipoptIntOption.max_iter 1000
ipoptStrOption.mu_strategy adaptive
ipoptNumOption.tol 1e-12
# Misc
ad.retape 0
# Definition file
# Dimensions
dim.state 4
dim.control 1
dim.boundaryconditions 4
dim.pathconstraints 0
dim.parameters 0
dim.constants 4
# Time interval
initial.time 0
final.time 20
# Constants
constant.0 0.7
constant.1 0.003
constant.2 0.3
constant.3 0.07
# Time discretisation
time.steps 2000
# Bounds for constraints (initial conditions for our case)
boundarycond.0.lowerbound 0.0458
boundarycond.0.upperbound 0.0458
boundarycond.1.lowerbound 0.249
boundarycond.1.upperbound 0.249
boundarycond.2.upperbound 0.25
boundarycond.2.upperbound 0.25
boundarycond.3.lowerbound 0.003
boundarycond.3.upperbound 0.003
# Bounds for variables (dynamical bound for each variable and control)
state.0.lowerbound 0
state.1.lowerbound 0
state.2.lowerbound 0
state.3.lowerbound 0
control.0.lowerbound 0.5524 #Optimal value calculated for our initial conditions
control.0.upperbound 0.5524
# Initialization for discretized problem (initial point for the optimization algorithm.. this we can touch only if the algorithm is not working, otherwise we can leave it like this)
state.0.init 0.1
state.1.init 0.1
state.2.init 0.1
state.3.init 0.1
control.0.init 0.5524
This is Ipopt version 3.13.3, running with linear solver mumps.
NOTE: Other linear solvers might be more efficient (see Ipopt documentation).
Number of nonzeros in equality constraint Jacobian...: 68003
Number of nonzeros in inequality constraint Jacobian.: 1
Number of nonzeros in Lagrangian Hessian.............: 52000
Total number of variables............................: 16004
variables with only lower bounds: 8004
variables with lower and upper bounds: 0
variables with only upper bounds: 0
Total number of equality constraints.................: 16003
Total number of inequality constraints...............: 1
inequality constraints with only lower bounds: 0
inequality constraints with lower and upper bounds: 0
inequality constraints with only upper bounds: 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
0 -1.0000000e-01 1.49e-01 2.00e+00 0.0 0.00e+00 - 0.00e+00 0.00e+00 0
1 -3.1931905e-01 8.67e-02 2.77e+00 -0.4 5.41e-01 - 3.53e-01 4.18e-01f 1
2 -3.1775095e-01 8.58e-02 1.34e+02 -6.8 1.53e-01 - 6.08e-01 1.03e-02h 1
3 -1.8448716e-01 5.86e-02 3.11e+02 -1.8 4.22e-01 - 1.00e+00 7.63e-01h 1
4 -6.5161268e-02 2.80e-02 1.44e+00 -2.2 1.19e-01 - 1.00e+00 1.00e+00h 1
5 -3.5807036e-02 1.98e-02 8.75e-01 -1.6 2.94e-02 - 1.00e+00 1.00e+00h 1
6 -3.1469765e-02 6.33e-03 4.20e+00 -1.7 2.06e-02 - 9.57e-01 1.00e+00f 1
7 -3.1177503e-02 2.83e-04 1.85e+00 -2.6 7.90e-03 - 9.98e-01 1.00e+00h 1
8 -3.1160454e-02 6.81e-07 5.64e-04 -4.2 4.74e-04 - 1.00e+00 1.00e+00h 1
9 -3.1160407e-02 8.12e-12 1.42e-03 -10.3 1.68e-06 - 9.99e-01 1.00e+00h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
10 -3.1161131e-02 5.51e-10 2.93e-11 -9.6 1.20e-05 - 1.00e+00 1.00e+00h 1
11 -3.1161147e-02 2.73e-13 5.92e-14 -12.3 2.67e-07 - 1.00e+00 1.00e+00h 1
12 -3.1161147e-02 9.02e-17 1.65e-15 -12.3 4.44e-11 - 1.00e+00 1.00e+00h 1
Number of Iterations....: 12
(scaled) (unscaled)
Objective...............: -3.1161146786448997e-02 -3.1161146786448997e-02
Dual infeasibility......: 1.6484280813939900e-15 1.6484280813939900e-15
Constraint violation....: 9.0205620750793969e-17 9.0205620750793969e-17
Complementarity.........: 5.0000136904343260e-13 5.0000136904343260e-13
Overall NLP error.......: 5.0000136904343260e-13 5.0000136904343260e-13
Number of objective function evaluations = 13
Number of objective gradient evaluations = 13
Number of equality constraint evaluations = 13
Number of inequality constraint evaluations = 13
Number of equality constraint Jacobian evaluations = 13
Number of inequality constraint Jacobian evaluations = 13
Number of Lagrangian Hessian evaluations = 12
Total CPU secs in IPOPT (w/o function evaluations) = 4.998
Total CPU secs in NLP function evaluations = 30.010
EXIT: Optimal Solution Found.
%% Cell type:markdown id: tags:
# Bistable gene-regulatory networks
## Model definition
We consider the piecewise linear system defined in Filippov sense as in [Augier and Yabo (2021)][id1]. In synthetic biology, this model can represent a _genetic toggle switch_, a synthetic regulatory network designed by [Gardner et al. (2000)][id2] in the bacteria E. coli through two genes lacI and tetR mutually repressing each other. Dynamics is given by
[id1]: https://hal.archives-ouvertes.fr/hal-03099681
[id2]: https://pubmed.ncbi.nlm.nih.gov/10659857/
$$
\left\{ \begin{array}{l}
\dot x_1 = -\gamma_1 x_1 + u(t) k_1 s^{-}(x_2, \theta_2), \\
\dot x_2 = -\gamma_1 x_2 + u(t)k_2 s^{-}(x_1, \theta_1),
\end{array}
\right.
$$
where for $j\in \{1,2\}$, $x_j\in \mathbb{R}$, and for $\theta\in \mathbb{R}$, $s^{-}(\cdot,\theta):\mathbb{R}\to \mathbb{R}$ is such that
$$
s^{-}(x,\theta)= \left\{ \begin{array}{ll}
1 & \textit{if } x < \theta, \\
0 & \textit{if } x > \theta,
\end{array} \right.
$$
and the control $u(\cdot) \in L^\infty ([0,t_f],[u_{\text{min}},u_{\text{max}}])$, with $0<u_{\text{min}}<1 \leq u_{\text{max}}$. Such system is piecewise linear in $\mathbb{R}^2$ and, in particular, is linear in the regular domains
$$
\begin{array}{ll}
&B_{00}=\left\{(x_1,x_2)\in \mathbb{R}^2 \mid 0<x_1<\theta_1, \ 0<x_2<\theta_2\right\},\\
& B_{01}=\left\{(x_1,x_2)\in \mathbb{R}^2 \mid 0<x_1<\theta_1, \ \theta_2<x_2<\frac{k_2}{\gamma_2}\right\},\\
& B_{10}=\left\{(x_1,x_2)\in \mathbb{R}^2 \mid \theta_1<x_1<\frac{k_1}{\gamma_1}, \ 0<x_2<\theta_2\right\},\\
&B_{11}=\left\{(x_1,x_2)\in \mathbb{R}^2 \mid \theta_1<x_1<\frac{k_1}{\gamma_1}, \ \theta_2<x_2<\frac{k_2}{\gamma_2}\right\}.
\end{array}
$$
The following image shows the dynamics for the open loop system (i.e. $u \equiv 1$):
<img src="openloop.png" width=436 height=337>
## Time-optimal control problem
The optimal control problem involves minimizing the time of a state transfer from $B_{10}$ to $B_{01}$, and so is defined for the state $x(t) = (x_1(t),x_2(t))$ as
$$
\left\{\begin{array}{l}
\textit{minimize} \ t_f \\
\dot x = f(x),\\
x(0) = x_0 \in B_{10}, \\ x_1(t_f)\in [0,\theta_1), \ x_2(t_f)=x_2^f > \theta_2, \\ u(\cdot) \in [u_{\text{min}},u_{\text{max}}].
\end{array}\right.
$$
with $f(x)$ being the right-hand side of the system previously introduced.
## Problem regularization
Bocop requires $s^-$ to be regularized to a smooth function. We then define, for $x\in \mathbb{R}$ and $k\in \mathbb{N}$, the Hill function
$$
\delta(x_i,\theta_i,k) = \frac{\theta_i^k}{x_i^k + \theta_i^k},
$$
which can approximate $s^-$ for large values of $k$ and, when $k \rightarrow \infty$, meets
$$
\lim_{k \rightarrow \infty} \delta(x_i,\theta_i,k) = \left\{ \begin{array}{ll}
1 & x_i < \theta_i, \\
0 & x_i > \theta_i, \\
1/2 & x_i = \theta_i.
\end{array} \right.
$$
System parameters are fixed to $\gamma_1 = 1.2$, $\gamma_2 = 2$, $\theta_1 = 0.6$, $\theta_2 = 0.4$ and $k_1 = k_2 = 1$, and control bounds are set to $u_{min}=0.5$ and $u_{max}=1.5$. Initial conditions are set to $x_1(0) = 0.8$, $x_2(0)=0.3$ and $x_2(t_f) = 0.7$. The parameter $k$ of the Hill function is set to $k = 200$.
<font color=white>[Thumbnail](thumbnail.png)</font>
%% Cell type:code id: tags:
``` python
!pygmentize problem.cpp
```
%% Output
// +++DRAFT+++ This class implements the OCP functions
// It derives from the generic class bocop3OCPBase
// OCP functions are defined with templates since they will be called
// from both the NLP solver (double arguments) and AD tool (ad_double arguments)
//#pragma once
#include <OCP.h>
template <typename Variable>
void OCP::finalCost(double initial_time, double final_time, const Variable *initial_state, const Variable *final_state, const Variable *parameters, const double *constants, Variable &final_cost)
{
// Minimize final time
final_cost = initial_state[2];
}
template <typename Variable>
void OCP::dynamics(double time, const Variable *state, const Variable *control, const Variable *parameters, const double *constants, Variable *state_dynamics)
{
Variable x1 = state[0];
Variable x2 = state[1];
Variable tf = state[2];
Variable u = control[0];
double g1 = constants[0];
double g2 = constants[1];
double t1 = constants[2];
double t2 = constants[3];
double k = constants[4];
Variable delta1 = pow(t2,k)/(pow(t2,k) + pow(x2,k)); // delta(x2,t2,k)
Variable delta2 = pow(t1,k)/(pow(t1,k) + pow(x1,k)); // delta(x1,t1,k)
state_dynamics[0] = -g1*x1 + u*delta1; // x1dot
state_dynamics[1] = -g2*x2 + u*delta2; // x2dot
state_dynamics[2] = 0; // tf (constant)
// Rescale the state
for (size_t i=0; i<2; i++)
state_dynamics[i] *= tf;
}
template <typename Variable>
void OCP::boundaryConditions(double initial_time, double final_time, const Variable *initial_state, const Variable *final_state, const Variable *parameters, const double *constants, Variable *boundary_conditions)
{
boundary_conditions[0] = initial_state[0]; // x1(0)
boundary_conditions[1] = initial_state[1]; // x2(0)
boundary_conditions[2] = final_state[0]; // x1(tf)
boundary_conditions[3] = final_state[1]; // x2(tf)
}
template <typename Variable>
void OCP::pathConstraints(double time, const Variable *state, const Variable *control, const Variable *parameters, const double *constants, Variable *path_constraints)
{}
void OCP::preProcessing()
{}
// ///////////////////////////////////////////////////////////////////
// explicit template instanciation for template functions, with double and double_ad 
// +++ could be in an included separate file ? 
// but needs to be done for aux functions too ? APPARENTLY NOT !
template void OCP::finalCost<double>(double initial_time, double final_time, const double *initial_state, const double *final_state, const double *parameters, const double *constants, double &final_cost);
template void OCP::dynamics<double>(double time, const double *state, const double *control, const double *parameters, const double *constants, double *state_dynamics);
template void OCP::boundaryConditions<double>(double initial_time, double final_time, const double *initial_state, const double *final_state, const double *parameters, const double *constants, double *boundary_conditions);
template void OCP::pathConstraints<double>(double time, const double *state, const double *control, const double *parameters, const double *constants, double *path_constraints);
template void OCP::finalCost<double_ad>(double initial_time, double final_time, const double_ad *initial_state, const double_ad *final_state, const double_ad *parameters, const double *constants, double_ad &final_cost);
template void OCP::dynamics<double_ad>(double time, const double_ad *state, const double_ad *control, const double_ad *parameters, const double *constants, double_ad *state_dynamics);
template void OCP::boundaryConditions<double_ad>(double initial_time, double final_time, const double_ad *initial_state, const double_ad *final_state, const double_ad *parameters, const double *constants, double_ad *boundary_conditions);
template void OCP::pathConstraints<double_ad>(double time, const double_ad *state, const double_ad *control, const double_ad *parameters, const double *constants, double_ad *path_constraints);
%% Cell type:code id: tags:
``` python
%matplotlib inline
import bocop
```
%% Cell type:code id: tags:
``` python
problem_path = "." # using local problem definition
bocop.build(problem_path)
bocop.build(problem_path, cmake_options = '-DCMAKE_CXX_COMPILER=g++')
```
%% Output
[EXEC] > ['cmake -DCMAKE_BUILD_TYPE=Debug -DPROBLEM_DIR=/opt/ct-gallery/bistable /opt/ct-gallery/env/lib/python3.7/site-packages/bocop']
[EXEC] > ['cmake -DCMAKE_BUILD_TYPE=Debug -DPROBLEM_DIR=/opt/ct-gallery/gallery/examples/bistable -DCPP_FILE=problem.cpp -DCMAKE_CXX_COMPILER=g++ /opt/ct-gallery/env/lib/python3.7/site-packages/bocop']
> -- The C compiler identification is GNU 7.5.0
> -- The CXX compiler identification is GNU 10.2.0
> -- Detecting C compiler ABI info
> -- Detecting C compiler ABI info - done
> -- Check for working C compiler: /opt/ct-gallery/env/bin/x86_64-conda-linux-gnu-cc - skipped
> -- Detecting C compile features
> -- Detecting C compile features - done
> -- Detecting CXX compiler ABI info
> -- Detecting CXX compiler ABI info - done
> -- Check for working CXX compiler: /usr/bin/c++ - skipped
> -- Check for working CXX compiler: /usr/bin/g++ - skipped
> -- Detecting CXX compile features
> -- Detecting CXX compile features - done
> -- Problem path: /opt/ct-gallery/bistable
> -- Problem path: /opt/ct-gallery/gallery/examples/bistable
> -- Using CPPAD found at /opt/ct-gallery/env/include/cppad/..
> -- Using IPOPT found at /opt/ct-gallery/env/lib/libipopt.so
> -- Found Python3: /opt/ct-gallery/env/bin/python3.7 (found version "3.7.8") found components: Interpreter Development Development.Module Development.Embed
> -- Found SWIG: /opt/ct-gallery/env/bin/swig (found suitable version "4.0.2", minimum required is "4")
> -- Build type: Debug
> -- Configuring done
> -- Generating done
> -- Build files have been written to: /opt/ct-gallery/bistable/build
[DONE] > ['cmake -DCMAKE_BUILD_TYPE=Debug -DPROBLEM_DIR=/opt/ct-gallery/bistable /opt/ct-gallery/env/lib/python3.7/site-packages/bocop']
> -- Build files have been written to: /opt/ct-gallery/gallery/examples/bistable/build
[DONE] > ['cmake -DCMAKE_BUILD_TYPE=Debug -DPROBLEM_DIR=/opt/ct-gallery/gallery/examples/bistable -DCPP_FILE=problem.cpp -DCMAKE_CXX_COMPILER=g++ /opt/ct-gallery/env/lib/python3.7/site-packages/bocop']
[EXEC] > make
> Scanning dependencies of target bocop
> [ 3%] Building CXX object src/CMakeFiles/bocop.dir/AD/dOCPCppAD.cpp.o
> [ 7%] Building CXX object src/CMakeFiles/bocop.dir/DOCP/dOCP.cpp.o
> [ 10%] Building CXX object src/CMakeFiles/bocop.dir/DOCP/dODE.cpp.o
> [ 14%] Building CXX object src/CMakeFiles/bocop.dir/DOCP/dControl.cpp.o
> [ 17%] Building CXX object src/CMakeFiles/bocop.dir/DOCP/dState.cpp.o
> [ 21%] Building CXX object src/CMakeFiles/bocop.dir/DOCP/solution.cpp.o
> [ 25%] Building CXX object src/CMakeFiles/bocop.dir/NLP/NLPSolverIpopt.cpp.o
> [ 28%] Building CXX object src/CMakeFiles/bocop.dir/OCP/OCP.cpp.o
> [ 32%] Building CXX object src/CMakeFiles/bocop.dir/tools/tools.cpp.o
> [ 35%] Building CXX object src/CMakeFiles/bocop.dir/tools/tools_interpolation.cpp.o
> [ 39%] Building CXX object src/CMakeFiles/bocop.dir/opt/ct-gallery/bistable/problem.cpp.o
> [ 42%] Linking CXX shared library /opt/ct-gallery/bistable/libbocop.so
> [ 39%] Building CXX object src/CMakeFiles/bocop.dir/opt/ct-gallery/gallery/examples/bistable/problem.cpp.o
> [ 42%] Linking CXX shared library /opt/ct-gallery/gallery/examples/bistable/libbocop.so
> [ 42%] Built target bocop
> Scanning dependencies of target bocopwrapper_swig_compilation
> [ 46%] Swig compile /opt/ct-gallery/env/lib/python3.7/site-packages/bocop/src/bocopwrapper.i for python
> Language subdirectory: python
> Search paths:
> ./
> /opt/ct-gallery/env/include/python3.7m/
> AD/
> DOCP/
> NLP/
> OCP/
> tools/
> /opt/ct-gallery/env/include/cppad/../
> /opt/ct-gallery/env/include/coin/
> /opt/ct-gallery/env/include/python3.7m/
> /opt/ct-gallery/env/lib/python3.7/site-packages/bocop/src/AD/
> /opt/ct-gallery/env/lib/python3.7/site-packages/bocop/src/DOCP/
> /opt/ct-gallery/env/lib/python3.7/site-packages/bocop/src/NLP/
> /opt/ct-gallery/env/lib/python3.7/site-packages/bocop/src/OCP/
> /opt/ct-gallery/env/lib/python3.7/site-packages/bocop/src/tools/
> ./swig_lib/python/
> /opt/ct-gallery/env/share/swig/4.0.2/python/
> ./swig_lib/
> /opt/ct-gallery/env/share/swig/4.0.2/
> Preprocessing...
> Starting language-specific parse...
> Processing types...
> C++ analysis...
> Processing nested classes...
> Generating wrappers...
> [ 46%] Built target bocopwrapper_swig_compilation
> Scanning dependencies of target bocopwrapper
> [ 50%] Building CXX object src/CMakeFiles/bocopwrapper.dir/CMakeFiles/bocopwrapper.dir/bocopwrapperPYTHON_wrap.cxx.o
> [ 53%] Linking CXX shared module /opt/ct-gallery/bistable/_bocopwrapper.so
> -- Moving python modules to /opt/ct-gallery/bistable
> [ 53%] Linking CXX shared module /opt/ct-gallery/gallery/examples/bistable/_bocopwrapper.so
> -- Moving python modules to /opt/ct-gallery/gallery/examples/bistable
> [ 53%] Built target bocopwrapper
> Scanning dependencies of target bocopApp
> [ 57%] Building CXX object src/CMakeFiles/bocopApp.dir/AD/dOCPCppAD.cpp.o
> [ 60%] Building CXX object src/CMakeFiles/bocopApp.dir/DOCP/dOCP.cpp.o
> [ 64%] Building CXX object src/CMakeFiles/bocopApp.dir/DOCP/dODE.cpp.o
> [ 67%] Building CXX object src/CMakeFiles/bocopApp.dir/DOCP/dControl.cpp.o
> [ 71%] Building CXX object src/CMakeFiles/bocopApp.dir/DOCP/dState.cpp.o
> [ 75%] Building CXX object src/CMakeFiles/bocopApp.dir/DOCP/solution.cpp.o
> [ 78%] Building CXX object src/CMakeFiles/bocopApp.dir/NLP/NLPSolverIpopt.cpp.o
> [ 82%] Building CXX object src/CMakeFiles/bocopApp.dir/OCP/OCP.cpp.o
> [ 85%] Building CXX object src/CMakeFiles/bocopApp.dir/tools/tools.cpp.o
> [ 89%] Building CXX object src/CMakeFiles/bocopApp.dir/tools/tools_interpolation.cpp.o
> [ 92%] Building CXX object src/CMakeFiles/bocopApp.dir/opt/ct-gallery/bistable/problem.cpp.o
> [ 92%] Building CXX object src/CMakeFiles/bocopApp.dir/opt/ct-gallery/gallery/examples/bistable/problem.cpp.o
> [ 96%] Building CXX object src/CMakeFiles/bocopApp.dir/main.cpp.o
> [100%] Linking CXX executable /opt/ct-gallery/bistable/bocopApp
> [100%] Linking CXX executable /opt/ct-gallery/gallery/examples/bistable/bocopApp
> [100%] Built target bocopApp
[DONE] > make
Done
0
%% Cell type:markdown id: tags:
State 0 and State 1 correspond to $x_1$ and $x_2$ respectively. State 2 represents the final time to minimize.
%% Cell type:code id: tags:
``` python
import matplotlib.pyplot as plt
bocop.run(problem_path)
plt.show()
```
%% Output
Done
%% Cell type:code id: tags:
``` python
import numpy as np
solution = bocop.readSolution(problem_path + "/problem.sol")
tf = solution.state[2][0]; t = solution.time_steps*tf
x1 = solution.state[0]; x2 = solution.state[1]
u = solution.control[0]
umin = 0.5
umax = 1.5
g1 = 1.4
g2 = 1.8
theta1 = 0.6
theta2 = 0.4
k1 = 1
k2 = 1
x1max = 0.9
x2max = 0.8
plt.figure(0, figsize=(10,4))
plt.subplots_adjust(hspace=0.4, wspace=0.4)
plt.subplot(121)
ax = plt.gca()
plt.xlabel('$x_1$'); plt.ylabel('$x_2$')
ax.set_ylim(0.,x1max); ax.set_xlim(0.,x2max)
xv = np.arange(0,theta1,0.01)
sumax = k2/g2*umax - (k2/g2*umax - theta2)*((k1/g1*umax - xv)/((k1/g1*umax - theta1)))**(g2/g1)
plt.plot(x1,x2,label='',zorder=3)
plt.plot(xv,sumax,'--', label='$(S_{u_{max}})$',linewidth=1, color='grey')
plt.legend()
plt.xlim([0,theta1*1.8])
ax.text(x1[0]+0.04, x2[0]-0.05, '$(x_1^0, x_2^0)$')
plt.scatter([x1[0],x1[-1]],[x2[0],x2[-1]], zorder=3, s=20)
plt.xticks([0, theta1], [0,'$\\theta_1$'])
plt.yticks([0, theta2, x2[-1]], [0,'$\\theta_2$','$x_2^f$'])
plt.grid(linestyle='dotted')
plt.subplot(122)
plt.plot(t[1:],u)
plt.xlabel('Time'); plt.ylabel('$u$')
plt.yticks([umin,1,umax], ['$u_{min}$', 1, '$u_{max}$'])
plt.grid(linestyle='dotted')
plt.show()
print("Final time is " + str(tf))
```
%% Output
Loading solution: ./problem.sol
Final time is 1.88756560904673
%% Cell type:code id: tags:
``` python
print("Bocop returns status {} with objective {:2.4g} and constraint violation {:2.4g}".format(solution.status,solution.objective,solution.constraints))
p0 = []
for i in range(solution.dim_state):
p0.append(solution.costate[i][0])
print("Costate at first time stage (t0+h/2): ",p0)
print("Multipliers for initial conditions: ",solution.boundarycond_multipliers[0:solution.dim_state])
```
%% Output
Bocop returns status 0 with objective 1.888 and constraint violation 1.486e-14
Costate at first time stage (t0+h/2): [-2.21241861069235, 0.680215961376846, 0.99503700918767]
Multipliers for initial conditions: [-2.18337787e+00 6.68757718e-01 1.38521948e-12]
......
......@@ -20,6 +20,7 @@ constant.3 0.4
constant.4 200
# Time discretisation
ode.discretization midpoint_implicit
time.steps 200
# Bounds for constraints
......@@ -43,3 +44,13 @@ state.0.init 0.1
state.1.init 0.1
state.2.init 3
control.0.init 0.5
# Ipopt
ipoptIntOption.print_level 5
ipoptIntOption.max_iter 1000
ipoptStrOption.mu_strategy adaptive
ipoptNumOption.tol 1e-12
# Misc
ad.retape 0
%% Cell type:markdown id: tags:
# Seasonal coffee leaf rust propagation
## Impulsive model definition
We consider the impulsive mathematical model for the propagation of the CLR in the coffee plantation. The latter is described by the dynamics of fungus during the production season (i.e. rainy season for some countries) through ordinary differential equations; and during the non-production period (i.e. dry season for some countries) represented by the impulsive effect. The model describes the evolution of susceptible branches $S$, infected branches $I$, urediniospores $U$ and berries $B$, and is given by dynamics
$$
\left\{
\begin{aligned}
&\dot{S}=\Lambda -\frac{\omega\nu U}{N}S -\mu S, \;t \neq nT; \\
&\dot{I}=\frac{\omega\nu U}{N}S -(\mu +d)I, \; t \neq nT; \\
&\dot{U}=\gamma I -(\nu+\mu_{U})U-v, \; t \neq nT; \\
&\dot B = \delta_S S + \delta_I I - \mu_B B, \; t \neq nT;\\
&S(nT^+)=\varphi_S S(nT); \\
&I(nT^+)=\varphi_I I(nT); \\
&U(nT^+)=\varphi_U U(nT);\\
&B(nT^+)=0;
\end{aligned}
\right.
$$
where $N(t)=S(t)+I(t)$ represents the total number of branches at time $t$, and $v$ is the inoculative biological control (e.g. a mycoparasite capable of consuming the urediniospores). Initial conditions are fixed $S(0)=S_0>0$, $I(0)=I_0 \geq 0$, $U(0)=U_0 \geq 0$, $B(0) = B_0 > 0$. The cost function to maximize is the profit for a time period $[0,kT]$, given by
$$
J(v) = \alpha_B \sum_{j=1}^k B(jT) - \alpha_v \int_0^{kT} v \, dt = \int_{0}^{kT} \Big( \alpha_B(\delta_S S + \delta_I I - \mu_B B) - \alpha_v v \Big) \, dt,
$$
where $\alpha_B$ and $\alpha_v$ are the market prices per coffee grain and per spore, respectively.
## Model in continuous form
The latter problem can be rewritten in the continuous form by defining $k$ set of states $x_n = (S_n,I_n,U_n,B_n)$ and $k$ controls $v_n(t)$, with dynamics defined for $t \in [0,T]$ as
$$
\left\{
\begin{aligned}
&\dot{S_n}=\Lambda -\omega\nu U_n\frac{S_n}{S_n + I_n} -\mu S_n, \\
&\dot{I_n}=\omega\nu U_n\frac{S_n}{S_n + I_n} -(\mu +d)I_n, \\
&\dot{U_n}=\gamma I_n -(\nu+\mu_{U})U_n - v_n, \\
&\dot B_n = \delta_S S_n + \delta_I I_n - \mu_B B_n,
\end{aligned}
\right.
$$
for $n = 1, \dots, k$, with boundary conditions
$$
\left\{
\begin{aligned}
& S_1(0)=S_0>0, \; I_1(0)=I_0 \geq 0, \\ & U_1(0)=U_0 \geq 0, \; B_1(0) = B_0 > 0, \\
& S_{n+1}(0)=\varphi_s S_n(T), \; I_{n+1}(0)=\varphi_I I_n(T), \\
& U_{n+1}(0)=\varphi_U U_n(T), \; B_{n+1}(0)= 0,
\end{aligned}
\right.
$$
for $n = 1, \dots, k-1$. Then, the cost function becomes
$$
J_c(v_1, \dots, v_k) = \alpha_B \sum_{n=0}^k B_n(T) - \alpha_v \sum_{n=0}^k \int_0^{T} v_n \, dt.
$$
<font color=white>[Thumbnail](thumbnail.png)</font>
%% Cell type:code id: tags:
``` python
!pygmentize problem.cpp
```
%% Output
// +++DRAFT+++ This class implements the OCP functions
// It derives from the generic class bocop3OCPBase
// OCP functions are defined with templates since they will be called
// from both the NLP solver (double arguments) and AD tool (ad_double arguments)
//#pragma once
#include <OCP.h>
template <typename Variable>
void OCP::finalCost(double initial_time, double final_time, const Variable *initial_state, const Variable *final_state, const Variable *parameters, const double *constants, Variable &final_cost)
{
int t;
final_cost = 0;
double prix = constants[13];
for( t = 0; t < 4; t = t + 1 ){
int q = t*5;
final_cost = final_cost - final_state[3 + q] + prix*final_state[4 + q];
}
}
template <typename Variable>
void OCP::dynamics(double time, const Variable *state, const Variable *control, const Variable *parameters, const double *constants, Variable *state_dynamics)
{
int t;
double Lam = constants[0];
double ome = constants[1];
double gam = constants[2];
double nu = constants[3];
double mu = constants[4];
double muU = constants[5];
double d = constants[6];
double deltaS = constants[10];
double deltaI = constants[11];
double muB = constants[12];
double ifcontrol = constants[14];
for( t = 0; t < 4; t = t + 1 ){
int q = t*5;
Variable S = state[0 + q];
Variable I = state[1 + q];
Variable U = state[2 + q];
Variable B = state[3 + q];
Variable v = control[t];
Variable lambda = (ome*nu*U)/(S+I);
state_dynamics[0 + q] = Lam - lambda*S - mu*S;
state_dynamics[1 + q] = lambda*S - (mu+d)*I;
state_dynamics[2 + q] = gam*I - (nu+muU)*U - v*ifcontrol;
state_dynamics[3 + q] = deltaS*S + deltaI*I - muB*B;
state_dynamics[4 + q] = v;
}
}
template <typename Variable>
void OCP::boundaryConditions(double initial_time, double final_time, const Variable *initial_state, const Variable *final_state, const Variable *parameters, const double *constants, Variable *boundary_conditions)
{
boundary_conditions[0] = initial_state[0];
boundary_conditions[1] = initial_state[1];
boundary_conditions[2] = initial_state[2];
boundary_conditions[3] = initial_state[3];
boundary_conditions[4] = initial_state[4];
double phiS = constants[7];
double phiI = constants[8];
double phiU = constants[9];
int t;
for( t = 0; t < 3; t = t + 1 ){
int q = t*5;
boundary_conditions[5 + q] = initial_state[5 + q] - final_state[0 + q]*phiS;
boundary_conditions[6 + q] = initial_state[6 + q] - final_state[1 + q]*phiI;
boundary_conditions[7 + q] = initial_state[7 + q] - final_state[2 + q]*phiU;
boundary_conditions[8 + q] = initial_state[8 + q] - final_state[3 + q]*0;
boundary_conditions[9 + q] = initial_state[9 + q] - final_state[4 + q];
}
}
template <typename Variable>
void OCP::pathConstraints(double time, const Variable *state, const Variable *control, const Variable *parameters, const double *constants, Variable *path_constraints)
{}
void OCP::preProcessing()
{}
// ///////////////////////////////////////////////////////////////////
// explicit template instanciation for template functions, with double and double_ad 
// +++ could be in an included separate file ? 
// but needs to be done for aux functions too ? APPARENTLY NOT !
template void OCP::finalCost<double>(double initial_time, double final_time, const double *initial_state, const double *final_state, const double *parameters, const double *constants, double &final_cost);
template void OCP::dynamics<double>(double time, const double *state, const double *control, const double *parameters, const double *constants, double *state_dynamics);
template void OCP::boundaryConditions<double>(double initial_time, double final_time, const double *initial_state, const double *final_state, const double *parameters, const double *constants, double *boundary_conditions);
template void OCP::pathConstraints<double>(double time, const double *state, const double *control, const double *parameters, const double *constants, double *path_constraints);
template void OCP::finalCost<double_ad>(double initial_time, double final_time, const double_ad *initial_state, const double_ad *final_state, const double_ad *parameters, const double *constants, double_ad &final_cost);
template void OCP::dynamics<double_ad>(double time, const double_ad *state, const double_ad *control, const double_ad *parameters, const double *constants, double_ad *state_dynamics);
template void OCP::boundaryConditions<double_ad>(double initial_time, double final_time, const double_ad *initial_state, const double_ad *final_state, const double_ad *parameters, const double *constants, double_ad *boundary_conditions);
template void OCP::pathConstraints<double_ad>(double time, const double_ad *state, const double_ad *control, const double_ad *parameters, const double *constants, double_ad *path_constraints);
%% Cell type:code id: tags:
``` python
%matplotlib inline
import bocop
import os,shutil
import matplotlib.pyplot as plt
```
%% Cell type:code id: tags:
``` python
problem_path = "." # using local problem definition
bocop.build(problem_path)
bocop.build(problem_path, cmake_options = '-DCMAKE_CXX_COMPILER=g++')
```
%% Output
[EXEC] > ['cmake -DCMAKE_BUILD_TYPE=Debug -DPROBLEM_DIR=/opt/ct-gallery/gallery/examples/clr /opt/ct-gallery/env/lib/python3.7/site-packages/bocop']
[EXEC] > ['cmake -DCMAKE_BUILD_TYPE=Debug -DPROBLEM_DIR=/opt/ct-gallery/gallery/examples/clr -DCPP_FILE=problem.cpp -DCMAKE_CXX_COMPILER=g++ /opt/ct-gallery/env/lib/python3.7/site-packages/bocop']
> -- The C compiler identification is GNU 7.5.0
> -- The CXX compiler identification is GNU 10.2.0
> -- Detecting C compiler ABI info
> -- Detecting C compiler ABI info - done
> -- Check for working C compiler: /opt/ct-gallery/env/bin/x86_64-conda-linux-gnu-cc - skipped
> -- Detecting C compile features
> -- Detecting C compile features - done
> -- Detecting CXX compiler ABI info
> -- Detecting CXX compiler ABI info - done
> -- Check for working CXX compiler: /usr/bin/c++ - skipped
> -- Check for working CXX compiler: /usr/bin/g++ - skipped
> -- Detecting CXX compile features
> -- Detecting CXX compile features - done
> -- Problem path: /opt/ct-gallery/gallery/examples/clr
> -- Using CPPAD found at /opt/ct-gallery/env/include/cppad/..
> -- Using IPOPT found at /opt/ct-gallery/env/lib/libipopt.so
> -- Found Python3: /opt/ct-gallery/env/bin/python3.7 (found version "3.7.8") found components: Interpreter Development Development.Module Development.Embed
> -- Found SWIG: /opt/ct-gallery/env/bin/swig (found suitable version "4.0.2", minimum required is "4")
> -- Build type: Debug
> -- Configuring done
> -- Generating done
> -- Build files have been written to: /opt/ct-gallery/gallery/examples/clr/build
[DONE] > ['cmake -DCMAKE_BUILD_TYPE=Debug -DPROBLEM_DIR=/opt/ct-gallery/gallery/examples/clr /opt/ct-gallery/env/lib/python3.7/site-packages/bocop']
[DONE] > ['cmake -DCMAKE_BUILD_TYPE=Debug -DPROBLEM_DIR=/opt/ct-gallery/gallery/examples/clr -DCPP_FILE=problem.cpp -DCMAKE_CXX_COMPILER=g++ /opt/ct-gallery/env/lib/python3.7/site-packages/bocop']
[EXEC] > make
> Scanning dependencies of target bocop
> [ 3%] Building CXX object src/CMakeFiles/bocop.dir/AD/dOCPCppAD.cpp.o
> [ 7%] Building CXX object src/CMakeFiles/bocop.dir/DOCP/dOCP.cpp.o
> [ 10%] Building CXX object src/CMakeFiles/bocop.dir/DOCP/dODE.cpp.o
> [ 14%] Building CXX object src/CMakeFiles/bocop.dir/DOCP/dControl.cpp.o
> [ 17%] Building CXX object src/CMakeFiles/bocop.dir/DOCP/dState.cpp.o
> [ 21%] Building CXX object src/CMakeFiles/bocop.dir/DOCP/solution.cpp.o
> [ 25%] Building CXX object src/CMakeFiles/bocop.dir/NLP/NLPSolverIpopt.cpp.o
> [ 28%] Building CXX object src/CMakeFiles/bocop.dir/OCP/OCP.cpp.o
> [ 32%] Building CXX object src/CMakeFiles/bocop.dir/tools/tools.cpp.o
> [ 35%] Building CXX object src/CMakeFiles/bocop.dir/tools/tools_interpolation.cpp.o
> [ 39%] Building CXX object src/CMakeFiles/bocop.dir/opt/ct-gallery/gallery/examples/clr/problem.cpp.o
> [ 42%] Linking CXX shared library /opt/ct-gallery/gallery/examples/clr/libbocop.so
> [ 42%] Built target bocop
> Scanning dependencies of target bocopwrapper_swig_compilation
> [ 46%] Swig compile /opt/ct-gallery/env/lib/python3.7/site-packages/bocop/src/bocopwrapper.i for python
> Language subdirectory: python
> Search paths:
> ./
> /opt/ct-gallery/env/include/python3.7m/
> AD/
> DOCP/
> NLP/
> OCP/
> tools/
> /opt/ct-gallery/env/include/cppad/../
> /opt/ct-gallery/env/include/coin/
> /opt/ct-gallery/env/include/python3.7m/
> /opt/ct-gallery/env/lib/python3.7/site-packages/bocop/src/AD/
> /opt/ct-gallery/env/lib/python3.7/site-packages/bocop/src/DOCP/
> /opt/ct-gallery/env/lib/python3.7/site-packages/bocop/src/NLP/
> /opt/ct-gallery/env/lib/python3.7/site-packages/bocop/src/OCP/
> /opt/ct-gallery/env/lib/python3.7/site-packages/bocop/src/tools/
> ./swig_lib/python/