Débutant

Introduction

Les maladies à transmission vectorielle sont causées par divers types d'agents pathogènes, notamment des virus et des bactéries, transmis par des vecteurs tels que les moustiques, les tiques, entre autres. Selon l'Organisation mondiale de la Santé, ces maladies entraînent plus de 700 000 décès chaque année.

Ross a été le premier à tenter de modéliser ce phénomène pour proposer des stratégies de lutte contre les épidémies de moaldies à transmission vectorielle. Le modèle qui porte son nom a été publié en 1911

Ross, R. (1911). Some Quantitative Studies in Epidemiology. 87, 466 - 467.

et sera à l'origine de nombreux autres modèles. Ross a envisagé un modèle simple du paludisme comportant humains et moustiques qui ne peuvent appartenir qu'à deux catégories : infectés et non infectés. Les naissances, compensent les décès, les nombres d'humains et de moustiques sont constants.

Dans les modèles inspirés par Ross, ceux, connus sous le nom de Ross-MacDonald, prennent comme variables les proportions d'infectés par rapport à la population totale pour s'affranchir de l'hypothèse des populations constantes, par exemple la proportion d'humains infectés, $h$ et la proportion de moustiques infectés, $m$.
Les populations d'hôtes et de vecteurs sont en mélange homogène. Chaque population est subdivisée en classes épidémiologiques, dans le cas le plus simple, deux classes, les individus sains (susceptibles) et infectieux (infectés). Les vecteurs, par exemple les moustiques, piquent à un taux $\alpha$ (nombre quotidien de piqûres par moustique sur humain, par exemple). Si $p_h$ désigne la probabilité de transmission de l'infection aux hôtes par piqûre, et $p_m$ la probabilité d'infection du moustique lors d'une piqûre sur un hôte infectieux, alors le taux d'infection des hôtes sains est donné par $\alpha p_h \xi m(1-h)$, tandis que le taux d'infection des vecteurs sains est donné par $\alpha p_m h (1-m)$. Ces formes fonctionnelles des taux d'infection sont caractéristiques de tous les modèles de type Ross-Macdonald. 

Les équations qui régissent l'évolution des proportions d'infectés dans les deux populations sont ainsi :

$$
\begin{array}{ll}
\frac{dm}{dt} &= \alpha p_m h (1-m) − \delta m \\
\frac{dh}{dt} &= \alpha p_h \xi m(1-h) −\gamma h 
\end{array}
$$

avec $\gamma$, le taux de guérison chez l'humain et $\delta$, le taux de mortalité (naturelle) des moustiques.

Le contrôle viable d'un modèle d'épidémiologie

Dans l'article

Michel de Lara, Lilian Sofia Sepulveda Salcedo. Viable Control of an Epidemiological Model. Mathematical Biosciences, 2016, 280, pp.24-37. ⟨10.1016/j.mbs.2016.07.010⟩. ⟨hal-01208484v2⟩

les auteurs proposent de s'intéresser aux stratégies qui visent à limiter le nombre d'humains infectés au moment du pic de l'épidémie.

Dynamiques et contraintes

Les auteurs considèrent une version contrôlée du modèle épidémiologique dynamique de Ross-Macdonald : la mortalité des vecteurs représente la variable de commande. 
En effet des actions de fumigation par exemple peuvent conduire à un taux de mortalité des moustiques bien plus élevé que le taux naturel $\delta$. Ainsi, sont définies les fonctions $u(.)$ qui contrôlent la population de moustiques en affectant leur mortalité en plus de leur mortalité naturelle. Ainsi, $u(t)\in [\delta;\bar{u}]$ pour tout $t$ avec $\bar{u}$ représentant l'efficacité maximale de ces méthodes de destruction des moustiques. 

La variable $\bar{H}\in [0;1]$ représente le taux maximal d'infectés tolérable dans la population humaine. Ainsi, éviter qu'au pic de l'épidémie, la proportion d'humains infectés soit supérieur à $\bar{H}$, c'est trouver des politiques de santé publique qui conduisent à ce que la proportion $h(t)$ d'humains infectés reste toujours inférieure au seuil maximal d'infection $\bar{H}$, c'est à dire qui respecte la contrainte (dans la terminologie de la théorie de la viabilité) :

$$
h(t)\leq \bar{H},\, \forall t\geq 0.
$$

La dynamique du système est ainsi décrite par :

$$
\left\{\begin{array}{ll}
\frac{dm}{dt} &= \alpha p_m h (1-m) − u m \\
\frac{dh}{dt} &= \alpha p_h \xi m(1-h) −\gamma h\\
u&\in   [\delta;\bar{u}]
\end{array}\right.
$$

Et l'ensemble des contraintes s'écrit :

$$
K(\bar{H}):=[0,1]\times [0;\bar{H}].
$$

Les proportions de moustiques infectés et d'humains infectés influencent sans doute la capacité à éviter la survenue d'un pic d'humains infectés supérieur à $\bar{H}$.
Résoudre le problème de la possibilité de préserver la population humaine d'un pic d'infectés trop élevé, c'est trouver l'ensemble des états $(m,h)$ à partir desquels il est possible de mettre en place des stratégies de contrôle de la population de moustiques qui contiennent la proportion d'infectés en dessous de ce seuil, c'est à dire le noyau de viabilité dans la terminologie de la théorie de la viabilité. 

Résultats

Théoriques

Pour des valeurs de paramètres $\alpha, \xi, \gamma, p_m$, $p_h$ et $\delta$ fixées, les valeurs des couples intensité maximale du pic, $\bar{H}$ et efficacité maximale de la fumigation $\bar{u}$ peuvent conduire à trois types de cas :

• le cas très favorable, le pic peut toujours être restreint à la valeur tolérable dès que celle-ci n'a pas été franchie, le noyau de viabilité est égale à l'ensemble de contraintes. Posons $A_h =\alpha p_h \xi $ et $A_m =\alpha p_m $, ceci se produit lorsque  : 
  $$
  \bar{H}\geq \frac{A_h}{A_h +\gamma}.
  $$
• le cas désespéré, dès lors qu'une infection survient, le pic d'humains infectés sera nécessairement supérieur à $\bar{H}$ même en opérant des fumigations d'efficacité maximale, le noyau de viabilité est réduit à l'origine. Ceci se produit lorsque  
  $$
  A-m(A-h + \gamma)\bar{H} + \gamma\bar{u} < A_mA_h.
  $$
• dans le dernier cas, le noyau est ni vide ni plein, la possibilité d'éviter un pic trop élevé dépendra de la situation initiale.
  Le noyau de viabilité peut être obtenu par le calcul d'une courbe intégrale.
  En effet, sa frontière est définie par une courbe fermée formée de sous-ensembles de la frontière de l'ensemble de contraintes $K$ et d'une courbe intégrale issue du point $(\bar{M},\bar{H})$ avec $\frac{dh}{dt}(\bar{M},\bar{H}) =0$ qui suit la dynamique inverse : 
  \begin{equation}
  \left\{
  \begin{array}{l}
  \frac{dm}{dt} &= -(\alpha p_m h (1-m) − u m) \\
  \frac{dh}{dt} &= -(\alpha p_h \xi m(1-h) −\gamma h)\\
  \end{array}
  \right.
  \end{equation}

 Numériques

Les auteurs ont estimé les valeurs des paramètres de la dynamique dans la cas de l'épidémie de dengue de 2013 à Cali e, Colombie : $\alpha =0.3600\;day^{-1}$, $p_m=0.2128$, $p_h = 0.1990$, $\xi= 1.0087 $ et $\delta = 0.0333 \;day^{-1}$, $\gamma =0.1 \;day^{-1}$.

Noyaux de viabilité

En prenant ces valeurs de paramètres, nous ne parvenons pas à reproduire les résultats numériques obtenus par les auteurs de l'article, sans doute une mauvaise interprétation de notre part à un moment donné.

$\bar{u} = 0.05$ $\bar{H} = 0.2$
La solution au problème de viabilité (3)-(4)-(8)-(10) dépend, en fin de compte, des conditions initiales relatives aux proportions de moustiques infectés et d'humains infectés.

La solution numérique H de l'équation différentielle (21a)–(21b) a la
propriété de prendre la valeur 0. Nous sommes dans le cas (21c) et nous obtenons
M∞ = 27 % et H(M∞) = 0. (26)
Ici, le seuil d'infection H = 1 % est si bas que les états initiaux présentant une forte proportion m
de moustiques infectés — à savoir supérieure à M∞ — ne sont pas viables. Les
valeurs M = 14 % et M∞ = 27 %, ainsi que la figure 2, ont des
implications pratiques pour une lutte viable contre la dengue. L'analyse du §3.3
suggère que, pour plafonner la proportion d'humains infectés au moment du pic, il
16
faut mesurer la proportion de moustiques infectés. Mesurer l'abondance des moustiques
est une tâche difficile, alors que la mesure de la proportion de moustiques infectés
peut être effectuée grâce à un échantillonnage approprié. Voici une conception possible
pour une politique viable.
• Surveillance sans fumigation. Lorsque la proportion d'humains infectés
est inférieure au seuil d'infection H = 1 % et que la proportion de
moustiques infectés est inférieure à la proportion M = 14 %, ne pas procéder à une fumigation.
• Surveillance avec fumigation (maximale). Lorsque la proportion de
moustiques infectés se situe entre les proportions M = 14 % et M∞ =
27 %, procéder à une fumigation à capacité maximale.
• Alerte. Lorsque la proportion de moustiques infectés dépasse M∞ = 27 %,
des mesures supplémentaires doivent être prises pour éviter un pic élevé
d'humains infectés.
17

Introduction

Le « problème de navigation de Zermelo » originel s'énonce comme suit : dans un plan illimité où la répartition du vent est définie par un champ de vecteurs dépendant de la position et du temps, un navire se déplace à vitesse constante par rapport à la masse d'air environnante. Comment le navire doit-il être dirigé pour rallier une destination donnée depuis un point de départ dans le temps le plus court ?

Ce problème peut-être modélisé sous la forme d'un système dynamique contrôlé dont on cherche la commande optimale, celle qui minimise le temps mis pour atteindre la cible.
Les variables d'état $x$ et $y$ désignent les coordonnées du navire en mouvement, tandis que $f(x, y)$ et $g(x, y)$ représentent les composantes de la vitesse du vent qui dépendent de la position. Deux variables de commande régissent l'évolution du navire : la vitesse $v$, dont la norme $||v||$ est bornée par un seuil dépendant de la position ($||v|| ≤ c(x,y)$), et la direction $u$, c'est-à-dire l'angle que forme le vecteur vitesse relative avec l'axe des $x$ qui est contrainte par les bornes $\alpha_{min}(x, y)$ et $\alpha_{max}(x, y)$.
L'évolution du navire est ainsi régie par le système dynamique contrôlé suivant :
\begin{equation}
\left\{
\begin{array}{l}
x'(t) &= f(x(t), y(t)) + v(t) \cos u(t) \\
y'(t) &= g(x(t), y(t)) + v(t) \sin u(t) \\
u(t) &\in [\alpha_{min}(x(t), y(t)),\alpha_{max}(x(t), y(t)] \text{ et } ||v(t)|| ≤ c(x(t), y(t)).
\end{array}
\right.
\end{equation}
 

Etant donné un point de départ $(x_0,y_0)$, il s'agit de trouver les fonctions de commande $u(t)$ et $v(t)$ qui permettent d'atteindre une cible en temps minimal. 

Le problème de Zermelo revisité avec les concepts de la théorie de la viabilité

Depuis

J.-P.AUBIN & H. FRANKOWSKA, « The viability kernel algorithm for computing value functions of infinite horizon optimal control problems », Journal of mathematical analysis and applications, vol. 201, 555–576, 1996.

des liens sont établis entre les problèmes de contrôle optimal et les noyaux de viabilité ou bassins de capture.

Le problème de temps minimal pour l'atteinte d'une cible en est un cas particulier et le problème de Zermelo a été souvent illustré. 

Dans leur livre

Aubin, J.-P., Bayen, A., & Saint-Pierre, P. Viability Theory: New Directions. Springer. 2011.

les auteurs revisitent le problème de Zermelo en remplaçant le "plan illimité" par un environnement $K\subset \mathbb{R}^2$ avec des obstacles, le « but » étant le port, il joue le rôle de la « cible » $C \subset K$ dans la terminologie de la théorie de la viabilité. Il s’agit bien entendu d’un problème plus complexe que le problème de Zermelo initial en raison de ces contraintes supplémentaires définies par l'ensemble $K$ dans lequel les trajectoires doivent rester. 

Associer à chaque point de l'environnement $K$ la valeur du temps minimal nécessaire pour atteindre la cible $C$ revient à calculer le bassin de capture dans la terminologie de la théorie de la viabilité pour le système initial augmenté d'une dimension, celle du coût associé au temps qui s'écoule. En effet, ce bassin de capture correspond à l'épigraphe de la fonction de temps minimal (i.e. l'ensemble des points situés au-dessus de sa courbe représentative) :

$$
t_{min}(K,C)(x_0,y_0) = \inf_{(x_0,y_0,z)\in Capt(K,C)}z
$$

De plus, l'ensemble des trajectoires issues de $(x_0,y_0)$ réalisant ce temps minimal d'atteinte de cible correspond à l'ensemble des trajectoires issues de $(x_0,y_0)$ contenues dans $Capt(K,C)$.

Résultats

Dynamique et ensemble de contraintes

Pour l'illustration numérique, nous choisissons comme zone de déplacement du bateau $K\subset [-6;1]\times[-2.5;2.5]$ privé de deux sous-ensembles représentant deux obstacles $0_1 =  [-4;3]\times[-0;2]$ et $0_2 =  [-2;-1]\times[-0.75;0.75]$ et comme cible $C = \{(x,y) | x^2+y^2\leq 0.04$. Ces sous-ensembles du plan sont dessinés dans le graphique ci-dessous :

La dynamique du courant est $f(x,y)=-0.25\tanh y $ et $g(x,y)=0$, sans contrainte sur le choix de l'angle de la vitesse relative $u(t) \in [0,2\pi]$ mais avec la norme de cette vitesse relative constante égale à 1, $|v(t)| =1$.  Ainsi,  l'évolution du navire est régie par l'équation :

\begin{equation}
\left\{
\begin{array}{l}
x'(t) &= \cos u(t) -0.25\tanh y(t)\\
y'(t) &= \sin u(t) \\
u(t) &\in [0,2\pi]
\end{array}
\right.
\end{equation}

Bassin de capture

Le bassin de capture est un sous-ensemble de dimension 3 (les 2 dimensions du problème initial + 1 dimension pour la fonction de coût associée au temps minimal).

Ce bassin de capture est épigraphique, c'est à dire qu'il se situe au-dessus de sa frontière inférieure dans la dimension de la fonction de coût. Nous reproduisons ci-dessous l'approximation de sa frontière inférieure calculée par le logiciel ViabLab à partir d'une grille de $701\times 501$ points dans l'espace initial de dimension 2 :

 
Les deux "trous" correspondent aux deux obstacles pour les points desquels la valeur du temps minimal pour atteindre la cible est $+\infty$.

Lignes de niveau de la fonction temps minimal

Ce résultat peut également être visualisé sous forme de lignes de niveaux dans l'environnement du problème initial : 

Trajectoires optimales

Comme précisé ci-dessus, les trajectoires optimales sont des sous-produits du noyau de viabilité. 
Ainsi, le logiciel ViabLab permet d'en calculer des approximations par des courbes linéaires par morceaux dont la précision dépend des approximations de pas de temps et d'espace conduisant à l'approximation du noyau de viabilité.

Ci-dessous sont représentées deux approximations de trajectoires optimales pour les points initiaux $(-4,-1,5)$ et $(-5,1)$ dans les deux représentations :

• à gauche dans l'espace 3D contenant la variable de coût (le temps qui s'écoule), les trajectoires sont issues des points $(-4,-1,5, t_{min}(4,-1,5))$ et $(-5,1,t_{min}(-5,1))$ et restent sur la frontière du bassin de capture
• à droite dans l'espace 2D du problème initial (projection des trajectoires de dimension 3)

Introduction

De nombreuses langues risquent de disparaître. Il est donc crucial de comprendre les mécanismes qui sous-tendent l'évolution du nombre de locuteurs et de déterminer s'il existe des mesures susceptibles d'aider à en préserver certaines.

Dans leur article,

Abrams, S. Strogatz, Modelling the dynamics of language death, Nature 424 (6951) (2003) 900

Abrams et Strogatz ont proposé un modèle mathématique pour étudier la concurrence entre les langues.
Ce modèle s'ajuste bien à plusieurs ensembles de données empiriques : il rend compte de manière satisfaisante des données historiques concernant le déclin du gallois,
du gaélique écossais, du quechua et d'autres langues menacées.

Dans le modèle d'Abrams et Strogatz ainsi que dans ceux qui en sont dérivés, on suppose que la taille de la population demeure constante. Par conséquent, les variables sont les proportions des différents groupes de locuteurs. Dans les modèles incluant le bilinguisme, la population se compose de trois groupes : les locuteurs monolingues de la langue A, les locuteurs monolingues de la langue B et les locuteurs bilingues AB ; le modèle est bidimensionnel, avec $\sigma_A$ représentant la proportion de locuteurs de la langue A et $\sigma_B$ la proportion de locuteurs de la langue B ($\sigma_{AB} = 1 - \sigma_A - \sigma_B$).

Au sein de toute sous-population linguistique, des forces et des influences amènent les membres d'un groupe à changer de langue au profit de l'autre. Dans le modèle d'Abrams et Strogatz, la vitesse à laquelle les locuteurs d'une langue apprennent la seconde langue dépend de l'attractivité de cette dernière. Dans leur définition la plus générale de l'attractivité, Abrams et Strogatz postulent qu'une langue est d'autant plus attractive qu'elle compte un grand nombre de locuteurs et qu'elle jouit d'un prestige élevé.

Dans leur article,

C. Bernard, S. Martin, Building strategies to ensure language coexistence in presence of bilingualism, Appl. Math. Comput. (2012), doi:10.1016/j.amc.2012.02.041

les auteurs modélisent uniquement les transitions des quatre types suivants : A → AB, AB → A, A → AB et AB → B ((les transitions A → B et B → A étant extrêmement rares en pratique). Ils supposent également une asymétrie entre monolingues et bilingues : A → AB (resp. B → AB) à un taux proportionnel à l'attractivité des locuteurs monolingues de A (resp. B) ; AB → A (resp. AB → B) à un taux proportionnel à l'attractivité de l'ensemble des locuteurs de A, y compris les bilingues (ainsi, certains bilingues peuvent devenir des locuteurs monolingues de A même si A ne compte aucun locuteur monolingue).
Par conséquent, le modèle bidimensionnel est défini par :

$$
\left\{\begin{array}{ll}
\sigma_A' &= (1-\sigma_A-\sigma_B)(1-\sigma_B)^as_A-\sigma_A\sigma_B^as_B\\
\sigma_B' &= (1-\sigma_A-\sigma_B)(1-\sigma_A)^as_B-\sigma_B\sigma_A^as_A
\end{array}\right.
$$

où $s_A$  (resp. $s_B$) désigne le prestige de la langue A (resp. B), et $a$ est un paramètre modélisant la manière dont l'attractivité d'une langue varie en fonction de la proportion de ses locuteurs.
Par commodité, les auteurs supposent que $s_A+s_B=1$ ce qui permet de substituer $s_A=s$ et $s_B=1-s$.

Le problème de la coexistence de deux langues

Bernard et Martin considèrent que le prestige, $s$, peut évoluer (sous l'effet de l'action publique, par exemple), mais que sa variation à chaque pas de temps est bornée. Ainsi, 

$$
s' = u\in U=[\underline{u};\bar{u}]
$$

Lorsque que $u>0$, le prestige de la langue A augmente et celui de B diminue, c'est l'inverse lorsque $u<0$.

La dynamique du système est ainsi décrite par :

$$
\left\{\begin{array}{ll}
\sigma_A' &= (1-\sigma_A-\sigma_B)(1-\sigma_B)^as-\sigma_A\sigma_B^a(1-s)\\
\sigma_B' &= (1-\sigma_A-\sigma_B)(1-\sigma_A)^a(1-s)-\sigma_B\sigma_A^as\\
s' &= u\in U=[\underline{u};\bar{u}]
\end{array}\right.
$$

Les auteurs proposent de rechercher des stratégies fondées sur les variations de prestige qui garantissent la coexistence des deux langues, c'est à dire qui maintiennent un niveau donné de locuteurs monolingues dans chacune des deux langues. L'ensemble des contraintes s'écrit donc :

$$
K:= [\underline{\sigma};1]\times[\underline{\sigma};1]\times [0;1]
$$

Résoudre ce problème, c'est trouver l'ensemble des états $(\sigma_A,sigma_B,s)$ à partir desquels il est possible de maintenir une certaine proportion des groupes monolingues en agissant sur la variation des prestiges relatifs des deux langues, c'est à dire le noyau de viabilité dans la terminologie de la théorie de la viabilité.

Résultats

Dans l'article cité ci-dessus, les auteurs ont établi une description analytique de la frontière du noyau de viabilité composée de courbes intégrales issues de points particuliers de la frontière de l'ensemble des contraintes.
Pour illustrer leur résultat, ils ont reproduit le noyau de viabilité pour les valeurs de paramètres $a=1.41$, $\underline{u}=-0.06$, $\bar{u}=0.06$, $\underline{\sigma}=0.1$.

Nous reproduisons ci-dessous 4 vues 3D de l'approximation du noyau de viabilité associé à ces valeurs de paramètres par le logiciel ViabLab à partir d'une grille de $201\times 201\times 201$ points (les lignes noires  représentent les arrêtes de l'ensemble de contraintes): 

Comme attendu, le noyau est invariant par la transformation $\sigma_A\rightarrow \sigma_B$, $\sigma_B\rightarrow \sigma_A$ et $s\rightarrow 1-s$.

Les situations non viables correspondent à des états à partir des desquels :

• même une augmentation continue d'intensité maximale du prestige de la langue A jusqu'à sa valeur maximale 1 ne permettra pas d'éviter que la proportion de ses locuteurs ne passe le seuil $\underline{\sigma}$,
• ou même une augmentation continue d'intensité maximale du prestige de la langue B jusqu'à sa valeur maximale 1 ne permettra pas d'éviter que la proportion de ses locuteurs ne passe le seuil $\underline{\sigma}$.

Introduction

L'un des enjeux majeurs de la gestion des prairies pâturées consiste à concilier des objectifs antagonistes : accroître les profits à court terme tout en préservant la viabilité à long terme du parcours. 
En effet, augmenter la charge pastorale (nombre d'animaux par hectare) permet d'accroître les gains à court terme, mais peut entraîner une dégradation pratiquement irréversible du milieu si la pression de pâturage qui en résulte est excessive. À l'inverse, maintenir une charge pastorale trop faible peut engendrer un manque à gagner considérable. 

Dans leur article 

Anderies, J., Janssen, M. & Walker, B. Grazing Management, Resilience, and the Dynamics of a Fire-driven Rangeland System. Ecosystems 5, 23–44 (2002). https://doi.org/10.1007/s10021-001-0053-9

les auteurs utilisent le modèle suivant pour rendre compte des interactions entre la pousse de l'herbe et la pression de pâturage : 

$$
\left\{\begin{array}{ll}
c' &= r_cs-\delta_cc\\
s'&= c(a_c+r_ss)(1-\frac{s}{s^*}-\alpha_{ws}(\frac{w}{w^*})^\beta)-\gamma_gs\\
w'&=r_ww(1-\frac{w}{w^*})
\end{array}\right.
$$

Concernant l'herbe, $c(t)$ et $s(t)$ représentent respectivement la biomasse de la couronne (racines plus points de croissance) et celle de la partie aérienne à l'instant $t$. Le paramètre $s^*$ correspond à la biomasse aérienne maximale par unité de surface, tandis que $\gamma_g\leq 1$ représente la fraction de la partie aérienne prélevée par la pression de pâturage exercée par l'ensemble des animaux présents. La couronne croît à un taux $r_cs$ en présence de parties aériennes et régresse à un taux $\delta_c$.
$a_c$ correspond au taux de croissance de la biomasse aérienne par unité de biomasse de la couronne, $r_s$ au taux de croissance de la biomasse aérienne par unité de biomasse de la couronne et par unité de biomasse aérienne.

La croissance de l'herbe est également affectée par la présence d'arbustes : $w(t)$ représente la biomasse des arbustes, $r_w$ le taux de croissance intrinsèque des arbustes, $w^*$ la biomasse maximale des arbustes par unité de surface, et $w_s(=w/w^*)$ modélise l'effet de compétition des arbustes sur l'herbe. Le paramètre $\beta$ mesure la non-linéarité relative de cet effet.

Le problème de viabilité du gestionnaire de pâtures

L'objectif de la gestion des parcours de pâturage consiste à identifier des stratégies d'intervention (modifications de la charge pastorale) permettant de concilier deux contraintes : nourrir suffisamment les animaux tout en préservant la qualité de l'herbe.

Dans l'article 

les auteurs étudient ce problème de viabilité en utilisant un modèle de pousse d'herbe simplifié en supposant que la biomasse de tiges maximale vaut 1 et que l'impact de la population d'arbustes est négligeable c'est à dire que : 

$$
s^* = 1 \text{ et } \alpha_{ws} = 0.$$

L'action possible sur le système est de faire varier la pression de pâturage et cette variation est supposée bornée, l'achat de nouveaux animaux  ou leur vente prend  un certain temps :

$$
\gamma'(t) = u(t) \in [u_{min},u_{max}].
$$

La dynamique du système est ainsi décrite par :

$$
\left\{\begin{array}{ll}
c' &= r_cs-\delta_cc\\
s'&= \frac{r_c}{\delta_c}s(a_c+r_ss)(1-s)-\gamma_gs\\
\gamma_g'&= u \in [u_{min},u_{max}]
\end{array}\right.
$$

La contrainte sur la qualité de l'herbe est modélisée par une valeur de la biomasse aérienne supérieure à un seuil, $s_{min}$, pour assurer la repousse. La contrainte sur la nourriture du bétail est modélisée par une pression de pâturage également supérieure à un seuil, $\underline{g}$, au cours du temps. L'ensemble des contraintes s'écrit donc :

$$
K(s_{min},\underline{g}):=R^+\times [s_{min},+\infty[\times  [\underline{g},1]
$$

Résoudre ce problème, c'est trouver l'ensemble des états $(c,s,\gamma_g)$ à partir desquels il est possible de nourrir le troupeau en conservant une herbe de qualité en agissant sur la variation du nombre d'animaux, c'est à dire le noyau de viabilité dans la terminologie de la théorie de la viabilité.

Résultats

Dans l'article cité ci-dessus les valeurs des paramètres sont les suivantes : $r_c=1$, $\delta_c=1$, $a_c=0.1$,  $r_s=3$,  $s^*=1$, $u_{min} = -0.05$, $u_{max} = 0.05$, $s_{min} = 0.1$, $\underline{g}=0.65$.

Nous reproduisons ci-dessous l'approximation du noyau de viabilité associé à ces valeurs de paramètres par le logiciel ViabLab à partir d'une grille de 201*201*201 points : 

Comme attendu, lorsque l'on part d'une pression de pâturage plus élevée, l'ensemble des combinaisons de biomasses tiges et couronnes pour lesquelles le système est viable est plus restreint.

Cette restriction est d'autant plus forte que la capacité de diminution de cette pression est faible, c'est à dire quand $u_{min}$ se rapproche de 0.

Introduction

L'un des enjeux majeurs de la gestion des prairies pâturées consiste à concilier des objectifs antagonistes : accroître les profits à court terme tout en préservant la viabilité à long terme du parcours. 
En effet, augmenter la charge pastorale (nombre d'animaux par hectare) permet d'accroître les gains à court terme, mais peut entraîner une dégradation pratiquement irréversible du milieu si la pression de pâturage qui en résulte est excessive. À l'inverse, maintenir une charge pastorale trop faible peut engendrer un manque à gagner considérable. 

Dans leur article 

Anderies, J., Janssen, M. & Walker, B. Grazing Management, Resilience, and the Dynamics of a Fire-driven Rangeland System. Ecosystems 5, 23–44 (2002). https://doi.org/10.1007/s10021-001-0053-9

les auteurs utilisent le modèle suivant pour rendre compte des interactions entre la pousse de l'herbe et la pression de pâturage : 

$$
\left\{\begin{array}{ll}
c' &= r_cs-\delta_cc\\
s'&= c(a_c+r_ss)(1-\frac{s}{s^*}-\alpha_{ws}(\frac{w}{w^*})^\beta)-\gamma_gs\\
w'&=r_ww(1-\frac{w}{w^*})
\end{array}\right.
$$

Concernant l'herbe, $c(t)$ et $s(t)$ représentent respectivement la biomasse de la couronne (racines plus points de croissance) et celle de la partie aérienne à l'instant $t$. Le paramètre $s^*$ correspond à la biomasse aérienne maximale par unité de surface, tandis que $\gamma_g\leq 1$ représente la fraction de la partie aérienne prélevée par la pression de pâturage exercée par l'ensemble des animaux présents. La couronne croît à un taux $r_cs$ en présence de parties aériennes et régresse à un taux $\delta_c$.
$a_c$ correspond au taux de croissance de la biomasse aérienne par unité de biomasse de la couronne, $r_s$ au taux de croissance de la biomasse aérienne par unité de biomasse de la couronne et par unité de biomasse aérienne.

La croissance de l'herbe est également affectée par la présence d'arbustes : $w(t)$ représente la biomasse des arbustes, $r_w$ le taux de croissance intrinsèque des arbustes, $w^*$ la biomasse maximale des arbustes par unité de surface, et $w_s(=w/w^*)$ modélise l'effet de compétition des arbustes sur l'herbe. Le paramètre $\beta$ mesure la non-linéarité relative de cet effet.

Le problème de viabilité du gestionnaire de pâtures

L'objectif de la gestion des parcours de pâturage consiste à identifier des stratégies d'intervention (modifications de la charge pastorale) permettant de concilier deux contraintes : nourrir suffisamment les animaux tout en préservant la qualité de l'herbe.

Dans l'article 

Martin, Sophie & Deffuant, Guillaume & Calabrese, Justin. (2011). Defining Resilience Mathematically: From Attractors To Viability. 10.1007/978-3-642-20423-4_2. 

les auteurs étudient ce problème de viabilité en utilisant un modèle de pousse d'herbe simplifié en supposant que la population d'arbuste est constante à la valeur d'équilibre et que la biomasse de la couronne prend instantanément la valeur d'équilibre associée à la biomasse aérienne : 

$$
w = w^* \text{ and  } c = \frac{r_c}{\delta_c}s.
$$

L'action possible sur le système est de faire varier la pression de pâturage et cette variation est supposée bornée, l'achat de nouveaux animaux  ou leur vente prend  un certain temps :

$$
\gamma'(t) = u(t) \in [u_{min},u_{max}].
$$

La dynamique du système est ainsi décrite par :

$$
\left\{\begin{array}{ll}
\gamma_g'&= u \in [u_{min},u_{max}]\\
s'&= \frac{r_c}{\delta_c}s(a_c+r_ss)(1-\frac{s}{s^*}-\alpha_{ws})-\gamma_gs.
\end{array}\right.
$$

La contrainte sur la qualité de l'herbe est modélisée par une valeur de la biomasse aérienne comprise entre deux bornes, la borne inférieure pour assurer le renouvellement, la borne supérieure pour préserver ses qualités nutritives. La contrainte sur la nourriture du bétail est modélisée par une pression de pâturage également supérieure à un seuil au cours du temps. L'ensemble des contraintes s'écrit donc :

$$
K(\underline{g},s_{min},s_{max}):=[\underline{g},1]\times [s_{min},s_{max}].
$$

Résoudre ce problème, c'est trouver l'ensemble des états $(s,\gamma_g)$ à partir desquels il est possible de nourrir le troupeau en conservant une herbe de qualité en agissant sur la variation du nombre d'animaux, c'est à dire le noyau de viabilité dans la terminologie de la théorie de la viabilité.

Résultats

Pour les valeurs des paramètres de l'article cité ci-dessus ($r_c=1$, $\delta_c=1$, $a_c=0.3$,  $r_s=3$,  $s^*=1$,  $\alpha_{ws}=0.5$, $u_{min} = -0.02$, $u_{max} = 0.02$, $\underline{g}=0$, $s_{min} = 0.1$, $s_{max} = 0.18$), le noyau de viabilité peut être obtenu par le calcul de deux courbes intégrales (les auteurs de l'article n'ont pas développé ce point).

En effet, sa frontière est définie par une courbe fermée formée de sous-ensembles de la frontière de l'ensemble de contraintes $K$ et de deux courbes intégrales : 

• l'une issue du point $(\gamma_g^1,s_{min}) $ avec $s'(\gamma_g^1,s_{min})) = 0$ qui suit la dynamique inverse : 
  \begin{equation}
  \left\{\begin{array}{ll}
  \gamma_g'&= -u_{min}\\
  s'&= -(\frac{r_c}{\delta_c}s(a_c+r_ss)(1-\frac{s}{s^*}-\alpha_{ws})-\gamma_gs)\\ 
  \end{array}\right.
  \end{equation}
• l'autre issue du point $(\gamma_g^2,s_{max}) $ avec $s'(\gamma_g^2,s_{max})) = 0$ qui suit la dynamique inverse : 
  \begin{equation}
  \left\{\begin{array}{ll}
  \gamma_g'&= -u_{max}\\
  s'&= -(\frac{r_c}{\delta_c}s(a_c+r_ss)(1-\frac{s}{s^*}-\alpha_{ws})-\gamma_gs)\\ 
  \end{array}\right.
  \end{equation}

La figure suivante montre la frontière du noyau de viabilité sous cette forme de courbe fermée ainsi que l'approximation obtenue avec le logiciel ViabLab sur une grille de 2001 * 2001 points. Le résultat obtenu est cohérent avec le graphique de l'article original et bien plus précis dans la partie de la frontière du noyau qui coïncide avec la frontière de l'ensemble des contraintes ainsi que dans celle définie par les courbes intégrales.

Voici un exemple dont on connaît le noyau de viabilité.

En effet, prenons une dynamique en deux dimensions dont l'orbite des points du cercle de centre l'origine du repère et de rayon $a$ est ce même cercle ; dont l'orbite des points de norme strictement inférieure à $a$ est une spirale qui tend vers l'origine et dont l'orbite des points de norme strictement supérieure à $a$ est une spirale qui diverge.

C'est le cas des solutions de l'équation $x'=f(x)$ avec 

$$
f(x_1,x_2):=(\frac{||x||-a}{||x||}x_1+x_2 ; -x_1+\frac{||x||-a}{||x||}x_2)
$$   

Prenons comme ensemble de contraintes, le carré  : 

$$
K := \{(x_1,x_2)\in R^2\,|\, |x_1|\leq c \text{ et } |x_2|\leq c \}.
$$

De manière évidente, l'ensemble des points dont sont issues des trajectoires qui restent dans $K$ indéfiniment  est le disque de centre l'origine du repère et de rayon $a$.

Cet exemple permet de tester la finesse d'approximations issues de logiciels génériques.

Ci-dessous la comparaison pour $a = 0.7$ et $c=1.5$ entre le résultat théorique (le disque dont la frontière est le cercle rouge) et le résultat approché obtenu avec le logiciel ViabLab sur une grille régulière de 2001*2001points (ensemble vert avec la frontière bleue) :