Resolver EDO usando una red neuronal informada por procesos físicos

Abrir script en vivo

En este ejemplo se muestra cómo entrenar una red neuronal informada por procesos físicos (PINN) para predecir las soluciones de una ecuación diferencial ordinaria (EDO).

No todas las ecuaciones diferenciales tienen una solución en forma cerrada. Existen muchos algoritmos numéricos tradicionales para hallar soluciones aproximadas a este tipo de ecuaciones. Una red neuronal informada por procesos físicos (PINN) [1] es una red neuronal que incorpora las leyes de la física en su estructura y el proceso de entrenamiento. Puede entrenar una red neuronal que genera la solución de una EDO que define un sistema físico. Este enfoque tiene varias ventajas, entre ellas, que proporciona soluciones aproximadas diferenciables en una forma analítica cerrada [2].

Este ejemplo muestra cómo:

generar datos de entrenamiento en el intervalo $x \in [0, 2]$ ;
definir una red neuronal que tome $x$ como entrada y devuelva como salida la solución aproximada de la EDO $\dot{y} = - 2 xy$ , evaluada en $x$ ;
entrenar la red con una función de pérdida personalizada;
comparar las predicciones de la red con la solución analítica.

Ecuaciones diferenciales ordinarias y función de pérdida

En este ejemplo, resolverá la EDO

$\dot{y} = - 2 xy,$

con la condición inicial $y (0) = 1$ . Esta EDO tiene la solución analítica

$y (x) = e^{- x^{2}} .$

Defina una función de pérdida personalizada que penalice las desviaciones respecto a la satisfacción de la EDO y la condición inicial. En este ejemplo, la función de pérdida es una suma ponderada de la pérdida de la EDO y la pérdida de la condición inicial:

$L_{θ} (x) = ‖ {y_{}^{˙}}_{θ} + 2 x y_{θ} ‖^{2} + k ‖ y_{θ} (0) - 1 ‖^{2}$

$θ$ son los parámetros de la red, $k$ es un coeficiente constante, $y_{θ}$ es la solución predicha por la red, y $\dot{y_{θ}}$ es la derivada de la solución predicha calculada mediante diferenciación automática. El término $‖ \dot{y_{θ}} + 2 {xy}_{θ} ‖^{2}$ es la pérdida de la EDO y cuantifica cuánto se desvía la solución predicha de satisfacer la definición de la EDO. El término $‖ y_{θ} (0) - 1 ‖^{2}$ es la pérdida de la condición inicial y cuantifica cuánto se desvía la solución predicha de satisfacer la definición de la condición inicial.

Generar datos de entrada y definir la red

Genere 10.000 puntos de datos de entrenamiento en el intervalo $x \in [0, 2]$ .

x = linspace(0,2,10000)';

Defina la red para resolver la EDO. Como la entrada es un número real $x \in R$ , especifique un tamaño de entrada de 1.

inputSize = 1;
layers = [
    featureInputLayer(inputSize,Normalization="none")
    fullyConnectedLayer(10)
    sigmoidLayer
    fullyConnectedLayer(1)
    sigmoidLayer];

Cree un objeto dlnetwork a partir del arreglo de capas.

net = dlnetwork(layers)

net = 
  dlnetwork with properties:

         Layers: [5×1 nnet.cnn.layer.Layer]
    Connections: [4×2 table]
     Learnables: [4×3 table]
          State: [0×3 table]
     InputNames: {'input'}
    OutputNames: {'layer_2'}
    Initialized: 1

  View summary with summary.

Definir la función de pérdida del modelo

Cree la función modelLoss, que aparece en la sección Función de pérdida del modelo del ejemplo, que toma como entradas una red neuronal, un minilote de datos de entrada y el coeficiente asociado a la pérdida de la condición inicial. Esta función devuelve la pérdida y los gradientes de la pérdida con respecto a los parámetros que se pueden aprender de la red neuronal.

Especificar las opciones de entrenamiento

Entrene con un tamaño de minilote de 100 durante 15 épocas.

numEpochs = 15;
miniBatchSize = 100;

Especifique las opciones para la optimización de SGDM. Especifique una tasa de aprendizaje de 0,5, un factor de caída de la tasa de aprendizaje de 0,5, un periodo de caída de la tasa de aprendizaje de 5 y un momento de 0,9.

initialLearnRate = 0.5;
learnRateDropFactor = 0.5;
learnRateDropPeriod = 5;
momentum = 0.9;

Especifique el coeficiente del término de condición inicial en la función de pérdida como 7. Este coeficiente especifica la contribución relativa de la condición inicial a la pérdida. Ajustar este parámetro puede ayudar a que el entrenamiento converja más rápidamente.

icCoeff = 7;

Entrenar un modelo

Para utilizar minilotes de datos durante el entrenamiento:

Cree un objeto arrayDatastore a partir de los datos de entrenamiento.
Cree un objeto minibatchqueue que tome el objeto arrayDatastore como entrada, especifique el tamaño del minilote para descartar minilotes parciales y que los datos de entrenamiento tengan el formato "BC" (lote, canal).

De forma predeterminada, el objeto minibatchqueue convierte los datos en objetos dlarray con el tipo subyacente single. También convierte cada salida en gpuArray si hay una GPU disponible. Utilizar una GPU requiere Parallel Computing Toolbox™ y un dispositivo GPU compatible. Para obtener información sobre los dispositivos compatibles, consulte GPU Computing Requirements (Parallel Computing Toolbox).

ads = arrayDatastore(x,IterationDimension=1);
mbq = minibatchqueue(ads, ...
    MiniBatchSize=miniBatchSize, ...
    PartialMiniBatch="discard", ...
    MiniBatchFormat="BC");

Inicialice el parámetro de velocidad para el solver SGDM.

velocity = [];

Calcule el número total de iteraciones para monitorizar el progreso del entrenamiento.

numObservationsTrain = numel(x);
numIterationsPerEpoch = floor(numObservationsTrain / miniBatchSize);
numIterations = numEpochs * numIterationsPerEpoch;

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="LogLoss", ...
    Info=["Epoch" "LearnRate"], ...
    XLabel="Iteration");

Entrene la red con un bucle de entrenamiento personalizado. Para cada época, cambie el orden de los datos y pase en bucle por minilotes de datos. Para cada minilote:

Evalúe la pérdida y los gradientes del modelo utilizando las funciones dlfeval y modelLoss.
Actualice los parámetros de red con la función sgdmupdate.
Muestre el progreso del entrenamiento.
Detenga el proceso si la propiedad Stop está establecida como true. El valor de la propiedad Stop del objeto TrainingProgressMonitor cambia a verdadero cuando hace clic en el botón Stop.

Cada época learnRateDropPeriod, multiplica la tasa de aprendizaje por learnRateDropFactor.

epoch = 0;
iteration = 0;
learnRate = initialLearnRate;
start = tic;

% Loop over epochs.
while epoch < numEpochs  && ~monitor.Stop
    epoch = epoch + 1;

    % Shuffle data.
    mbq.shuffle

    % Loop over mini-batches.
    while hasdata(mbq) && ~monitor.Stop

        iteration = iteration + 1;

        % Read mini-batch of data.
        X = next(mbq);

        % Evaluate the model gradients and loss using dlfeval and the modelLoss function.
        [loss,gradients] = dlfeval(@modelLoss, net, X, icCoeff);

        % Update network parameters using the SGDM optimizer.
        [net,velocity] = sgdmupdate(net,gradients,velocity,learnRate,momentum);

        % Update the training progress monitor.
        recordMetrics(monitor,iteration,LogLoss=log(loss));
        updateInfo(monitor,Epoch=epoch,LearnRate=learnRate);
        monitor.Progress = 100 * iteration/numIterations;

    end
    
    % Reduce the learning rate.
    if mod(epoch,learnRateDropPeriod)==0
        learnRate = learnRate*learnRateDropFactor;
    end
end

Probar un modelo

Pruebe la precisión de la red comparando sus predicciones con la solución analítica.

Genere datos de prueba en el intervalo $x \in [0, 4]$ para ver si la red es capaz de extrapolar fuera del intervalo de entrenamiento $x \in [0, 2]$ .

xTest = linspace(0,4,1000)';

Realice predicciones con la función minibatchpredict. De forma predeterminada, la función minibatchpredict usa una GPU en caso de que esté disponible. Para utilizar una GPU se requiere una licencia de Parallel Computing Toolbox™ y un dispositivo GPU compatible. Para obtener información sobre los dispositivos compatibles, consulte GPU Computing Requirements (Parallel Computing Toolbox). De lo contrario, la función usa la CPU. Para especificar el entorno de ejecución, use la opción ExecutionEnvironment.

yModel = minibatchpredict(net,xTest);

Evalúe la solución analítica.

yAnalytic = exp(-xTest.^2);

Compare la solución analítica y la predicción del modelo representándolas en una escala logarítmica.

figure;
plot(xTest,yAnalytic,"-")
hold on
plot(xTest,yModel,"--")
legend("Analytic","Model")
xlabel("x")
ylabel("y (log scale)")
set(gca,YScale="log")

El modelo se aproxima con precisión a la solución analítica en el intervalo de entrenamiento $x \in [0, 2]$ y extrapola en el intervalo $x \in (2, 4]$ con menor precisión.

Calcule el error relativo cuadrático medio en el intervalo de entrenamiento $x \in [0, 2]$ .

yModelTrain = yModel(1:500);
yAnalyticTrain = yAnalytic(1:500);

errorTrain = mean(((yModelTrain-yAnalyticTrain)./yAnalyticTrain).^2)

errorTrain = single

0.0029

Calcule el error relativo cuadrático medio en el intervalo extrapolado $x \in (2, 4]$ .

yModelExtra = yModel(501:1000);
yAnalyticExtra = yAnalytic(501:1000);

errorExtra = mean(((yModelExtra-yAnalyticExtra)./yAnalyticExtra).^2)

errorExtra = single

6.9634e+05

Observe que el error relativo cuadrático medio es mayor en el intervalo extrapolado que en el intervalo de entrenamiento.

Función de pérdida del modelo

La función modelLoss toma como entradas un objeto dlnetwork net, un minilote de datos de entrada X y el coeficiente asociado a la pérdida de la condición inicial icCoeff. Esta función devuelve los gradientes de la pérdida con respecto a los parámetros que se pueden aprender en net y la pérdida correspondiente. La pérdida se define como una suma ponderada de la pérdida de la EDO y la pérdida de la condición inicial. La evaluación de esta pérdida requiere derivadas de segundo orden. Para habilitar la diferenciación automática de segundo orden, utilice la función dlgradient y establezca el argumento de nombre-valor EnableHigherDerivatives en true.

function [loss,gradients] = modelLoss(net, X, icCoeff)
y = forward(net,X);

% Evaluate the gradient of y with respect to x. 
% Since another derivative will be taken, set EnableHigherDerivatives to true.
dy = dlgradient(sum(y,"all"),X,EnableHigherDerivatives=true);

% Define ODE loss.
eq = dy + 2*y.*X;

% Define initial condition loss.
ic = forward(net,dlarray(0,"CB")) - 1;

% Specify the loss as a weighted sum of the ODE loss and the initial condition loss.
loss = mean(eq.^2,"all") + icCoeff * ic.^2;

% Evaluate model gradients.
gradients = dlgradient(loss, net.Learnables);

end

Referencias

Maziar Raissi, Paris Perdikaris y George Em Karniadakis, Physics Informed Deep Learning (Part I): Data-driven Solutions of Nonlinear Partial Differential Equations https://arxiv.org/abs/1711.10561
Lagaris, I. E., A. Likas, and D. I. Fotiadis. "Artificial Neural Networks for Solving Ordinary and Partial Differential Equations". IEEE Transactions on Neural Networks 9, n.º 5 (septiembre 1998): 987–1000. https://doi.org/10.1109/72.712178.

Consulte también

dlarray | dlfeval | dlgradient