Modelado de sistemas dinámicos usando EDO neuronales

En este ejemplo se muestra cómo entrenar una red neuronal con ecuaciones diferenciales ordinarias (EDO) para aprender la dinámica de un sistema físico.

Las EDO neuronales [1] son operaciones de deep learning definidas por la solución de una EDO. Más concretamente, una EDO neuronal es una operación que se puede utilizar en cualquier arquitectura y, dada una entrada, define su salida como la solución numérica de la EDO

$$y^{\prime }=f(t,y,\theta )$$

para el horizonte temporal $$(t_0,t_1)$$ y la condición inicial $$y(t_0)=y_0$$. El lado derecho $$f(t,y,\theta )$$ de la EDO depende de un conjunto de parámetros entrenables $$\theta$$, que el modelo aprende durante el proceso de entrenamiento. En este ejemplo, $$f(t,y,\theta )$$ se modela con una función de modelo que contiene operaciones totalmente conectadas y activaciones no lineales. La condición inicial $$y_0$$ es la entrada de toda la arquitectura, como en el caso del ejemplo, o la salida de una operación anterior.

En este ejemplo se muestra cómo entrenar una red neuronal con EDO neuronales para aprender la dinámica $$x$$ de un sistema físico dado, descrito por la siguiente EDO:

$$x^{\prime }=Ax$$,

donde $$A$$ es una matriz de 2 por 2.

La red neuronal de este ejemplo toma como entrada una condición inicial y calcula la solución de la EDO a través del modelo de EDO neuronal aprendido.

La operación con la EDO neuronal, dada una condición inicial, genera como salida la solución de un modelo de EDO. En este ejemplo, especifique un bloque con una capa totalmente conectada, una capa tanh y otra capa totalmente conectada como el modelo de EDO.

En este ejemplo, la EDO que define el modelo se resuelve numéricamente con el par Runge-Kutta explícito (4,5) de Dormand y Prince [2]. El paso hacia atrás utiliza la diferenciación automática para aprender los parámetros entrenables $$\theta$$ mediante la retropropagación por cada operación del solver EDO.

La función aprendida $$f(t,y,\theta )$$ se utiliza como el lado derecho para calcular la solución del mismo modelo para condiciones iniciales adicionales.

Sintetizar datos de dinámicas objetivo

Defina la dinámica objetivo como un modelo de EDO lineal $$x^{\prime }=Ax$$, con x0 como la condición inicial, y calcule su solución numérica xTrain con ode45 en el intervalo de tiempo [0 15]. Para calcular un dato de validación preciso, establezca la tolerancia relativa del solver numérico ode45 en ${10}^{-7}$. Después, puede utilizar el valor de xTrain como dato de validación para aprender una dinámica aproximada con un modelo de EDO neuronal.

x0 = [2; 0];
A = [-0.1 -1; 1 -0.1];
trueModel = @(t,y) A*y;

numTimeSteps = 2000;
T = 15;
odeOptions = odeset(RelTol=1.e-7);
t = linspace(0, T, numTimeSteps);
[~, xTrain] = ode45(trueModel, t, x0, odeOptions);
xTrain = xTrain';

Visualice los datos de entrenamiento en una gráfica.

figure
plot(xTrain(1,:),xTrain(2,:))
title("Ground Truth Dynamics") 
xlabel("x(1)") 
ylabel("x(2)")
grid on

Definir e inicializar parámetros de modelos

La función de modelo consta de una sola llamada a dlode45 para resolver la EDO definida por la dinámica aproximada $$f(t,y,\theta )$$ durante 40 unidades de tiempo.

neuralOdeTimesteps = 40;
dt = t(2);
timesteps = (0:neuralOdeTimesteps)*dt;

Defina los parámetros que se pueden aprender para usarlos en la llamada a dlode45 y recopílelos en la variable neuralOdeParameters. La función initializeGlorot toma como entrada el tamaño de los parámetros que se pueden aprender sz y el número de salidas y de entradas de las operaciones totalmente conectadas, y devuelve un objeto dlarray con el tipo subyacente single con valores establecidos usando la inicialización Glorot. La función initializeZeros toma como entrada el tamaño de los parámetros que se pueden aprender y devuelve los parámetros como un objeto dlarray con el tipo subyacente single. Las funciones de ejemplo de inicialización están incluidas en este ejemplo como archivos de soporte. Para acceder a estas funciones, abra este ejemplo como un script en vivo. Para obtener más información sobre cómo inicializar parámetros que se pueden aprender para funciones de modelo, consulte Initialize Learnable Parameters for Model Function.

Inicialice la estructura de los parámetros.

neuralOdeParameters = struct;

Inicialice los parámetros para las operaciones totalmente conectadas en el modelo de EDO. La primera operación totalmente conectada toma como entrada un vector de tamaño stateSize y aumenta su longitud hasta hiddenSize. A la inversa, la segunda operación totalmente conectada toma como entrada un vector de longitud hiddenSize y disminuye su longitud hasta stateSize.

stateSize = size(xTrain,1);
hiddenSize = 20;

neuralOdeParameters.fc1 = struct;
sz = [hiddenSize stateSize];
neuralOdeParameters.fc1.Weights = initializeGlorot(sz, hiddenSize, stateSize);
neuralOdeParameters.fc1.Bias = initializeZeros([hiddenSize 1]);

neuralOdeParameters.fc2 = struct;
sz = [stateSize hiddenSize];
neuralOdeParameters.fc2.Weights = initializeGlorot(sz, stateSize, hiddenSize);
neuralOdeParameters.fc2.Bias = initializeZeros([stateSize 1]);

Muestre los parámetros del modelo que se pueden aprender.

neuralOdeParameters.fc1

ans = struct with fields:
    Weights: [20×2 dlarray]
       Bias: [20×1 dlarray]

neuralOdeParameters.fc2

ans = struct with fields:
    Weights: [2×20 dlarray]
       Bias: [2×1 dlarray]

Definir modelos de EDO neuronales

Cree la función odeModel, enumerada en la sección Modelo de EDO del ejemplo, que toma como entrada la entrada de tiempo (sin utilizar), la solución correspondiente y los parámetros de la función de EDO. La función aplica una operación totalmente conectada, una operación tanh, y otra operación totalmente conectada a los datos de entrada usando los pesos y sesgos dados por los parámetros.

Definir la función de modelo

Cree la función model, enumerada en la sección Función de modelo del ejemplo, que calcula las salidas del modelo de deep learning. La función model toma como entrada los parámetros del modelo y los datos de entrada. La función genera como salida la solución de la EDO neuronal.

Definir la función de pérdida del modelo

Cree la función modelLoss, enumerada en la sección Función de pérdida del modelo del ejemplo, que toma como entrada los parámetros del modelo, un minilote de datos de entrada con los objetivos correspondientes, y devuelve la pérdida y los gradientes de la pérdida con respecto a los parámetros que se pueden aprender.

Especificar las opciones de entrenamiento

Especifique las opciones para la optimización de Adam.

gradDecay = 0.9;
sqGradDecay = 0.999;
learnRate = 0.002;

Entrene durante 1200 iteraciones con un minilote de tamaño 200.

numIter = 1200;
miniBatchSize = 200;

Cada 50 iteraciones, resuelva la dinámica aprendida y muéstrela frente a la validación en un diagrama de fases para mostrar la ruta de entrenamiento.

plotFrequency = 50;

Entrenar modelos usando un bucle de entrenamiento personalizado

Inicialice los parámetros averageGrad y averageSqGrad para el solver Adam.

averageGrad = [];
averageSqGrad = [];

Inicialice el objeto TrainingProgressMonitor. Dado que el cronómetro empieza cuando crea el objeto de monitorización, asegúrese de crear el objeto cerca del bucle de entrenamiento.

monitor = trainingProgressMonitor(Metrics="Loss",Info=["Iteration","LearnRate"],XLabel="Iteration");

Entrene la red con un bucle de entrenamiento personalizado.

Para cada iteración:

• Cree un minilote de datos a partir de los datos sintetizados con la función createMiniBatch, enumerada en la sección Función de creación de minilotes del ejemplo.

• Evalúe la pérdida y los gradientes del modelo utilizando la función dlfeval y la función modelLoss, enumerada en la sección Función de pérdida del modelo del ejemplo.

• Actualice los parámetros del modelo con la función adamupdate.

• Actualice la gráfica del progreso del entrenamiento.

numTrainingTimesteps = numTimeSteps;
trainingTimesteps = 1:numTrainingTimesteps;
plottingTimesteps = 2:numTimeSteps;

iteration = 0;

while iteration < numIter && ~monitor.Stop
    iteration = iteration + 1;

    % Create batch
    [X, targets] = createMiniBatch(numTrainingTimesteps, neuralOdeTimesteps, miniBatchSize, xTrain);

    % Evaluate network and compute loss and gradients
    [loss,gradients] = dlfeval(@modelLoss,timesteps,X,neuralOdeParameters,targets);

    % Update network
    [neuralOdeParameters,averageGrad,averageSqGrad] = adamupdate(neuralOdeParameters,gradients,averageGrad,averageSqGrad,iteration,...
        learnRate,gradDecay,sqGradDecay);

    % Plot loss
    recordMetrics(monitor,iteration,Loss=loss);

    % Plot predicted vs. real dynamics
    if mod(iteration,plotFrequency) == 0  || iteration == 1

        % Use ode45 to compute the solution 
        y = dlode45(@odeModel,t,dlarray(x0),neuralOdeParameters,DataFormat="CB");

        plot(xTrain(1,plottingTimesteps),xTrain(2,plottingTimesteps),"r--")

        hold on
        plot(y(1,:),y(2,:),"b-")
        hold off

        xlabel("x(1)")
        ylabel("x(2)")
        title("Predicted vs. Real Dynamics")
        legend("Training Ground Truth", "Predicted")

        drawnow
    end
    updateInfo(monitor,Iteration=iteration,LearnRate=learnRate);
    monitor.Progress = 100*iteration/numIter;
end

Evaluar modelos

Utilice el modelo para calcular soluciones aproximadas con distintas condiciones iniciales.

Defina cuatro condiciones iniciales nuevas distintas de las utilizadas para entrenar el modelo.

tPred = t;

x0Pred1 = sqrt([2;2]);
x0Pred2 = [-1;-1.5];
x0Pred3 = [0;2];
x0Pred4 = [-2;0];

Resuelva numéricamente la dinámica verdadera de la EDO con ode45 para las cuatro condiciones iniciales nuevas.

[~, xTrue1] = ode45(trueModel, tPred, x0Pred1, odeOptions);
[~, xTrue2] = ode45(trueModel, tPred, x0Pred2, odeOptions);
[~, xTrue3] = ode45(trueModel, tPred, x0Pred3, odeOptions);
[~, xTrue4] = ode45(trueModel, tPred, x0Pred4, odeOptions);

Resuelva numéricamente la EDO con la dinámica de EDO neuronales aprendida.

xPred1 = dlode45(@odeModel,tPred,dlarray(x0Pred1),neuralOdeParameters,DataFormat="CB");
xPred2 = dlode45(@odeModel,tPred,dlarray(x0Pred2),neuralOdeParameters,DataFormat="CB");
xPred3 = dlode45(@odeModel,tPred,dlarray(x0Pred3),neuralOdeParameters,DataFormat="CB");
xPred4 = dlode45(@odeModel,tPred,dlarray(x0Pred4),neuralOdeParameters,DataFormat="CB");

Visualizar predicciones

Visualice las soluciones predichas para distintas condiciones iniciales frente a las soluciones de la validación con la función plotTrueAndPredictedSolutions, enumerada en la sección Función de representación de soluciones verdaderas y predichas del ejemplo.

figure
subplot(2,2,1)
plotTrueAndPredictedSolutions(xTrue1, xPred1);
subplot(2,2,2)
plotTrueAndPredictedSolutions(xTrue2, xPred2);
subplot(2,2,3)
plotTrueAndPredictedSolutions(xTrue3, xPred3);
subplot(2,2,4)
plotTrueAndPredictedSolutions(xTrue4, xPred4);

Funciones de ayuda

Función de modelo

La función model, que define la red neuronal utilizada para hacer predicciones, consta de una única llamada a una EDO neuronal. Para cada observación, esta función toma un vector de longitud stateSize, que se utiliza como condición inicial para resolver numéricamente la EDO con la función odeModel, que representa el lado derecho que se puede aprender $$f(t,y,\theta )$$ de la EDO que se desea resolver como el lado derecho y un vector de puntos de tiempo tspan que define la hora a la que se genera como salida la solución numérica. La función usa el vector tspan para cada observación, independientemente de la condición inicial, ya que el sistema aprendido es autónomo. Es decir, la función odeModel no depende de manera explícita del tiempo.

function X = model(tspan,X0,neuralOdeParameters)

X = dlode45(@odeModel,tspan,X0,neuralOdeParameters,DataFormat="CB");

end

Modelo de EDO

La función odeModel es el lado derecho que se puede aprender utilizado en la llamada a dlode45. Toma como entrada un vector de tamaño stateSize, lo amplía hasta que tenga la longitud hiddenSize y aplica una función de no linealidad tanh. Después, la función comprime el vector de nuevo para que tenga la longitud stateSize.

function y = odeModel(~,y,theta)

y = tanh(theta.fc1.Weights*y + theta.fc1.Bias);
y = theta.fc2.Weights*y + theta.fc2.Bias;

end

Función de pérdida del modelo

Esta función toma como entradas un vector tspan, un conjunto de condiciones iniciales X0, los parámetros que se pueden aprender neuralOdeParameters y secuencias objetivo targets. Calcula las predicciones con la función model y las compara con las secuencias objetivo dadas. Por último, calcula la pérdida y el gradiente de la pérdida con respecto a los parámetros que se pueden aprender de la EDO neuronal.

function [loss,gradients] = modelLoss(tspan,X0,neuralOdeParameters,targets)

% Compute predictions.
X = model(tspan,X0,neuralOdeParameters);

% Compute L1 loss.
loss = l1loss(X,targets,NormalizationFactor="all-elements",DataFormat="CBT");

% Compute gradients.
gradients = dlgradient(loss,neuralOdeParameters);

end

Función de creación de minilotes

La función createMiniBatch crea un lote de observaciones de la dinámica objetivo. Toma como entrada el número total de unidades de tiempo de los datos de validación numTimesteps, el número de unidades de tiempo consecutivas que se desea devolver para cada observación numTimesPerObs, el número de observaciones miniBatchSize y los datos de validación X.

function [x0, targets] = createMiniBatch(numTimesteps,numTimesPerObs,miniBatchSize,X)

% Create batches of trajectories.
s = randperm(numTimesteps - numTimesPerObs, miniBatchSize);

x0 = dlarray(X(:, s));
targets = zeros([size(X,1) miniBatchSize numTimesPerObs]);

for i = 1:miniBatchSize
    targets(:, i, 1:numTimesPerObs) = X(:, s(i) + 1:(s(i) + numTimesPerObs));
end

end

Función de representación de soluciones verdaderas y predichas

La función plotTrueAndPredictedSolutions toma como entrada la solución verdadera xTrue, la solución aproximada xPred calculada con el modelo de EDO neuronal aprendido y la condición inicial correspondiente x0Str. Calcula el error entre las soluciones verdadera y predicha y lo representa en un diagrama de fases.

function plotTrueAndPredictedSolutions(xTrue,xPred)

xPred = squeeze(xPred)';

err = mean(abs(xTrue(2:end,:) - xPred), "all");

plot(xTrue(:,1),xTrue(:,2),"r--",xPred(:,1),xPred(:,2),"b-",LineWidth=1)

title("Absolute Error = " + num2str(err,"%.4f"))
xlabel("x(1)")
ylabel("x(2)")

xlim([-2 3])
ylim([-2 3])

legend("Ground Truth","Predicted")

end

[1] Chen, Ricky T. Q., Yulia Rubanova, Jesse Bettencourt y David Duvenaud. “Neural Ordinary Differential Equations.” Preimpresión, presentado el 13 de diciembre de 2019. https://arxiv.org/abs/1806.07366.

[2] Shampine, Lawrence F. y Mark W. Reichelt. “The MATLAB ODE Suite.” SIAM Journal on Scientific Computing 18, n.º 1 (enero 1997): 1–22. https://doi.org/10.1137/S1064827594276424.

Copyright 2021–2023 The MathWorks, Inc.

Consulte también

dlode45 | dlarray | dlgradient | dlfeval | adamupdate

Temas

• Train Neural ODE Network
• Train Latent ODE Network with Irregularly Sampled Time-Series Data
• Specify Training Options in Custom Training Loop
• Entrenar una red con un bucle de entrenamiento personalizado
• Define Custom Training Loops, Loss Functions, and Networks
• List of Functions with dlarray Support