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

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 parcial (EDP).

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. Por ejemplo, puede entrenar una red neuronal que genera la solución de una EDP que define un sistema físico.

En este ejemplo se entrena una PINN que toma muestras $$(x,t)$$ como entrada y genera $$u(x,t)$$, donde $$u$$ es la solución de la ecuación de Burgers.

$$\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}-\frac{0.01}{\pi }\frac{{\partial }^{2}u}{\partial x^2}=0,$$

con $$u(x,0)=-sin(\pi x)$$ como condición inicial y $$u(-1,t)=0$$ y $$u(1,t)=0$$ como condiciones límite.

Este diagrama ilustra el flujo de datos a través de la PINN.

Entrenar este modelo no requiere recopilar datos con antelación. Puede generar datos usando la definición de la EPD y las restricciones.

Generar datos de entrenamiento

Entrenar el modelo requiere un conjunto de datos de puntos de colocación que imponga las condiciones límite y las condiciones iniciales, y que cumpla la ecuación de Burgers.

Seleccione 25 puntos de tiempo equidistantes para imponer cada una de las condiciones límite $$u(-1,t)=0$$ y $$u(1,t)=0$$.

numBoundaryConditionPoints = [25 25];

x0BC1 = -1*ones(1,numBoundaryConditionPoints(1));
x0BC2 = ones(1,numBoundaryConditionPoints(2));

t0BC1 = linspace(0,1,numBoundaryConditionPoints(1));
t0BC2 = linspace(0,1,numBoundaryConditionPoints(2));

u0BC1 = zeros(1,numBoundaryConditionPoints(1));
u0BC2 = zeros(1,numBoundaryConditionPoints(2));

Seleccione 50 puntos espaciales equidistantes para imponer la condición inicial $$u(x,0)=-sin(\pi x)$$.

numInitialConditionPoints  = 50;

x0IC = linspace(-1,1,numInitialConditionPoints);
t0IC = zeros(1,numInitialConditionPoints);
u0IC = -sin(pi*x0IC);

Agrupe los datos para las condiciones inicial y límite.

X0 = [x0IC x0BC1 x0BC2];
T0 = [t0IC t0BC1 t0BC2];
U0 = [u0IC u0BC1 u0BC2];

Muestree uniformemente 10.000 puntos $$(t,x)\in (0,1)\times (-1,1)$$ para imponer que la salida de la red cumpla la ecuación de Burgers.

numInternalCollocationPoints = 10000;

points = rand(numInternalCollocationPoints,2);

dataX = 2*points(:,1)-1;
dataT = points(:,2);

Definir la arquitectura de una red neuronal

Defina una arquitectura de red neuronal de perceptrón multicapa.

• Para las capas totalmente conectadas repetidas, especifique un tamaño de salida de 20.

• Para la capa totalmente conectada final, especifique un tamaño de salida de 1.

Especifique los hiperparámetros de la red.

numBlocks = 8;
fcOutputSize = 20;

Cree la red neuronal. Para construir fácilmente una red neuronal con bloques de capas que se repiten, puede repetir bloques de capas usando la función repmat.

fcBlock = [
    fullyConnectedLayer(fcOutputSize)
    tanhLayer];

layers = [
    featureInputLayer(2)
    repmat(fcBlock,[numBlocks 1])
    fullyConnectedLayer(1)]

layers = 
  18×1 Layer array with layers:

     1   ''   Feature Input     2 features
     2   ''   Fully Connected   20 fully connected layer
     3   ''   Tanh              Hyperbolic tangent
     4   ''   Fully Connected   20 fully connected layer
     5   ''   Tanh              Hyperbolic tangent
     6   ''   Fully Connected   20 fully connected layer
     7   ''   Tanh              Hyperbolic tangent
     8   ''   Fully Connected   20 fully connected layer
     9   ''   Tanh              Hyperbolic tangent
    10   ''   Fully Connected   20 fully connected layer
    11   ''   Tanh              Hyperbolic tangent
    12   ''   Fully Connected   20 fully connected layer
    13   ''   Tanh              Hyperbolic tangent
    14   ''   Fully Connected   20 fully connected layer
    15   ''   Tanh              Hyperbolic tangent
    16   ''   Fully Connected   20 fully connected layer
    17   ''   Tanh              Hyperbolic tangent
    18   ''   Fully Connected   1 fully connected layer

Convierta el arreglo de capas en un objeto dlnetwork.

net = dlnetwork(layers)

net = 
  dlnetwork with properties:

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

  View summary with summary.

Entrenar una PINN puede dar como resultado una mayor precisión cuando los parámetros que se pueden aprender tienen el tipo de datos doble. Convierta los parámetros de aprendizaje de la red a doble con la función dlupdate. Tenga en cuenta que no todas las redes neuronales son compatibles con los parámetros de aprendizaje de tipo doble, por ejemplo, las redes que utilizan optimizaciones GPU que dependen de parámetros de aprendizaje con el tipo simple.

net = dlupdate(@double,net);

Definir la función de pérdida del modelo

Cree la función modelLoss que toma como entrada la red neuronal, las entradas de la red y las condiciones inicial y límite, y devuelve los gradientes de la pérdida con respecto a los parámetros que se pueden aprender y la correspondiente pérdida.

El modelo se entrena imponiendo que, dada una entrada $$(x,t)$$, la salida de la red $$u(x,t)$$ cumpla la ecuación de Burgers, las condiciones límite y la condición inicial. En particular, dos cantidades contribuyen a minimizar la pérdida:

$$\text{loss}={\text{MSE}}_f+{\text{MSE}}_u$$,

donde $${\text{MSE}}_f=\frac{1}{N_f}\sum _{i=1}^{N_f}{|f(x_f^i,t_f^i)|}^2$$, $${\text{MSE}}_u=\frac{1}{N_u}\sum _{i=1}^{N_u}{|u(x_u^i,t_u^i)-u^i|}^2$$, la función $$f$$ es el lado izquierdo de la ecuación de Burgers, $$\{x_u^i,t_u^i{\}}_{i=1}^{N_u}$$ corresponde a los puntos de colocación en el límite del dominio computacional y tiene en cuenta la condición límite e inicial, y $$\{x_f^i,t_f^i{\}}_{i=1}^{N_f}$$ son puntos del interior del dominio.

function [loss,gradients] = modelLoss(net,X,T,X0,T0,U0)

% Make predictions with the initial conditions.
XT = cat(1,X,T);
U = forward(net,XT);

% Calculate derivatives with respect to X and T.
X = stripdims(X);
T = stripdims(T);
U = stripdims(U);
Ux = dljacobian(U,X,1);
Ut = dljacobian(U,T,1);

% Calculate second-order derivatives with respect to X.
Uxx = dldivergence(Ux,X,1);

% Calculate mseF. Enforce Burger's equation.
f = Ut + U.*Ux - (0.01./pi).*Uxx;
mseF = mean(f.^2);

% Calculate mseU. Enforce initial and boundary conditions.
XT0 = cat(1,X0,T0);
U0Pred = forward(net,XT0);
mseU = l2loss(U0Pred,U0);

% Calculated loss to be minimized by combining errors.
loss = mseF + mseU;

% Calculate gradients with respect to the learnable parameters.
gradients = dlgradient(loss,net.Learnables);

end

Especificar las opciones de entrenamiento

Especifique las opciones de entrenamiento:

• Entrenar usando el solver L-BFGS y usar las opciones predeterminadas para el estado del solver L-BFGS.

• Entrenar durante 1500 iteraciones.

• Detener el entrenamiento antes de tiempo cuando la norma de los gradientes o los pasos sea menor que $$10^{-5}$$.

solverState = lbfgsState;
maxIterations = 1500;
gradientTolerance = 1e-5;
stepTolerance = 1e-5;

Entrenar redes neuronales

Convierta los datos de entrenamiento a objetos dlarray. Especifique que las entradas X y T tengan el formato "BC" (lote, canal) y que las condiciones iniciales tengan el formato "CB" (canal, lote).

X = dlarray(dataX,"BC");
T = dlarray(dataT,"BC");
X0 = dlarray(X0,"CB");
T0 = dlarray(T0,"CB");
U0 = dlarray(U0,"CB");

Acelere la función de pérdida usando la función dlaccelerate. Para obtener más información sobre la aceleración de funciones, consulte Deep Learning Function Acceleration for Custom Training Loops.

accfun = dlaccelerate(@modelLoss);

Cree un identificador de función que contenga la función de pérdida para el paso de actualización L-BFGS. Para evaluar la función dlgradient dentro de la función modelLoss usando diferenciación automática, utilice la función dlfeval.

lossFcn = @(net) dlfeval(accfun,net,X,T,X0,T0,U0);

Inicialice el objeto TrainingProgressMonitor. En cada iteración, represente la pérdida y monitorice la norma de los gradientes y pasos. 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="TrainingLoss", ...
    Info=["Iteration" "GradientsNorm" "StepNorm"], ...
    XLabel="Iteration");

Entrene la red con un bucle de entrenamiento personalizado usando el conjunto de datos completo en cada iteración. Para cada iteración:

• Actualice los parámetros de red que se pueden aprender y el estado del solver con la función lbfgsupdate.

• Actualice la monitorización del progreso del entrenamiento.

• Detenga el entrenamiento antes de tiempo si la norma de los gradientes es menor que las tolerancias especificadas o si la búsqueda lineal falla.

iteration = 0;
while iteration < maxIterations && ~monitor.Stop
    iteration = iteration + 1;
    [net, solverState] = lbfgsupdate(net,lossFcn,solverState);

    updateInfo(monitor, ...
        Iteration=iteration, ...
        GradientsNorm=solverState.GradientsNorm, ...
        StepNorm=solverState.StepNorm);

    recordMetrics(monitor,iteration,TrainingLoss=solverState.Loss);

    monitor.Progress = 100*iteration/maxIterations;

    if solverState.GradientsNorm < gradientTolerance || ...
            solverState.StepNorm < stepTolerance || ...
            solverState.LineSearchStatus == "failed"
        break
    end

end

Evaluar la precisión del modelo

Para valores de $$t$$ en 0.25, 0.5, 0.75 y 1, compare los valores predichos del modelo de deep learning con las soluciones verdaderas de la ecuación de Burgers.

Establezca los valores objetivo de tiempo en los que probar el modelo. Para cada valor de tiempo, calcule la solución en 1.001 puntos equidistantes en el intervalo [-1,1].

tTest = [0.25 0.5 0.75 1];
numObservationsTest = numel(tTest);

szXTest = 1001;
XTest = linspace(-1,1,szXTest);
XTest = dlarray(XTest,"CB");

Pruebe el modelo. Para cada una de las entradas de prueba, prediga las soluciones de la EDP usando la PINN y compárelas con las soluciones proporcionadas por la función solveBurgers, enumerada en la sección Función para resolver la ecuación de Burgers del ejemplo. Para acceder a esta función, abra el ejemplo como un script en vivo. Evalúe la precisión calculando el error relativo entre las predicciones y los objetivos.

UPred = zeros(numObservationsTest,szXTest);
UTest = zeros(numObservationsTest,szXTest);

for i = 1:numObservationsTest
    t = tTest(i);

    TTest = repmat(t,[1 szXTest]);
    TTest = dlarray(TTest,"CB");
    XTTest = cat(1,XTest,TTest);

    UPred(i,:) = forward(net,XTTest);
    UTest(i,:) = solveBurgers(extractdata(XTest),t,0.01/pi);
end

err = norm(UPred - UTest) / norm(UTest)

err = 
0.0237

Visualice las predicciones de prueba en una gráfica.

figure
tiledlayout("flow")

for i = 1:numel(tTest)
    nexttile
    
    plot(XTest,UPred(i,:),"-",LineWidth=2);

    hold on
    plot(XTest, UTest(i,:),"--",LineWidth=2)
    hold off

    ylim([-1.1, 1.1])
    xlabel("x")
    ylabel("u(x," + t + ")")
end

legend(["Prediction" "Target"])

Función para resolver la ecuación de Burgers

La función solveBurgers devuelve la solución verdadera de la ecuación de Burgers en los valores de tiempo t según lo descrito en [2].

function U = solveBurgers(X,t,nu)

% Define functions.
f = @(y) exp(-cos(pi*y)/(2*pi*nu));
g = @(y) exp(-(y.^2)/(4*nu*t));

% Initialize solutions.
U = zeros(size(X));

% Loop over x values.
for i = 1:numel(X)
    x = X(i);

    % Calculate the solutions using the integral function. The boundary
    % conditions in x = -1 and x = 1 are known, so leave 0 as they are
    % given by initialization of U.
    if abs(x) ~= 1
        fun = @(eta) sin(pi*(x-eta)) .* f(x-eta) .* g(eta);
        uxt = -integral(fun,-inf,inf);
        fun = @(eta) f(x-eta) .* g(eta);
        U(i) = uxt / integral(fun,-inf,inf);
    end
end

end

Referencias

1. 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

2. C. Basdevant, M. Deville, P. Haldenwang, J. Lacroix, J. Ouazzani, R. Peyret, P. Orlandi, A. Patera, Spectral and finite difference solutions of the Burgers equation, Computers & fluids 14 (1986) 23–41.

Consulte también

dlarray | dlfeval | dlgradient

Temas

• Define Custom Training Loops, Loss Functions, and Networks
• Make Predictions Using Model Function
• Specify Training Options in Custom Training Loop