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

В этом примере показано, как обучить Физику информированную нейронную сеть (PINN) численно вычислять решение уравнения Бургера при помощи ограниченной памяти BFGS (LBFGS) алгоритм.

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

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

ut+uux-0.01π2ux2=0,

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

Вместо обучения сеть с помощью trainNetwork функция или использование пользовательского учебного цикла, который обновляет параграницы с помощью sgdmupdate или подобные функции, этот пример оценивает настраиваемые параметры при помощи fmincon функция (требует Optimization Toolbox™). fmincon функция находит минимум ограниченных нелинейных многомерных функций.

Пример обучает модель путем осуществления этого, учитывая вход (x,t), выход сети u(x,t) выполняет уравнение Бургера, граничные условия и начальное условие. Чтобы обучить модель, экс-клен использует ограниченную память BFGS (LBFGS) алгоритм, который является приближенным методом ньютона, который аппроксимирует алгоритм Бройдена Флетчера Голдфарба Шэнно.

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

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

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

Выберите 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];

Выберите 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).

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

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

Алгоритм в этом примере требует, чтобы настраиваемые параметры были на первом уровне stucture, не используйте вложенные структуры на этом шаге. fmincon функция требует, чтобы learnable, чтобы быть удвоилось.

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

numLayers = 9;
numNeurons = 20;

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

parameters = struct;

sz = [numNeurons 2];
parameters.fc1_Weights = initializeHe(sz,2,'double');
parameters.fc1_Bias = initializeZeros([numNeurons 1],'double');

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

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

    sz = [numNeurons numNeurons];
    numIn = numNeurons;
    parameters.(name + "_Weights") = initializeHe(sz,numIn,'double');
    parameters.(name + "_Bias") = initializeZeros([numNeurons 1],'double');
end

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

sz = [1 numNeurons];
numIn = numNeurons;
parameters.("fc" + numLayers + "_Weights") = initializeHe(sz,numIn,'double');
parameters.("fc" + numLayers + "_Bias") = initializeZeros([1 1],'double');

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

parameters
parameters = struct with fields:
    fc1_Weights: [20×2 dlarray]
       fc1_Bias: [20×1 dlarray]
    fc2_Weights: [20×20 dlarray]
       fc2_Bias: [20×1 dlarray]
    fc3_Weights: [20×20 dlarray]
       fc3_Bias: [20×1 dlarray]
    fc4_Weights: [20×20 dlarray]
       fc4_Bias: [20×1 dlarray]
    fc5_Weights: [20×20 dlarray]
       fc5_Bias: [20×1 dlarray]
    fc6_Weights: [20×20 dlarray]
       fc6_Bias: [20×1 dlarray]
    fc7_Weights: [20×20 dlarray]
       fc7_Bias: [20×1 dlarray]
    fc8_Weights: [20×20 dlarray]
       fc8_Bias: [20×1 dlarray]
    fc9_Weights: [1×20 dlarray]
       fc9_Bias: [1×1 dlarray]

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

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

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

Задайте fmincon Целевая функция

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

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

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

  • Оптимизируйте использование fmincon optmizer с алгоритмом LBFGS для не больше, чем 7 500 итераций и вычислений функции.

  • Оцените с допуском оптимальности 1e-5.

  • Предоставьте градиенты для алгоритма.

options = optimoptions('fmincon', ...
    'HessianApproximation','lbfgs', ...
    'MaxIterations',7500, ...
    'MaxFunctionEvaluations',7500, ...
    'OptimalityTolerance',1e-5, ...
    'SpecifyObjectiveGradient',true);

Обучите сеть Используя fmincon

Обучите сеть с помощью fmincon функция.

fmincon функция требует, чтобы настраиваемые параметры были заданы как вектор. Преобразуйте параметры в вектор с помощью paramsStructToVector функция (присоединенный к этому примеру как вспомогательный файл). Чтобы преобразовать назад в структуру параметров, также возвратите названия параметра и размеры.

[parametersV,parameterNames,parameterSizes] = parameterStructToVector(parameters);
parametersV = extractdata(parametersV);

Преобразуйте обучающие данные в dlarray объекты с форматом 'CB' (образуйте канал, пакет).

dlX = dlarray(dataX','CB');
dlT = dlarray(dataT','CB');
dlX0 = dlarray(X0,'CB');
dlT0 = dlarray(T0,'CB');
dlU0 = dlarray(U0,'CB');

Создайте указатель на функцию с одним входом, который задает целевую функцию.

objFun = @(parameters) objectiveFunction(parameters,dlX,dlT,dlX0,dlT0,dlU0,parameterNames,parameterSizes);

Обновите настраиваемые параметры с помощью fmincon функция. В зависимости от количества итераций это может требовать времени к запущенному. Чтобы включить подробный многословный выход, установите 'Display' опция оптимизации к 'iter-detailed'.

parametersV = fmincon(objFun,parametersV,[],[],[],[],[],[],[],options);
Solver stopped prematurely.

fmincon stopped because it exceeded the function evaluation limit,
options.MaxFunctionEvaluations = 7.500000e+03.

Для предсказания преобразуйте вектор из параметров к структуре с помощью parameterVectorToStruct функция (присоединенный к этому примеру как вспомогательный файл).

parameters = parameterVectorToStruct(parametersV,parameterNames,parameterSizes);

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

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

Установите целевые времена тестировать модель в. В течение каждого раза вычислите решение в 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

fmincon Целевая функция

objectiveFunction функция задает целевую функцию, которая будет использоваться алгоритмом LBFGS.

Этот objectiveFunction funtion берет в качестве входа, вектора из настраиваемых параметров parametersV, сетевые входные параметры, dlX и dlT, начальные и граничные условия dlX0, dlT0, и dlU0, и имена и размеры настраиваемых параметров parameterNames и parameterSizes, соответственно, и возвращает потерю, которая будет минимизирована fmincon функциональный loss и вектор, содержащий градиенты потери относительно настраиваемых параметров gradientsV.

function [loss,gradientsV] = objectiveFunction(parametersV,dlX,dlT,dlX0,dlT0,dlU0,parameterNames,parameterSizes)

% Convert parameters to structure of dlarray objects.
parametersV = dlarray(parametersV);
parameters = parameterVectorToStruct(parametersV,parameterNames,parameterSizes);

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

% Return loss and gradients for fmincon.
gradientsV = parameterStructToVector(gradients);
gradientsV = extractdata(gradientsV);
loss = extractdata(loss);

end

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

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

loss=MSEf+MSEu,

где MSEf=1Nfi=1Nf|f(xfi,tfi)|2 и MSEu=1Nui=1Nu|u(xui,tui)-ui|2.

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

Вычисление MSEf требует производных ut,ux,2ux2 из выхода 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))/2;

% 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

Ссылки

  1. Maziar Raissi, Париж Пердикэрис и Джордж Эм Карниадакис, Физика Информированное Глубокое обучение (Первая часть): управляемые данными Решения Нелинейных Дифференциальных уравнений с частными производными https://arxiv.org/abs/1711.10561

  2. К. Бэсдевэнт, М. Девилл, П. Холденванг, Ж. Лакруа, Дж. Оуэззэни, Р. Пеирет, П. Орланди, А. Патера, решения для Спектральной и конечной разности уравнения Burgers, Компьютеров & жидкостей 14 (1986) 23–41.

Смотрите также

| |

Похожие темы