Introduction
Ce tutoriel illustre l'ajustement d'une donnée initiale paramétrée d'un modèle de réaction diffusion. Les deux paramètres (p1,p2) représentent le centre de la distribution de la donnée initiale. Il s'agit ici de déterminer les valeurs du couple (p1,p2) minimisant l'erreur quadratique entre la simulation et les données.
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 donc d'une étude portant sur 1 espèce avec une donnée initiale paramétrée par 2 paramètres:
import mse
MY_OUTPUT_DIR=os.getenv("MSE_OUTPUT_DIR", "/home/mseuser/MSE/demoIC/")
mse.setStringParam(mse.MSE_KEY_OUTPUT_DIR,MY_OUTPUT_DIR)
mse.setIntParam(mse.MSE_KEY_MODEL_NBSPECIES,1)
mse.setIntParam(mse.MSE_KEY_MODEL_NBPARAMS,2)
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 30 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,40.0)
mse.setIntParam(mse.MSE_KEY_GEOM_SYM,0)
mse.setIntParam(mse.MSE_KEY_GEOM_REVEDGE,1)
mse.buildGeom()
mse.buildDirectories()
os.system("cd "+MY_OUTPUT_DIR+";/opt/freefem++/src/nw/FreeFem++-nw "+MY_OUTPUT_DIR+"/buildXY.edp")
Description du modèle de réaction-diffusion
Le modèle étudié ici est de la forme:
Le coefficient de diffusion homogène est 0.02:
Le terme de réaction est 5u(1-u):
L'état initial est défini de la façon suivante:
mse.setdpiDiffusionxExp(0,0,0,"0")
mse.setdpiDiffusionxExp(1,0,0,"0")
mse.setDiffusionyExp(0,0,"0.02")
mse.setdpiDiffusionyExp(0,0,0,"0")
mse.setdpiDiffusionyExp(1,0,0,"0")
mse.buildDiffusion()
mse.setdyfxExp(0,0,"-10*a + 5")
mse.setdPiFj(0,0,"0")
mse.setdPiFj(1,0,"0")
mse.buildModel()
mse.buildModelParams()
mse.buildWVarf()
mse.printTSOneStepW()
mse.setIntParam(mse.MSE_KEY_FORCE_INITIAL_STATE,1)
mse.setInitialState(0,0,"exp(-40*((x-p1)*(x-p1)+(y-p2)*(y-p2)))")
mse.setdPICiFj(0,0,"(-80*p1 + 80*x)*exp(-40*(-p1 + x)^2 - 40*(-p2 + y)^2)")
mse.setdPICiFj(1,0,"(-80*p2 + 80*y)*exp(-40*(-p1 + x)^2 - 40*(-p2 + y)^2)")
mse.printInitialState()
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.05. La non-linéarité sera gérée par un algorithme de Newton-Raphson et une tolérance absolue de 0.02.
mse.setRealParam(mse.MSE_KEY_SIM_MINSTEP,0.05)
mse.setIntParam(mse.MSE_KEY_SIM_MAXMULTSTEP,1)
mse.setRealParam(mse.MSE_KEY_SIM_NEWTON_CRIT,0.01)
mse.setIntParam(mse.MSE_KEY_SIM_NEWTON_MIN_IT,1)
mse.setIntParam(mse.MSE_KEY_SIM_NEWTON_MAX_IT,5)
mse.setIntParam(mse.MSE_KEY_SIM_FORCE_POSITIV,1)
mse.buildTS()
Définition de la fonction objectif à minimiser
Il s'agit ici de minimiser la différence quadratique entre la simulation et les données pour quelques points de l'espace et du temps.
La fonction de coût associé au paramètre (p1,p2) est:
Avec u la simulation pour les paramètre (p1,p2) et uref la simulation servant de référence.
Cela se spécifie de la façon suivante dans l'interface:
mse.buildContrast()
Description des données
Pour cela, il faut créer les fichiers suivants dans le répertoire de l'étude initialement spécifié "/home/mseuser/demoIC/PROCESSUS/DATA".
Ici, nous avons deux dates et trois lieux d'observations:
Toutes les couples dates X lieux sont présents dans les données:
Pour cet exemple, les valeurs ont été calculées à partir d'une simulation avec les paramètres (p1=-0.25, p2=0.3).
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:
mse.setBfgsLowerBound(0,-1.0)
mse.setBfgsLowerBound(1,-1.0)
mse.setBfgsUpperBound(0,1.0)
mse.setBfgsUpperBound(1,1.0)
mse.setBfgsInitialState(0,0.5)
mse.setBfgsInitialState(1,0.5)
mse.setBfgsTol(1e-3)
mse.buildBfgs()
Résultat de l'optimisation
Il faut 21 appels au simulateur pour converger vers la valeur (p1=-2.499, p2=0.3). Le graphique suivant représente la trajectoire des itérations de l'algorithme d'optimisation.
