Синтезируйте данные целевой динамики

Задайте целевую динамику как линейную модель ODE $$x^{\prime }=Ax$$, с x0 как его начальное условие, и вычисляют его числовое решение xTrain с ode45 во временном интервале [0 15]. Чтобы вычислить точные достоверные данные, установите относительную погрешность ode45 числовой решатель к ${10}^{-7}$. Позже, вы используете значение xTrain как достоверные данные для изучения аппроксимированной динамики с нейронной моделью ODE.

x0 = [2; 0];
A = [-0.1 -1; 1 -0.1];
trueModel = @(t,y) A*y;

numTimeSteps = 2000;
T = 15;
odeOptions = odeset(RelTol=1.e-7);
t = linspace(0, T, numTimeSteps);
[~, xTrain] = ode45(trueModel, t, x0, odeOptions);
xTrain = xTrain';

Визуализируйте обучающие данные в графике.

plot(xTrain(1,:),xTrain(2,:))
legend("Ground truth dynamics") 
xlabel("x(1)") 
ylabel("x(2)")
grid on

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

Функция модели состоит из одного вызова dlode45 решать ОДУ, заданный аппроксимированной динамикой $$f(t,y,\theta )$$ для 40 временных шагов.

neuralOdeTimesteps = 40;
dt = t(2);
timesteps = (0:neuralOdeTimesteps)*dt;

Задайте настраиваемые параметры, чтобы использовать в вызове dlode45 и соберите их в переменной neuralOdeParameters. Функциональный initializeGlorot берет в качестве входа размер настраиваемых параметров sz и количество выходных параметров и количество входных параметров полностью связанных операций, и возвращают dlarray объект с базовым типом 'single' с использованием набора значений инициализация Glorot. Функциональный initializeZeros берет в качестве входа размер настраиваемых параметров и возвращает параметры как dlarray объект с базовым типом 'single'. Функции, взятые в качестве примера, инициализации присоединены к этому примеру как к вспомогательным файлам. Чтобы получить доступ к этим функциям, откройте этот пример как live скрипт. Для получения дополнительной информации об инициализации настраиваемых параметров для функций модели, смотрите, Инициализируют Настраиваемые параметры для Функции Модели.

Инициализируйте структуру параметров.

neuralOdeParameters = struct;

Инициализируйте параметры для полностью связанных операций в модели ODE. Первая полностью связанная операция берет в качестве входа вектор из размера stateSize и увеличивает его длину до hiddenSize. С другой стороны вторая полностью связанная операция берет в качестве входа вектор из длины hiddenSize и уменьшает его длину к stateSize.

stateSize  = size(xTrain,1);
hiddenSize = 20;

neuralOdeParameters.fc1 = struct;
sz = [hiddenSize stateSize];
neuralOdeParameters.fc1.Weights = initializeGlorot(sz, hiddenSize, stateSize);
neuralOdeParameters.fc1.Bias    = initializeZeros([hiddenSize 1]);

neuralOdeParameters.fc2 = struct;
sz = [stateSize hiddenSize];
neuralOdeParameters.fc2.Weights = initializeGlorot(sz, stateSize, hiddenSize);
neuralOdeParameters.fc2.Bias    = initializeZeros([stateSize 1]);

Отобразите настраиваемые параметры модели.

neuralOdeParameters.fc1

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

neuralOdeParameters.fc2

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

Задайте нейронную модель ОДУ

Создайте функциональный odeModel, перечисленный в разделе ODE Model примера, который берет в качестве входа, который время ввело (неиспользованный), соответствующее решение и параметры функции ОДУ. Функция применяет полностью связанную операцию, tanh операцию и другую полностью связанную операцию к входным данным с помощью весов и смещает данный параметрами.

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

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

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

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

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

Задайте опции для оптимизации Адама.

gradDecay = 0.9;
sqGradDecay = 0.999;
learnRate = 0.002;

Обучайтесь для 1 200 итераций с "мини-пакетным размером" 200.

numIter = 1200;
miniBatchSize = 200;

Каждые 50 итераций, решите изученную динамику и отобразите их против основной истины в схеме фазы, чтобы показать учебный путь.

plotFrequency = 50;

Обучите модель Используя пользовательский учебный цикл

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

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

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

figure(2)
lineLossTrain = animatedline(Color=[0.85 0.325 0.098]);
ylim([0 inf])
xlabel("Iteration")
ylabel("Loss")
grid on

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

averageGrad   = [];
averageSqGrad = [];

Инициализируйте lossHistory массив, чтобы записать эволюцию учебной потери.

lossHistory = [];
numTrainingTimesteps = numTimeSteps;
trainingTimesteps = 1:numTrainingTimesteps;
plottingTimesteps = 2:numTimeSteps;

start = tic;

for iter=1:numIter
    
    % Create batch 
    [dlx0, targets] = createMiniBatch(numTrainingTimesteps, neuralOdeTimesteps, miniBatchSize, xTrain);

    % Evaluate network and compute gradients
    [grads,loss] = dlfeval(@modelGradients,timesteps,dlx0,neuralOdeParameters,targets);
    
    % Update network 
    [neuralOdeParameters,averageGrad,averageSqGrad] = adamupdate(neuralOdeParameters,grads,averageGrad,averageSqGrad,iter,...
        learnRate,gradDecay,sqGradDecay);
    
    % Plot loss
    currentLoss = double(extractdata(loss));
    
    figure(2)
    addpoints(lineLossTrain, iter, currentLoss);
    drawnow
    
    % Plot predicted vs. real dynamics
    if mod(iter,plotFrequency) == 0  || iter == 1
        figure(3)
        clf

        % Use ode45 to compute the solution 
        y = dlode45(@odeModel,t,dlarray(x0),neuralOdeParameters,DataFormat="CB");
        y = extractdata(y);
        
        plot(xTrain(1,plottingTimesteps),xTrain(2,plottingTimesteps),"r--")
        hold on
        plot(y(1,:),y(2,:),"b-")
        xlabel("x(1)")
        ylabel("x(2)")
        hold off
        D = duration(0,0,toc(start),Format="hh:mm:ss");
        title("Iter = " + iter + ", loss = " + num2str(currentLoss) + ", Elapsed: " + string(D))
        legend("Training ground truth", "Predicted")
    end
end

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

Используйте модель, чтобы вычислить аппроксимированные решения с различными начальными условиями.

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

tPred = t;

x0Pred1 = sqrt([2;2]);
x0Pred2 = [-1;-1.5];
x0Pred3 = [0;2];
x0Pred4 = [-2;0];

Численно решите ОДУ истинная динамика с ode45 для четырех новых начальных условий.

[~, xTrue1] = ode45(trueModel, tPred, x0Pred1, odeOptions);

[~, xTrue2] = ode45(trueModel, tPred, x0Pred2, odeOptions);

[~, xTrue3] = ode45(trueModel, tPred, x0Pred3, odeOptions);

[~, xTrue4] = ode45(trueModel, tPred, x0Pred4, odeOptions);

Численно решите ОДУ с изученной нейронной динамикой ОДУ.

xPred1 = dlode45(@odeModel,tPred,dlarray(x0Pred1),neuralOdeParameters,DataFormat="CB");
xPred1 = extractdata(squeeze(xPred1))';

xPred2 = dlode45(@odeModel,tPred,dlarray(x0Pred2),neuralOdeParameters,DataFormat="CB");
xPred2 = extractdata(squeeze(xPred2))';

xPred3 = dlode45(@odeModel,tPred,dlarray(x0Pred3),neuralOdeParameters,DataFormat="CB");
xPred3 = extractdata(squeeze(xPred3))';

xPred4 = dlode45(@odeModel,tPred,dlarray(x0Pred4),neuralOdeParameters,DataFormat="CB");
xPred4 = extractdata(squeeze(xPred4))';

Визуализируйте предсказания

Визуализируйте предсказанные решения для различных начальных условий против решений для основной истины с функциональным plotTrueAndPredictedSolutions, перечисленный в разделе Plot True и Predicted Solutions Function примера.

subplot(2,2,1)
plotTrueAndPredictedSolutions(xTrue1, xPred1);

subplot(2,2,2)
plotTrueAndPredictedSolutions(xTrue2, xPred2);

subplot(2,2,3)
plotTrueAndPredictedSolutions(xTrue3, xPred3);

subplot(2,2,4)
plotTrueAndPredictedSolutions(xTrue4, xPred4);

Функции помощника

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

model функция, которая задает нейронную сеть, используемую, чтобы сделать предсказания, состоит из одного нейронного вызова ОДУ. Для каждого наблюдения эта функция берет вектор из длины stateSize, который используется в качестве начального условия для решения численно ОДУ с функциональным odeModel, который представляет learnable правую сторону $$f(t,y,\theta )$$ из ОДУ, которое будет решено как правая сторона и вектор из моментов времени tspan определение времени, в которое выводится числовое решение. Функция использует векторный tspan для каждого наблюдения, независимо от начального условия, поскольку изученная система автономна. Таким образом, odeModel функция явным образом не зависит вовремя.

function X = model(tspan,X0,neuralOdeParameters)

X = dlode45(@odeModel,tspan,X0,neuralOdeParameters,DataFormat="CB");

end

Модель ОДУ

odeModel функция является learnable правой стороной, используемой в вызове dlode45. Это берет в качестве входа вектор из размера stateSize, увеличивает его так, чтобы это имело длину hiddenSize, и применяет функцию нелинейности tanh. Затем функция сжимает вектор снова, чтобы иметь длину stateSize.

function y = odeModel(~,y,theta)

y = tanh(theta.fc1.Weights*y + theta.fc1.Bias);
y = theta.fc2.Weights*y + theta.fc2.Bias;

end

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

Эта функция берет в качестве входных параметров векторный tspan, набор начальных условий dlX0, настраиваемые параметры neuralOdeParameters, и целевые последовательности targets. Это вычисляет предсказания с model функция, и сравнивает их с данными целевыми последовательностями. Наконец, это вычисляет градиент относительно настраиваемых параметров нейронного ОДУ.

function [gradients,loss] = modelGradients(tspan,dlX0,neuralOdeParameters,targets)

% Compute predictions.
dlX = model(tspan,dlX0,neuralOdeParameters);

% Compute L1 loss.
loss = l1loss(dlX,targets,NormalizationFactor="all-elements",DataFormat="CBT");

% Compute gradients.
gradients = dlgradient(loss,neuralOdeParameters);
end

Создайте функцию мини-пакетов

createMiniBatch функция создает пакет наблюдений за целевой динамикой. Это берет в качестве входа общее количество временных шагов достоверных данных numTimesteps, количество последовательных временных шагов, которые будут возвращены для каждого наблюдения numTimesPerObs, количество наблюдений miniBatchSize, и достоверные данные X.

function [x0, targets] = createMiniBatch(numTimesteps,numTimesPerObs,miniBatchSize,X)

% Create batches of trajectories.
s = randperm(numTimesteps - numTimesPerObs, miniBatchSize);

x0 = dlarray(X(:, s));
targets = zeros([size(X,1) miniBatchSize numTimesPerObs]);

for i = 1:miniBatchSize
    targets(:, i, 1:numTimesPerObs) = X(:, s(i) + 1:(s(i) + numTimesPerObs));
end

end

Постройте истинную и предсказанную функцию решений

plotTrueAndPredictedSolutions функционируйте берет в качестве входа истинное решение xTrue, аппроксимированное решение xPred вычисленный с изученной нейронной моделью ODE и соответствующим начальным условием x0Str. Это вычисляет ошибку между истинными и предсказанными решениями и строит его в схеме фазы.

function plotTrueAndPredictedSolutions(xTrue,xPred)

err = mean(abs(xTrue(2:end,:) - xPred), "all");

plot(xTrue(:,1),xTrue(:,2),"r--",xPred(:,1),xPred(:,2),"b-",LineWidth=1)

title("Absolute error = " + num2str(err,'%.4f') )
xlabel("x(1)")
ylabel("x(2)")

xlim([-2 3])
ylim([-2 3])

legend("Ground truth","Predicted",Location="best")
end

[1] Чен, Рики Т. К., Юлия Рубанова, Джесси Бетанкур и Давид Дювено. “Нейронные Обыкновенные дифференциальные уравнения”. Предварительно распечатайте, представленный 13 декабря 2019. https://arxiv.org/abs/1806.07366.

[2] Шемпин, Лоуренс Ф. и Марк В. Рейчелт. “Пакет ODE MATLAB”. SIAM Journal на Научных вычислениях 18, № 1 (январь 1997): 1–22. https://doi.org/10.1137/S1064827594276424.