Tutoriel - Estimation de la condition initiale

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:

IHM demo IC general

import os
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:

tuto demo IC geometrie

MY_GEOM_DIR=os.getenv("MSE_GEOM_DIR","/home/usermse/MSE/GEODATA/CARREm1p1/")
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:

im

Le coefficient de diffusion homogène est 0.02:

tuto demo IC

Le terme de réaction est 5u(1-u):

tuto reaction

L'état initial est défini de la façon suivante:

tuto demoIC def IC

mse.setDiffusionxExp(0,0,"0.02")
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.setfxExp(0,"5*a*(1-a)")
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.initInitialState()
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.

demoIC simulateur

 

mse.setRealParam(mse.MSE_KEY_SIM_TEND,1.1)
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:

demo IC cout

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:

demoIC contraste

mse.setIntParam(mse.MSE_KEY_CONTRAST_TYPE,6)
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:

demoIC optim

mse.initBfgsStruct()
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.

demoIC traj optim