ENSTBB-IPB 2018-2019

Séance 3 : Régression Non Linéaire avec R

Exercice 0.1 Modèle exponentiel
Procédé en batch dans un fermenteur
L’objectif d’un ”procédé en batch” de génie fermentaire est de déterminer les caractéristiques cinétiques d’un microor-
ganisme en particulier son taux de croissance µ. Pour cela, on va ensemencer un bioréacteur avec (entre autre) une
certaine concentration de microorganismes et de substrat dont on va mesurer l’évolution au cours du temps. Des mesures
de densité optique seront régulièrement effectuées à l’aide d’un spectrophotomètre. Les échantillons prélevés seront dilués
afin de rester dans la zone de linéarité du spectrophotomètre, zone pour laquelle la densité optique est proportionnelle à
la concentration. Après différentes calibrations, on a obtenu les concentrations suivantes
Temps ti (h) Biomasse xi (g/L)
0 0,11
1 0,15
2 0,19
3 0,26
4 0,33
5 0,47
6 0,63
7 0,95
8 1,36
On veut modéliser l’évolution de la concentration x sous les conditions ”batch” par l’équation différentielle suivante :
dx
(
= µx
dt
x(0) = x0

où µ le taux de croissance de la biomasse est supposé constant au cours du temps et x0 est la concentration initiale (ce
sont des constantes strictement positives).
1. Lire le fichier de données chemostat.csv dans la data frame chemo. Tracer le nuage de points expérimentaux (ti , xi ).
> chemo <- read.csv2("chemostat.csv")
> attach(chemo)# pour simplifier les noms des variables
> plot(ti,xi,main="Croissance de la biomasse", xlab="temps (h)",ylab="biomasse g/l")

Croissance de la biomasse
1.4
1.0
biomasse g/l

0.6
0.2

0 2 4 6 8

temps (h)

dx
2. Résoudre l’équation différentielle (où µ est une constante): = µx(t) x(0) = x0 .
dt
x(t) = x0 exp(µt)
3. On veut déterminer le taux de croissance µ à partir des données (ti , xi ). Afin de ne pas donner une importance
exagérée à x0 , on considèrera que x0 est aussi un paramètre à déterminer (sinon la courbe passera exactement par
ce point (t0 = 0, x0 = 0.11) alors qu’elle ne passera qu’à proximité des autres).Quelle fonction S veut-on alors
minimiser ? Ce problème est-il linéaire par rapport au(x) paramètre(s) ?
On cherche le minimum de la fonction S définie par
n
X
S(x0 , µ) = (x0 exp(µti ) − xi )2
i=1

Le problème n’est pas linéaire: la recherche des points critiques ne conduit pas à un système linéaire. On va devoir
utiliser la fonction nls pour le résoudre. Mais cette fonction utilise un procédé itératif et nécessite donc d’avoir des
valeurs initiales des paramètres pour être mise en oeuvre.
4. Linéarisation du problème: recherche des valeurs initiales pour (x0 , µ).
(a) On définit z(t) = ln(x(t)). Montrer que z(t) = a1 + a2 t (ie z suit un modèle linéaire par rapport à t, une
droite). Exprimer µ et x0 en fonction de a1 et a2 .

z(t) = ln(x0 exp(µt)) = ln(x0 ) + µt

On a donc a1 = ln(x0 ) et a2 = µ et par conséquent x0 = exp(a1 ) et µ = a2 . On va maintenant chercher les
paramètres (a1 , a2 ) qui minimisent la fonction
n
X
Sl (a1 , a2 ) = (a1 + a2 ti − ln(xi ))2
i=1

avec la fonction lm puisque nous nous sommes ramenés à un problème linéaire.

(b) Déterminer les paramètres (a1 , a2 ) puis en déduire (x0 , µ).
> lin=lm(log(xi)~ti,data=chemo)
> summary(lin)
Call:
lm(formula = log(xi) ~ ti, data = chemo)

Residuals:
Min 1Q Median 3Q Max
-0.08847 -0.04460 -0.01712 0.05198 0.08861

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -2.259256 0.039685 -56.93 1.35e-10 ***
ti 0.309766 0.008336 37.16 2.66e-09 ***
---
Signif. codes:
0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.06457 on 7 degrees of freedom

Multiple R-squared: 0.995, Adjusted R-squared: 0.9942
F-statistic: 1381 on 1 and 7 DF, p-value: 2.656e-09
> a1=coef(lin)[1]
> a2=coef(lin)[2]
> x0=exp(a1)
> mu=a2
> print(c(a1,a2,x0,mu))
(Intercept) ti (Intercept) ti
-2.2592562 0.3097660 0.1044281 0.3097660
Affichage et contrôle des résultats
 > eq = paste0("Equation: ln(x) = ", round(a1,4)," + ", round(a2,4),"*t" )
> plot(ti,log(xi),main="Logarithme de la courbe croissance",sub=eq,xlab="temps (h)",ylab="log(biomas
> abline(lin,col="red" ,lwd = 2)

Logarithme de la courbe croissance

0.0
log(biomasse)

−1.0
−2.0

0 2 4 6 8

temps (h)
Equation: ln(x) = −2.2593 + 0.3098*t
On dispose donc de valeurs approchées de (x0 , µ) obtenues en linéarisant le problème. On revient à notre
problème initial.
5. Utilisation de la fonction nls: recherche des paramètres (x0 , µ).

(a) Définir une liste contenant les valeurs x0 et µ déterminées précedemment.

> initialisation=list(a=x0,b=mu)
Cette liste contient les valeurs de départ que nous allons utiliser avec la fonction nls pour déterminer plus
précisément (x0 , µ).
(b) Déterminer une estimation de (x0 , µ) par la méthode des moindres carrés à l’aide de la fonction nls (Nonlinear
Least Squares) de R studio. Afficher les résultats de nls avec la commande summary. Que signifie ”Number of
iterations to convergence ” ?
> modele <- nls(xi~ a*exp(b*ti), data = chemo, start = initialisation)
> summary(modele)
Formula: xi ~ a * exp(b * ti)

Parameters:
Estimate Std. Error t value Pr(>|t|)
a 0.085647 0.006086 14.07 2.17e-06 ***
b 0.343939 0.009952 34.56 4.40e-09 ***
---
Signif. codes:
0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.02583 on 7 degrees of freedom

Number of iterations to convergence: 4

Achieved convergence tolerance: 6.075e-06
Pour utiliser nls, il faut donner des valeurs de départ (start) pour les paramètres, valeurs contenues dans
initialisation. Le procédé est itératif: ”Number of iterations to convergence ” donne le nombre d’itérations
nécessaires pour converger vers une valeur acceptable.
(c) Tracer le nuage de points (ti , xi ) et la courbe théorique obtenue ainsi que son équation. Penser à contrôler les
résidus.
> x0=coef(modele)[1];mu=coef(modele)[2]
> eq = paste0("Equation: x = (", round(x0,4)," )*exp(", round(mu,4),"*t)" )# sous titre equation
> titre=paste0("détermination taux de croissance mu=",round(mu,4))#titre du graphique
> plot(ti,xi,main=titre,sub=eq,xlab="temps (h)",ylab="biomasse (g/L)")#tracé du nuage de points
> curve(x0*exp(mu*x),from=0,to=10,col="red" ,lwd = 2,add=TRUE)# ajout courbe prédite

détermination taux de croissance mu=0.3439

1.4
biomasse (g/L)

1.0
0.6
0.2

0 2 4 6 8

temps (h)
Equation: x = (0.0856 )*exp(0.3439*t) Contrôle des
résidus: on ne peut pas utiliser plot(modele) car cette fonction n’est implémentée que suite à l’utilisation de lm
(uniquement pour les modèles linéaires). Mais on peut tracer les résidus en fonction des valeurs prédites ou
en fonction du temps, vérifier la normalité des résidus,...

> layout(matrix(1:4,2,2))# fenetre graphique coupée en 4

> plot(ti,residuals(modele),main="residus en fonction du temps",xlab="temps (h)",ylab="résidus")
> plot(fitted(modele),residuals(modele),main="residus en fonction des valeurs prédites",xlab="valeur
> # manque de points mais les residus augmentent-ils avec t ou x=>hétéroscédasticité ?
> # residus normaux ? si les residus suivent une loi normale, les points doivent ^ etre alignés.
> qqnorm(residuals(modele))
> # si la droite passe par l'origine, les residus sont bien de moyenne nulle
> qqline(residuals(modele))
> # residus studentisés: non implémentés avec nls

residus en fonction du temps Normal Q−Q Plot

Sample Quantiles
résidus

−0.04

−0.04

0 2 4 6 8 −1.5 −1.0 −0.5 0.0 0.5 1.0 1.5

temps (h) Theoretical Quantiles

residus en fonction des valeurs prédites

résidus

−0.04

0.2 0.4 0.6 0.8 1.0 1.2

valeurs prédites
On ne peut pas
utiliser plot(modele) lorsque l’on a utilisé nls, ni rstudent. Beaucoup de fonctions ne sont pas implémentées
lorsque le problème est non linéaire. On remarque que les résidus ne semblent pas aléatoirement répartis (il
 serait souhaitable d’avoir davantage de points).
(d) Prédire les valeurs de x pour t = 0.5 à t = 7, 5h par pas de 1.
> newt<-seq(0.5,9.5,1)
> pc=predict(modele,data.frame(ti= newt), level = 0.95, interval = "confidence")
> print(pc)
[1] 0.1017181 0.1434726 0.2023671 0.2854373 0.4026072
[6] 0.5678746 0.8009830 1.1297807 1.5935475 2.2476872
La fonction predict donne uniquement les valeurs prédites. Les intervalles de confiance et de prédiction ne sont
pas implémentés en non linéaire (cf help).
> layout(1)
> plot( pc ~ newt, type='p',pch=3,main="Prédiction de la croissance de la biomasse", xlab="temps (h)

Prédiction de la croissance de la biomasse

2.0
biomasse g/l

1.5
1.0
0.5

2 4 6 8

temps (h)

Exercice 0.2 Détermination du Kla d’un fermenteur

La capacité d’oxygénation d’un fermenteur est déterminée en suivant la cinétique de transfert d’oxygène en absence
de microorganismes. On rappelle que la concentration en oxygène dissous, O2L , dans un fermenteur en absence de
microorganismes est donnée par l’équation différentielle
dO2L ∗
= Kla(O2L − O2L ) O2L (0) = O2L0
dt
∗
où Kla est le coefficient de transfert (gaz-liquide) de l’oxygène, supposé constant, O2L la concentration saturante d’oxygène
∗
liquide à 30o C. On suppose que O2L = 100. A l’aide de données expérimentales ci-dessous, on veut déterminer Kla et
O2L0 .

ti temps (s) 10 20 30 40 50 60 70 100 130

O2Li (mgl−1 ) 27 46 60 67,5 79 85 86 94 98
On veut déterminer les paramètres Kla et O2L0 .
∗
1. Résoudre l’équation différentielle vérifiée par O2L . (On montrera que O2L (t) = (O2L (0) − O2L )e−Klat + O2L
∗
).
2. Pour déterminer Kla et O2L0 , quelle fonction S1 cherche-t-on à minimiser ? On cherche à minimiser la fonction
n
X
S(Kla, O2L0 ) = (O2L (ti ) − O2Li )2
i=1

3. Déterminer des valeurs initiales Kla et O2L0 en linéarisant le problème.

 > # Création des données
> t <- c(10,20,30,40,50,60,70,100,130) # abscisse, variable explicative
> CL <- c(27,46,60,67.5,79,85,86,94,98) # ordonnée, variable expliquée
> # ******************************
> # détermination des valeurs initiales $Kla$ et $O_{2L_0}$ par linéarisation et régression
> donnees <- data.frame(t, CL)
> CS<-100;z=log(CS-CL)
> modlin <- lm(z~t, data = donnees)
> C00 <- -exp(coef(modlin)[1])+CS; kla0 <- -coef(modlin)[2]
> eq = paste0("Equation: ln(CS-CL)=", round(coef(modlin)[1],4),"+ ", -round(kla0,4),"*t" )# sous titre
> titre="Détermination Kla par linéarisation";#titre du graphique
> #tracé du nuage de points
> plot(t, z,main=titre,sub=eq,xlab="temps (sec)",ylab="ln(CS-CL)")# nuage
> abline(modlin,col="red")
>

Détermination Kla par linéarisation

4.0
3.0
ln(CS−CL)

2.0
1.0

20 40 60 80 100 120

temps (sec)
Equation: ln(CS−CL)=4.5802+ −0.0292*t

4. Déterminer Kla et O2L0 à l’aide de la fonction nls (Nonlinear Least Squares) de R studio.
> modele <- nls(CL~(C0-CS)*exp(-kla*t)+CS, data = donnees, start = list(C0 = C00, kla = kla0))
> b0 <- coef(modele)[1]
> kla <- coef(modele)[2]
> eq = paste0("Equation: CL =(", round(b0,4)," -100)*exp(-", round(kla,4),"*t)+100" )# sous titre
> titre="Détermination Kla";#titre du graphique
> #tracé du nuage de points
> plot(t, CL,main=titre,sub=eq,xlab="temps (sec)",ylab="Concentration Oxygène (mg/l)")# nuage
> curve((b0-CS)*exp(-kla*x)+CS,from=0,to=150,col="red" ,lwd = 2,add=TRUE)# ajout courbe prédite
 Détermination Kla
Concentration Oxygène (mg/l)

90
70
50
30

20 40 60 80 100 120

temps (sec)
Equation: CL =(2.3525 −100)*exp(−0.0293*t)+100

5. Etudier la qualité de la prédiction.

> layout(matrix(1:4,2,2))# fenetre graphique coupée en 4
> plot(donnees$t,residuals(modele))# résidus en fonction du temps
> plot(fitted(modele),residuals(modele))
> qqnorm(residuals(modele))
> qqline(residuals(modele))

Normal Q−Q Plot

residuals(modele)

Sample Quantiles
1

1
−2

−2

20 40 60 80 100 120 −1.5 −1.0 −0.5 0.0 0.5 1.0 1.5

donnees$t Theoretical Quantiles

residuals(modele)

1
−2

30 40 50 60 70 80 90 100

fitted(modele)