Notebook source code:
examples/clr/clr.ipynb
Run the notebook yourself on binder
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
\[\begin{split}\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.\end{split}\]
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
\[\begin{split}\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.\end{split}\]
for \(n = 1, \dots, k\), with boundary conditions
\[\begin{split}\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.\end{split}\]
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.\]
[1]:
!pygmentize problem.cpp
// +++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);
[2]:
%matplotlib inline
import bocop
import os,shutil
import matplotlib.pyplot as plt
[3]:
problem_path = "." # using local problem definition
bocop.build(problem_path, cmake_options = '-DCMAKE_CXX_COMPILER=g++')
[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/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 -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/
> /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/gallery/examples/clr/_bocopwrapper.so
> -- Moving python modules to /opt/ct-gallery/gallery/examples/clr
> [ 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/gallery/examples/clr/problem.cpp.o
> [ 96%] Building CXX object src/CMakeFiles/bocopApp.dir/main.cpp.o
> [100%] Linking CXX executable /opt/ct-gallery/gallery/examples/clr/bocopApp
> [100%] Built target bocopApp
[DONE] > make
Done
[3]:
0
[4]:
files = ['control', 'nocontrol']
for file in files:
shutil.copyfile(problem_path + "/" + file + '.def', problem_path + "/" + 'problem.def')
bocop.run(problem_path, graph=0)
os.rename("problem.sol",file + '.sol')
Executing bocop ...
Done
Executing bocop ...
Done
[5]:
import numpy as np
from matplotlib import ticker
formatter = ticker.ScalarFormatter(useMathText=True)
formatter.set_scientific(True)
formatter.set_powerlimits((-1,1))
label = ['Optimal control','No treatment']
files = ['control', 'nocontrol']
v = [[],[]]
q = 4 # Number of periods simulated
states = 5 # Number of states
names = ['Biological control $v$','$S$','$I$','$U$','$B$','$c$']
plt.figure(figsize=(6,7))
plt.subplot(321)
plt.subplots_adjust(hspace=0.7, wspace=0.3)
for o in range(0,len(files)):
solution = bocop.readSolution(problem_path + "/" + files[o] + ".sol")
t = solution.time_steps; period = t[-1]; tt = []
# Construct the real time interval
for i in range(0,q):
tt.append(t + i*period)
t = np.hstack(tt)
t0 = t[0]; t1 = t[-1]
tmp = []
for i in range(0,q):
tmp.append(solution.control[i])
tmp.append(0)
v[o].append(np.hstack(tmp))
for i in range(0,q):
tmp = []
for p in range(0,states-1):
tmp.append(solution.state[p*states + i])
v[o].append(np.hstack(tmp))
for i in range(0,5):
if i == 0:
plt.subplot(311) # Control plot is bigger
else:
plt.subplot(322 + i)
ax = plt.gca()
plt.title(names[i])
plt.xlabel('Time (Day)')
ax.set_xlim(t0,t1)
plt.plot([t0,t1], [np.max(v[o][0]), np.max(v[o][0])], linestyle=':', color='grey', linewidth=1)
plt.plot([t0,t1], [np.min(v[o][0]), np.min(v[o][0])], linestyle=':', color='grey', linewidth=1)
for o in range(0,len(files)):
curve = v[o][i]
plt.plot(t,curve,linewidth=1,label=label[o])
for j in range(0,q):
plt.axvline(x=period*j, linestyle=':', color='grey', linewidth = 1)
lim = ax.get_ylim()
ax.set_ylim((0,lim[1]))
ax.yaxis.set_major_formatter(formatter)
if i == 0: plt.legend(bbox_to_anchor=(0, 1.3, 1, 0), loc="lower center", ncol=2)
plt.show()
Loading solution: ./control.sol
Loading solution: ./nocontrol.sol
