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

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

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

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

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

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

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

numLayers = 9;
numNeurons = 20;

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

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

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

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]

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

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

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

Определите fmincon Целевая функция

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

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

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

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

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

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

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

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

The 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)$$ выполняет уравнение Бургера, граничные условия и интиальное условие. В частности, две величины способствуют минимизации потерь:

$$\text{loss}={\text{MSE}}_f+{\text{MSE}}_u$$,

где $${\text{MSE}}_f=\frac{1}{N_f}\sum _{i=1}^{N_f}{|f(x_f^i,t_f^i)|}^2$$ и $${\text{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}$$ являются точками внутри области.

Вычисление $${\text{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))/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

Ссылки