물리정보 신경망을 사용하여 ODE 해 구하기

이 예제에서는 상미분 방정식(ODE)의 해를 예측하는 물리정보 신경망(PINN)을 훈련시키는 방법을 보여줍니다.

모든 미분 방정식이 닫힌 형식의 해를 갖는 것은 아닙니다. 이러한 방정식 유형의 근사해를 구하기 위해 다양한 전통적 수치 알고리즘을 사용할 수 있습니다. 물리정보 신경망(PINN)[1]은 신경망 자체의 구조와 훈련 과정에 물리 법칙을 통합시킨 신경망입니다. 물리 시스템을 정의하는 ODE의 해를 출력하는 신경망을 훈련시킬 수 있습니다. 이 접근법에는 몇 가지 이점이 있으며, 그중 하나는 닫힌 해석적 형태의 미분 가능한 근사해를 제공한다는 점입니다[2].

이 예제에서는 다음을 수행하는 방법을 보여줍니다.

$x \in [0, 2]$ 범위에서 훈련 데이터 생성.
$x$ 를 입력값으로 받고 $x$ 에서 계산된 ODE $\dot{y} = - 2 xy$ 의 근사해를 반환하는 신경망 정의.
사용자 지정 손실 함수를 사용하여 신경망 훈련.
신경망 예측값과 해석적 해 비교.

ODE 및 손실 함수

이 예제에서는 다음 ODE를 풉니다.

$\dot{y} = - 2 xy,$

여기서 초기 조건은 $y (0) = 1$ 입니다. 이 ODE의 해석적 해는 다음과 같습니다.

$y (x) = e^{- x^{2}} .$

ODE와 초기 조건을 충족하는 값에서 벗어난 편차에 벌점을 적용하는 사용자 지정 손실 함수를 정의합니다. 이 예제에서 손실 함수는 ODE 손실과 초기 조건 손실의 가중합입니다.

$L_{θ} (x) = ‖ {y_{}^{˙}}_{θ} + 2 x y_{θ} ‖^{2} + k ‖ y_{θ} (0) - 1 ‖^{2}$

$θ$ 는 신경망 파라미터이고, $k$ 는 상수 계수이며, $y_{θ}$ 는 신경망에서 예측된 해이고, $\dot{y_{θ}}$ 는 자동 미분을 사용하여 계산된 예측 해의 도함수입니다. 항 $‖ \dot{y_{θ}} + 2 {xy}_{θ} ‖^{2}$ 은 ODE 손실로, 예측된 해가 ODE 정의를 충족하는 값에서 얼마나 벗어나 있는지를 정량화합니다. 항 $‖ y_{θ} (0) - 1 ‖^{2}$ 는 초기 조건 손실로, 예측된 해가 초기 조건을 충족하는 데 얼마나 벗어나 있는지를 정량화합니다.

입력 데이터 생성 및 신경망 정의하기

$x \in [0, 2]$ 범위에서 10,000개의 훈련 데이터 점을 생성합니다.

x = linspace(0,2,10000)';

ODE를 푸는 신경망을 정의합니다. 입력값이 실수 $x \in R$ 이므로 입력 크기를 1로 지정합니다.

inputSize = 1;
layers = [
    featureInputLayer(inputSize,Normalization="none")
    fullyConnectedLayer(10)
    sigmoidLayer
    fullyConnectedLayer(1)
    sigmoidLayer];

계층 배열에서 dlnetwork 객체를 만듭니다.

net = dlnetwork(layers)

net = 
  dlnetwork with properties:

         Layers: [5×1 nnet.cnn.layer.Layer]
    Connections: [4×2 table]
     Learnables: [4×3 table]
          State: [0×3 table]
     InputNames: {'input'}
    OutputNames: {'layer_2'}
    Initialized: 1

  View summary with summary.

모델 손실 함수 정의하기

이 예제의 모델 손실 함수 섹션에 나와 있는, 신경망, 입력 데이터로 구성된 미니 배치, 초기 조건 손실과 관련된 계수를 입력값으로 받는 함수 modelLoss를 만듭니다. 이 함수는 손실과 신경망의 학습 가능한 파라미터에 대한 손실 기울기를 반환합니다.

훈련 옵션 지정하기

크기가 100인 미니 배치로 Epoch 15회만큼 훈련시킵니다.

numEpochs = 15;
miniBatchSize = 100;

SGDM 최적화에 대한 옵션을 지정합니다. 학습률을 0.5로, 학습률 감소 인자를 0.5로, 학습률 감소 주기를 5로, 모멘텀을 0.9로 지정합니다.

initialLearnRate = 0.5;
learnRateDropFactor = 0.5;
learnRateDropPeriod = 5;
momentum = 0.9;

손실 함수의 초기 조건 항 계수를 7로 지정합니다. 이 계수는 손실에 대한 초기 조건의 상대적 기여도를 지정합니다. 이 파라미터를 수정하여 훈련 수렴의 속도를 높일 수 있습니다.

icCoeff = 7;

모델 훈련시키기

훈련 중에 데이터의 미니 배치를 사용하려면 다음을 수행하십시오.

훈련 데이터에서 arrayDatastore 객체를 만듭니다.
arrayDatastore 객체를 입력값으로 받고, 미니 배치 크기를 지정하여 일부 미니 배치를 버리도록 설정하며, 훈련 데이터 형식 "BC"(배치, 채널)를 갖는 minibatchqueue 객체를 만듭니다.

기본적으로 minibatchqueue 객체는 기본 유형 single을 사용하여 데이터를 dlarray 객체로 변환합니다. 또한 사용 가능한 GPU가 있으면 각 출력값을 gpuArray로 변환합니다. GPU를 사용하려면 Parallel Computing Toolbox™와 지원되는 GPU 장치가 필요합니다. 지원되는 장치에 대한 자세한 내용은 GPU 연산 요구 사항 (Parallel Computing Toolbox) 항목을 참조하십시오.

ads = arrayDatastore(x,IterationDimension=1);
mbq = minibatchqueue(ads, ...
    MiniBatchSize=miniBatchSize, ...
    PartialMiniBatch="discard", ...
    MiniBatchFormat="BC");

SGDM 솔버에 대한 속도 파라미터를 초기화합니다.

velocity = [];

훈련 진행 상황 모니터의 총 반복 횟수를 계산합니다.

numObservationsTrain = numel(x);
numIterationsPerEpoch = floor(numObservationsTrain / miniBatchSize);
numIterations = numEpochs * numIterationsPerEpoch;

TrainingProgressMonitor 객체를 초기화합니다. monitor 객체를 생성할 때 타이머가 시작되므로 이 객체를 훈련 루프와 가깝게 생성해야 합니다.

monitor = trainingProgressMonitor( ...
    Metrics="LogLoss", ...
    Info=["Epoch" "LearnRate"], ...
    XLabel="Iteration");

사용자 지정 훈련 루프를 사용하여 신경망을 훈련시킵니다. 각 Epoch마다 데이터를 섞고 루프를 사용해 데이터의 미니 배치를 순회합니다. 각 미니 배치에 대해 다음을 수행합니다.

dlfeval 함수와 modelLoss 함수를 사용하여 모델 손실과 모델 기울기를 평가합니다.
sgdmupdate 함수를 사용하여 신경망 파라미터를 업데이트합니다.
훈련 과정을 표시합니다.
Stop 속성이 true이면 중지합니다. 중지 버튼을 클릭하면 TrainingProgressMonitor 객체의 Stop 속성값이 true로 변경됩니다.

learnRateDropPeriod Epoch마다 학습률에 learnRateDropFactor를 곱합니다.

epoch = 0;
iteration = 0;
learnRate = initialLearnRate;
start = tic;

% Loop over epochs.
while epoch < numEpochs  && ~monitor.Stop
    epoch = epoch + 1;

    % Shuffle data.
    mbq.shuffle

    % Loop over mini-batches.
    while hasdata(mbq) && ~monitor.Stop

        iteration = iteration + 1;

        % Read mini-batch of data.
        X = next(mbq);

        % Evaluate the model gradients and loss using dlfeval and the modelLoss function.
        [loss,gradients] = dlfeval(@modelLoss, net, X, icCoeff);

        % Update network parameters using the SGDM optimizer.
        [net,velocity] = sgdmupdate(net,gradients,velocity,learnRate,momentum);

        % Update the training progress monitor.
        recordMetrics(monitor,iteration,LogLoss=log(loss));
        updateInfo(monitor,Epoch=epoch,LearnRate=learnRate);
        monitor.Progress = 100 * iteration/numIterations;

    end
    
    % Reduce the learning rate.
    if mod(epoch,learnRateDropPeriod)==0
        learnRate = learnRate*learnRateDropFactor;
    end
end

모델 테스트하기

예측값과 해석적 해를 비교하여 신경망의 정확도를 테스트합니다.

신경망이 $x \in [0, 2]$ 훈련 범위를 벗어나 외삽 가능한지 보기 위해 범위 $x \in [0, 4]$ 에서 테스트 데이터를 생성합니다.

xTest = linspace(0,4,1000)';

minibatchpredict 함수를 사용하여 예측을 수행합니다. 기본적으로 minibatchpredict 함수는 GPU를 사용할 수 있으면 GPU를 사용합니다. GPU를 사용하려면 Parallel Computing Toolbox™ 라이선스와 지원되는 GPU 장치가 필요합니다. 지원되는 장치에 대한 자세한 내용은 GPU 연산 요구 사항 (Parallel Computing Toolbox) 항목을 참조하십시오. GPU를 사용할 수 없는 경우, 함수는 CPU를 사용합니다. 실행 환경을 지정하려면 ExecutionEnvironment 옵션을 사용하십시오.

yModel = minibatchpredict(net,xTest);

해석적 해를 구합니다.

yAnalytic = exp(-xTest.^2);

해석적 해와 모델 예측값을 로그 스케일로 플로팅하여 서로 비교합니다.

figure;
plot(xTest,yAnalytic,"-")
hold on
plot(xTest,yModel,"--")
legend("Analytic","Model")
xlabel("x")
ylabel("y (log scale)")
set(gca,YScale="log")

모델은 훈련 범위 $x \in [0, 2]$ 에서는 해석적 해를 정확히 근사하고, $x \in (2, 4]$ 범위에서는 상대적으로 낮은 정확도로 외삽합니다.

훈련 범위 $x \in [0, 2]$ 에서 평균제곱 상대 오차를 계산합니다.

yModelTrain = yModel(1:500);
yAnalyticTrain = yAnalytic(1:500);

errorTrain = mean(((yModelTrain-yAnalyticTrain)./yAnalyticTrain).^2)

errorTrain = single

0.0029

외삽된 범위 $x \in (2, 4]$ 에서 평균제곱 상대 오차를 계산합니다.

yModelExtra = yModel(501:1000);
yAnalyticExtra = yAnalytic(501:1000);

errorExtra = mean(((yModelExtra-yAnalyticExtra)./yAnalyticExtra).^2)

errorExtra = single

6.9634e+05

평균제곱 상대 오차가 훈련 범위에서보다 외삽된 범위에서 더 높게 나타납니다.

모델 손실 함수

modelLoss 함수는 dlnetwork 객체 net, 입력 데이터 X로 구성된 미니 배치 및 초기 조건 손실과 관련된 계수 icCoeff를 입력값으로 받습니다. 이 함수는 net의 학습 가능한 파라미터에 대한 손실 기울기와 해당 손실을 반환합니다. 손실은 ODE 손실과 초기 조건 손실의 가중합으로 정의됩니다. 이 손실을 평가하려면 2계 도함수가 필요합니다. 2계 자동 미분을 활성화하려면, 함수 dlgradient를 사용하고 EnableHigherDerivatives 이름-값 인수를 true로 설정하십시오.

function [loss,gradients] = modelLoss(net, X, icCoeff)
y = forward(net,X);

% Evaluate the gradient of y with respect to x. 
% Since another derivative will be taken, set EnableHigherDerivatives to true.
dy = dlgradient(sum(y,"all"),X,EnableHigherDerivatives=true);

% Define ODE loss.
eq = dy + 2*y.*X;

% Define initial condition loss.
ic = forward(net,dlarray(0,"CB")) - 1;

% Specify the loss as a weighted sum of the ODE loss and the initial condition loss.
loss = mean(eq.^2,"all") + icCoeff * ic.^2;

% Evaluate model gradients.
gradients = dlgradient(loss, net.Learnables);

end

참고 문헌

Maziar Raissi, Paris Perdikaris, and George Em Karniadakis, Physics Informed Deep Learning (Part I): Data-driven Solutions of Nonlinear Partial Differential Equations https://arxiv.org/abs/1711.10561
Lagaris, I. E., A. Likas, and D. I. Fotiadis. “Artificial Neural Networks for Solving Ordinary and Partial Differential Equations.” IEEE Transactions on Neural Networks 9, no. 5 (September 1998): 987–1000. https://doi.org/10.1109/72.712178.

참고 항목

dlarray | dlfeval | dlgradient