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

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

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

Учитывая вычислительную область $[- 1, 1] \times [0, 1]$ , эти примеры используют физическую информированную нейронную сеть (PINN) [1] и обучают многослойную нейронную сеть перцептрона, которая берёт выборки $(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]);

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

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

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

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

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

numLayers = 9;
numNeurons = 20;

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

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

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

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]

Задайте модели и функции градиентов модели

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

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

Настройка опций обучения

Обучите модель для 3000 эпох с мини-пакетом размером 1000.

numEpochs = 3000;
miniBatchSize = 1000;

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

executionEnvironment = "auto";

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

initialLearnRate = 0.01;
decayRate = 0.005;

Обучите сеть

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

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

Форматируйте данные с метками размерности 'BC' (пакет, канал). По умолчанию в minibatchqueue объект преобразует данные в dlarray объекты с базовым типом single.
Train на графическом процессоре в соответствии со значением 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}$ ошибка.

Установите целевое время, чтобы протестировать модель в. Для каждого времени вычислите решение с равными интервалами 1001 точек в область значений [-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')

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

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

The 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 между каждым.

Функция model принимает как вход параметры модели 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, Paris Perdikaris, and George Em Karniadakis, Physics Informed Deep Learning (Part I): управляемые данными решения нелинейных дифференциальных уравнений с частными производными https://arxiv.org/abs/1711.10561
C. Basdevant, M. Deville, P. Haldenwang, J. Lacroix, J. Ouazzani, R. Peyret, P. Orlandi, A. Patera, Спектральные и конечные дифференциальные решения уравнения Бургерса, компьютеры и жидкости 14 (1986) 23-41.

См. также

dlarray | dlfeval | dlgradient | minibatchqueue

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