Parameter Identification with PINN (Physics Informed NN)
43 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
Hello,
I have a system of PDEs simulating a heat transfer system, and I have experimental data of the behaviour of the system. I am trying to use a PINN for the identification of 4 unkwnon parameters.
I am trying first to validate the model with known parameters, but I am unable to reach any reasonable parameter value. The learning rate is too slow, any recommendations?
Thanks for the help!
To run the code you will need to download the following .mat: https://cyclomed-my.sharepoint.com/:u:/g/personal/administracion_cyclomed_onmicrosoft_com/ETlg6E6520JDnBoIOv35Bq0B0o6ZCXvAvvqG3Gn1Ht3ssA?e=hKclKH
For context this is a simplified version of the PDE:
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1696651/image.png)
% Identify parameters of single blow equation
% 1) Fluid (u1) --> B*du1/dt + du1/dx = 1/Pe*d2u1/dx^2 + NTU*(u2 - u1) + NTU_W(u3 - u1)
% 2) Regenerator (u2) -->du2/dx = K_R*d2u2/dx^2 - NTU*(u2 - u1)
% 3) Wall (u3) --> du3/dx = R_TC*K_W*d2u3/dx^2 - R_TC*NTU_W*(u3 - u1) - NTU_W(u3 - u3_init)
% Load Experiment data
load("Test 11.mat");
%% Neural Network Definition
batchSize = 500;
dimension = 2; % X = (t,x)
% Set up neural net.
hiddenSize = 100;
net = [
featureInputLayer(dimension)
fullyConnectedLayer(hiddenSize)
tanhLayer
fullyConnectedLayer(hiddenSize)
tanhLayer
fullyConnectedLayer(hiddenSize)
tanhLayer
fullyConnectedLayer(3)];
% Create struct of parameters
parameters.net = dlnetwork(net);
% Unknown PDE parameters
parameters.Pe = dlarray(1000);
parameters.NTU = dlarray(200);
parameters.NTU_W = dlarray(1);
parameters.NTU_EXT = dlarray(50);
% Known PDE parameters
constants.B = B;
constants.K_R = K_R;
constants.K_W = K_W;
constants.R_TC = R_TC;
constants.R_Q = R_Q;
constants.T_W = u3(1,1);
%% Experiment data -- Real PDE data
% Define the size of the original matrix
[rows, cols] = size(u1);
% Generate the grid for the original matrix
[X, Y] = meshgrid(1:cols, 1:rows);
% Generate the grid for the interpolated matrix
[XI, YI] = meshgrid(linspace(1, cols, batchSize), linspace(1, rows, batchSize));
% Interpolate the matrix along each dimension
UTest1 = (interp2(X, Y, u1, XI, YI));
UTest2 = (interp2(X, Y, u2, XI, YI));
UTest3 = (interp2(X, Y, u3, XI, YI));
% Select random points
[rows, columns] = size(UTest1);
X = ceil(rand(batchSize,1)*columns);
Y = ceil(rand(batchSize,1)*rows);
for k = 1 : length(X)
row = Y(k);
col = X(k);
UTest.u1(k) = UTest1(row, col);
UTest.u2(k) = UTest2(row, col);
UTest.u3(k) = UTest3(row, col);
end
% Transfor to dlarray
x = dlarray((X./batchSize)', "CB");
t = dlarray((Y.*t(end)./batchSize)', "CB");
X = [t;x];
% Experiment data
UTest.u1 = dlarray(UTest.u1, "CB");
UTest.u2 = dlarray(UTest.u2, "CB");
UTest.u3 = dlarray(UTest.u3, "CB");
%% Solve inverse PINN
% Training loop
avgG = [];
avgSqG = [];
maxIters = 10000000;
lossFcn = dlaccelerate(@modelLoss);
for iter = 1:maxIters
[loss,gradients] = dlfeval(lossFcn,X,parameters,constants,UTest);
[parameters,avgG,avgSqG] = adamupdate(parameters,gradients,avgG,avgSqG,iter);
if mod(iter, 1000) == 0
fprintf("Iteration : %d, Loss : %.4f \n",iter, extractdata(loss));
fprintf("Iteration: %d, Predicted Pe: %.3f, Actual Pe: %.3f, Predicted NTU: %.3f, Actual NTU: %.3f\n", ...
iter, extractdata(parameters.Pe), Pe, extractdata(parameters.NTU), NTU);
fprintf("Iteration: %d, Predicted NTU_W: %.3f, Actual NTU_W: %.3f, Predicted NTU_EXT: %.3f, Actual NTU_EXT: %.3f\n", ...
iter, extractdata(parameters.NTU_W), NTU_W, extractdata(parameters.NTU_EXT), NTU_Q);
end
end
%% PINN definition
function [loss,gradients] = modelLoss(X,parameters,constants,UTest)
% X = (t,x).
% PDE
u = predict(parameters.net, X);
u1 = u(1,:);
u2 = u(2,:);
u3 = u(3,:);
% First partial derivatives
du1 = dlgradient(sum(u1,2), X, RetainData=true, EnableHigherDerivatives=true);
du2 = dlgradient(sum(u2,2), X, RetainData=true, EnableHigherDerivatives=true);
du3 = dlgradient(sum(u3,2), X, RetainData=true, EnableHigherDerivatives=true);
du1dt = du1(1,:);
du1dx = du1(2,:);
du2dt = du2(1,:);
du2dx = du2(2,:);
du3dt = du3(1,:);
du3dx = du3(2,:);
% Second partial derivatives
d2u1 = dlgradient(sum(du1dx,2),X,RetainData=true);
d2u2 = dlgradient(sum(du2dx,2),X,RetainData=true);
d2u3 = dlgradient(sum(du3dx,2),X,RetainData=true);
d2u1dx2 = d2u1(2,:);
d2u2dx2 = d2u2(2,:);
d2u3dx2 = d2u3(2,:);
% PDE residual
odeResidual = (constants.B*du1dt + du1dx - 1/parameters.Pe*d2u1dx2 - parameters.NTU*(u2 - u1) - parameters.NTU_W*(u3 - u1)).^2;
odeResidual = odeResidual + (du2dt - constants.K_R*d2u2dx2 + parameters.NTU*(u2 - u1)).^2;
odeResidual = odeResidual + (du3dt - constants.K_W*d2u3dx2 + constants.R_TC*parameters.NTU_W*(u3 - u1) + constants.R_Q*parameters.NTU_EXT*(u3 - constants.T_W)).^2;
% Compute the mean square error of the ODE residual
pdeLoss = mean(odeResidual,"all");
% Compute the L2 difference between the predicted xpred and the true x.
reconstructionLoss = l2loss(u1,UTest.u1) + l2loss(u2,UTest.u2) + l2loss(u3,UTest.u3);
loss = pdeLoss + reconstructionLoss;
[gradients.net, gradients.Pe, gradients.NTU, gradients.NTU_W, gradients.NTU_EXT] = dlgradient(pdeLoss,parameters.net.Learnables,parameters.Pe, parameters.NTU, parameters.NTU_W, parameters.NTU_EXT);
end
1 Kommentar
José Eduardo
am 5 Jun. 2024
Hello Carlos,
I'm not an expert in matlab, but supposing everything is working as it should I would recommend you three things:
1) Be more selective of your starting guesses in your 4 parameters, by more selective I mean, choose parameters that are at maximum learning_rate*epochs/4 of distance from the real parameters.
2) Use a regularization lambda for each of yout losses so that they be in the same power of 10 when added up, and try to let the data loss a bit more than the other losses.
3) Check the maximum frequency of your output solution, if it is too high, the PINN method will have a huge difficulty to converge. From my expierence, for higher frequencies you need to lower the learning rates and increase the epochs as well as hope it can fastly converge the output function in the beginning so that the parameters (scale factors of the PDE) are the only unknonws.
I hope I can help, good luck
José
Antworten (1)
yin
am 19 Jul. 2024
Have you solved this problem? Can you share how you achieved it?
1 Kommentar
José Eduardo Alves
am 19 Jul. 2024
I have a solution for another problem if you want you can check my github, but it is in python for a simple case study of a forced mass-spring system.
Siehe auch
Kategorien
Mehr zu General PDEs finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!