Modélisation d’un pendule double par IA

7.8. Modélisation d’un pendule double par IA#

Marc Buffat, dpt mécanique, Lyon 1

inspiré de Deep Learning with Real Data and the Chaotic Double Pendulum on https://cloud4scieng.org

Pendule double

from scipy.integrate import odeint 
import time
import math
import numpy as np
import pylab as py
from math import *
import numpy as np
from scipy.optimize import newton
from matplotlib import animation, rc
from matplotlib import pyplot as plt
import pandas as pd
import torch
import torch.nn as nn
from torch.autograd import Variable
from sklearn.preprocessing import MinMaxScaler
from IPython.display import HTML

7.8.1. Problème#

Expérience du double pendule

from IPython.display import YouTubeVideo
YouTubeVideo('DVr7KnpDAOc', width=800, height=500)

7.8.2. Modélisation mécanique#

On doit donc définir:

  1. les parametres du problème: longueur et masse des 2 pendules \(l_1,m_1,l_2,m_2\)

  2. les 2 ddl du problème: \(\theta_1(t)\) et \(\theta_2(t)\)

Système conservatif: équation de Lagrange

(7.1)#\[\begin{align} \alpha \left( \beta \omega_{1}^{2} \sin{(\theta_{1})} + \beta \frac{d^{2}\theta_{1}}{d t^{2}} + \omega_{1}^{2} \sin{(\theta_{1})} + \frac{d^{2}\theta_{1}}{d t^{2}} \right) + \sin{(\theta_{1} - \theta_{2})} \left(\frac{d\theta_{2}}{d t}\right)^{2} + \cos{(\theta_{1} - \theta_{2})} \frac{d^{2}}{d t^{2}} \theta_{2} = 0\\ \alpha \omega_{1}^{2} \sin{\theta_{2}} - \alpha \sin{(\theta_{1} - \theta_{2})} \left(\frac{d \theta_{1}}{d t}\right)^{2} + \alpha \cos{(\theta_{1}- \theta_{2})} \frac{d^{2} \theta_{1}}{d t^{2}} + \frac{d^{2} \theta_{2}}{d t^{2}} = 0 \end{align}\]

7.8.2.1. Simulation numérique#

Pour résoudre numériquement ce système, nous devons le transformer en un système d’équations différentielles du premier ordre.

On choisit comme vecteur d’état \(U = [\theta_1, \dot{\theta_1},\theta_2, \dot{\theta_2}]\), et on écrit le système sous la forme $\( \dot{U} = \mathbf{F} U\)$

7.8.2.2. Modèle numérique#

  • \(u_0 =[\theta_1(0),0.,\theta_2(0),0.]\) conditions initiales (dim=4)

vecteur d’état \(V ~=~ [\theta_1, \dot{\theta_1}, \theta_2, \dot{\theta_2} ]\)

  • u[0] = angle du premier pendule

  • u[1] = vitesse angulaire du premier pendule

  • u[2] = angle du second pendule

  • u[3] = vitesse angulaire du second pendule

Le vector \(V ~=~ [\theta_1, u_1, \theta_2, u_2 ]\) contient les 2 angles \(\theta_1\) et \(\theta_2\) et les vitesses angulaires \(u_1\) and \(u_2\)

\[ u_1 = \frac{d \theta_1}{dt} \mbox{ et } u_2 = \frac{d \theta_2}{dt} \]

En posant \(c=cos(\theta_1 - \theta_2 )\) et \(s = sin(\theta_1 - \theta_2 )\)

\[ \frac{d u_1}{dt} ~=~ \frac{m_2 g sin(\theta_2)c - m_2 s ( L_1c u_1^2 + L_2 u_2^2 )- (m_1+m_2)g sin(\theta_1)}{ L_1 (m_1+m_2s^2)} \]
\[ \frac{d u_2}{dt} ~=~ \frac{(m_1+m_2)(L_1u_1^2s - gsin(u_2) + gsin(u_1)c) + m_2L_2u_2^2sc)}{L_2 (m_1+m_2s^2)} \]

7.8.2.3. Paramètres de l’étude#

  • paramètres \(m_1,m_2,L_1,L_2,g\)

  • calcul de la solution en N temps équi-répartis

m1 = 2                 # mass of pendulum 1 (in kg)
m2 = 1                 # mass of pendulum 2 (in kg)
L1 = 1.4               # length of pendulum 1 (in meter)
L2 = 1                 # length of pendulum 2 (in meter)
g = 9.8                # gravitatioanl acceleration constant (m/s^2)

tfinal = 25.0       # Final time. Simulation time = 0 to tfinal.
Nt = 751
t = np.linspace(0, tfinal, Nt)

7.8.2.4. Résolution numérique#

# Differential equations describing the system: function F(u)
def double_pendulum(u,t,m1,m2,L1,L2,g):
    # du = derivatives
    # u = variables
    # p = parameters
    # t = time variable
    return

# calcul solution pour une valeur initiale theta1,theta2 (en degré)
def calcul_sol(theta1, theta2):
    tfinal = 25.0       # Final time. Simulation time = 0 to tfinal.
    Nt = 751
    sol = None
    return sol

7.8.2.5. Calcul de la solution#

# tracer d'une solution
def plot_sol(sol):
    global t,L1,L2
    theta1 = sol[:,0]
    theta2 = sol[:,2]
    x1 =  L1*np.sin(theta1);          # First Pendulum
    y1 = -L1*np.cos(theta1);
    x2 = x1 + L2*np.sin(theta2);     # Second Pendulum
    y2 = y1 - L2*np.cos(theta2);
    # tracer
    plt.figure(figsize=(12,6))
    plt.subplot(1,2,1)
    plt.plot(t,theta1,label="$\\theta_1$")
    plt.plot(t,theta2,label="$\\theta_2$")
    plt.legend()
    plt.xlabel('t')
    plt.subplot(1,2,2)
    plt.plot(x1,y1,label="masse1")
    plt.plot(x2,y2,label="masse2")
    plt.legend()
    plt.xlabel('x')
    plt.axis('equal');
    return

# solution cas  stable theta1=-20 deg. theta2 = 35 deg.
### BEGIN SOLUTION
### END SOLUTION

# solution cas plus instable theta=10 deg. theta2=80 deg.
### BEGIN SOLUTION
### END SOLUTION

7.8.3. Objectif#

Peut on prédire la trajectoire du système sans faire la simulation, en utilisant un apprentissage à partir d’une base de données issue de la simulation.

7.8.3.1. prédiction#

En se donnant la valeur de solution sur 10 valeurs en temps consécutives, peut on prédire l’évolution de la solution sur des temps longs: p.e. sur les 100 pas en temps suivants

7.8.3.2. base de données#

Pour créer la base de données, on découpe la simulation avec une fenêtre glissante de 10 valeurs en temps consécutives pour prédire la valeur de la solution à l’instant t

\[ sol(t) = \mathbf{F}\left(sol(t-\Delta t),sol(t-2\Delta t),..sol(t-10\Delta t)\right) \]

7.8.4. Base de données d’apprentissage#

  • train_seq: données (data) tableau dim (10,4)

  • train_label: résultat dim (4)

création de la BD avec une technique de fenêtre glissante et d’adimensionalisation

# creation du training set
# taille de la fenetre
train_window = 10

def create_inout_sequences(input_data, tw):
    inout_seq = []

# normalisation de la base de données (avec scikit learn MinMaxScaler )
def preprocess_sol(sol, train_window):
    return

# solution numérique pour calculer la base de donnees
sol = calcul_sol(-20., 35.)
# generation de la base de données
train_inout_seq = preprocess_sol(sol, train_window)
print("BD:",len(train_inout_seq))
print("valeurs BD:",train_inout_seq[0])

7.8.5. Modèle LSTM#

réseau de neuronnes récurrents: LSTM « neural network » nn de pyTorch

paramétres

  • taille des données élémentaires input_size= 4 : 4 données / temps = \(\theta_1, \theta_2, \dot{\theta_1},\dot{\theta_2} \)

  • nbre de couches de neuronnes cachés (depliement de la récurrence en temps) hidden_size=100 pas en temps

  • optimisation: méthode d’optimisation par gradient, algorithme stochastique de gradient d’ADAM

  • utilisation des états précédents (états cn (long-term memory) et états cachés hn (short-term memory))

     lstm = nn.LSTM(input_size, hidden_size, num_LSTM=1)
 
 output, (hn,cn) = lstm(input, (h0,c0))

attention pytorch a sa propre version des tableaux (torch tensor)

  • conversion données: tableau 1D 100 -> sortie dim 4

       Linear(100, 4)

apprentissage

  • à partir de la base de données initiale: simulation (751,4)

  • fenêtrage: 741 fenêtres glissantes (10,5) -> (4)

  • on choisit au hasard 500 fenêtres pour faire l’apprentissage ( pour 1 passage)

  • on réitère le processus complet epochs fois. epochs est un hyperparamètre qui définit le nombre de fois que l’algorithme d’apprentissage parcours l’ensemble des données d’entraînement. Ce paramétre peut etre très grand !

# model LSTM
class MyLSTM(nn.Module):
    def __init__(self):
        """initialisation """
        super().__init__()
        # reseau LSTM
        self.lstm = nn.LSTM(4, 100)
        # transformation lineaire des données 
        self.linear = nn.Linear(100, 4)
        # h0,c0 valeurs initiales des etats internes (c0) et des etats cachés (h0) 
        self.c_h = (torch.zeros(1,1,100),
                    torch.zeros(1,1,100))
        return
    def forward(self, x):
        """prediction """
        h, self.c_h= self.lstm(x.view(len(x) ,1, -1), self.c_h)
        predictions = self.linear(h.view(len(x), -1))
        return predictions[-1]

7.8.5.1. Apprentissage#

# apprentissage (False): attention très tres long !!!
model = MyLSTM()
use_saved_model = True
use_saved_model = False   # recalcul
#
if use_saved_model == False:
    loss_function = nn.MSELoss()
    hidden_layer_size = 100
    optimizer = torch.optim.Adam(model.parameters(), lr=0.0005)
    epochs = 2002 #
    epochs = 52   # calcul + simple
    model_file = 'lstm-pend1'
    print("recalcul model:{} nbepochs={}".format(model_file, epochs))

# lecture du model sauvegardé
if use_saved_model:
    model = MyLSTM()
    sol = calcul_sol(-20., 35.)
    model.load_state_dict(torch.load('lstm-pend2'))  #this model is in the stable region
    #sol = calcul_sol(-50.,140.)
    #model.load_state_dict(torch.load('lstm-pend3'))  #this model is from chaotic region

7.8.5.2. Prédiction du model#

  • à partir de données data, on utilise les train_window=10 données initiales pour prédire la solution aux steps=100 pas en temps suivants

  • test sur la base de données d’apprentissage avec steps=100 et steps=200

  • test sur une autre base de données

# calcul de la prediction sur steps pas a partir de data
def OneStep(data, steps = 100):
    print('data set length =', len(data), "pred=",steps)
    train_window = 10
    # mise a l'echelle
    scaler = MinMaxScaler(feature_range=(-1, 1))
    train_data_normalized = scaler.fit_transform(data.reshape(-1, 4))
    test_inputs = train_data_normalized[0:train_window].reshape(train_window,4).tolist()
    s2 = train_data_normalized.reshape(len(data),4).tolist()
    realdata = data
    model.eval()
    preds = test_inputs.copy()
    t2 = test_inputs
    hidden_layer_size = 100
    x = 0
    #fut_pred = len(data) - train_window
    fut_pred = steps
    for i in range(fut_pred):
        seq = torch.FloatTensor(np.array(t2[i:]))
        # etat interne et cachés
        model.c_h = (torch.zeros(1, 1, hidden_layer_size),
                     torch.zeros(1, 1, hidden_layer_size))
        x = model(seq)
        preds.append(x.detach().numpy())
        t2.append(x.detach().numpy())
    actual_predictions = scaler.inverse_transform(np.array(preds).reshape(-1,4))
    # estimation erreur
    err = 0.0
    cnt = steps
    maxerr = 0.0
    maxloc = 0
    # estimation erreur sur theta2 
    for i in range(cnt):
        er =np.abs(realdata[i,2]-actual_predictions[i,2])
        err += er
        if er > maxerr:
            maxerr = er
            maxloc = i
    print("mean error =", err/cnt)
    print('maxerr =', maxerr, ' at ', maxloc)
    return actual_predictions

## calcul de la prediction a partir de la bd de test
### BEGIN SOLUTION
### END SOLUTION

# tracer de la solution
def plot_traj(data,predictions,steps=100):
    # the following will plot the lower mass path for steps using 
    # the actual ODE solver and the predicitons
    global t
    u0 = data[:,0]     # theta_1 
    u1 = data[:,1]     # omega 1
    u2 = data[:,2]     # theta_2 
    u3 = data[:,3]     # omega_2 
    up0 = predictions[:,0]     # theta_1 
    up1 = predictions[:,1]     # omega 1
    up2 = predictions[:,2]     # theta_2 
    up3 = predictions[:,3]     # omega_2 
    x1 = L1*np.sin(u0);          # First Pendulum
    y1 = -L1*np.cos(u0);
    x2 = x1 + L2*np.sin(u2);     # Second Pendulum
    y2 = y1 - L2*np.cos(u2);
    xp1 = L1*np.sin(up0);          # First Pendulum
    yp1 = -L1*np.cos(up0);
    xp2 = xp1 + L2*np.sin(up2);     # Second Pendulum
    yp2 = yp1 - L2*np.cos(up2);
    #
    plt.figure(figsize=(12,6))
    plt.subplot(1,2,1)
    plt.plot(x2[0:steps], y2[0:steps], color='r',label="test")
    plt.plot(xp2[0:steps],yp2[0:steps] , color='g',label="pred.")
    plt.plot([x2[0]],[y2[0]],'ok')
    plt.legend()
    plt.axis('equal')
    plt.title("trajectoire masse 2")
    plt.subplot(1,2,2)
    plt.plot(t[:steps],u2[:steps]-up2[:steps])
    plt.xlabel('t')
    plt.title("erreur sur $\\theta_2(t)$ ");
    return

7.8.5.3. test sur base de donnée de test#

# prediction 1 sur la base de données d'apprentissage steps=100
### BEGIN SOLUTION
### END SOLUTION

# prediction 2 sur un temps plus long: steps=200
### BEGIN SOLUTION
### END SOLUTION

7.8.5.4. prédiction d’une nouvelle simulation#

## calcul pour theta1=-25 deg. , theta2 = 30 deg.
# base entrainement
#sol1 =  calcul_sol(-20., 35.)
sol1 = None
### BEGIN SOLUTION
# base entrainement
#sol1 =  calcul_sol(-20., 35.)
### END SOLUTION

## calcul pour theta1=-30 deg. , theta2 = 40 deg.
sol1 = None
### BEGIN SOLUTION
### END SOLUTION

7.8.6. Conclusion#

Nous remarquons que l’apprentissage LSTM avec une unique solution permettait au LSTM de générer des solutions raisonnablement bonnes à partir de points proches du point initial de la solution.

7.8.7. Amélioration: Physics-Informed Neural Networks (PINN)#

  1. contrainte sur l’évolutionn du système :

  • système conservatif: Ec + Ep = cste

\[ E(\theta_1, \theta_2, \dot{\theta_1}, \dot{\theta_2}) = cste\]

  • minimisation sous contrainte

  1. Réseaux de neuronnes inspiré par la physique !!!

c’est nouvelle classe de réseaux de neurones qui hybride apprentissage automatique et lois physiques permettant une « meilleur » prédiction des réseaux de neuronnes (Deep Learning) pour des applications physiques gérées par des équations (ODE ou PDE):

7.8.8. FIN#
Contenu