Introduction
Ce tutoriel illustre l'ajustement d'une diffusion hétérogène et d'un modèle proie prédateur. Les paramètres p1 et p2 caractérisent la diffusion, p3 représente le taux de prédation, p4 est la mortalité du prédateur.
La suite de cette page décrit comment mener cette étude dans MSE, d'une part via l'interface, d'autre part via un script python.
Description de l'étude
Il s'agit d'une étude portant sur 2 espèces et 4 paramètres:
import mse
MY_OUTPUT_DIR=os.getenv("MSE_OUTPUT_DIR", "/home/mseuser/MSE/demo2s4p/")
mse.setStringParam(mse.MSE_KEY_OUTPUT_DIR,MY_OUTPUT_DIR)
mse.setIntParam(mse.MSE_KEY_MODEL_NBSPECIES,2)
mse.setIntParam(mse.MSE_KEY_MODEL_NBPARAMS,4)
mse.initGradEdp()
Description de la zone d'étude
La zone d'étude est le carré [-1,1]2. La triangulation sera créée à partir de 45 points par unité de longueur des arêtes:
mse.setStringParam(mse.MSE_KEY_GEOM_DIR,MY_GEOM_DIR)
mse.setRealParam(mse.MSE_KEY_GEOM_NPERLENGTH,45.0)
mse.setIntParam(mse.MSE_KEY_GEOM_SYM,0)
mse.setIntParam(mse.MSE_KEY_GEOM_REVEDGE,1)
mse.buildGeom()
mse.buildDirectories()
Description du modèle proie-prédateur
Le script définissant la donnée initiale est defICGen.edp représentée sur la figure suivante :
La figure suivante représente la fonction H :
La covariable H a été créée à partir d'une image comme décrit dans le tutoriel 'importer une image'.
Ce modèle est décrit dans l'interface de la façon suivante :
mse.setdyfxExp(0,0,"-H*b*p3 - 2*a + 1")
mse.setdyfxExp(1,0,"-H*a*p3")
mse.setdPiFj(0,0,"0")
mse.setdPiFj(1,0,"0")
mse.setdPiFj(2,0,"-H*a*b")
mse.setdPiFj(3,0,"0")
mse.setfxExp(1,"H*p3*a*b-p4*b")
mse.setdyfxExp(0,1,"H*b*p3")
mse.setdyfxExp(1,1,"H*a*p3 - p4")
mse.setdPiFj(0,1,"0")
mse.setdPiFj(1,1,"0")
mse.setdPiFj(2,1,"H*a*b")
mse.setdPiFj(3,1,"-b")
mse.buildModel()
mse.buildModelParams()
mse.setDiffusionxExp(0,0,"0.1")
mse.setdpiDiffusionxExp(0,0,0,"0")
mse.setdpiDiffusionxExp(1,0,0,"0")
mse.setdpiDiffusionxExp(2,0,0,"0")
mse.setdpiDiffusionxExp(3,0,0,"0")
mse.setDiffusionyExp(0,0,"0.1")
mse.setdpiDiffusionyExp(0,0,0,"0")
mse.setdpiDiffusionyExp(1,0,0,"0")
mse.setdpiDiffusionyExp(2,0,0,"0")
mse.setdpiDiffusionyExp(3,0,0,"0")
mse.setDiffusionxExp(1,0,"p1*H+p2*(1-H)")
mse.setdpiDiffusionxExp(0,1,0,"H")
mse.setdpiDiffusionxExp(1,1,0,"-H + 1")
mse.setdpiDiffusionxExp(2,1,0,"0")
mse.setdpiDiffusionxExp(3,1,0,"0")
mse.setDiffusionyExp(1,0,"p1*H+p2*(1-H)")
mse.setdpiDiffusionyExp(0,1,0,"H")
mse.setdpiDiffusionyExp(1,1,0,"-H + 1")
mse.setdpiDiffusionyExp(2,1,0,"0")
mse.setdpiDiffusionyExp(3,1,0,"0")
mse.buildDiffusion()
mse.buildWVarf()
mse.printTSOneStepW()
Définition des paramètres du simulateur
On choisit de simuler sur la période de temps [0,1], avec un pas de temps fixe égale à 0.1. La non-linéarité sera gérée par un algorithme de Newton-Raphson et une tolérance absolue de 0.02.
Définition de la fonction objectif à minimiser
Il s'agit ici de minimiser l'intégrale de la différence quadratique entre la simulation et les données pour l'état final de la simulation.
La fonction de coût associé aux paramètres (p1,p2,p3,p4) est :
La date d'observation est précisée dans un fichier dataDates.txt qu'il faudra créer dans le répertoire de l'étude initialement spécifié "/home/mseuser/DEMO2s4p/PROCESSUS/DATA".
Cela se spécifie de la façon suivante dans l'interface:
mse.buildContrast()
Description des données
La simulation de référence a été calculée à partir d'une simulation avec les paramètres (p1=2, p2=5,p3=10,p4=2).
Description du problème d'optimisation
Il s'agit d'un problème d'optimisation sous contraintes. Pour le résoudre, un algorithme de Quasi-Newton BFGS est couplé au simulateur. Il s'agit d'un algorithme itératif pour lequel le point de départ est précisé ainsi que les conditions de convergence.
Cela se spécifie de la façon suivante dans l'interface:
mse.buildContrast()
Résultat de l'optimisation
Le graphique suivant montre l'évolution de la fonction objectif au cours des itérations de l'algorithme d'optimisation. On constate que le résultat dépend du point de départ et de la tolérance.
A la fin de l'optimisation, on obtient quelques renseignements sur la fonction objectif:
Quelques images de simulation du modèle optimal
La figure ci-dessous représente la population de proie pour t=0, t=0.5 et t=1:
La figure ci-dessous représente la population du prédateur pour t=0, t=0.5 et t=1:
