Решите дифференциальные уравнения с частными производными Используя глубокое обучение

В этом примере показано, как решить уравнение Бургера с помощью глубокого обучения.

Уравнение Бургера является дифференциальным уравнением с частными производными (PDE), которое возникает в различных областях прикладной математики. В частности, гидроаэромеханика, нелинейная акустика, газовая динамика и потоки трафика.

Учитывая вычислительную область $[- 1, 1] \times [0, 1]$ , это примеры используют физику сообщила нейронной сети (PINN) [1] и обучают многоуровневую perceptron нейронную сеть, которая берет выборки $(x, t)$ как введено, где $x \in [- 1, 1]$ пространственная переменная, и $t \in [0, 1]$ переменная времени и возвращается $u (x, t)$ , где u является решением уравнения Бургера:

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

с $u (x, t = 0) = - s i n (π x)$ как начальное условие, и $u (x = - 1, t) = 0$ и $u (x = 1, t) = 0$ как граничные условия.

Пример обучает модель путем осуществления этого, учитывая вход $(x, t)$ , выход сети $u (x, t)$ выполняет уравнение Бургера, граничные условия и начальное условие.

Обучение эта модель не требует данных о сборе заранее. Можно сгенерировать данные с помощью определения УЧП и ограничений.

Сгенерируйте обучающие данные

Обучение модель требует набора данных узлов коллокации, которые осуществляют граничные условия, осуществите начальные условия и выполните уравнение Бургера.

Выберите 25 равномерно распределенных моментов времени, чтобы осуществить каждое из граничных условий $u (x = - 1, t) = 0$ и $u (x = 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));

Выберите 50 равномерно распределенных пространственных точек, чтобы осуществить начальное условие $u (x, t = 0) = - s i n (π x)$ .

numInitialConditionPoints  = 50;

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

Группа вместе данные для начальных и граничных условий.

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

Выберите 10 000 точек, чтобы осуществить выход сети, чтобы выполнить уравнение Бургера.

numInternalCollocationPoints = 10000;

pointSet = sobolset(2);
points = net(pointSet,numInternalCollocationPoints);

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

Создайте datastore массивов, содержащий обучающие данные.

ds = arrayDatastore([dataX dataT]);

Задайте модель глубокого обучения

Задайте многоуровневую perceptron архитектуру с 9 полностью операции подключения с 20 скрытыми нейронами. Первое полностью операция connect имеет два входных канала, соответствующие входным параметрам $x$ и $t$ . Последнее полностью операция connect имеет ту выход $u (x, t)$ .

Задайте и инициализируйте параметры модели

Задайте параметры для каждой из операций и включайте их в struct. Используйте формат parameters.OperationName.ParameterName где parameters struct, OperationName имя операции (например, "fc1") и ParameterName имя параметра (например, "Веса").

Задайте количество слоев и количество нейронов для каждого слоя.

numLayers = 9;
numNeurons = 20;

Инициализируйте параметры для первого полностью операция connect. Первое полностью операция connect имеет два входных канала.

parameters = struct;

sz = [numNeurons 2];
parameters.fc1.Weights = initializeHe(sz,2);
parameters.fc1.Bias = initializeZeros([numNeurons 1]);

Инициализируйте параметры для каждого остающегося промежуточного звена, полностью соединяют операции.

for layerNumber=2:numLayers-1
    name = "fc"+layerNumber;

    sz = [numNeurons numNeurons];
    numIn = numNeurons;
    parameters.(name).Weights = initializeHe(sz,numIn);
    parameters.(name).Bias = initializeZeros([numNeurons 1]);
end

Инициализируйте параметры для финала полностью операция connect. Финал полностью операция connect имеет ту выходной канал.

sz = [1 numNeurons];
numIn = numNeurons;
parameters.("fc" + numLayers).Weights = initializeHe(sz,numIn);
parameters.("fc" + numLayers).Bias = initializeZeros([1 1]);

Просмотрите сетевые параметры.

parameters

parameters = struct with fields:
    fc1: [1×1 struct]
    fc2: [1×1 struct]
    fc3: [1×1 struct]
    fc4: [1×1 struct]
    fc5: [1×1 struct]
    fc6: [1×1 struct]
    fc7: [1×1 struct]
    fc8: [1×1 struct]
    fc9: [1×1 struct]

Просмотрите параметры первого полносвязного слоя.

parameters.fc1

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

Модель Define и функции градиентов модели

Создайте функциональный model, перечисленный в разделе Model Function в конце примера, который вычисляет выходные параметры модели глубокого обучения. Функциональный model берет в качестве входа параметры модели и сетевые входные параметры, и возвращает выходной параметр модели.

Создайте функциональный modelGradients, перечисленный в разделе Model Gradients Function в конце примера, который берет в качестве входа параметры модели, сетевые входные параметры и начальные и граничные условия, и возвращает градиенты потери относительно настраиваемых параметров и соответствующей потери.

Задайте опции обучения

Обучите модель в течение 3 000 эпох с мини-пакетным размером 1 000.

numEpochs = 3000;
miniBatchSize = 1000;

Чтобы обучаться на графическом процессоре, если вы доступны, задайте среду выполнения "auto". Используя графический процессор требует Parallel Computing Toolbox™ и поддерживаемого устройства графического процессора. Для получения информации о поддерживаемых устройствах смотрите Поддержку графического процессора Релизом (Parallel Computing Toolbox) (Parallel Computing Toolbox).

executionEnvironment = "auto";

Задайте опции оптимизации ADAM.

initialLearnRate = 0.01;
decayRate = 0.005;

Обучение сети

Обучите сеть с помощью пользовательского учебного цикла.

Создайте minibatchqueue возразите, что процессы и управляют мини-пакетами данных во время обучения. Для каждого мини-пакета:

Формат данные с размерностью маркирует 'BC' (пакет, канал). По умолчанию, minibatchqueue объект преобразует данные в dlarray объекты с базовым типом single.
Обучайтесь на графическом процессоре согласно значению executionEnvironment переменная. По умолчанию, minibatchqueue объект преобразует каждый выход в gpuArray если графический процессор доступен.

mbq = minibatchqueue(ds, ...
    'MiniBatchSize',miniBatchSize, ...
    'MiniBatchFormat','BC', ...
    'OutputEnvironment',executionEnvironment);

Преобразуйте начальные и граничные условия в dlarray. Для точек входных данных задайте формат с размерностями 'CB' (образуйте канал, пакет).

dlX0 = dlarray(X0,'CB');
dlT0 = dlarray(T0,'CB');
dlU0 = dlarray(U0);

Если обучение с помощью графического процессора, преобразуйте начальную букву и условия к gpuArray.

if (executionEnvironment == "auto" && canUseGPU) || (executionEnvironment == "gpu")
    dlX0 = gpuArray(dlX0);
    dlT0 = gpuArray(dlT0);
    dlU0 = gpuArray(dlU0);
end

Инициализируйте параметры для решателя Адама.

averageGrad = [];
averageSqGrad = [];

Ускорьте функцию градиентов модели использование dlaccelerate функция. Чтобы узнать больше, смотрите, Ускоряют Пользовательские Учебные Функции Цикла.

accfun = dlaccelerate(@modelGradients);

Инициализируйте график процесса обучения.

figure
C = colororder;
lineLoss = animatedline('Color',C(2,:));
ylim([0 inf])
xlabel("Iteration")
ylabel("Loss")
grid on

Обучите сеть.

Для каждой итерации:

Считайте мини-пакет данных из мини-пакетной очереди
Оцените градиенты модели и потерю с помощью ускоренных градиентов модели и dlfeval функции.
Обновите скорость обучения.
Обновите настраиваемые параметры с помощью adamupdate функция.

В конце каждой эпохи обновите учебный график со значениями потерь.

start = tic;

iteration = 0;

for epoch = 1:numEpochs
    reset(mbq);

    while hasdata(mbq)
        iteration = iteration + 1;

        dlXT = next(mbq);
        dlX = dlXT(1,:);
        dlT = dlXT(2,:);

        % Evaluate the model gradients and loss using dlfeval and the
        % modelGradients function.
        [gradients,loss] = dlfeval(accfun,parameters,dlX,dlT,dlX0,dlT0,dlU0);

        % Update learning rate.
        learningRate = initialLearnRate / (1+decayRate*iteration);

        % Update the network parameters using the adamupdate function.
        [parameters,averageGrad,averageSqGrad] = adamupdate(parameters,gradients,averageGrad, ...
            averageSqGrad,iteration,learningRate);
    end

    % Plot training progress.
    loss = double(gather(extractdata(loss)));
    addpoints(lineLoss,iteration, loss);

    D = duration(0,0,toc(start),'Format','hh:mm:ss');
    title("Epoch: " + epoch + ", Elapsed: " + string(D) + ", Loss: " + loss)
    drawnow
end

Проверяйте эффективность ускоренной функции путем проверки частоты успешных обращений и уровня заполняемости.

accfun

accfun = 
  AcceleratedFunction with properties:

          Function: @modelGradients
           Enabled: 1
         CacheSize: 50
           HitRate: 99.9984
         Occupancy: 2
         CheckMode: 'none'
    CheckTolerance: 1.0000e-04

Оцените точность модели

Для значений $t$ в 0,25, 0.5, 0.75, и 1, сравнивают ожидаемые значения модели глубокого обучения с истинными решениями уравнения Бургера с помощью $l^{2}$ ошибка.

Установите целевые времена тестировать модель в. В течение каждого раза вычислите решение в 1 001 равномерно распределенной точке в области значений [-1,1].

tTest = [0.25 0.5 0.75 1];
numPredictions = 1001;
XTest = linspace(-1,1,numPredictions);

figure

for i=1:numel(tTest)
    t = tTest(i);
    TTest = t*ones(1,numPredictions);

    % Make predictions.
    dlXTest = dlarray(XTest,'CB');
    dlTTest = dlarray(TTest,'CB');
    dlUPred = model(parameters,dlXTest,dlTTest);

    % Calcualte true values.
    UTest = solveBurgers(XTest,t,0.01/pi);

    % Calculate error.
    err = norm(extractdata(dlUPred) - UTest) / norm(UTest);

    % Plot predictions.
    subplot(2,2,i)
    plot(XTest,extractdata(dlUPred),'-','LineWidth',2);
    ylim([-1.1, 1.1])

    % Plot true values.
    hold on
    plot(XTest, UTest, '--','LineWidth',2)
    hold off

    title("t = " + t + ", Error = " + gather(err));
end

subplot(2,2,2)
legend('Predicted','True')

Графики показывают, как близко предсказания к истинным значениям.

Решите функцию уравнения бургера

solveBurgers функция возвращает истинное решение уравнения Бургера во времена t как обрисовано в общих чертах в [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

Функция градиентов модели

Модель обучена путем осуществления этого, учитывая вход $(x, t)$ выход сети $u (x, t)$ выполняет уравнение Бургера, граничные условия и начальное условие. В частности, два количества способствуют потере, которая будет минимизирована:

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

где ${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}$ .

Здесь, ${x_{u}^{i}, t_{u}^{i}}_{i = 1}^{N_{u}}$ соответствуйте узлам коллокации на контуре вычислительной области и счета и граничное и начальное условие. ${x_{f}^{i}, t_{f}^{i}}_{i = 1}^{N_{f}}$ точки во внутренней части области.

Вычисление ${MSE}_{f}$ требует производных $\frac{\partial u}{\partial t}, \frac{\partial u}{\partial x}, \frac{\partial^{2} u}{\partial x^{2}}$ из выхода $u$ из модели.

Функциональный modelGradients берет в качестве входа, параметры модели parameters, сеть вводит dlX и dlT, начальные и граничные условия dlX0, dlT0, и dlU0, и возвращает градиенты потери относительно настраиваемых параметров и соответствующей потери.

function [gradients,loss] = modelGradients(parameters,dlX,dlT,dlX0,dlT0,dlU0)

% Make predictions with the initial conditions.
U = model(parameters,dlX,dlT);

% Calculate derivatives with respect to X and T.
gradientsU = dlgradient(sum(U,'all'),{dlX,dlT},'EnableHigherDerivatives',true);
Ux = gradientsU{1};
Ut = gradientsU{2};

% Calculate second-order derivatives with respect to X.
Uxx = dlgradient(sum(Ux,'all'),dlX,'EnableHigherDerivatives',true);

% Calculate lossF. Enforce Burger's equation.
f = Ut + U.*Ux - (0.01./pi).*Uxx;
zeroTarget = zeros(size(f), 'like', f);
lossF = mse(f, zeroTarget);

% Calculate lossU. Enforce initial and boundary conditions.
dlU0Pred = model(parameters,dlX0,dlT0);
lossU = mse(dlU0Pred, dlU0);

% Combine losses.
loss = lossF + lossU;

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

end

Функция модели

Модель, обученная в этом примере, состоит из ряда полностью операций подключения с tanh операцией между каждым.

Функция модели берет в качестве входа параметры модели parameters и сеть вводит dlX и dlT, и возвращает выход dlU модели.

function dlU = model(parameters,dlX,dlT)

dlXT = [dlX;dlT];
numLayers = numel(fieldnames(parameters));

% First fully connect operation.
weights = parameters.fc1.Weights;
bias = parameters.fc1.Bias;
dlU = fullyconnect(dlXT,weights,bias);

% tanh and fully connect operations for remaining layers.
for i=2:numLayers
    name = "fc" + i;

    dlU = tanh(dlU);

    weights = parameters.(name).Weights;
    bias = parameters.(name).Bias;
    dlU = fullyconnect(dlU, weights, bias);
end

end

Ссылки

Maziar Raissi, Париж Пердикэрис и Джордж Эм Карниадакис, Физика Информированное Глубокое обучение (Первая часть): управляемые данными Решения Нелинейных Дифференциальных уравнений с частными производными https://arxiv.org/abs/1711.10561
К. Бэсдевэнт, М. Девилл, П. Холденванг, Ж. Лакруа, Дж. Оуэззэни, Р. Пеирет, П. Орланди, А. Патера, решения для Спектральной и конечной разности уравнения Burgers, Компьютеров & жидкостей 14 (1986) 23–41.

Документация