L'accumulation de nutriments (comme le phosphore ou l'azote) dans l'eau d'un lac peut amener un changement d'état qui s'accompagne de prolifération d'algues, de dégradation de la qualité de l'eau et de la biodiversité, éventuellement de bloom bactérien : c'est l'eutrophisation. Le problème du lac et des exploitations riveraines consiste à déterminer s'il est possible de concilier la pratique d'une activité qui apporte des nutriments et la conservation du lac dans un état souhaitable (oligotrophe, par opposition à eutrophe).
Ce problème est décrit en détail dans :
S. Martin. The cost of restoration as a way of defining resilience: a viability approach applied to a model of lake eutrophication. Ecol. Soc. http://www.ecologyandsociety.org/vol9/iss2/art8
On souhaite que les apports $L$ soient supérieurs à un seuil $L_{min}$, pour tenir compte des besoins de l'activité des exploitations riveraines ; et que la concentration du phosphore total $P$ soit inférieure à un seuil $P_{max}$ pour conserver le lac oligotrophe. Ces états souhaitables constituent l'ensemble de contraintes $K=[L_{min}, +\infty[ \times [0,P_{max}]$.
L'évolution de la concentration du phosphore total dans le lac est modélisée par une pseudo-sygmoïde :
$$\frac {dP} {dt}=-bP(t)+L(t)+r\frac {P(t)^{q}} {m^{q} + P(t)^{q}} \qquad$$
On suppose que l'évolution des apports de phosphore peut être contrôlée (par des unités de dépollution, la mise en place de zones humides, le changement de pratique agricoles ou industrielles, etc.) et on modélise ces contrôles par une grandeur unique $u\in U=[u_min,u_max]$. La dynamique des apports est modélisée par :
$$\frac {dL} {dt}=u \in [- u_{min}, u_{max}] \qquad$$
Le problème de viabilité est donc défini par :
\begin{equation} \label{eq:systemlac}
(P)\left\{
\begin{array}{l}
\frac{dL}{dt}=u\in U=\left[ u_{min},u_{max}\right] \\
\frac{dP}{dt}=-b P(t) + L(t) +r\frac{P(t)^{q}}{m^{q} + P(t)^{q}} \\
\left(L(t),P(t)\right) \in K=[L_{min}, +\infty[ \times [0,P_{max}]
\end{array}
\right.
\end{equation}
Le noyau de viabilité peut être obtenu par le calcul d'une courbe intégrale (voir page 22 de https://arxiv.org/pdf/2107.02684 ). La figure suivante montre le résultat pour les paramètres suivants : $b=1,95$ an$^{-1}$ ; $q=1,9$ ; $m=19,44\ \mu gl^{-1}$; $r=72,22\ \mu gl^{-1}$ an$^{-1}$; $L_{min}=1,25\ \mu gl^{-1}$ ; $P_{max}=17,39\ \mu gl^{-1}$ ; $|u_{min}|=u_{max}=3,15$.
Une approximation peut être calculée directement avec des logiciels de calcul de noyau comma ViabLab.
Avec une précédente version de VIABLAB (2.1), voici un extrait du code.
/*! \file lac_classiqueb07r1.h
string prefix="lac_classique_b08r1_test" ;
/***********************************************************
* Some model specific parameters can be defined here
*************************************************************/
int axisAcc = 5000 ; //nb de points par axe
double umax = 0.09;
double m_lac = 1.0 ;
double q = 8 ; // change sediment_p definition if q changes
double b = 0.8;
double r = 1.0;
double L_min = 0.1;
double L_max = 1.0;
double P_min = 0.0;
double P_max = 1.4;
////////////////////////////////////
/// variable d'états
////////////////////////////////////
/*! \var dim
* \brief State dimension
*/
const int dim=2;
/*! \var dicret_type
* \brief the type of discretization scheme used for the dynamics
* EE or 1 = Euler Explicit scheme
* EI or 2 = Euler Implicit
* RK2I or 3 = RK2 Implicit (RK2 for -F)
* RK2E or 4 = RK2 Explicit (RK2 for F)
*/
const int discret_type=RK4E;
const int dynType=CC;
double STATE_MIN[dim]={L_min,P_min}; // L,P
double STATE_MAX[dim]={L_max,P_max}; // L,P
unsigned long long int nbPointsState[dim] = {axisAcc,axisAcc};
/* *****************************************
* Definition of constraints and target
*************************************************** */
inline double constraintsX( double * x ){
return 1.0;
}
/* *****************************************
* Definition of constraints and target
*************************************************** */
const int dimc=1;
double CONTROL_MIN[dimc]={-umax};
double CONTROL_MAX[dimc]={umax};
unsigned long long int nbPointsControl[dimc] = {3};
inline double constraintsXU( double * x, double * u ) {
return 1.0;
}
bool globalDeltaT=false;
/*
* Direction le long de laquelle l'ensemble sera représenté par des suites de bits
*/
unsigned long long int dirTramage =1;
int gridMethod=BS;
int setType=VIAB;
int refine= 0;
int periodic[dim]={0};
int computeSet=1;
int saveBoundary=0;
/*!
* \var computeLC : indicates which method to use to compute Lipschitz constant
* 0= analytic global value used
* 1= local calculation using the jacobian matrix
* 2= local calculation using the finite differences
*/
const int computeLC=2;
/*!
* \var computeM : indicates which method to use to compute the bound for dynamics
* 0= analytic global value used
* 1= local calculation using the localDynBounds function
* 2= local calculation using explicit maximization on controls
*/
const int computeM=2;
// ///////////////////////////////////
//Fonctions utiles pour la dynamique
// //////////////////////////////////
double sediment_p(double p) {return p*p*p*p*p*p*p*p/(m_lac*m_lac*m_lac*m_lac*m_lac*m_lac*m_lac*m_lac + p*p*p*p*p*p*p*p) ;} ;
////////////////////////////////////
/// fonction Dynamique
////////////////////////////////////
void dynamics(double * x, double *u, double * image)
{
image[0]= u[0];
image[1]= x[0] - b* x[1] + r* sediment_p(x[1]) ;
}
Le contrôle est discrétisé sur 3 valeurs étant donné les spécificités du problème. Les valeurs des paramètres sont : $b=0.8$ an$^{-1}$ ; $q=8$ ; $m=1\ \mu gl^{-1}$; $r=1\ \mu gl^{-1}$ an$^{-1}$; $L_{min}=0.1\ \mu gl^{-1}$ ; $P_{max}=1.4\ \mu gl^{-1}$ ; $|u_{min}|=u_{max}=0.09$.
La discrétisation est faite sur 5000 points / axe. L'approximation calculée par ViabLab est montrée sur la figure suivante. En rouge les seuils des contraintes. On voit que l'approximation du noyau est bien faite par l'extérieur.
