Projet Machin Learning — Couplage Modèle-Données & Réseaux de neurones

req_pkgs <- c(
  "tidyverse",  
  "scales",
  "patchwork",
  "nloptr",     
  "RSQLite",    
  "DBI",
  "progress",
  "neuralnet", 
  "Metrics"     
)

library(tidyverse)
library(scales)
library(patchwork)
library(nloptr)
library(DBI)
library(RSQLite)
library(progress)
library(neuralnet)
library(Metrics)

set.seed(123) 

# Paramètres globaux
N <- 40     # horizon temporel (n = 0..N)
M <- 100    # nombre de courbes

1 1. Introduction

1.1 Exercice 1.1 — Data-Model Coupling

1.1.1 1) Principe global

Le couplage modèle–données consiste à :

• définir un modèle mathématique (souvent déterministe) avec un petit nombre de paramètres,
• disposer de données observées (souvent bruitées),
• ajuster (calibrer/estimer) les paramètres du modèle en minimisant une mesure d’écart entre sorties du modèle et données.

Autrement dit, on utilise les données pour « verrouiller » les paramètres d’un modèle explicatif.

1.1.2 2) Intérêts et limites

Intérêts : - interprétabilité (peu de paramètres, sens physique/statistique), - extrapolation possible si la dynamique du modèle est pertinente, - besoin de moins de données qu’un modèle purement ML dans certains cas.

Limites : - biais de modèle : si le modèle est mal choisi, même un ajustement parfait est impossible, - paramètres parfois très sensibles → instabilité, non-identifiabilité, - optimisation potentiellement difficile (objectif non convexe, bruit, minima locaux).

2 2. Construction et analyse d’un modèle simple

On considère le modèle (Ex. 2.1) :

\[ U_{n+1} = a\,U_n^{c} + b,\quad n\in\{0,\dots,N-1\},\quad U_0 = U(0). \]

2.1 Exercice 2.1 — Construction du modèle

# 2.1) Fonction Generate_Output
# format de stockage choisi pour la série temporelle :
#un tibble avec colonnes: n, U

Generate_Output <- function(U0, N, a, b, c) {
  # Sécurité / validation légère
  stopifnot(is.numeric(U0), is.numeric(N), is.numeric(a), is.numeric(b), is.numeric(c))
  stopifnot(N >= 1)

  U <- numeric(N + 1)
  U[1] <- U0

  for (n in 1:N) {
    # U[n+1] = a * U[n]^c + b
    U[n + 1] <- a * (U[n]^c) + b
  }

  tibble(n = 0:N, U = U)
}

# Test rapide
Generate_Output(U0 = 1, N = 10, a = 3, b = 2, c = 0.5) %>% print(n = 12)

# A tibble: 11 × 2
       n     U
   <int> <dbl>
 1     0  1   
 2     1  5   
 3     2  8.71
 4     3 10.9 
 5     4 11.9 
 6     5 12.3 
 7     6 12.5 
 8     7 12.6 
 9     8 12.7 
10     9 12.7 
11    10 12.7 

2.1.1 4) Générer 6 courbes (N=40, a∈[2,5], b∈[1,5], c∈[0.1,0.9])

#choix de 6 triplets $(a,b,c)$ variés dans les plages demandées
param_grid_6 <- tribble(
  ~id, ~a,  ~b,  ~c,
  "A", 2.0, 1.0, 0.1,
  "B", 2.0, 5.0, 0.9,
  "C", 3.5, 1.0, 0.5,
  "D", 3.5, 5.0, 0.5,
  "E", 5.0, 1.0, 0.9,
  "F", 5.0, 5.0, 0.1
)

U0_demo <- 1

curves_6 <- param_grid_6 %>%
  mutate(curve = pmap(list(U0 = U0_demo, N = N, a = a, b = b, c = c), Generate_Output)) %>%
  unnest(curve)

ggplot(curves_6, aes(n, U, color = id)) +
  geom_line(linewidth = 1) +
  labs(
    title = "6 courbes du modèle U_{n+1} = a U_n^c + b",
    subtitle = paste0("U0 = ", U0_demo, ", N = ", N),
    x = "n", y = "U_n",
    color = "Courbe"
  ) +
  theme_minimal()

2.1.2 5) Influence qualitative des paramètres (commentaire)

• a : facteur multiplicatif sur le terme (U_n^c).
  Augmenter a accroît à la fois la croissance et l’amplitude des trajectoires.

• b : terme additif constant.
  Augmenter b translate la dynamique vers le haut (effet d’offset).

• c : paramètre contrôlant la non-linéarité de la dynamique.
  Lorsque (0 < c < 1), l’effet de (U_n) est atténué : les grandes valeurs sont comprimées, ce qui tend à ralentir la croissance.
  À l’inverse, lorsque (c > 1), les grandes valeurs de (U_n) sont amplifiées, ce qui peut conduire à une croissance plus rapide, voire à des comportements explosifs selon les valeurs de (a).

Phénomènes modélisables : dynamiques de croissance et relaxation non linéaires, processus où l’état futur dépend d’une puissance de l’état courant (avec forçage constant).

2.2 Exercice 2.2 — Sensibilité des paramètres

2.2.1 1) Notion

La sensibilité d’un paramètre mesure l’impact d’une petite variation de ce paramètre sur la sortie du modèle (ici la courbe (U_n)).

2.2.2 2–4) Approche + indicateur

On construit un indicateur simple : pour chaque paramètre \(p\in\{a,b,c\}\), on compare la courbe de référence à la courbe obtenue après une perturbation relative de \(+\delta\) (ex : 1%), puis on mesure une distance relative.

Ici on utilise une distance relative type :

\[ S(p)=\frac{1}{N+1}\sum_{n=0}^{N}\frac{|U_n^{(p+\delta)}-U_n^{ref}|}{|U_n^{ref}|+\varepsilon}. \]

# 2.2) Indicateur de sensibilité
rel_curve_diff <- function(U_ref, U_new, eps = 1e-8) {
  # U_ref et U_new : vecteurs numériques de même longueur
  stopifnot(length(U_ref) == length(U_new))
  mean(abs(U_new - U_ref) / (abs(U_ref) + eps))
}

sensitivity_indicator <- function(U0, N, a, b, c, delta = 0.01) {
  ref <- Generate_Output(U0, N, a, b, c)$U

  Ua <- Generate_Output(U0, N, a * (1 + delta), b, c)$U
  Ub <- Generate_Output(U0, N, a, b * (1 + delta), c)$U
  Uc <- Generate_Output(U0, N, a, b, c * (1 + delta))$U

  tibble(
    param = c("a", "b", "c"),
    S = c(
      rel_curve_diff(ref, Ua),
      rel_curve_diff(ref, Ub),
      rel_curve_diff(ref, Uc)
    )
  )
}

sens_demo <- sensitivity_indicator(U0 = 1, N = N, a = 3.5, b = 2.5, c = 0.5, delta = 0.01)
sens_demo

# A tibble: 3 × 2
  param       S
  <chr>   <dbl>
1 a     0.0141 
2 b     0.00263
3 c     0.0196 

ggplot(sens_demo, aes(param, S)) +
  geom_col() +
  labs(title = "Indicateur de sensibilité (perturbation +1%)", x = "Paramètre", y = "S(p)") +
  theme_minimal()

2.2.3 5) Conclusion (à partir des résultats obtenus)

Les valeurs obtenues (perturbation relative +1%) sont :

• \(S(a)=\text{0.01408}\)
• \(S(b)=\text{0.002626}\)
• \(S(c)=\text{0.01958}\)

On observe donc l’ordre de sensibilité : \[ c \;>\; a \;>>\; b. \]

Interprétation : - Le paramètre \(c\) (exposant) est le plus influent : une variation très faible de \(c\) modifie la dynamique (non-linéarité) à tous les instants, ce qui amplifie l’écart cumulé. - Le paramètre \(a\) a un effet important (gain multiplicatif), mais un peu moindre que \(c\) ici. - Le paramètre \(b\) est le moins sensible dans cette configuration : il agit comme un décalage additif, donc l’impact relatif (normalisé par \(|U_n|\)) est plus faible.

Conséquence pour la calibration : on s’attend à ce que l’estimation de \(c\) soit la plus délicate (instabilité / variance plus forte si les données sont bruitées), puis \(a\), tandis que \(b\) est en général plus robuste.

2.3 Exercice 2.3 — Lois de probabilité pour les paramètres

2.3.1 1) Inputs / outputs

• Inputs : \(U_0\), \(N\), et les paramètres \((a,b,c)\).
• Outputs : la série \((U_n)_{n=0,\dots,N}\) (la courbe).

2.3.2 2–4) Choix de lois

Choix raisonnable (simple, borné) :

• \(a\sim\mathcal{U}(2,5)\)
• \(b\sim\mathcal{U}(1,5)\)
• \(c\sim\mathcal{U}(0.1,0.9)\)
• \(U_0\sim\mathcal{U}(0.5,2)\) (exemple)

On pourra ajuster ces lois si les sorties explosent (voir Ex 3.2.1).

# 2.3) générateurs de paramètres
sample_params <- function(M) {
  tibble(
    a = runif(M, 2, 5),
    b = runif(M, 1, 5),
    c = runif(M, 0.1, 0.9),
    U0 = runif(M, 0.5, 2.0)
  )
}

params_demo <- sample_params(10000)

params_demo %>%
  pivot_longer(everything()) %>%
  ggplot(aes(value)) +
  geom_histogram(bins = 50) +
  facet_wrap(~name, scales = "free") +
  labs(title = "Histogrammes des lois choisies pour a, b, c, U0") +
  theme_minimal()

3 3. Tests de simulation

3.1 Exercice 3.1 — Notion de bruit

• Le bruit représente les erreurs de mesure, variabilité instrumentale, variabilité non modélisée, etc.
• Problèmes : sur-apprentissage, instabilité de calibration, biais des paramètres si le bruit n’est pas correctement pris en compte.
• Loi classique : bruit additif gaussien centré, indépendant :
  \(\varepsilon \sim \mathcal{N}(0,\sigma^2)\).

3.2 Exercice 3.2 — Génération d’une base d’apprentissage

3.2.1 1) Vérifier la stabilité relative des sorties

On définit un critère de “stabilité” pratique : pas d’explosion numérique (pas de Inf/NaN) et valeurs max raisonnables.

# 3.2) Génération de base de courbes + stabilité

is_stable_curve <- function(curve_U, max_allowed = 1e6) {
  all(is.finite(curve_U)) && max(abs(curve_U)) < max_allowed
}

generate_learning_base <- function(M, N, sigma = 0.2, max_allowed = 1e6, max_tries = 20000) {
  # Génère M courbes stables par rejet si nécessaire
  pb <- progress_bar$new(total = M, format = "  génération [:bar] :percent eta: :eta")

  keep <- vector("list", M)
  meta <- vector("list", M)

  i <- 1
  tries <- 0
  while (i <= M && tries < max_tries) {
    tries <- tries + 1
    p <- sample_params(1)
    cur <- Generate_Output(p$U0, N, p$a, p$b, p$c)

    if (is_stable_curve(cur$U, max_allowed = max_allowed)) {
      keep[[i]] <- cur$U
      meta[[i]] <- p
      pb$tick()
      i <- i + 1
    }
  }
  if (i <= M) stop("Impossible de générer M courbes stables avec les lois actuelles. Ajuste les lois ou le critère.")

  meta_df <- bind_rows(meta) %>% mutate(curve_id = row_number())
  output_curves <- do.call(rbind, keep) %>% as.matrix()
  colnames(output_curves) <- paste0("U_", 0:N)

  # Bruit
  noise <- matrix(rnorm(M * (N + 1), mean = 0, sd = sigma), nrow = M, ncol = N + 1)
  main_learning_base <- output_curves + noise

  list(
    meta = meta_df,
    Output_Curves = output_curves,
    Noise = noise,
    Main_Learning_Base = main_learning_base
  )
}

# Paramètres de bruit 
sigma <- 0.2

db <- generate_learning_base(M = M, N = N, sigma = sigma, max_allowed = 1e6)

dim(db$Main_Learning_Base)

[1] 100  41

3.2.2 2–5) Résultats + aperçu

# Aperçu de quelques courbes bruitées
set.seed(123)
sample_ids <- sample(1:M, 6)

df_plot <- tibble(curve_id = rep(sample_ids, each = N + 1),
                  n = rep(0:N, times = length(sample_ids)),
                  U = as.vector(t(db$Main_Learning_Base[sample_ids, ])))

ggplot(df_plot, aes(n, U, group = curve_id, color = factor(curve_id))) +
  geom_line(linewidth = 1) +
  labs(title = "Exemples de courbes bruitées (Main_Learning_Base)", color = "curve_id") +
  theme_minimal()

3.2.3 6) Construire une base SQL

On crée une base SQLite locale contenant : - une table meta (a,b,c,U0,curve_id), - une table curves_wide (curve_id + colonnes U_0..U_N bruitées), - éventuellement une table curves_long (format long utile pour ML).

# 3.2.6) Base SQL SQLite
sqlite_path <- "learning_base.sqlite"
if (file.exists(sqlite_path)) file.remove(sqlite_path)

[1] TRUE

con <- dbConnect(RSQLite::SQLite(), sqlite_path)

# meta
dbWriteTable(con, "meta", db$meta)

# courbes bruitées au format wide
curves_wide <- as_tibble(db$Main_Learning_Base) %>%
  mutate(curve_id = 1:M, .before = 1)

dbWriteTable(con, "curves_wide", curves_wide)

# format long 
curves_long <- curves_wide %>%
  pivot_longer(cols = starts_with("U_"), names_to = "n", values_to = "U") %>%
  mutate(n = as.integer(str_remove(n, "U_")))

dbWriteTable(con, "curves_long", curves_long)

dbDisconnect(con)

sqlite_path

[1] "learning_base.sqlite"

3.2.4 7) Split train/test (0.7M / 0.3M)

# 3.2.7) Split Train/Test
set.seed(123)
train_ids <- sample(1:M, size = floor(0.7 * M), replace = FALSE)
test_ids  <- setdiff(1:M, train_ids)

Main_Training_Base <- db$Main_Learning_Base[train_ids, , drop = FALSE]
Main_Test_Base     <- db$Main_Learning_Base[test_ids, , drop = FALSE]

meta_train <- db$meta %>% filter(curve_id %in% train_ids)
meta_test  <- db$meta %>% filter(curve_id %in% test_ids)

c(n_train = nrow(Main_Training_Base), n_test = nrow(Main_Test_Base))

n_train  n_test 
     70      30 

# ici on affiche des effectifs (nombre de courbes), nes pourcentages sont :
tibble(
  n_train = nrow(Main_Training_Base),
  n_test = nrow(Main_Test_Base),
  pct_train = nrow(Main_Training_Base) / M,
  pct_test = nrow(Main_Test_Base) / M
)

# A tibble: 1 × 4
  n_train n_test pct_train pct_test
    <int>  <int>     <dbl>    <dbl>
1      70     30       0.7      0.3

Remarque : ce découpage est un simple split train/test (hold-out), ce n’est pas une validation croisée.

3.3 Exercice 3.3 — Apprentissage des paramètres (calibration)

3.3.1 1) Intérêt

• Retrouver (a,b,c) à partir de données bruitées permet d’obtenir un modèle ajusté, interprétable, capable de prédire/reproduire des séries similaires.

3.3.2 2) Indicateur de différence relative entre deux courbes

On réutilise rel_curve_diff.

3.3.3 3–5) Fonctions RelDiff et f_obj

# 3.3) Fonctions d'écart et objectif d'optimisation

RelDiff <- function(a_t, b_t, c_t, U0, curve_obs, N) {
  # curve_obs : vecteur de longueur N+1 (observé/bruité)
  curve_pred <- Generate_Output(U0, N, a_t, b_t, c_t)$U

  # pour sécurité choisir `optim()` / L-BFGS-B exige une valeur FINIE pour la fonction objectif.
  # Si la simulation explose (Inf/NaN), on renvoie une pénalité finie.
  if (any(!is.finite(curve_pred))) return(1e9)

  d <- rel_curve_diff(curve_obs, curve_pred)
  if (!is.finite(d)) return(1e9)

  d
}

# objectif : moyenne des différences relatives sur tout le training set
f_obj <- function(theta, U0_list, curves_obs, N) {
  # theta = c(a, b, c)
  a_t <- theta[1]; b_t <- theta[2]; c_t <- theta[3]

  #pénalités simples si hors bornes (sécurité, en plus des bornes de l'algo)
  if (a_t < 0 || b_t < 0 || c_t <= 0) return(1e9)

  diffs <- map2_dbl(
    .x = U0_list,
    .y = seq_len(nrow(curves_obs)),
    .f = function(U0, i) RelDiff(a_t, b_t, c_t, U0, curves_obs[i, ], N)
  )
  if (any(!is.finite(diffs))) return(1e9)
  mean(diffs)
}

#test de l'objectif sur un theta arbitraire
U0_train <- meta_train$U0
test_theta <- c(3, 2, 0.5)
f_obj(test_theta, U0_train, Main_Training_Base, N)

[1] 0.5395052

3.3.4 4) Interprétation de mapply

Predicted_Curves <- mapply(Generate_Output, All_U0, N, a_t, b_t, c_t)
Cela génère une courbe prédite par valeur de U0, en utilisant les paramètres (a_t,b_t,c_t) identiques pour toutes les courbes.

3.3.5 6) But de f_obj

• f_obj est une fonction objectif (“objective function”) : on veut trouver \((a,b,c)\) minimisant l’erreur moyenne sur la base d’apprentissage.

3.3.6 7) DIRECT : principe

DIRECT (DIviding RECTangles) est une méthode déterministe de recherche globale (sans gradient) : elle découpe l’espace des paramètres en rectangles, évalue l’objectif, puis raffine les zones prometteuses. Utile pour des fonctions non convexes et/ou bruitées.

3.3.7 8) Ajustement par DIRECT (via nloptr)

On borne les paramètres dans des plages raisonnables :

• \(a\in[0,10]\)
• \(b\in[0,10]\)
• \(c\in[0.01,2]\)

# 3.3.8) Ajustement par DIRECT (nloptr)
lower <- c(a = 0, b = 0, c = 0.01)
upper <- c(a = 10, b = 10, c = 2)

direct_res <- nloptr(
  x0 = c(1, 1, 0.5),
  eval_f = function(theta) f_obj(theta, U0_train, Main_Training_Base, N),
  lb = lower,
  ub = upper,
  opts = list(
    algorithm = "NLOPT_GN_DIRECT_L",
    maxeval = 500,     # augmente si besoin
    print_level = 0
  )
)

direct_res$solution

[1] 2.6063947 3.8272452 0.3698743

direct_res$objective

[1] 0.5037756

Résultat obtenu (DIRECT, NLOPT_GN_DIRECT_L, maxeval=500)
- Solution : \(a=\text{2.6064}\), \(b=\text{3.8272}\), \(c=\text{0.36987}\)
- Valeur de l’objectif : \(f_{obj}=\text{0.503776}\)

Commentaire : DIRECT explore globalement l’espace des paramètres. Ici il trouve une solution donnant une erreur relative moyenne d’environ 0.50 sur le training. La valeur de l’objectif est cohérente avec celle observée ensuite par descente de gradient (objectif très proche), ce qui suggère que l’objectif possède au moins une “vallée” de solutions proches en qualité d’ajustement.

3.4 Exercice 3.4 — Stabilité des paramètres (bootstrap)

3.4.1 1) Table des 50 courbes tirées (exemple d’un tirage)

On veut construire une table contenant les 50 courbes tirées avec remise et les infos associées. Ci-dessous, on montre cela pour un tirage (itération 1) : identifiants, paramètres vrais \((a,b,c,U_0)\) et un extrait de la courbe observée bruitée.

set.seed(123)
idx_demo <- sample(seq_len(nrow(Main_Training_Base)), size = 50, replace = TRUE)

# table infos (sur les courbes tirées)
demo_info <- meta_train %>%
  mutate(row_in_train = row_number()) %>%
  filter(row_in_train %in% idx_demo) %>%
  select(curve_id, U0, a, b, c)

# Table courbes (format large, 50 x (N+1))
demo_curves <- as_tibble(Main_Training_Base[idx_demo, , drop = FALSE]) %>%
  mutate(sample_row = row_number(), .before = 1)

demo_info %>% slice_head(n = 10)

# A tibble: 10 × 5
   curve_id    U0     a     b     c
      <int> <dbl> <dbl> <dbl> <dbl>
 1        7 1.25   2.39  2.92 0.191
 2        9 1.58   3.38  4.63 0.542
 3       12 1.06   3.01  2.50 0.169
 4       14 1.89   2.44  4.35 0.544
 5       18 1.28   3.10  3.95 0.513
 6       21 1.53   2.74  2.75 0.702
 7       22 1.45   2.23  3.99 0.503
 8       23 1.96   4.39  3.14 0.695
 9       30 1.14   2.80  3.08 0.159
10       32 0.951  3.80  4.78 0.869

#petit extrait pour vérification (premières colonnes de 3 courbes)
demo_curves %>%
  slice_head(n = 3) %>%
  select(sample_row, U_0 = 1, U_1 = 2, U_2 = 3, U_3 = 4, U_4 = 5)

# A tibble: 3 × 5
    U_0   U_1   U_2   U_3   U_4
  <int> <dbl> <dbl> <dbl> <dbl>
1     1 0.501  5.39 10.4  13.6 
2     2 2.04   6.48  8.65  9.67
3     3 1.59   9.14 15.9  19.4 

On répète 10 fois : - tirer 50 courbes (avec remise) depuis le training set, - ajuster \((a,b,c)\) par DIRECT, - stocker les triplets.

# 3.4) Bootstrap de calibration
bootstrap_fit_direct <- function(B = 10, sample_size = 50, maxeval = 400) {
  res <- vector("list", B)
  for (k in 1:B) {
    idx <- sample(seq_len(nrow(Main_Training_Base)), size = sample_size, replace = TRUE)
    curves_k <- Main_Training_Base[idx, , drop = FALSE]
    U0_k <- U0_train[idx]

    fit_k <- nloptr(
      x0 = c(1, 1, 0.5),
      eval_f = function(theta) f_obj(theta, U0_k, curves_k, N),
      lb = lower,
      ub = upper,
      opts = list(
        algorithm = "NLOPT_GN_DIRECT_L",
        maxeval = maxeval,
        print_level = 0
      )
    )
    res[[k]] <- tibble(
      iter = k,
      a = fit_k$solution[1],
      b = fit_k$solution[2],
      c = fit_k$solution[3],
      obj = fit_k$objective
    )
  }
  bind_rows(res)
}

boot_tbl <- bootstrap_fit_direct(B = 10, sample_size = 50, maxeval = 400)
boot_tbl

# A tibble: 10 × 5
    iter     a     b     c   obj
   <int> <dbl> <dbl> <dbl> <dbl>
 1     1  1.66  4.44 0.428 0.493
 2     2  2.90  3.75 0.342 0.556
 3     3  3.81  2.53 0.342 0.532
 4     4  1.67  4.83 0.419 0.491
 5     5  2.74  3.72 0.366 0.491
 6     6  4.01  2.38 0.342 0.561
 7     7  3.04  3.50 0.342 0.543
 8     8  2.76  3.94 0.333 0.481
 9     9  2.90  3.70 0.317 0.499
10    10  2.90  3.19 0.309 0.472

# Statistiques demandées : moyenne, variance relative, écart-type relatif
summary_boot <- boot_tbl %>%
  summarise(
    a_mean = mean(a), a_var_rel = var(a) / (mean(a)^2), a_sd_rel = sd(a) / abs(mean(a)),
    b_mean = mean(b), b_var_rel = var(b) / (mean(b)^2), b_sd_rel = sd(b) / abs(mean(b)),
    c_mean = mean(c), c_var_rel = var(c) / (mean(c)^2), c_sd_rel = sd(c) / abs(mean(c))
  )

summary_boot

# A tibble: 1 × 9
  a_mean a_var_rel a_sd_rel b_mean b_var_rel b_sd_rel c_mean c_var_rel c_sd_rel
   <dbl>     <dbl>    <dbl>  <dbl>     <dbl>    <dbl>  <dbl>     <dbl>    <dbl>
1   2.84    0.0708    0.266   3.60    0.0446    0.211  0.354    0.0126    0.112

3.4.2 5) Description des résultats et intérêt (à partir des sorties)

Les 10 ajustements DIRECT (sur 50 courbes tirées avec remise) donnent :

• Moyennes : \(\bar a=\text{2.84}\), \(\bar b=\text{3.6}\), \(\bar c=\text{0.354}\).
• Écart-type relatif :
  • \(\mathrm{sd}_{rel}(a)=\text{0.266}\) (≈ 26.6%)
    
  • \(\mathrm{sd}_{rel}(b)=\text{0.211}\) (≈ 21.1%)
    
  • \(\mathrm{sd}_{rel}(c)=\text{0.112}\) (≈ 11.2%)

Lecture : - Le paramètre \(c\) apparaît le plus stable (dispersion relative la plus faible). - \(a\) est le plus variable parmi les trois (dispersion relative la plus forte), puis \(b\).

Intérêt de l’approche bootstrap : - Elle donne une mesure empirique de l’incertitude (variabilité) des paramètres estimés due au bruit et au fait qu’on n’observe qu’un sous-échantillon de courbes. - Elle permet de choisir une valeur “robuste” (ici la moyenne) pour figer \[(a,b,c)\] ensuite (Ex 3.5).

3.5 Exercice 3.5 — Précision du modèle ajusté

3.5.1 1) Choix d’indicateurs

On utilise : - RMSE (erreur quadratique moyenne) : pénalise fortement les gros écarts, - MAE (erreur absolue moyenne) : robuste, - Différence relative moyenne (déjà définie).

# 3.5) Indicateurs d'ajustement
predict_curves_matrix <- function(U0_list, N, a, b, c) {
  # Renvoie une matrice (len(U0_list) x (N+1)) des courbes prédites
  preds <- map(U0_list, ~Generate_Output(.x, N, a, b, c)$U)
  do.call(rbind, preds)
}

curve_metrics <- function(curves_obs, curves_pred) {
  stopifnot(all(dim(curves_obs) == dim(curves_pred)))

  rmse_i <- apply((curves_obs - curves_pred)^2, 1, mean) %>% sqrt()
  mae_i  <- apply(abs(curves_obs - curves_pred), 1, mean)
  rdiff_i <- map2_dbl(
    split(curves_obs, row(curves_obs)),
    split(curves_pred, row(curves_pred)),
    ~rel_curve_diff(.x, .y)
  )

  tibble(
    RMSE_mean = mean(rmse_i),
    MAE_mean = mean(mae_i),
    RelDiff_mean = mean(rdiff_i)
  )
}

#les paramètres fixés, on prend la moyenne bootstrap
a_hat <- summary_boot$a_mean
b_hat <- summary_boot$b_mean
c_hat <- summary_boot$c_mean

pred_train <- predict_curves_matrix(U0_train, N, a_hat, b_hat, c_hat)
metrics_train <- curve_metrics(Main_Training_Base, pred_train)
metrics_train

# A tibble: 1 × 3
  RMSE_mean MAE_mean RelDiff_mean
      <dbl>    <dbl>        <dbl>
1     1145.     943.        0.504

U0_test <- meta_test$U0
pred_test <- predict_curves_matrix(U0_test, N, a_hat, b_hat, c_hat)
metrics_test <- curve_metrics(Main_Test_Base, pred_test)
metrics_test

# A tibble: 1 × 3
  RMSE_mean MAE_mean RelDiff_mean
      <dbl>    <dbl>        <dbl>
1      492.     407.        0.459

# Visualiser quelques fits (train + test)
plot_fit_examples <- function(curves_obs, U0_list, a, b, c, title_prefix, n_examples = 4) {
  ids <- sample(seq_len(nrow(curves_obs)), n_examples)
  obs <- curves_obs[ids, , drop = FALSE]
  pred <- predict_curves_matrix(U0_list[ids], N, a, b, c)

  df <- tibble(
    curve = rep(seq_len(n_examples), each = N + 1),
    n = rep(0:N, times = n_examples),
    obs = as.vector(t(obs)),
    pred = as.vector(t(pred))
  ) %>%
    pivot_longer(cols = c(obs, pred), names_to = "type", values_to = "U")

  ggplot(df, aes(n, U, linetype = type)) +
    geom_line(linewidth = 1) +
    facet_wrap(~curve, scales = "free_y") +
    labs(title = paste0(title_prefix, " — exemples (obs vs pred)"), linetype = "") +
    theme_minimal()
}

set.seed(123)

p1 <- plot_fit_examples(Main_Training_Base, U0_train, a_hat, b_hat, c_hat, "Training")
p2 <- plot_fit_examples(Main_Test_Base, U0_test, a_hat, b_hat, c_hat, "Test")
p1 / p2

3.5.2 4) Description des résultats (à partir des sorties)

En fixant \((a,b,c)\) aux moyennes bootstrap (\(\bar a=\text{2.84}\), \(\bar b=\text{3.6}\), \(\bar c=\text{0.354}\)) :

• Training : RMSE ≈ 1145, MAE ≈ 943, RelDiff ≈ 0.504
  
• Test : RMSE ≈ 492, MAE ≈ 407, RelDiff ≈ 0.459

Interprétation : - L’erreur relative moyenne est de l’ordre de 0.46–0.50, ce qui est cohérent avec la valeur de l’objectif obtenue lors des optimisations (≈ 0.50). - Les erreurs absolues (RMSE/MAE) sont plus grandes sur le training que sur le test : cela arrive si les courbes du training contiennent des amplitudes plus élevées (RMSE/MAE dépendent fortement de l’échelle). - Le fait que la RelDiff soit légèrement plus faible sur le test suggère qu’il n’y a pas de sur-ajustement évident ; on reste dans un régime de généralisation comparable.

3.6 Exercice 3.6 — Comparaison avec un algorithme de descente de gradient

3.6.1 1) Algorithme choisi

On utilise optim() avec méthode L-BFGS-B : - c’est une méthode quasi-Newton (utilise un gradient approché), - gère des bornes sur les paramètres, - dépend de l’initialisation et peut converger vers des minima locaux

3.6.2 2) Ajuster avec initialisation à 0

# 3.6) Optim gradient (L-BFGS-B)
gd_fit <- optim(
  par = c(0, 0, 0.01), # on évite c=0 exact
  fn = function(theta) f_obj(theta, U0_train, Main_Training_Base, N),
  method = "L-BFGS-B",
  lower = lower,
  upper = upper,
  control = list(maxit = 200)
)

gd_fit$par

[1] 3.1254426 3.2106642 0.3333344

gd_fit$value

[1] 0.5037524

3.6.3 3) 5 initialisations différentes

inits <- tribble(
  ~a0, ~b0, ~c0,
  0,   0,   0.01,
  1,   1,   0.5,
  5,   5,   0.5,
  2,   8,   1.2,
  8,   2,   0.2
)

gd_runs <- pmap_dfr(inits, function(a0, b0, c0) {

  fit <- tryCatch(
    optim(
      par = c(a0, b0, c0),
      fn = function(theta) f_obj(theta, U0_train, Main_Training_Base, N),
      method = "L-BFGS-B",
      lower = lower,
      upper = upper,
      control = list(maxit = 300)
    ),
    error = function(e) NULL
  )

  #si optim échoue (fn non finie,...), on renvoie une ligne pénalisée
  if (is.null(fit) || !is.finite(fit$value)) {
    return(tibble(
      a0 = a0, b0 = b0, c0 = c0,
      a = NA_real_, b = NA_real_, c = NA_real_,
      obj = 1e9, conv = NA_integer_
    ))
  }

  tibble(
    a0 = a0, b0 = b0, c0 = c0,
    a = fit$par[1], b = fit$par[2], c = fit$par[3],
    obj = fit$value, conv = fit$convergence
  )
})

gd_runs

# A tibble: 5 × 8
     a0    b0    c0     a     b     c     obj  conv
  <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>   <dbl> <int>
1     0     0  0.01  3.13 3.21  0.333 5.04e-1     0
2     1     1  0.5   3.12 3.21  0.334 5.04e-1    52
3     5     5  0.5   2.56 3.77  0.383 5.04e-1     0
4     2     8  1.2   2    8     1.2   1   e+9     0
5     8     2  0.2   6.05 0.245 0.205 5.04e-1    52

# A noter au niveau robustesse des statistiques
#certaines initialisations peuvent échouer (objectif pénalisé / convergence anormale)
# pour éviter de biaiser la synthèse, on calcule aussi des statistiques sur les runs valides
gd_runs_valid <- gd_runs %>% filter(is.finite(obj), obj < 1e8, conv == 0)
gd_runs_valid

# A tibble: 2 × 8
     a0    b0    c0     a     b     c   obj  conv
  <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <int>
1     0     0  0.01  3.13  3.21 0.333 0.504     0
2     5     5  0.5   2.56  3.77 0.383 0.504     0

3.6.4 4) Moyenne, variance relative, écart-type relatif

# Synthèse "brute" (sur les 5 tests, comme demandé)
gd_summary_all <- gd_runs %>%
  summarise(
    a_mean = mean(a, na.rm = TRUE),
    a_var_rel = var(a, na.rm = TRUE) / (mean(a, na.rm = TRUE)^2),
    a_sd_rel  = sd(a, na.rm = TRUE) / abs(mean(a, na.rm = TRUE)),

    b_mean = mean(b, na.rm = TRUE),
    b_var_rel = var(b, na.rm = TRUE) / (mean(b, na.rm = TRUE)^2),
    b_sd_rel  = sd(b, na.rm = TRUE) / abs(mean(b, na.rm = TRUE)),

    c_mean = mean(c, na.rm = TRUE),
    c_var_rel = var(c, na.rm = TRUE) / (mean(c, na.rm = TRUE)^2),
    c_sd_rel  = sd(c, na.rm = TRUE) / abs(mean(c, na.rm = TRUE))
  )

gd_summary_all

# A tibble: 1 × 9
  a_mean a_var_rel a_sd_rel b_mean b_var_rel b_sd_rel c_mean c_var_rel c_sd_rel
   <dbl>     <dbl>    <dbl>  <dbl>     <dbl>    <dbl>  <dbl>     <dbl>    <dbl>
1   3.37     0.217    0.465   3.69     0.568    0.754  0.491     0.670    0.818

# Synthèse "valide", on retire les runs pénalisés / non convergés)
gd_summary_valid <- gd_runs_valid %>%
  summarise(
    a_mean = mean(a), a_var_rel = var(a) / (mean(a)^2), a_sd_rel = sd(a) / abs(mean(a)),
    b_mean = mean(b), b_var_rel = var(b) / (mean(b)^2), b_sd_rel = sd(b) / abs(mean(b)),
    c_mean = mean(c), c_var_rel = var(c) / (mean(c)^2), c_sd_rel = sd(c) / abs(mean(c))
  )

gd_summary_valid

# A tibble: 1 × 9
  a_mean a_var_rel a_sd_rel b_mean b_var_rel b_sd_rel c_mean c_var_rel c_sd_rel
   <dbl>     <dbl>    <dbl>  <dbl>     <dbl>    <dbl>  <dbl>     <dbl>    <dbl>
1   2.84    0.0201    0.142   3.49    0.0130    0.114  0.358   0.00959   0.0979

Il est important de noter que certaines initialisations peuvent conduire à une évaluation pénalisée ou à une convergence non standard. Ces runs sont exclus de l’analyse afin d’éviter de biaiser les statistiques et la comparaison des performances.

On a donc calculer : - une synthèse “brute” sur les 5 runs (conforme au tableau demandé), - puis une synthèse “sur runs valides” (pour juger l’algorithme).

3.6.5 5) Précision avec paramètres moyens

# on utilise les paramètres moyens issus des runs "valides" pour évaluer correctement le gradient
a_gd <- gd_summary_valid$a_mean
b_gd <- gd_summary_valid$b_mean
c_gd <- gd_summary_valid$c_mean

pred_train_gd <- predict_curves_matrix(U0_train, N, a_gd, b_gd, c_gd)
pred_test_gd  <- predict_curves_matrix(U0_test,  N, a_gd, b_gd, c_gd)

metrics_train_gd <- curve_metrics(Main_Training_Base, pred_train_gd)
metrics_test_gd  <- curve_metrics(Main_Test_Base, pred_test_gd)

metrics_train_gd

# A tibble: 1 × 3
  RMSE_mean MAE_mean RelDiff_mean
      <dbl>    <dbl>        <dbl>
1     1145.     943.        0.504

metrics_test_gd

# A tibble: 1 × 3
  RMSE_mean MAE_mean RelDiff_mean
      <dbl>    <dbl>        <dbl>
1      492.     407.        0.458

3.6.6 6) Comparaison DIRECT vs gradient

DIRECT (NLOPT_GN_DIRECT_L) :
- Solution : \((a,b,c)=(\text{2.6064},\;\text{3.8272},\;\text{0.36987})\)
- Objectif : 0.503776

Gradient (L-BFGS-B) :
- Initialisation à 0 : \((a,b,c)=(\text{3.1254},\;\text{3.2107},\;\text{0.33333})\), objectif 0.503752
- Sur 5 initialisations : certaines initialisations échouent (pénalité 1e9) ou convergent de façon non standard ; on interprète donc surtout les runs valides.

Lecture : plusieurs triplets peuvent donner une valeur d’objectif très proche (vallée de minima), donc le problème peut être partiellement non-identifiable. DIRECT est plus robuste à l’initialisation ; le gradient est plus rapide mais plus sensible.

3.7 Exercice 3.7 — Conclusion de l’approche

1. On construit un modèle, on définit des lois sur les paramètres, on simule une base bruitée, puis on calibre (a,b,c) en minimisant un écart entre prédictions et observations.
   
2. Comparaison à une approche “classique” : ici on apprend des paramètres du modèle (approche mécanistique), plutôt que d’apprendre directement une fonction prédictive sans structure (ML).

4 4. Approche Machine Learning

4.1 Exercice 4.1 — Principe

• ML : apprendre une relation entrée→sortie à partir des données, en optimisant une loss.
• Intérêts : flexible, capture non-linéarités complexes.
• Limites : moins interprétable, demande souvent plus de données, risque d’overfitting.

Réseau de neurones : composition de couches linéaires + non-linéarités, apprentissage par rétropropagation. Le nombre de neurones et le nombre de couches d’un réseau de neurones déterminent sa capacité de représentation. Un réseau comportant trop peu de neurones ou de couches peut conduire à un sous-apprentissage, car il n’est pas suffisamment expressif pour capturer les relations non linéaires présentes dans les données. À l’inverse, un réseau trop profond ou trop large peut provoquer un sur-apprentissage, en s’ajustant excessivement aux données d’apprentissage, et entraîner une dégradation des performances sur les données de test, ainsi qu’un coût de calcul plus élevé. Un compromis doit donc être trouvé, généralement en analysant l’écart de performance entre les jeux d’apprentissage et de test.

4.2 Exercice 4.2 — Ajustement par réseau de neurones

4.2.1 1) Construire la table Learn_Base (format long)

la table demandée est de colonnes (n_curve, n, U0, TS) où TS = valeur de la série à l’instant n.

# 4.2.1) Learn_Base à partir de Main_Learning_Base (bruité)
Learn_Base <- tibble(
  n_curve = rep(1:M, each = N),
  n = rep(1:N, times = M),
  U0 = rep(db$meta$U0, each = N),
  TS = as.vector(t(db$Main_Learning_Base[, 2:(N + 1)])) # TS_{n,curve} = U_n, n=1..N
)

Learn_Base %>% slice_head(n = 10)

# A tibble: 10 × 4
   n_curve     n    U0    TS
     <int> <int> <dbl> <dbl>
 1       1     1  1.21  7.66
 2       1     2  1.21 11.2 
 3       1     3  1.21 12.7 
 4       1     4  1.21 12.2 
 5       1     5  1.21 12.2 
 6       1     6  1.21 12.5 
 7       1     7  1.21 12.3 
 8       1     8  1.21 12.5 
 9       1     9  1.21 12.6 
10       1    10  1.21 12.4 

4.2.2 2) Split Train_Base_1 / Test_Base_1

L’enoncé ici inverse les proportion (0.3M “first”, 0.7M “last”).
On suit exactement : - Train_Base_1 = 0.3M premières courbes (donc 0.3M*(N) lignes), - Test_Base_1 = 0.7M dernières courbes.

# Ce split est DIFFERENT du split 70/30 précédent.
M_train1 <- floor(0.3 * M)
train1_curves <- 1:M_train1
test1_curves  <- (M_train1 + 1):M

Train_Base_1 <- Learn_Base %>% filter(n_curve %in% train1_curves)
Test_Base_1  <- Learn_Base %>% filter(n_curve %in% test1_curves)

c(nrow_train1 = nrow(Train_Base_1), nrow_test1 = nrow(Test_Base_1))

nrow_train1  nrow_test1 
       1200        2800 

4.2.3 3) Standardisation (fonction S)

On code la fonction telle que décrite : on centre-réduit par la moyenne de Train et l’amplitude (range) de Train.

# 4.2.3) Fonction de standardisation
S <- function(TrainB, TestB) {
  TrainB_S <- matrix(0, nrow(TrainB), 4)
  TestB_S  <- matrix(0, nrow(TestB), 4)

  TrainB_S[,1] <- TrainB[,1]
  TrainB_S[,2] <- TrainB[,2]
  TrainB_S[,3] <- TrainB[,3]
  TrainB_S[,4] <- TrainB[,4]

  TestB_S[,1] <- TestB[,1]
  TestB_S[,2] <- TestB[,2]
  TestB_S[,3] <- TestB[,3]
  TestB_S[,4] <- TestB[,4]

  for (i in 2:4) {
    TrainB_S[,i] <- (TrainB[,i] - mean(TrainB[,i])) / diff(range(TrainB[,i]))
    TestB_S[,i]  <- (TestB[,i]  - mean(TrainB[,i])) / diff(range(TrainB[,i]))
  }

  list(TrainB_S, TestB_S)
}

# Application
Train_mat <- as.matrix(Train_Base_1 %>% select(n_curve, n, U0, TS))
Test_mat  <- as.matrix(Test_Base_1  %>% select(n_curve, n, U0, TS))

Learn_Base_S <- S(Train_mat, Test_mat)
Train_Base_1_S <- as.data.frame(Learn_Base_S[[1]])
Test_Base_1_S  <- as.data.frame(Learn_Base_S[[2]])

colnames(Train_Base_1_S) <- c("n_curve","n","U0","TS")
colnames(Test_Base_1_S)  <- c("n_curve","n","U0","TS")

head(Train_Base_1_S)

  n_curve          n          U0          TS
1       1 -0.5000000 -0.01765081 -0.02915330
2       1 -0.4743590 -0.01765081 -0.02905984
3       1 -0.4487179 -0.01765081 -0.02902096
4       1 -0.4230769 -0.01765081 -0.02903227
5       1 -0.3974359 -0.01765081 -0.02903360
6       1 -0.3717949 -0.01765081 -0.02902549

4.2.4 4) Entraînement d’un réseau de neurones (3 couches cachées: 3,6,3)

On modélise TS ~ n + U0 (deux entrées).
On utilise neuralnet::neuralnet.

# 4.2.4) Entraînement réseau de neurones
set.seed(123)

# neuralnet préfère des variables numériques simples
Train_nn <- Train_Base_1_S %>% mutate(n = as.numeric(n), U0 = as.numeric(U0), TS = as.numeric(TS))
Test_nn  <- Test_Base_1_S  %>% mutate(n = as.numeric(n), U0 = as.numeric(U0), TS = as.numeric(TS))

nn_time <- system.time({
  nn_fit <- neuralnet(
    formula = TS ~ n + U0,
    data = Train_nn,
    hidden = c(3, 6, 3),
    linear.output = TRUE,      # sortie continue
    lifesign = "minimal",
    stepmax = 1e6
  )
})

nn_time

utilisateur     système      écoulé 
       1.50        0.20        2.03 

4.2.5 5) Diagramme + nombre de paramètres

# Diagramme
plot(nn_fit)

Nombre de paramètres (poids + biais) pour une architecture 2-3-6-3-1 : - couche 1 : (2 + 1) * 3 = 9 - couche 2 : (3 + 1) * 6 = 24 - couche 3 : (6 + 1) * 3 = 21 - sortie : (3 + 1) * 1 = 4
Total = 58 paramètres.

Comparaison : - Couplage modèle-données : 3 paramètres \((a,b,c)\), - Réseau : 58 paramètres → plus flexible mais plus “boîte noire”.

4.3 Exercice 4.3 — Étude du modèle ajusté

4.3.1 1) Courbes observées vs prédites (train/test)

On prédit TS point par point, puis on reconstitue des courbes par n_curve.

# 4.3.1) Prédictions NN + figures

nn_pred_train <- neuralnet::compute(nn_fit, Train_nn %>% select(n, U0))$net.result %>% as.vector()
nn_pred_test  <- neuralnet::compute(nn_fit, Test_nn  %>% select(n, U0))$net.result %>% as.vector()

Train_nn_pred <- Train_nn %>% mutate(TS_pred = nn_pred_train)
Test_nn_pred  <- Test_nn  %>% mutate(TS_pred = nn_pred_test)

plot_nn_curves <- function(df, title_prefix, n_curves = 4) {
  ids <- sample(unique(df$n_curve), n_curves)
  dsub <- df %>% filter(n_curve %in% ids) %>%
    select(n_curve, n, TS, TS_pred) %>%
    pivot_longer(cols = c(TS, TS_pred), names_to = "type", values_to = "value")

  ggplot(dsub, aes(n, value, linetype = type)) +
    geom_line(linewidth = 1) +
    facet_wrap(~n_curve, scales = "free_y") +
    labs(title = paste0(title_prefix, " — courbes (obs vs NN)"), linetype = "") +
    theme_minimal()
}

set.seed(123)

plot_nn_curves(Train_nn_pred, "Train_Base_1_S")

plot_nn_curves(Test_nn_pred, "Test_Base_1_S")

4.3.2 2) Indicateurs de précision (mêmes que Ex 3.5)

Ici on compare TS point à point (sur la base long), ce qui est cohérent car les indicateurs précédents étaient par courbe.

# 4.3.2) Indicateurs pour NN (mêmes indicateurs que Ex 3.5, calculés par courbe)
nn_curve_metrics <- function(df, eps = 1e-8) {
  per_curve <- df %>%
    group_by(n_curve) %>%
    summarise(
      RMSE = sqrt(mean((TS - TS_pred)^2)),
      MAE  = mean(abs(TS - TS_pred)),
      RelDiff = mean(abs(TS - TS_pred) / (abs(TS) + eps)),
      .groups = "drop"
    )

  tibble(
    RMSE_mean = mean(per_curve$RMSE),
    MAE_mean = mean(per_curve$MAE),
    RelDiff_mean = mean(per_curve$RelDiff)
  )
}

m_train_nn <- nn_curve_metrics(Train_nn_pred)
m_test_nn  <- nn_curve_metrics(Test_nn_pred)

m_train_nn

# A tibble: 1 × 3
  RMSE_mean MAE_mean RelDiff_mean
      <dbl>    <dbl>        <dbl>
1    0.0203   0.0176        0.536

m_test_nn

# A tibble: 1 × 3
  RMSE_mean MAE_mean RelDiff_mean
      <dbl>    <dbl>        <dbl>
1     0.105   0.0896         2.80

4.3.3 Résultats et commentaires

Sur les données standardisées (cohérent car le réseau a été entraîné sur cette échelle) :

• Train_Base_1_S : RMSE_mean = 0.0203, MAE_mean = 0.0176, RelDiff_mean = 0.5355
• Test_Base_1_S : RMSE_mean = 0.1051, MAE_mean = 0.0896, RelDiff_mean = 2.8042

Commentaire : - Les erreurs sont faibles sur le training : le réseau approxime très bien la relation \((n,U_0)\mapsto TS\) sur l’échantillon d’apprentissage. - Les erreurs augmentent sur le test (ce qui est normal), mais restent modérées : le réseau généralise globalement. - Le saut train→test peut s’expliquer par le split “0.3M first / 0.7M last” qui n’est pas aléatoire : la distribution des courbes en test peut différer de celle du train.

4.3.4 3) Indicateurs sur Train_Base_2_S / Test_Base_2_S (split 70/30 aléatoire)

On définit ici une Base 2 cohérente avec le split aléatoire 70/30 utilisé dans l’approche de couplage (Ex 3.2.7) :
les courbes d’identifiants train_ids sont dans Train_Base_2, et celles d’identifiants test_ids dans Test_Base_2.

# 4.3.3) Base 2 : split 70/30 (mêmes IDs que l'approche couplage)

# 1) construction des bases, format long, non standardisé
Train_Base_2 <- Learn_Base %>% filter(n_curve %in% train_ids)
Test_Base_2  <- Learn_Base %>% filter(n_curve %in% test_ids)

# 2) standardisation (Train_Base_2 -> référence)
Train2_mat <- as.matrix(Train_Base_2 %>% select(n_curve, n, U0, TS))
Test2_mat  <- as.matrix(Test_Base_2  %>% select(n_curve, n, U0, TS))

Learn_Base_2_S <- S(Train2_mat, Test2_mat)
Train_Base_2_S <- as.data.frame(Learn_Base_2_S[[1]])
Test_Base_2_S  <- as.data.frame(Learn_Base_2_S[[2]])

colnames(Train_Base_2_S) <- c("n_curve","n","U0","TS")
colnames(Test_Base_2_S)  <- c("n_curve","n","U0","TS")

# 3) entrainement d'un NN sur Base 2 
Train2_nn <- Train_Base_2_S %>% mutate(n = as.numeric(n), U0 = as.numeric(U0), TS = as.numeric(TS))
Test2_nn  <- Test_Base_2_S  %>% mutate(n = as.numeric(n), U0 = as.numeric(U0), TS = as.numeric(TS))

set.seed(123)
nn_time2 <- system.time({
  nn_fit2 <- neuralnet(
    formula = TS ~ n + U0,
    data = Train2_nn,
    hidden = c(3, 6, 3),
    linear.output = TRUE,
    lifesign = "minimal",
    stepmax = 1e6
  )
})

# 4) prédictions (échelle standardisée)
nn_pred_train2 <- neuralnet::compute(nn_fit2, Train2_nn %>% select(n, U0))$net.result %>% as.vector()
nn_pred_test2  <- neuralnet::compute(nn_fit2, Test2_nn  %>% select(n, U0))$net.result %>% as.vector()

Train2_pred_S <- Train2_nn %>% mutate(TS_pred = nn_pred_train2)
Test2_pred_S  <- Test2_nn  %>% mutate(TS_pred = nn_pred_test2)

# indicateurs (standardisés, comme demandé : Base_2_S)
m_train2_nn_S <- nn_curve_metrics(Train2_pred_S)
m_test2_nn_S  <- nn_curve_metrics(Test2_pred_S)

m_train2_nn_S

# A tibble: 1 × 3
  RMSE_mean MAE_mean RelDiff_mean
      <dbl>    <dbl>        <dbl>
1    0.0270   0.0247        0.671

m_test2_nn_S

# A tibble: 1 × 3
  RMSE_mean MAE_mean RelDiff_mean
      <dbl>    <dbl>        <dbl>
1    0.0241   0.0218        0.511

# 5)pour comparaison avec l’approche couplage : retour à l’échelle originale (déstandardisation de TS)
TS_mean2 <- mean(Train2_mat[,4])
TS_rng2  <- diff(range(Train2_mat[,4]))

Train2_pred_orig <- Train_Base_2 %>%
  mutate(TS_pred = nn_pred_train2 * TS_rng2 + TS_mean2)

Test2_pred_orig <- Test_Base_2 %>%
  mutate(TS_pred = nn_pred_test2 * TS_rng2 + TS_mean2)

m_train2_nn_orig <- nn_curve_metrics(Train2_pred_orig)
m_test2_nn_orig  <- nn_curve_metrics(Test2_pred_orig)

m_train2_nn_orig

# A tibble: 1 × 3
  RMSE_mean MAE_mean RelDiff_mean
      <dbl>    <dbl>        <dbl>
1     1018.     930.         33.8

m_test2_nn_orig

# A tibble: 1 × 3
  RMSE_mean MAE_mean RelDiff_mean
      <dbl>    <dbl>        <dbl>
1      907.     821.         40.1

Interprétation (Base 2) : - Ici, l’apprentissage se fait sur 70% des courbes et le split est aléatoire : train et test ont donc des distributions plus proches qu’avec Base 1 (split « first/last »). - On s’attend en général à une meilleure généralisation (écart train→test plus faible) qu’avec Base 1.

Les indicateurs calculés sur la Base 2 montrent des performances globalement comparables à celles obtenues sur la Base 1, avec toutefois une dégradation plus marquée sur le jeu de test.
Cette différence s’explique par le caractère non déterministe de l’apprentissage du réseau de neurones et par la dépendance de ses performances au choix du jeu d’apprentissage.
En comparaison, l’approche de couplage modèle–données présente des résultats plus stables, car elle repose sur l’estimation d’un nombre réduit de paramètres interprétables et sur une structure de modèle imposée.

Comparaison avec l’approche couplage (Ex 3.5) :
Les indicateurs m_train2_nn_orig / m_test2_nn_orig sont calculés en échelle originale, donc directement comparables à metrics_train / metrics_test.
Les différences de performance se justifient notamment par : - le couplage ajuste un seul triplet \((a,b,c)\) commun à toutes les courbes (modèle parcimonieux mais biais de modèle si les courbes sont hétérogènes), - le NN apprend une relation empirique \((n,U_0)\mapsto U_n\) (plus flexible, mais il « moyenne » implicitement l’effet de la variabilité des paramètres \((a,b,c)\) entre courbes).

4.4 Exercice 4.4 — Conclusion et temps de calcul

# 4.4) Temps de calcul pour comparaison
time_direct <- system.time({
  tmp <- nloptr(
    x0 = c(1, 1, 0.5),
    eval_f = function(theta) f_obj(theta, U0_train, Main_Training_Base, N),
    lb = lower,
    ub = upper,
    opts = list(algorithm = "NLOPT_GN_DIRECT_L", maxeval = 500, print_level = 0)
  )
})

time_gd <- system.time({
  tmp <- tryCatch(
    optim(
      par = c(1, 1, 0.5),
      fn = function(theta) f_obj(theta, U0_train, Main_Training_Base, N),
      method = "L-BFGS-B",
      lower = lower,
      upper = upper,
      control = list(maxit = 300)
    ),
    error = function(e) NULL
  )
})

tibble(
  DIRECT_sec = time_direct["elapsed"],
  GD_sec = time_gd["elapsed"],
  NN_sec = nn_time["elapsed"]
)

# A tibble: 1 × 3
  DIRECT_sec GD_sec NN_sec
       <dbl>  <dbl>  <dbl>
1       19.5   21.4   2.03

Conclusion attendue :

• Couplage modèle–données : très interprétable, peu de paramètres, calibration parfois coûteuse si objectif non convexe ; performances liées à la pertinence du modèle.
• ML (NN) : flexible, plus de paramètres donc risque d’overfit, besoin de standardisation/validation ; peut mieux approximer si la loi vraie est complexe.

4.4.1 Interprétation des temps

Temps mesurés sur la machine utilisée :

• DIRECT ≈ 19.5 s
• Gradient (L-BFGS-B) ≈ 21.38 s
• Réseau de neurones (apprentissage) ≈ 2.03 s

Commentaire : - DIRECT et le gradient coûtent du temps car chaque évaluation de l’objectif nécessite de simuler beaucoup de courbes (ici 70) sur \(N=40\). - Le réseau de neurones est beaucoup plus rapide à entraîner ici (base petite, architecture 3-6-3), et une fois entraîné la prédiction est quasi instantanée. - En revanche, le couplage modèle–données fournit des paramètres interprétables (3 paramètres) alors que le réseau utilise 58 paramètres (moins interprétable).