Parameter Identification with PINN (Physics Informed NN)

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:

% 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

[José Eduardo]

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é

Answer 1

[yin]

Have you solved this problem? Can you share how you achieved it?

[José Eduardo Alves]

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.

Check here. You can change the branch for other additions in the problem.