Lösen partieller Differenzialgleichungen mithilfe physikbasierter neuronaler Netze

In diesem Beispiel verwendet:

Dieses Beispiel zeigt, wie Sie ein physikbasiertes neuronales Netz (PINN, physics-informed neural network) zur Vorhersage der Lösungen einer partiellen Differenzialgleichung trainieren können.

Ein physikbasiertes neuronales Netz (PINN) [1] ist ein neuronales Netz, in dessen Struktur und Trainingsprozess physikalische Gesetze integriert sind. Beispielsweise können Sie ein neuronales Netz trainieren, das die Lösung einer partiellen Differenzialgleichung ausgibt, die ein physikalisches System definiert.

In diesem Beispiel wird ein PINN trainiert, das die Beispieldatensätze $(x, t)$ als Eingaben und Ausgaben $u (x, t)$ verwendet, wobei $u$ die Lösung der folgenden Burgersgleichung ist:

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

wobei $u (x, 0) = - s i n (π x)$ die Anfangsbedingung ist und $u (- 1, t) = 0$ und $u (1, t) = 0$ die Grenzbedingungen darstellen.

Dieses Diagramm zeigt den Datenfluss durch das PINN.

Diagram of the data flow of the neural network. The input is x_1, x_2, ... x_N, and t. The output is the PDE solution u(x,t).

Für das Training dieses Modells müssen keine Daten im Voraus gesammelt werden. Sie können Daten mithilfe der Definition der partiellen Differenzialgleichung und den Randbedingungen generieren.

Generieren von Trainingsdaten

Zum Training des Modells ist ein Datensatz mit Kollokationspunkten erforderlich, die die Grenzbedingungen und Anfangsbedingungen erzwingt und die Burgersgleichung erfüllen.

Wählen Sie 25 Zeitpunkte in gleichmäßigen Abständen aus, um die Grenzbedingungen $u (- 1, t) = 0$ und $u (1, t) = 0$ zu erzwingen.

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));

Wählen Sie 50 Raumpunkte in gleichmäßigen Abständen aus, um die Anfangsbedingung $u (x, 0) = - s i n (π x)$ zu erzwingen.

numInitialConditionPoints  = 50;

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

Gruppieren Sie die Daten für die Anfangs- und Grenzbedingungen.

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

Tasten Sie gleichmäßig 10.000 Punkte $(t, x) \in (0, 1) \times (- 1, 1)$ ab, um zu erzwingen, dass die Ausgabe des Netzes die Burgersgleichung erfüllt.

numInternalCollocationPoints = 10000;

points = rand(numInternalCollocationPoints,2);

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

Definieren der Architektur des neuronalen Netzes

Definieren Sie eine mehrschichtige neuronale Perzeptron-Netzarchitektur.

Diagram of the neual network architecture. The layers are connected in series. The input is passed through 8 fully connected layers that are each proceeded by a tanh layer. The output of the final tanh layer is connected to a fully connected layer, which outputs the network predictions.

Geben Sie für die wiederholten vollständig vernetzten Schichten eine Ausgabegröße von 20 an.
Geben Sie für die letzte vollständig vernetzte Schicht eine Ausgabegröße von 1 an.

Geben Sie die Netz-Hyperparameter an.

numBlocks = 8;
fcOutputSize = 20;

Erstellen Sie das neuronale Netz. Um ein neuronales Netz mit wiederholten Schichtblöcken einfach zu konstruieren, können Sie mithilfe der Funktion repmat Schichtblöcke wiederholen.

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   Fully connected layer with output size 20
     3   ''   Tanh              Hyperbolic tangent
     4   ''   Fully Connected   Fully connected layer with output size 20
     5   ''   Tanh              Hyperbolic tangent
     6   ''   Fully Connected   Fully connected layer with output size 20
     7   ''   Tanh              Hyperbolic tangent
     8   ''   Fully Connected   Fully connected layer with output size 20
     9   ''   Tanh              Hyperbolic tangent
    10   ''   Fully Connected   Fully connected layer with output size 20
    11   ''   Tanh              Hyperbolic tangent
    12   ''   Fully Connected   Fully connected layer with output size 20
    13   ''   Tanh              Hyperbolic tangent
    14   ''   Fully Connected   Fully connected layer with output size 20
    15   ''   Tanh              Hyperbolic tangent
    16   ''   Fully Connected   Fully connected layer with output size 20
    17   ''   Tanh              Hyperbolic tangent
    18   ''   Fully Connected   Fully connected layer with output size 1

Konvertieren Sie das Schichtarray zu einem dlnetwork-Objekt.

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.

Das Training eines PINNs kann zu besserer Genauigkeit führen, wenn die lernbaren Parameter den Datentyp „Double“ aufweisen. Konvertieren Sie die lernbaren Parameter des Netzes mithilfe der Funktion dlupdate zu „Double“. Beachten Sie, dass nicht alle neuronalen Netze lernbare Elemente vom Typ „Double“ unterstützen; beispielsweise Netzwerke, die eine GPU-Optimierung verwenden, die lernbare Elemente vom Typ „Single“ benötigt.

net = dlupdate(@double,net);

Definieren der Modellverlustfunktion

Erstellen Sie die Funktion modelLoss, die das neuronale Netz, die Netzeingänge und die Anfangs- und Grenzbedingungen als Eingaben verwendet und die Verlustgradienten der lernbaren Parameter und den entsprechenden Verlust zurückgibt.

Das Modell wird trainiert, indem erzwungen wird, dass für die Eingabe $(x, t)$ die Ausgabe des Netzes $u (x, t)$ die Burgersgleichung, die Grenzbedingungen und die Anfangsbedingung erfüllt. Insbesondere zwei Größen tragen zu einem minimalen Verlust bei:

$loss = {MSE}_{f} + {MSE}_{u}$ ,

wobei ${MSE}_{f} = \frac{1}{N_{f}} \sum_{i = 1}^{N_{f}} {| f (x_{f}^{i}, t_{f}^{i}) |}^{2}$ , ${MSE}_{u} = \frac{1}{N_{u}} \sum_{i = 1}^{N_{u}} {| u (x_{u}^{i}, t_{u}^{i}) - u^{i} |}^{2}$ , die Funktion $f$ die linke Seite der Burgersgleichung darstellt und ${x_{u}^{i}, t_{u}^{i}}_{i = 1}^{N_{u}}$ den Kollokationspunkten am Rand der rechnerischen Domäne entspricht und sowohl Rand- als auch Anfangsbedingung berücksichtigt und ${x_{f}^{i}, t_{f}^{i}}_{i = 1}^{N_{f}}$ Punkte im Inneren der Domäne sind.

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

Festlegen von Trainingsoptionen

Legen Sie die Trainingsoptionen fest:

Trainieren Sie mit dem L-BFGS-Solver und verwenden Sie die Standardoptionen für den L-BFGS-Solverzustand.
Trainieren Sie über 1500 Iterationen hinweg
Stoppen Sie das Training frühzeitig, wenn die Norm der Gradienten oder der Schritte kleiner als $1 0^{- 5}$ ist.

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

Trainieren von neuronalen Netzen

Konvertieren Sie die Trainingsdaten zu dlarray-Objekten. Legen Sie fest, dass die Eingaben X und T das Format "BC" (Batch, Kanal) aufweisen und die Anfangsbedingungen das Format "CB" (Kanal, Batch) aufweisen.

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

Beschleunigen Sie die Verlustfunktion mithilfe der Funktion dlaccelerate. Weitere Informationen über die Funktionsbeschleunigung finden Sie unter Deep Learning Function Acceleration.

accfun = dlaccelerate(@modelLoss);

Erstellen Sie einen Function Handle mit der Verlustfunktion für den L-BFGS-Aktualisierungsschritt. Um die dlgradient-Funktion innerhalb der modelLoss-Funktion mit automatischer Differenzierung auszuwerten, verwenden Sie die dlfeval-Funktion.

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

Initialisieren Sie das TrainingProgressMonitor-Objekt. Plotten Sie bei den Verlust bei jeder Iteration und überwachen Sie die Norm der Gradienten und Schritte. Da der Timer startet, wenn Sie das Monitor-Objekt erstellen, achten Sie darauf, das Objekt nah an der Trainingsschleife zu erstellen.

monitor = trainingProgressMonitor( ...
    Metrics="TrainingLoss", ...
    Info=["Iteration" "GradientsNorm" "StepNorm"], ...
    XLabel="Iteration");

Trainieren Sie das Netz mit einer benutzerdefinierten Trainingsschleife mithilfe des vollständigen Datensatzes bei jeder Iteration. Bei jeder Iteration:

Aktualisieren Sie die lernbaren Netzparameter und den Solver-Zustand mithilfe der Funktion lbfgsupdate.
Aktualisieren Sie den Trainingsfortschritts-Monitor.
Stoppen Sie das Training frühzeitig, wenn die Norm der Gradienten unterhalb der angegebenen Toleranzen liegt oder die Zeilensuche fehlschlägt.

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

Evaluieren der Modellgenauigkeit

Vergleichen Sie für die Werte 0,25, 0,5, 0,75 und 1 von $t$ die vorhersagten Werte des Deep-Learning-Modells mit den tatsächlichen Lösungen der Burgersgleichung.

Legen Sie die Zielzeitpunkte fest, zu denen das Modell getestet werden soll. Berechnen Sie für jeden Zeitpunkt die Lösung für 1001 gleichmäßig verteilte Punkte im Bereich [-1,1].

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

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

Testen Sie das Modell. Sagen Sie für jede Testeingabe die Lösungen der partiellen Differenzialgleichung mithilfe des PINN vorher und vergleichen Sie die Ergenisse mit den von der Funktion solveBurgers errechneten Lösungen, aufgeführt im Abschnitt Funktion zum Lösen der Burgersgleichung im Beispiel. Um auf diese Funktion zuzugreifen, öffnen Sie das Beispiel als Live-Skript. Beurteilen Sie die Genauigkeit, indem Sie den relativen Fehler zwischen den Vorhersagen und Zielen berechnen.

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

Visualisieren Sie die Testvorhersagen in einem Diagramm.

figure
tiledlayout("flow")

for i = 1:numel(tTest)
    t = tTest(i);
    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"])

Funktion zum Lösen der Burgersgleichung

Die Funktion solveBurgers gibt die tatsächliche Lösung der Burgersgleichung zum Zeitpunkt t gemäß [2] zurück.

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

Referenzen

Maziar Raissi, Paris Perdikaris, and George Em Karniadakis, Physics Informed Deep Learning (Part I): Data-driven Solutions of Nonlinear Partial Differential Equations https://arxiv.org/abs/1711.10561
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.

Siehe auch

dlarray | dlfeval | dlgradient