TP: Physics Informed Neural Networks (PINNs) - suite

Dans ce TP nous commencerons par illustrer l'assimilation de données et l'estimation de paramètre sur notre EDO $$ y'(t) + a \, y(t) = f(t), $$ avec $y(0) = 0$, et $$ f(t) = \begin{cases} e^{-bt}\sin(t) & \text{si } t \geq 0, \\ 0 & \text{si } t < 0, \end{cases} $$ qui admet pour solution explicite: $$ y(t) = \begin{cases} \frac{1}{(a-b)^2 + 1} \big( ((a-b) \sin(t) - \cos(t)) e^{-bt} + e^{-at} \big) & \text{si } t \geq 0, \\ 0 & \text{si } t < 0. \end{cases} $$

Commençons par illustrer numériquement la solution de cette équation.

import numpy as np
import matplotlib.pyplot as plt

import torch
import torch.nn as nn
import torch.optim as optim
from torch.autograd import grad

# paramètres de l'EDO
a = 0.1
b = 0.1

# grille de points auxquelles on veut évaluer la solution y
t = np.linspace(-5., 30., 200)

# définition de la solution explicite
sol_expl = lambda t: ( ( (a-b) * np.sin(t) - np.cos(t) ) * np.exp(-b*t) + np.exp(-a*t) ) / ((a-b)**2 + 1) *(t>0)

plt.plot(t, sol_expl(t), label='solution y(t)')
plt.plot(t, np.sin(t)*np.exp(-b*t)*(t>0) , label='source f(t)')

plt.legend(fontsize=14)
plt.xlabel('t', fontsize=18)
plt.title('explicite solution and source f', fontsize=18)
plt.show()

I) Assimilation de données avec un PINN

Nous allons maintenant entrainer un PINNs sur notre EDO, mais rajoutons à la loss des observations de la solution de cette EDO.

# On génère des observations de la solution explicite (données artificielles)
# qu'on va ajouter à la loss (L_{data}) pour l'entrainement

# on définit la grille de temps
n_obs = 10 # nombre d'observation

t_obs = t[-1]*np.random.rand(n_obs)

# on définit les observations
y_obs = sol_expl(t_obs) + 0.05*np.random.randn(n_obs) # on rajoute du bruit

# Illustration des données artificielle par rapport à la solution explicite
plt.plot(t, sol_expl(t), label="ground truth")
plt.plot(t_obs, y_obs, 'o', label="data")

plt.xlabel("t")
plt.legend(fontsize=14)
plt.xlabel('t', fontsize=18)
plt.ylabel('y(t)', fontsize=18)
plt.title('observations of the ode', fontsize=18)
plt.show()

Question:

• Entrainer un réseau de neurones uniquement à l'aide des données observées;
• Entrainer un PINN incluant les données observées;
• Comparer, à l'aide d'une illustration graphique, les capacités de généralisation des deux réseaux sur l'intervalle de temps $[-5, 30]$. On pourra ajouter à ce graphique la solution explicite.

On commencera par définir:

• Une classe de réseau de neurones;
• Une fonction d'entrainement qu'à partir des données;
• Une fonction d'entrainement d'un PINN faisant intervenir les données.

On reprend la classe de neurones et la fonction d'entrainement simple des TP précédents.

class MLP(nn.Module):
    def __init__(self, layers, activations):
        """
        Parameters:
        - layers: List[int], the number of neurons in each layer (including input and output layers).
        - activations: List[nn.Module], the activation functions for each layer.
          The length of activations should be len(layers) - 2 (one for each hidden layer).
        """
        super(MLP, self).__init__()
        
        if len(activations) != len(layers) - 2:
            raise ValueError("The number of activation functions must match the number of hidden layers.")

        # Build the network
        self.layers = nn.Sequential()
        for i in range(len(layers) - 1):
            self.layers.add_module(f"linear_{i}", nn.Linear(layers[i], layers[i + 1]))
            if i < len(activations):  # Apply activation to all but the last layer
                self.layers.add_module(f"activation_{i}", activations[i])

    def forward(self, x):
        return self.layers(x)

# fonction d'entrainement
def train(model, optimizer, loss_fct, t_data, y_data, epochs=1000):
    
    """
    Trains a given model on specific data for a defined number of epochs.

    Parameters:
    - model : nn.Module
        The neural network model to be trained.
    - optimizer : torch.optim.Optimizer
        The optimizer used to update the model's weights.
    - loss_fct : 
        The loss function to measure the error between the model's predictions and the target data.
    - t_data : torch.Tensor
        Input data (time variables).
    - y_data : torch.Tensor
        Target data corresponding to the input data.
    - epochs : int, optional (default 1000)
        The number of epochs (iterations over the entire dataset) for training.

    """
    
    for epoch in range(epochs):

        # reset gradients for the optimizer
        optimizer.zero_grad()

        # compare the model to the data
        y_pred = model(t_data)
        loss = loss_fct(y_pred, y_data)

        # optimisation step
        loss.backward()
        optimizer.step()
        
        if epoch % 1000 == 0:
            print(f"Epoch {epoch}, Loss: {loss.item()}")

    print(f"Epoch {epoch}, Loss: {loss.item()}")

On définit maintenant la fonction d'entrainement qui utilise à la fois le modèle mathématique et les données synthétiques. Pour cela on reprend la fonction d'entrainement d'un PINN dans laquelle on rajoute le terme de loss correspondant aux données.

def train_pinn_data(model, optimizer, loss_fct, t_data, t_collo, y0, y_data, a, f, epochs=1000):

    """
    Trains a Physics-Informed Neural Network (PINN) using both data and physics-based constraints.

    Parameters:
        - model (torch.nn.Module): The neural network model that approximates the solution to the ODE.
        - optimizer (torch.optim.Optimizer): The optimizer used to update the model's weights during training.
        - loss_fct: The loss function to compute the error between the predicted and true data values.
        - t_data (torch.Tensor): Input time values corresponding to observed data points for training.
        - t_collo (torch.Tensor): Input time values used to enforce the physics-based ODE constraints.
        - y0 (torch.Tensor): Initial condition of the ODE (not explicitly used in this implementation).
        - y_data (torch.Tensor): Observed data values corresponding to the input `t1` for supervised learning.
        - a (float):  Coefficient of the linear term in the ODE y' + a y - f = 0.
        - f (torch.Tensor): The forcing term in the ODE.
        - epochs (int, optional): Number of training iterations (default: 1000).
    """
    
    for epoch in range(epochs):

        # reset gradients for the optimizer
        optimizer.zero_grad()

        # compare the model with respect to the data
        y_pred = model(t_data)
        loss_data = loss_fct(y_data, y_pred) 

        ### compute the loss involving the ODE terms

        # evaluate the model on the ODE grid
        y = model(t_collo)

        # evaluate the derivatives involved in the ode
        dydt = grad(y, t_collo, grad_outputs=torch.ones_like(y), create_graph=True)[0]
        residual = dydt + a * y - f
        
        # evaluate the residual, that is the mismatch with the ode
        loss_ode = residual.square().sum()

        #evaluate the mismatch with the initial condition
        loss_ic = (model(torch.tensor([0.]))-y0)**2

        # compute the loss
        loss =  loss_ode + loss_ic + loss_data

        # optimization step
        loss.backward()
        optimizer.step()
        if epoch % 1000 == 0:
            print(f"Epoch {epoch}, Loss: {loss.item()}")
            
    print(f"Epoch {epoch}, Loss: {loss.item()}")

On va maintenant entrainer les réseaux. Pour cela, mettons les données dans le bon format.

# Conversion des données "artificielles" en torch.tensor
t_data = torch.tensor(t_obs, dtype=torch.float32).reshape(-1,1) 
y_data = torch.tensor(y_obs, dtype=torch.float32).reshape(-1,1) 

# On définit les points de collocation
t_collo = torch.linspace(t[0], t[-1], 100, dtype=torch.float32).reshape(-1,1).requires_grad_()
# on va calculer y' par rapport à ces valeurs

# On définit le terme source
with torch.no_grad():
    f = np.sin(t_collo) * np.exp(-b*t_collo) *(t_collo>0)

Définissons les deux réseaux qu'on va entrainer.

# architecture du réseau
layers = [1, 50, 50, 50, 1]  # Input, hidden layers, output
activations = [nn.Tanh(), nn.Tanh(), nn.Tanh()]  # Different activations for each hidden layer

# condition initiale
y0 = 0.

# on définit 2 réseaux : un PINN et un entrainé juste sur les données
model_std = MLP(layers, activations)
model_pinn = MLP(layers, activations)

# loss quadratique pour L_{data}
loss_fct = torch.nn.MSELoss(reduction='sum')

# méthode d'optimisation pour le réseau entrainé que sur des données
optimizer_std = optim.Adam(model_std.parameters(), lr=0.005)

# entrainement du réseau qu'à partir des données
train(model_std, optimizer_std, loss_fct, t_data, y_data, 40000)

Epoch 0, Loss: 1.5798189640045166
Epoch 1000, Loss: 0.05634869635105133
Epoch 2000, Loss: 0.05616194009780884
Epoch 3000, Loss: 0.056059204041957855
Epoch 4000, Loss: 0.05612415075302124
Epoch 5000, Loss: 0.05582982301712036
Epoch 6000, Loss: 0.054485760629177094
Epoch 7000, Loss: 0.05363202840089798
Epoch 8000, Loss: 0.008142339065670967
Epoch 9000, Loss: 0.006374146323651075
Epoch 10000, Loss: 0.0026368407998234034
Epoch 11000, Loss: 0.007776322774589062
Epoch 12000, Loss: 0.014769842848181725
Epoch 13000, Loss: 0.009971349500119686
Epoch 14000, Loss: 0.0026500301901251078
Epoch 15000, Loss: 0.004339136183261871
Epoch 16000, Loss: 0.0026113337371498346
Epoch 17000, Loss: 0.008989427238702774
Epoch 18000, Loss: 0.00212152604945004
Epoch 19000, Loss: 0.003635855857282877
Epoch 20000, Loss: 0.006638651248067617
Epoch 21000, Loss: 0.0023539222311228514
Epoch 22000, Loss: 0.002346342895179987
Epoch 23000, Loss: 0.014741810038685799
Epoch 24000, Loss: 0.052617717534303665
Epoch 25000, Loss: 0.011784191243350506
Epoch 26000, Loss: 0.005956662353128195
Epoch 27000, Loss: 0.00416391808539629
Epoch 28000, Loss: 0.011775477789342403
Epoch 29000, Loss: 0.002432651584967971
Epoch 30000, Loss: 0.006452460307627916
Epoch 31000, Loss: 0.01971275731921196
Epoch 32000, Loss: 0.0022261415142565966
Epoch 33000, Loss: 0.0060270847752690315
Epoch 34000, Loss: 0.0078819515183568
Epoch 35000, Loss: 0.01083344779908657
Epoch 36000, Loss: 0.007879367098212242
Epoch 37000, Loss: 0.001906237448565662
Epoch 38000, Loss: 0.007879149168729782
Epoch 39000, Loss: 0.007876201532781124
Epoch 39999, Loss: 0.00802676659077406

# méthode d'optimisation pour le PINN
optimizer_pinn = optim.Adam(model_pinn.parameters(), lr=0.005)

# entrainement du pinn avec les données
train_pinn_data(model_pinn, optimizer_pinn, loss_fct, t_data, t_collo, y0, y_data, a, f, 20000)

Epoch 0, Loss: 8.657977104187012
Epoch 1000, Loss: 0.18530361354351044
Epoch 2000, Loss: 0.031076328828930855
Epoch 3000, Loss: 0.01294733863323927
Epoch 4000, Loss: 0.014712825417518616
Epoch 5000, Loss: 0.014337877742946148
Epoch 6000, Loss: 0.0142395393922925
Epoch 7000, Loss: 0.01214692648500204
Epoch 8000, Loss: 0.014813145622611046
Epoch 9000, Loss: 0.016587508842349052
Epoch 10000, Loss: 0.017921241000294685
Epoch 11000, Loss: 0.02227081172168255
Epoch 12000, Loss: 0.016530705615878105
Epoch 13000, Loss: 0.014213841408491135
Epoch 14000, Loss: 0.01361708901822567
Epoch 15000, Loss: 0.01533418707549572
Epoch 16000, Loss: 0.016133449971675873
Epoch 17000, Loss: 0.012564685195684433
Epoch 18000, Loss: 0.012125207111239433
Epoch 19000, Loss: 0.013291883282363415
Epoch 19999, Loss: 0.014044342562556267

Procédons à la comparaison de la solution explicite et des deux réseaux entrainés sur la grille de temps fine t.

# évaluation du model sur une grille de temps fine
t_val = torch.tensor(t, dtype=torch.float32).reshape(-1,1)

# comparaison de la solution explicite avec les prédictions des réseaux de neurones
plt.plot(t, sol_expl(t), label="ground truth")
plt.plot(t_obs, y_obs, 'o', label="data")
plt.plot(t_val, model_pinn(t_val).detach(), label="pinn predictions")
plt.plot(t_val, model_std(t_val).detach(), label="non pinn predictions") # on utilise .detach()"

plt.legend(fontsize=14)
plt.xlabel('t', fontsize=18)
plt.ylabel('y(t)', fontsize=18)
plt.title('fit of the ode data with a pinn', fontsize=18)
plt.show()

Le PINN est capable de représenter les données tout en respectant le modèle mathématique. Le réseau juste entrainé sur les données représente bien les données, mais fera des prédictions complétement erronées si on l'évalue à des temps hors du jeu de données.

II) Assimilation de données et problème inverse à l'aide d'un PINN

Supposons que :

• nous ayons des observations bruitées : $$ \mathcal{D} = \{(t^{data}_1, y^{data}_1), \dots, (t^{data}_n, y^{data}_n)\}, $$ telles que $$ y^{data}_j = y(t^{data}_j) + \epsilon_j, \qquad j = 1, \dots, n, $$
• nous ne connaissons pas $a$.

Dans cette section, pour illustrer la capacité des PINNs à estimer des parametres dans une EDO, on se donne un $a$, et on va générer de fausse données. Ensuite, on va faire comme si on ne connaissait pas $a$, et on va essayer de retrouver sa valeur.

Le problème inverse consiste ici à estimer $a$ à partir des données et de l'entrainement d'un réseau de neurones (problème d'optimisation). En d'autres termes, ce paramètre sera aussi appris.

Pour ce faire, on va redéfinir la classe python de réseau de neurones pour y incorporer ce paramètre.

Commençons par générer de fausses données, pour lesquelles $a=0.1$ et essayons de retrouver cette valeurs en entrainant un PINNs avec une loss dépendant de l'EDO et des (fausses) données.

# on génère les fausses données sur l'intervalle à partitr du temps t=0
n_obs = 30  # nombre d'observation

t_obs = t[-1]*np.random.rand(n_obs)
t_obs = np.sort(t_obs)

y_obs = sol_expl(t_obs) + 0.05*np.random.randn(n_obs) # on rajoute du bruit

# illustration des données par rapport à la solution explicite 
# pour a=0.1

plt.plot(t, sol_expl(t), label="ground truth")
plt.plot(t_obs, y_obs, 'o', label="data")
plt.xlabel("t")
plt.legend(fontsize=14)
plt.xlabel('t', fontsize=18)
plt.ylabel('y(t)', fontsize=18)
plt.title('observations of the ode', fontsize=18)
plt.show()

Question: A l'aide d'un PINN estimer le paramètre $a$

• On commencera par définir une classe de réseaux de neurones appropriée faisant intervenir le paramètre à apprendre.
• On définira une fonction d'entrainement du PINN adaptée au problème, et faisant intervenir les données.
• On comparera sur un graphique:
  • les données;
  • les prédictions du réseau de neurones après entrainement sur la grille t;
  • la solution explicite avec le paramètre $a$ estimé.

Définissons la classe en ajoutant le paramètre à estimer.

class PINN_IP(nn.Module):
    def __init__(self, layers, activations):
        """
        Parameters:
        - layers: List[int], the number of neurons in each layer (including input and output layers).
        - activations: List[nn.Module], the activation functions for each layer.
          The length of activations should be len(layers) - 2 (one for each hidden layer).
        """
        super().__init__()
        
        if len(activations) != len(layers) - 2:
            raise ValueError("The number of activation functions must match the number of hidden layers.")

        # Build the network
        self.layers = nn.Sequential()
        for i in range(len(layers) - 1):
            self.layers.add_module(f"linear_{i}", nn.Linear(layers[i], layers[i + 1]))
            if i < len(activations):  # Apply activation to all but the last layer
                self.layers.add_module(f"activation_{i}", activations[i])

        self.a = nn.Parameter(torch.tensor(1., requires_grad=True))  # add the ODE parameter to estimate (to learn) 

    def forward(self, x):
        return self.layers(x)

La fonction d'entrainement est similaire à celle des PINNs. La seule modification est sur comment on fait intervenir le paramètre $a$ de l'EDO.

def train_IP(model, optimizer, loss_fct, t_data, t_collo, y0, y_data, f, epochs=1000):
    
    """
    Trains a Physics-Informed Neural Network (PINN) using both data and physics-based constraints.

    Parameters:
        - model (torch.nn.Module): The neural network model that approximates the solution to the ODE.
        - optimizer (torch.optim.Optimizer): The optimizer used to update the model's weights during training.
        - loss_fct: The loss function to compute the error between the predicted and true data values.
        - t_data (torch.Tensor): Input time values corresponding to observed data points for training.
        - t_collo (torch.Tensor): Input time values used to enforce the physics-based ODE constraints.
        - y0 (torch.Tensor): Initial condition of the ODE (not explicitly used in this implementation).
        - y_data (torch.Tensor): Observed data values corresponding to the input `t1` for supervised learning.
        - f (torch.Tensor): The forcing term in the ODE.
        - epochs (int, optional): Number of training iterations (default: 1000).
    """
    
    for epoch in range(epochs):

        # reset gradients for the optimizer
        optimizer.zero_grad()

        # compare the model with respect to the data
        y_pred = model(t_data)
        loss_data = loss_fct(y_data, y_pred) 

        ### compute the loss involving the ODE terms

        # evaluate the model on the ODE grid
        y = model(t_collo)

        # evaluate the derivatives involved in the ode
        dydt = grad(y, t_collo, grad_outputs=torch.ones_like(y), create_graph=True)[0]
        residual = dydt + model.a * y - f          # the parameter a plays a role here
        
        # evaluate the residual, that is the mismatch with the ode
        loss_ode = residual.square().sum()

        #evaluate the mismatch with the initial condition
        loss_ic = (model(torch.tensor([0.]))-y0)**2

        # compute the loss
        loss =  loss_ode + loss_ic + loss_data

        # optimization step
        loss.backward()
        optimizer.step()
        if epoch % 1000 == 0:
            print(f"Epoch {epoch}, Loss: {loss.item()}")
            
    print(f"Epoch {epoch}, Loss: {loss.item()}")

Mettons en forme les données.

# conversion des données en tenseur torch
t_data = torch.tensor(t_obs, dtype=torch.float32).reshape(-1,1) 
y_data = torch.tensor(y_obs, dtype=torch.float32).reshape(-1,1)

# on définit les points de collocation
t_collo = torch.linspace(t[0], t[-1], 100, dtype=torch.float32).reshape(-1,1).requires_grad_()
# on va calculer des dérivées par rapport à ces valeurs

# on définit le terme sour de l'EDO
with torch.no_grad():
    f = np.sin(t_collo) * np.exp(-b*t_collo) *(t_collo>0)

# loss quadratique pour L_{data}
loss_fct = torch.nn.MSELoss(reduction='sum')

On définit le réseau de neurones.

# on définit l'architecture du réseau
layers = [1, 40, 40, 40, 1]  # Input, hidden layers, output
activations = [nn.Tanh(), nn.Tanh(), nn.Tanh()]  # Different activations for each hidden layer

# condition initiale de l'EDO
y0 = 0.

model_pinn = PINN_IP(layers, activations)
print(model_pinn)

PINN_IP(
  (layers): Sequential(
    (linear_0): Linear(in_features=1, out_features=40, bias=True)
    (activation_0): Tanh()
    (linear_1): Linear(in_features=40, out_features=40, bias=True)
    (activation_1): Tanh()
    (linear_2): Linear(in_features=40, out_features=40, bias=True)
    (activation_2): Tanh()
    (linear_3): Linear(in_features=40, out_features=1, bias=True)
  )
)

# méthode d'optimisation pour le PINN
optimizer_pinn = optim.Adam(model_pinn.parameters(), lr=0.005)

train_IP(model_pinn, optimizer_pinn, loss_fct, t_data, t_collo, y0, y_data, f, epochs=20000)

Epoch 0, Loss: 17.530567169189453
Epoch 1000, Loss: 1.290029525756836
Epoch 2000, Loss: 0.3914357125759125
Epoch 3000, Loss: 0.31573277711868286
Epoch 4000, Loss: 0.2900945544242859
Epoch 5000, Loss: 0.43634235858917236
Epoch 6000, Loss: 0.1264444887638092
Epoch 7000, Loss: 0.10852628946304321
Epoch 8000, Loss: 0.08978502452373505
Epoch 9000, Loss: 0.07505929470062256
Epoch 10000, Loss: 0.09754650294780731
Epoch 11000, Loss: 0.09586494415998459
Epoch 12000, Loss: 0.07240636646747589
Epoch 13000, Loss: 0.06465081870555878
Epoch 14000, Loss: 0.06450511515140533
Epoch 15000, Loss: 0.06816313415765762
Epoch 16000, Loss: 0.06450008600950241
Epoch 17000, Loss: 0.05646022781729698
Epoch 18000, Loss: 0.054925281554460526
Epoch 19000, Loss: 0.05465039610862732
Epoch 19999, Loss: 0.05068134516477585

On obtient pour estimation de $a$

model_pinn.a.item()

0.09289278835058212

Cette estimation est très bonne sachant que la vraie valeur est $a=0.1$.

sol_estim = lambda a, t: ( ( (a-b) * np.sin(t) - np.cos(t) ) * np.exp(-b*t) + np.exp(-a*t) ) / ((a-b)**2 + 1) *(t>0)

# on définit une grille fine pour l'évaluation du réseau
t_val = torch.tensor(t, dtype=torch.float32).reshape(-1,1)

# comparaison de la solution explicite avec les prédictions des réseaux de neurones
# après avoir estimé $a$
plt.plot(t, sol_estim(model_pinn.a.item(), t), label="ground truth")
plt.plot(t_obs, y_obs, 'o', label="data")
plt.plot(t_val, model_pinn(t_val).detach(), label="predictions")

plt.xlabel("t")
plt.legend(fontsize=14)
plt.xlabel('t', fontsize=18)
plt.ylabel('y(t)', fontsize=18)
plt.title(f'fit of the ode, estimation a = {model_pinn.a.item():.3f}', fontsize=18)
plt.show()

III) Modèle de Lotka Volterra

Le modèle de Lotka-Volterra est un système d'EDO décrivant une dynamique de populations \textit{proies-prédateurs} : \begin{align*} x'(t) &= a \, x(t) - b \, x(t)y(t),\\ y'(t) &= c \, x(t)y(t) - d \, y(t), \end{align*} où :

• $x(t)$ représente la population des proies au temps $t$,
• $y(t)$ représente la population des prédateurs au temps $t$,
• $a > 0$ est le taux de reproduction des proies,
• $b > 0$ est le taux de mortalité due à la prédation,
• $c > 0$ est le taux de reproduction des prédateurs en fonction des proies rencontrées,
• $d > 0$ est le taux de mortalité des prédateurs.

Commençons par voir comment résoudre cette EDO avec odeint

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt
import pandas as pd

import torch
import torch.nn as nn
import torch.optim as optim
from torch.autograd import grad

def lotka_volterra(u, t, a, b, c, d):
    x, y = u
    dxdt = a*x - b*x*y
    dydt = -c*y + d*x*y
    return [dxdt, dydt]

# Parameters
a = 3.    
b = 1.    
c = 1.    
d = 1.   

# Initial conditions
x0 = 1.   # initial prey population
y0 = 1.    # initial predator population
initial_state = [x0, y0]

# Time vector
t_start = 0.
t_end = 10.
n_obs = 1000
t = np.linspace(t_start, t_end, n_obs)

# Solve ODE system
sol = odeint(lotka_volterra, [x0, y0], t, args=(a, b, c, d))

plt.figure(figsize=(15, 6))
fig, axes = plt.subplots(1, 2, figsize=(15, 6))  

axes[0].plot(t, sol[:,0], label='Prey (x)')
axes[0].plot(t, sol[:,1], label='Predator (y)')
axes[0].set_title('Lotka-Volterra model simulation, population evolution', fontsize=18)
axes[0].set_xlabel('t', fontsize=18)
axes[0].set_ylabel('Population density', fontsize=18)
axes[0].legend(fontsize=18)
axes[0].grid(True)

axes[1].plot(sol[:,0], sol[:,1], label='')
axes[1].set_title('Lotka-Volterra model simulation, population evolution', fontsize=18)
axes[1].set_xlabel('Prey (x)', fontsize=18)
axes[1].set_ylabel('Predator (y)', fontsize=18)
axes[1].grid(True)
plt.tight_layout()
plt.show()

<Figure size 1080x432 with 0 Axes>

Dans le fichier prey_predator_data.csv, on a des données sur l'évolution des populations de proie et prédateur. Commençons par visualiser ces données.

df = pd.read_csv("prey_predator_data.csv")

t_data = torch.tensor(df['Time'].values, dtype=torch.float32).reshape(-1, 1)#.requires_grad_()
u_data = torch.tensor(df[['Prey', 'Predator']].values, dtype=torch.float32)

plt.figure(figsize=(15, 6))
fig, axes = plt.subplots(1, 2, figsize=(15, 6)) 

with torch.no_grad():
    axes[0].plot(t_data, u_data[:,0], "o", label='Prey (x) - data')
    axes[0].plot(t_data, u_data[:,1], "o", label='Predator (y) - data')

axes[0].set_title('data prey predator', fontsize=18)
axes[0].set_xlabel('t', fontsize=18)
axes[0].set_ylabel('data population density', fontsize=18)
axes[0].legend(fontsize=12)
axes[0].grid(True)

with torch.no_grad():
    axes[1].plot(u_data[:,0], u_data[:,1], "o", label='data')

axes[1].set_title('data vs Lotka-Volterra model', fontsize=18)
axes[1].set_xlabel('Prey (x)', fontsize=18)
axes[1].set_ylabel('Predator (y)', fontsize=18)
axes[1].legend(fontsize=18)
axes[1].grid(True)

plt.tight_layout()

plt.show()

<Figure size 1080x432 with 0 Axes>

Question: A l'aide d'un PINN estimer les paramètres biologiques $a,b,c,d$

• On commencera par définir une classe de réseaux de neurones appropriée faisant intervenir les 4 paramètres à apprendre.
• On définira une fonction d'entrainement adaptée au problème.
• On comparera sur un graphique:
  • les données;
  • les prédictions du réseau de neurones après entrainement;
  • la simulation du modèle de Lotka-Volterra avec odeint et les paramètres estimés.

On reprend notre classe de réseau de neurones à laquelle on ajoute les paramètres à estimer.

class PINN_LV_IP(nn.Module):
    def __init__(self, layers, activations):
        """
        Parameters:
        - layers: List[int], the number of neurons in each layer (including input and output layers).
        - activations: List[nn.Module], the activation functions for each layer.
          The length of activations should be len(layers) - 2 (one for each hidden layer).
        """
        super().__init__()
        
        if len(activations) != len(layers) - 2:
            raise ValueError("The number of activation functions must match the number of hidden layers.")

        # Build the network
        self.layers = nn.Sequential()
        for i in range(len(layers) - 1):
            self.layers.add_module(f"linear_{i}", nn.Linear(layers[i], layers[i + 1]))
            if i < len(activations):  # Apply activation to all but the last layer
                self.layers.add_module(f"activation_{i}", activations[i])

        # ODE parameters to estimate (to learn) 
        self.a = nn.Parameter(torch.tensor(1., requires_grad=True))
        self.b = nn.Parameter(torch.tensor(1., requires_grad=True))
        self.c = nn.Parameter(torch.tensor(1., requires_grad=True))
        self.d = nn.Parameter(torch.tensor(1., requires_grad=True))

    def forward(self, x):
        return self.layers(x)

On définit la fonction d'entrainement. Attention, on a ici une EDO en dimension 2, il faut faire attention à la taille des tenseurs torch.

def train_LV_IP(model, optimizer, loss_fct, t_data, t_collo, u0, u_data, epochs=1000):
    
    for epoch in range(epochs):

        # reset gradients for the optimizer
        optimizer.zero_grad()

        # compare the model with respect to the data
        u_pred_data = model(t_data)
        loss_data = loss_fct(u_pred_data, u_data)

        ### compute the loss involving the ODE terms

        # evaluate the model on the ODE grid
        u_pred = model(t_collo)
        x, y = u_pred[:,0:1], u_pred[:,1:2]

        # evaluate the derivatives involved in the ode
        dxdt = grad(x, t_collo, grad_outputs=torch.ones_like(x), create_graph=True)[0]
        dydt = grad(y, t_collo, grad_outputs=torch.ones_like(y), create_graph=True)[0]

        # evaluate the residual, that is the mismatch with the ode
        dx = model.a*x - model.b*x*y
        dy = -model.c*y + model.d*x*y
    
        loss_ode = (dxdt - dx).square().sum() + (dydt - dy).square().sum()

        #evaluate the mismatch with the initial condition
        u_pred_ic = model(torch.tensor([0.], dtype=torch.float32))
        loss_ic = (u_pred_ic - u0).square().sum()

        # définition de la loss
        loss = loss_ode + loss_data + loss_ic

        # optimization step
        loss.backward()
        optimizer.step()
    
        if epoch%1000 == 0:
            print(f"Epoch {epoch}, Loss: {loss.item()}")
            
    print(f"Epoch {epoch}, Loss: {loss.item()}")

# on définit les points de colocations
t_collo = torch.linspace(t_start, t_end, 100).reshape(-1,1).requires_grad_()
# on va calculer des dérivées par rapport à ces valeurs

# loss quadratique pour L_{data}
loss_fct = torch.nn.MSELoss(reduction='sum')

# condition initiale de l'EDO
u0 = torch.tensor([x0, y0], dtype=torch.float32)

# on définit l'architecture du réseau
layers = [1, 32, 64, 32, 2]  # Input, hidden layers, output
activations = [nn.Tanh(), nn.Tanh(), nn.Tanh()]  # Different activations for each hidden layer

model_lv_ip = PINN_LV_IP(layers, activations)
print(model_lv_ip)

PINN_LV_IP(
  (layers): Sequential(
    (linear_0): Linear(in_features=1, out_features=32, bias=True)
    (activation_0): Tanh()
    (linear_1): Linear(in_features=32, out_features=64, bias=True)
    (activation_1): Tanh()
    (linear_2): Linear(in_features=64, out_features=32, bias=True)
    (activation_2): Tanh()
    (linear_3): Linear(in_features=32, out_features=2, bias=True)
  )
)

learning_rate = 0.001 

optimizer_lv_ip = torch.optim.Adam(model_lv_ip.parameters(), lr=learning_rate)

train_LV_IP(model_lv_ip, optimizer_lv_ip, loss_fct, t_data, t_collo, u0, u_data, epochs=30000)

Epoch 0, Loss: 283.7341003417969
Epoch 1000, Loss: 74.58319091796875
Epoch 2000, Loss: 44.38751220703125
Epoch 3000, Loss: 12.934341430664062
Epoch 4000, Loss: 5.593912601470947
Epoch 5000, Loss: 4.114262580871582
Epoch 6000, Loss: 3.4386048316955566
Epoch 7000, Loss: 3.176316976547241
Epoch 8000, Loss: 2.61506724357605
Epoch 9000, Loss: 2.3250067234039307
Epoch 10000, Loss: 2.071730136871338
Epoch 11000, Loss: 1.8399827480316162
Epoch 12000, Loss: 1.7000688314437866
Epoch 13000, Loss: 1.5772970914840698
Epoch 14000, Loss: 1.5243178606033325
Epoch 15000, Loss: 1.5147485733032227
Epoch 16000, Loss: 1.6741383075714111
Epoch 17000, Loss: 1.3543474674224854
Epoch 18000, Loss: 1.3197346925735474
Epoch 19000, Loss: 1.4140877723693848
Epoch 20000, Loss: 1.286712646484375
Epoch 21000, Loss: 1.4145796298980713
Epoch 22000, Loss: 1.271097183227539
Epoch 23000, Loss: 1.4051527976989746
Epoch 24000, Loss: 1.4629477262496948
Epoch 25000, Loss: 1.5788968801498413
Epoch 26000, Loss: 1.244924783706665
Epoch 27000, Loss: 1.2394773960113525
Epoch 28000, Loss: 1.244489073753357
Epoch 29000, Loss: 1.4729983806610107
Epoch 29999, Loss: 1.551688551902771

print(f'estimation de a {model_lv_ip.a.item():.2f}')
print(f'estimation de b {model_lv_ip.b.item():.2f}')
print(f'estimation de c {model_lv_ip.c.item():.2f}')
print(f'estimation de d {model_lv_ip.d.item():.2f}')

estimation de a 4.40
estimation de b 1.60
estimation de c 1.39
estimation de d 1.45

# resolution de l'ODE (modèle de Lotka-Volterra) avec les paramètres estimés
params = (model_lv_ip.a.item(), model_lv_ip.b.item(), model_lv_ip.c.item(), model_lv_ip.d.item())
sol = odeint(lotka_volterra, [x0, y0], t, args=params)

t_test = torch.tensor(t, dtype=torch.float32).reshape(-1,1)
u = model_lv_ip(t_test).detach().numpy()

# Extract prey and predator populations
x_pred = u[:, 0]
y_pred = u[:, 1]


plt.figure(figsize=(15, 6))
fig, axes = plt.subplots(1, 2, figsize=(15, 6)) 

axes[0].plot(t, sol[:,0], label='prey (x) ground truth ')
axes[0].plot(t, sol[:,1], label='predator (y) ground truth')
axes[0].plot(t, x_pred, label='prey (x) prediction')
axes[0].plot(t, y_pred, label='predator (y) prediction')

with torch.no_grad():
    axes[0].plot(t_data, u_data[:,0], "o", label='prey (x) - data')
    axes[0].plot(t_data, u_data[:,1], "o", label='Ppredator (y) - data')

axes[0].set_title('data vs Lotka-Volterra model', fontsize=18)
axes[0].set_xlabel('t', fontsize=18)
axes[0].set_ylabel('Population density', fontsize=18)
axes[0].legend(fontsize=12)
axes[0].grid(True)

axes[1].plot(sol[:,0], sol[:,1], label='ground truth')
axes[1].plot(x_pred, y_pred, label='predictions')

with torch.no_grad():
    axes[1].plot(u_data[:,0], u_data[:,1], "o", label='data')

axes[1].set_title('data vs Lotka-Volterra model', fontsize=18)
axes[1].set_xlabel('Prey (x)', fontsize=18)
axes[1].set_ylabel('Predator (y)', fontsize=18)
axes[1].legend(fontsize=18)
axes[1].grid(True)

plt.tight_layout()

plt.show()

<Figure size 1080x432 with 0 Axes>