# *************************************************************** # ETAPE d'ANALYSE de SENSIBILITE du modèle de fermentation * # *************************************************************** rm(list=ls(all = TRUE)) set.seed(2) # chargement du package AS du logiciel R library(deSolve); library(sensitivity); ## Modele de fermenteur en réaction batch # --------------------------------------- fermenteur <- function(t, x, parms) { with(parms, { X=x[1];N=x[2];E=x[3];S=x[4]; # calcul de mu1(S) mu1 = mu1max*N/(KN+N) # calcul de mu2(E,S) mu2 = mu2max*S/(KS+S)*KE/(KE+E) # second membre de l'équation de B dX = mu1*X # second membre de l'équation de N dN = -k1*mu1*X # second membre de l'équation de E dE = mu2*X # second membre de l'équation de S dS = -k2*mu2*X return(list(c(dX,dN,dE,dS))) }) } ## Simulation du modèle et comparaison aux données réelles # ------------------------------------------------------- # parametres du modele parms0 <- list( k1 = 0.1, # - k2 = 1.9, # - mu1max = 1.2, # 1/h mu2max = 1.7, # 1/h KN = 1.6, # g/L KE = 10.0, # g/L KS = 0.01 # g/L ) # vecteur de temps times <- seq(from = 0, to = 180, by = 0.1) # condition initiale x0 <- c(0.01, 0.425, 0, 200) # intégration numérique out <- ode(x0, times, fermenteur, parms0) # on renomme les colonnes de la matrice de sortie colnames(out) <- c("time","X", "N", "E", "S") ## Chargement des données data <- read.csv("batch.txt") data[,"N"]=data[,"N"]/1000 data[,"dCO2dt"]=data[,"dCO2dt"]/100 # trace des solutions simulees et des données par(mfrow=c(2,2)) plot(out[,"time"],out[,"X"],xlab="time",ylab="",type = 'l',col="red",ylim=c(0,7.5)) points(data[,"time"],data[,"X"],xlab="time",col="red") lines(out[,"time"],out[,"N"]*10,xlab="time",type = 'l',col="green") points(data[,"time"],data[,"N"]*10,xlab="time",col="green") legend(50, 5, legend=c("Biomasse X (simulée)","Biomasse X (données)", "Azote N (x10) (simulée)", "Azote N (x10) (données)"), col=c("red", "red","green","green"),pch=c(NA,1,NA,1),lty=c(1,NA,1,NA), cex=0.8) plot(out[,"time"],out[,"S"],xlab="time",ylab="",type = 'l',col="black",ylim=c(0,210)) points(data[,"time"],data[,"S"],xlab="time",col="black") lines(out[,"time"],out[,"E"],xlab="time",type = 'l',col="blue") points(data[,"time"],data[,"E"],xlab="time",col="blue") legend(60,200, legend=c("Sucre S (simulée)","Sucre S (données)", "Ethanol E (simulée)","Ethanol E (données)"), col=c("black","black", "blue","blue"),pch=c(NA,1,NA,1),lty=c(1,NA,1,NA), cex=0.8) dCO2 <- as.numeric(parms0["mu2max"])*out[,"S"]/(as.numeric(parms0["KS"])+out[,"S"])*as.numeric(parms0["KE"])/(as.numeric(parms0["KE"])+out[,"E"])*out[,"X"] plot(out[,"time"],dCO2,xlab="time",ylab="",type = 'l',col="magenta",ylim=c(0,3.0)) points(data[,"time"],data[,"dCO2dt"],xlab="time",col="magenta") legend(60,2, legend=c("dCO2/dt (simulée)", "dCOE/dt (données)"), col=c("magenta", "magenta"), pch=c(NA,1),lty=c(1,NA), cex=0.8) ## Analyse de sensibilité # ------------------------ # fonction qui calcule les sorties d'intéret du systeme FermenteurAS <- function(parametre){ # on va regarder plusieurs sorties d'interêt: # 1. la valeur maximale de biomasse (levures) atteinte # 2. la valeur maximale d'ethanol atteinte # 3. la valeur maximale de vitesse de dégagement de CO2 atteinte # 4. la vitesse maximale de croissance des levures # 5. la vitesse maximale de production d'ethanol # 6. le temps que dure la fermentation (temps auquel le sucre est épuisé) # il y a donc 6 sorties et 7 parametres # parametre est une matrice dont chaque ligne contient un jeu de paramètres # c'est la matrice du plan d'experiences # Cette matrice a donc 6 colonnes, 1 par parametres et autant de ligne que # de jeux de paramètres à tester nb_param <- 7 nb_sorties <-6 # initialisation pour la sauvegarde sorties <- matrix(0, nrow=nrow(parametre), ncol=nb_sorties) # vecteur temps pour la simulation times <- seq(from = 0, to = 180, by = 0.1) # boucle des scenarios du plan d'experiences for (i in 1:nrow(parametre)) { # parametres du modele, mis a jour en fonction du scenario parms <- list( k1=parametre[i,1], k2=parametre[i,2], mu1max=parametre[i,3], mu2max=parametre[i,4], KN=parametre[i,5], KE=parametre[i,6], KS=parametre[i,7] ) # condition initiale x0 <- c(0.01, 0.425, 0, 200) # intégration numérique out <- ode(x0, times, fermenteur, parms) # on renomme les colonnes de la matrice de sortie colnames(out) <- c("time","X", "N", "E", "S") # calcul de la vitesse de dégagement du CO2 dCO2 <- as.numeric(parms["mu2max"])*out[,"S"]/(as.numeric(parms["KS"])+out[,"S"])*as.numeric(parms["KE"])/(as.numeric(parms["KE"])+out[,"E"])*out[,"X"] # vitesse de croissance des levures dX = as.numeric(parms["mu1max"])*out[,"N"]/(as.numeric(parms["KN"])+out[,"N"])*out[,"X"] # vitesse de production de l'ethanol dE = as.numeric(parms["mu2max"])*out[,"S"]/(as.numeric(parms["KS"])+out[,"S"])*as.numeric(parms["KE"])/(as.numeric(parms["KE"])+out[,"E"])*out[,"X"] # calcul des sorties d'interet # 1. valeur maximale de biomasse (levures) atteinte sorties[i,1] <- max(out[,"X"]) # 2. valeur maximale d'ethanol atteinte sorties[i,2] <- max(out[,"E"]) # 3. la valeur maximale de vitesse de dégagement de CO2 atteinte sorties[i,3] <- max(dCO2) # 4. vitesse maximale de croissance des levures sorties[i,4] <- max(dX) # 5. vitesse maximale de production d'ethanol sorties[i,5] <- max(dE) # 6. temps que dure la fermentation (temps auquel le sucre est épuisé) sorties[i,6] <- out[which.min(abs(out[,"S"]-10)),"time"] }# fin boucle scenarios du plan d'experiences return(sorties) } # fin fonction du modele nb_param <- 7 nb_sorties <-6 # DEFINITION DES GAMMES DE VALEURS pour chaque parametre # valeur nominale ValNominale <- unlist(parms0) # borne inferieure binfparam <- ValNominale*0.4 # borne superieure bsupparam <- ValNominale*1.6 # noms des parametres Paramname <- c("k1","k2","mu1max","mu2max","KN","KE","KS") # noms des sorties Sortiesname <- c("max X","max E","max dCO2/dt", "max dX/dt", "max dE/dt", "temps de fermentation") ### METHODE 1 : OAT standard # Plan d'experience PlanExp <- matrix(0, nrow=2*nb_param+1, ncol=nb_param) PlanExp[1,]=ValNominale # modification des valeurs parametre par parametre for (i in 1:nb_param){ PlanExp[2*i,]=ValNominale PlanExp[2*i+1,]=ValNominale PlanExp[2*i,i]=binfparam[i] PlanExp[2*i+1,i]=bsupparam[i] } # verification de l'echantillonage des parametres plot.new() par(mfrow=c(2,2)) for (i in 1:3){ plot(PlanExp[,2*(i-1)+1],PlanExp[,2*i],xlab=Paramname[2*(i-1)+1], ylab=Paramname[2*i]) # correspond au plan d'experience } plot(PlanExp[,nb_param-1],PlanExp[,nb_param],xlab=Paramname[nb_param-1], ylab=Paramname[nb_param]) # correspond au plan d'experience title("Plan d'experiences OAT standard", outer=TRUE) # simulations des scenarios SortiesOAT <- FermenteurAS(PlanExp); # graphiques des sorties contre les valeurs des parametres #sortie 1 par(mfrow=c(2,3),ask=F) for(i in 1:nb_param){ plot(PlanExp[,i], SortiesOAT[,1], xlab=Paramname[i], ylab=Sortiesname[1]) points(PlanExp[1,i],SortiesOAT[1,1],pch=21,col="red") } # sortie 2 par(mfrow=c(2,3),ask=F) for(i in 1:nb_param){ plot(PlanExp[,i], SortiesOAT[,2], xlab=Paramname[i], ylab=Sortiesname[2]) points(PlanExp[1,i],SortiesOAT[1,2],pch=21,col="red") } # Indices de sensibilite et d'elasticite # initialisation a 0 SensibP <- matrix(rep(0,nb_sorties*nb_param),nrow=nb_param,ncol=nb_sorties); # effet absolu d'une augmentation des parametres SensibM <- matrix(rep(0,nb_sorties*nb_param),nrow=nb_param,ncol=nb_sorties); # effet absolu d'une diminution des parametres ElastP <- matrix(rep(0,nb_sorties*nb_param),nrow=nb_param,ncol=nb_sorties); # effet relatif d'une diminution des parametres ElastM <- matrix(rep(0,nb_sorties*nb_param),nrow=nb_param,ncol=nb_sorties); # effet relatif d'une diminution des parametres # calcul des indices for (i in 1:nb_param){ # parametres for (j in 1:nb_sorties){ # sorties SensibP[i,j] <- (SortiesOAT[(2*i),j] - SortiesOAT[1,j])/(PlanExp[(2*i),i] - PlanExp[1,i]); # 1ere ligne de PAR = reference : parametres aux valeurs nominales SensibM[i,j] <- - (SortiesOAT[(2*i+1),j] - SortiesOAT[1,j])/(PlanExp[(2*i+1),i] - PlanExp[1,i]); # a mettre en <0 pour les graphes qui suivent ElastP[i,j] <- PlanExp[1,i]/SortiesOAT[1,j] * (SortiesOAT[(2*i),j] - SortiesOAT[1,j]) /(PlanExp[(2*i),i] - PlanExp[1,i]); ElastM[i,j] <- - PlanExp[1,i]/SortiesOAT[1,j] * (SortiesOAT[(2*i+1),j] - SortiesOAT[1,j]) /(PlanExp[(2*i+1),i] - PlanExp[1,i]); # a mettre en <0 pour les graphes qui suivent } } # trace de l'indice de sensibilite pour chaque sortie par(mfrow=c(2,2),ask=F) for (j in 1:nb_sorties) { xval <- c(); yval <- c(); for (i in 1:nb_param) { xval <- c(xval,i,i); yval <- c(yval,SensibP[i,j],SensibM[i,j]); } plot(xval, yval, type = "h", lwd=2, xaxt="n", xlab="", ylab="sensibilite", main=Sortiesname[j]) axis(1,at=c(1:nb_param),labels=Paramname) abline(0,0) } # trace de l'indice d'elasticite pour chaque sortie for (j in 1:nb_sorties) { xval <- c(); yval <- c(); for (i in 1:nb_param) { xval <- c(xval,i,i); yval <- c(yval,ElastP[i,j],ElastM[i,j]); } plot(xval, yval, type = "h", lwd=2, xaxt="n", xlab="", ylab="elasticite", main=Sortiesname[j]) axis(1,at=c(1:nb_param),labels=Paramname) abline(0,0) } ### METHODE 2 : MORRIS ?morris # ouvre l'aide sur la fonction morris # Plan d'experience (independant du modele) etudeMorris <- morris( model = FermenteurAS, factors = Paramname, r = 10, design = list(type = "oat", levels = 6, grid.jump = 3), scale=T, binf=binfparam, bsup=bsupparam) # levels : nb de valeurs par parametres # resume summary(etudeMorris) # verification de l'echantillonage des parametres par(mfrow=c(2,2)) for (i in 1:3){ plot(etudeMorris$X[,2*(i-1)+1],etudeMorris$X[,2*i],xlab=Paramname[2*(i-1)+1], ylab=Paramname[2*i]) # correspond au plan d'experience } plot(etudeMorris$X[,nb_param-1],etudeMorris$X[,nb_param],xlab=Paramname[nb_param-1], ylab=Paramname[nb_param]) # correspond au plan d'experience title("Plan d'experiences Morris", outer=TRUE) # sorties grapiques: Morris scatter plot mu* vs sigma op = par(mfrow=c(1,2)) for(i in 1:nb_sorties){ plot(etudeMorris,y_col=i) title(main=Sortiesname[i]) }