# ***************************************************************
#      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])
}