# ***************************************************
#      ETAPE de CALIBRATION: modèle de bioreacteur  *
# ***************************************************

rm(list=ls(all = TRUE)) 
set.seed(2)

## Packages a installer/charger
library(deSolve)
library(stats4)

## Chargement et affichage des donnees
data <- read.csv("batch.txt")
data[,"N"]=data[,"N"]/1000
data[,"dCO2dt"]=data[,"dCO2dt"]/100
Nbdata <- length(data[,"time"]) # nombre de données

## Tracé de N en fonction de B et de S en fonction de E
## -----------------------------------------------------
par(mfrow=c(1,2))
plot(data[,"X"],data[,"N"],xlab="X",ylab = "N",col="red")
plot(data[,"E"],data[,"S"],xlab="E",ylab = "S",col="red")

# Identification des paramètres k1 et k2 du modèle de fermentation en réacteur batch
# ----------------------------------------------------------------------------------

# identification de k1
# --------------------
fit <- lm(-data[,"N"] ~ data[,"X"])
c1_est <- fit$coefficients[1]
k1_est <- fit$coefficients[2]

par(mfrow=c(1,2))
plot(data[,"X"],data[,"N"],xlab="X",ylab = "N",col="red")
lines(data[,"X"],-k1_est*data[,"X"]-c1_est,col="red")
sprintf('estimation de k1 = %e',as.numeric(k1_est))

# identification de k2
# --------------------
fit <- lm(-data[,"S"] ~ data[,"E"])
c2_est <- fit$coefficients[1]
k2_est <- fit$coefficients[2]

plot(data[,"E"],data[,"S"],xlab="time",ylab = "X",col="red")
lines(data[,"E"],-k2_est*data[,"E"]-c2_est,col="red")
sprintf('estimation de k2 = %e',as.numeric(k2_est))

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

# Fonction cout à minimiser
cout <- function(parms,cas_identif,beta,data,k2_est) {
  
  # affichage des paramètres qu'on cherche à identifier
  print(parms)
  
  # récupération des valeurs des paramètres
  k1val <- as.numeric(parms["k1"])
  if (cas_identif == 1){
    k2val <- as.numeric(parms["k2"])}
  else
  {k2val <- k2_est}
  mu1maxval <- as.numeric(parms["mu1max"])
  mu2maxval <- as.numeric(parms["mu2max"])
  KNval <- as.numeric(parms["KN"])
  KEval <- as.numeric(parms["KE"])
  KSval <- as.numeric(parms["KS"])
  X0val <- as.numeric(parms["X0"])
  N0val <- as.numeric(parms["N0"])
  E0val <- as.numeric(parms["E0"])
  S0val <- as.numeric(parms["S0"])
  
  # stockage des valeurs des paramètres nécessaires pour la simulation du modèle
  param <- list(k1 = k1val, k2 = k2val, mu1max = mu1maxval,
                mu2max = mu2maxval, KN = KNval, KE = KEval,
                KS = KSval)
  
  times <- seq(from = 0, to = 90, by = 0.1)
  # condition initiale
  x0  <- c(X0val,N0val,E0val,S0val)
  # intégration numérique du modèle
  out <- ode(x0, data[,"time"], fermenteur, param)
  # on donne des noms aux colonnes
  colnames(out) <- c("time","X", "N", "E","S")

  # Calcul de la fonction coût à minimiser donnée par la somme
  # des carrés des écarts
  SumSquares <- as.numeric(beta["betaX"])*sum(((data[1:Nbdata-1,"X"]-out[1:Nbdata-1,"X"])**2
                           +(data[2:Nbdata,"X"]-out[2:Nbdata,"X"])**2)/2
                          *(data[2:Nbdata,"time"]-data[1:Nbdata-1,"time"]))
  SumSquares <- SumSquares+as.numeric(beta["betaN"])*sum(((data[1:Nbdata-1,"N"]-out[1:Nbdata-1,"N"])**2
                           +(data[2:Nbdata,"N"]-out[2:Nbdata,"N"])**2)/2
                          *(data[2:Nbdata,"time"]-data[1:Nbdata-1,"time"]))
  SumSquares <- SumSquares+as.numeric(beta["betaE"])*sum(((data[1:Nbdata-1,"E"]-out[1:Nbdata-1,"E"])**2
                           +(data[2:Nbdata,"E"]-out[2:Nbdata,"E"])**2)/2
                          *(data[2:Nbdata,"time"]-data[1:Nbdata-1,"time"]))
  SumSquares <- SumSquares+as.numeric(beta["betaS"])*sum(((data[1:Nbdata-1,"S"]-out[1:Nbdata-1,"S"])**2
                           +(data[2:Nbdata,"S"]-out[2:Nbdata,"S"])**2)/2
                          *(data[2:Nbdata,"time"]-data[1:Nbdata-1,"time"]))
  
  return(SumSquares)
}

# coefficients de pondération beta
betaXval <- 1/sum(((data[1:Nbdata-1,"X"])**2+(data[2:Nbdata,"X"])**2)/2
               *(data[2:Nbdata,"time"]-data[1:Nbdata-1,"time"]))
betaNval <- 1/sum(((data[1:Nbdata-1,"N"])**2+(data[2:Nbdata,"N"])**2)/2
               *(data[2:Nbdata,"time"]-data[1:Nbdata-1,"time"]))
betaEval <- 1/sum(((data[1:Nbdata-1,"E"])**2+(data[2:Nbdata,"E"])**2)/2
               *(data[2:Nbdata,"time"]-data[1:Nbdata-1,"time"]))
betaSval <- 1/sum(((data[1:Nbdata-1,"S"])**2+(data[2:Nbdata,"S"])**2)/2
               *(data[2:Nbdata,"time"]-data[1:Nbdata-1,"time"]))
beta <- c(betaX=betaXval, betaN=betaNval,betaE=betaEval,betaS=betaSval)

# Cas 1: on identifie tous les paramètres du modèle + les conditions initiales
cas_identif <-1
# Jeu de parametres de départ
StartList <- c(k1 = 0.01, k2 = 2.0, mu1max = 1.2, mu2max = 1.2,
               KN = 1.6, KE = 12., KS = 0.03, X0=data[1,"X"], N0=data[1,"N"],
               E0=data[1,"E"], S0=data[1,"S"])

Estim1 <- optim(StartList, ## Param de départ
               cout, ## fonction à optimiser
               gr=NULL,
               cas_identif,beta,data,k2_est,
               method= "L-BFGS-B", lower=c(0,0,0,0,0,0,0,0,0,0,0)
               ## algo choisi (CG ) gradient conjugué algo possible : 
               #method = c("Nelder-Mead", "BFGS", "CG", "L-BFGS-B", "SANN", "Brent"),
)

# Cas 2: on identifie tous les paramètres du modèle sauf k2 + les conditions initiales
cas_identif <-2
# Jeu de parametres de départ
StartList <- c(k1 = 0.01, mu1max = 1.2, mu2max = 1.2,
               KN = 1.6, KE = 12., KS = 0.03, X0=data[1,"X"], N0=data[1,"N"],
               E0=data[1,"E"], S0=data[1,"S"])

Estim2 <- optim(StartList, ## Param de départ
               cout, ## fonction à optimiser
               gr=NULL,
               cas_identif,beta,data,k2_est,
               method= "L-BFGS-B", lower=c(0,0,0,0,0,0,0,0,0,0)
               ## algo choisi (CG ) gradient conjugué algo possible : 
               #method = c("Nelder-Mead", "BFGS", "CG", "L-BFGS-B", "SANN", "Brent"),
)


# simulation du modèle identifié et comparaison aux données
# ---------------------------------------------------------
# vecteur de temps
times <- seq(from = 0, to = 90, by = 0.1)

# ** modéle identifié 1
# parametres du modele
parms1 <- list(
  k1 = Estim1$par["k1"],
  k2 = Estim1$par["k2"],
  mu1max = Estim1$par["mu1max"],
  mu2max = Estim1$par["mu2max"],
  KN = Estim1$par["KN"],
  KE = Estim1$par["KE"],
  KS = Estim1$par["KS"]
)

# condition initiale
x0  <- c(Estim1$par["X0"],Estim1$par["N0"],Estim1$par["E0"],Estim1$par["S0"])

# intégration numérique
out1 <- ode(x0, times, fermenteur, parms1)
# on renomme les colonnes de la matrice de sortie
colnames(out1) <- c("time","X", "N", "E", "S")

# ** modéle identifié 2
# parametres du modele
parms2 <- list(
  k1 = Estim2$par["k1"],
  k2 = k2_est,
  mu1max = Estim2$par["mu1max"],
  mu2max = Estim2$par["mu2max"],
  KN = Estim2$par["KN"],
  KE = Estim2$par["KE"],
  KS = Estim2$par["KS"]
)

# condition initiale
x0  <- c(Estim2$par["X0"],Estim2$par["N0"],Estim2$par["E0"],Estim2$par["S0"])

# intégration numérique
out2 <- ode(x0, times, fermenteur, parms2)
# on renomme les colonnes de la matrice de sortie
colnames(out2) <- c("time","X", "N", "E", "S")

# trace des solutions simulees et des données
par(mfrow=c(3,2))
plot(out1[,"time"],out1[,"X"],xlab="time",ylab="biomasse",type = 'l',col="magenta")
lines(out2[,"time"],out2[,"X"],xlab="time",type = 'l',col="blue")
points(data[,"time"],data[,"X"],xlab="time",col="red")
legend(50, 5, legend=c("identif 1","identif 2", "données"),
       col=c("magenta", "blue","red"),pch=c(NA,NA,1),lty=c(1,1,NA), cex=0.8)

plot(out1[,"time"],out1[,"N"],xlab="time",ylab="azote",type = 'l',col="magenta")
lines(out2[,"time"],out2[,"N"],xlab="time",type = 'l',col="blue")
points(data[,"time"],data[,"N"],xlab="time",col="red")

plot(out1[,"time"],out1[,"E"],xlab="time",ylab="ethanol",type = 'l',col="magenta")
lines(out2[,"time"],out2[,"E"],xlab="time",type = 'l',col="blue")
points(data[,"time"],data[,"E"],xlab="time",col="red")

plot(out1[,"time"],out1[,"S"],xlab="time",ylab="sucre",type = 'l',col="magenta")
lines(out2[,"time"],out2[,"S"],xlab="time",type = 'l',col="blue")
points(data[,"time"],data[,"S"],xlab="time",col="red")

dCO2_1 <- as.numeric(parms1["mu2max"])*out1[,"S"]/(as.numeric(parms1["KS"])+out1[,"S"])*as.numeric(parms1["KE"])/(as.numeric(parms1["KE"])+out1[,"E"])*out1[,"X"]
dCO2_2 <- as.numeric(parms2["mu2max"])*out2[,"S"]/(as.numeric(parms2["KS"])+out2[,"S"])*as.numeric(parms2["KE"])/(as.numeric(parms2["KE"])+out2[,"E"])*out2[,"X"]

plot(out1[,"time"],dCO2_1,xlab="time",ylab="dCO2/dt",type = 'l',col="magenta")
lines(out2[,"time"],dCO2_2,xlab="time",type = 'l',col="blue")
points(data[,"time"],data[,"dCO2dt"],xlab="time",col="red")