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

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

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

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

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

с $u (x, t = 0) = -$ sin (øx) в качестве начального $условия, и u ($ x = - $1, t) = 0$ и u (x = 1, t) = 0 в качестве граничных условий.

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

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

Создание данных обучения

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

Выберите 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) = -$ sin (ø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];

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

numInternalCollocationPoints = 10000;

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

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

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

ds = arrayDatastore([dataX dataT]);

Определение модели глубокого обучения

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

Определение и инициализация параметров модели

Определите параметры для каждой операции и включите их в структуру. Использовать формат parameters.OperationName.ParameterName где parameters - структура, OperationName - имя операции (например, «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 принимает в качестве входных данных параметры модели и сетевые входные данные и возвращает выходные данные модели.

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

Укажите параметры обучения

Тренируйте модель на 3000 эпох с размером мини-партии 1000.

numEpochs = 3000;
miniBatchSize = 1000;

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

executionEnvironment = "auto";

Укажите параметры оптимизации ADAM.

initialLearnRate = 0.01;
decayRate = 0.005;

Железнодорожная сеть

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

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

Форматирование данных с метками размеров 'BC' (партия, канал). По умолчанию minibatchqueue объект преобразует данные в dlarray объекты с базовым типом single.
Поезд на GPU в соответствии со значением 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

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

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 сравните предсказанные значения модели глубокого обучения с истинными решениями уравнения Бургера, используя $^{}$ ошибку l2.

Задайте целевое время для тестирования модели в. Для каждого случая вычислите решение в 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')

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

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

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 =_{} {MSEf +}_{}$ MSEu,

где $_{} \frac{}{_{}}_{}^{_{}} {MSEf=1Nf∑i=1Nf'f (_{}^{} xfi,_{}^{} tfi}^{)}$ | 2 $и_{} \frac{}{_{}}_{}^{_{}} {_{}^{MSEu=1Nu∑i=1Nu'u} (_{}^{} xui,^{}}^{tui}$ ) -ui | 2.

Здесь ${_{}^{xui},_{}^{} {tui}_{} i}^{_{=}}$ 1Nu соответствуют точкам коллокации на границе вычислительной области и учитывают как границу, так и начальное ${условие}_{.}^{} {_{}^{xfi},_{}^{{tfi}_{}}}$ i = 1Nf - точки во внутренней части области.

Вычисление $_{MSEf}$ требует $\frac{}{} \frac{}{} \frac{^{}}{^{∂u∂t,∂u∂x,∂2u∂x2}}$ производных выходных данных $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

Ссылки

Мазиар Райсси, Париж Пердикарис, и Джордж Эм Карниадакис, физика информированного глубокого обучения (часть I): управляемые данными решения нелинейных дифференциальных уравнений в частных производных https://arxiv.org/abs/1711.10561
C. Basdevant, M. Deville, P. Haldenwang, J. Lacroix, J. Ouazzani, R. Peiret, P. Orlandi, A. Patera, Спектральные и конечные дифференциальные решения уравнения Бургерса, Computers & Fluids 14 (1986) 23-41.

См. также

dlarray | dlfeval | dlgradient | minibatchqueue

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