Building a continuous-time state-space model for the heat equation

Hello,

I'm struggeling a litle bit with a thermal model based on the two and thre dimensional heat equation.

My aim is to estimate the HCT with the aid of "experimental data" and a grey-box model.

The thermal BC's are quite simple, and illustrated below.

I generate dummy data via this matalb script:

% generate_dummy_tempdata_and_model_fixed.m
% Generate Dummy-data and prepare (2D/3D)
clc, close all, clear all
%% === Parameter ===
modellTyp = "3D";                 % "2D" or "3D"
nx = 4; ny = 4; nz = 4;           % Node count
dx = 0.02;dy = 0.03; dz = 0.04;   % Node distance [m]
times = 0:1:30;                   % Time range [s]
T0 = 300;                         % Initial temperature [K]
dt = times(2) - times(1);         % Time step [s]
%% === 1) Generate dummy-data and save as CSV ===
data = [];
switch modellTyp
    case "2D"
        for iy = 0:ny-1
            for ix = 0:nx-1
                x = ix * dx;
                y = iy * dy;
                Tcurve = T0 + 10 * sin(2*pi*times/60 + ix*0.3 + iy*0.5);
                data = [data; x, y, Tcurve];
            end
        end
        % Header ass string array
        header = ["x","y", "t" + string(times)];
        csvFile = 'example_tempdata2D.csv';
    case "3D"
        for iz = 0:nz-1
            for iy = 0:ny-1
                for ix = 0:nx-1
                    x = ix * dx;
                    y = iy * dy;
                    z = iz * dz;
                    Tcurve = T0 + 10 * sin(2*pi*times/60 + ix*0.2 + iy*0.3 + iz*0.4);
                    data = [data; x, y, z, Tcurve];
                end
            end
        end
        % Header as string array
        header = ["x","y","z", "t" + string(times)];
        csvFile = 'example_tempdata3D.csv';
end
% Bulid table and save
Ttbl_out = array2table(data, 'VariableNames', header);
writetable(Ttbl_out, csvFile);
disp(['Dummy-data saved: ', csvFile]);

The core of the simulation is the transformation of the 2d and 3d diffusion equation into an ode.

I do it with this script:

function [A,B,C,D] = fHeat2D3D(theta, Ts)
% fHeat2D3D Continuous-time grey-box state-space for 2D/3D heat diffusion
% Signature for idgrey: fHeat2D3D(theta, Ts)
% Inputs:
%   theta = [Lambda, rho, cp, hBot1, ..., hBotN]
%   Ts    = sample time (ignored in continuous-time model)
% Globals: nx, ny, nz, dx, dy, dz, modellTyp, h_top

global nx ny nz dx dy dz modellTyp h_top

% Extract physical parameters
Lambda = theta(1);
rho    = theta(2);
cp     = theta(3);

% Determine expected number of bottom convection coefficients
if strcmp(modellTyp,'2D')
    expectedNbot = nx;
else
    expectedNbot = nx * nz;
end
% Validate theta length
if length(theta) ~= 3 + expectedNbot
    error('Theta must have length 3+%d (got %d)', expectedNbot, length(theta));
end
% Extract bottom convection coefficients
hBot = theta(4:end);

% Number of states
if strcmp(modellTyp,'2D')
    N = nx * ny;
else
    N = nx * ny * nz;
end

% Build A matrix (2D/3D Laplacian)
A = zeros(N);
dists = [dx, dx, dy, dy];
for idx = 1:N
    tmp = idx - 1;
    ix  = mod(tmp, nx) + 1;
    iy  = floor(tmp / nx) + 1;
    if strcmp(modellTyp,'2D')
        iz = 1;
    else
        iz = floor(tmp / (nx * ny)) + 1;
    end
    diagSum = 0;
    nbr = [-1,0; 1,0; 0,-1; 0,1];
    for k = 1:4
        x2 = ix + nbr(k,1);
        y2 = iy + nbr(k,2);
        if x2>=1 && x2<=nx && y2>=1 && y2<=ny
            z2 = iz;
            idx2 = x2 + (y2-1)*nx + (z2-1)*nx*ny;
            K = Lambda / (rho * cp * dists(k)^2);
            A(idx, idx2) = K;
            diagSum = diagSum + K;
        end
    end
    A(idx, idx) = -diagSum;
end

% Build B matrix
% Identify side (Dirichlet) nodes
side = false(N,1);
for idx = 1:N
    tmp = idx - 1;
    x   = mod(tmp, nx) + 1;
    y   = floor(tmp / nx) + 1;
    if strcmp(modellTyp,'2D')
        side(idx) = (x == 1) || (x == nx);
    else
        z = floor(tmp / (nx * ny)) + 1;
        side(idx) = (x==1)||(x==nx)||(y==1)||(y==ny);
    end
end
Nside = sum(side);
B = zeros(N, Nside + 1 + expectedNbot);

% Side BC inputs
sidx = find(side);
for i = 1:Nside
    B(sidx(i), i) = 1;
end
% Top convection BC
if strcmp(modellTyp,'2D')
    topIdx = find(~side & floor((0:N-1)/nx)==0);
else
    topIdx = find(~side & floor((0:N-1)/(nx*ny))==nz-1);
end
for i = 1:numel(topIdx)
    B(topIdx(i), Nside+1) = h_top;
end
% Bottom convection BC
if strcmp(modellTyp,'2D')
    botIdx = find(~side & floor((0:N-1)/nx)==ny-1);
else
    botIdx = find(~side & floor((0:N-1)/(nx*ny))==0);
end
for i = 1:numel(botIdx)
    B(botIdx(i), Nside+1+i) = hBot(i);
end

% Outputs and feedthrough
C = eye(N);
D = zeros(N, size(B,2));
end

Finally the ode ("function [A,B,C,D] = fHeat2D3D(theta, Ts)") is used as continuous-time state-space model to estimate the system parameter HTC(bottom), heat conductivity, density and specific heat capacity via idgrey and greyest.

% greybox_heat_diffusion.m
% Parametric continuous-time grey-box model for 2D/3D heat conduction

clear; clc; close all;
global nx ny nz dx dy dz modellTyp h_top t_env_top

%% === User setup ===
modellTyp     = "2D";            % "2D" or "3D"
nx = 3; ny = 3; nz = 3;          % Count of nodes in x, y, z direction(ignore nz for 2D)
dx = 0.02; dy = 0.03; dz = 0.04; % Distance between nodes [m]

dt = 1;                           % Time step [s]
times = 0:dt:60;                  % Time range [s]
csv2D = 'example_tempdata2D.csv';
csv3D = 'example_tempdata3D.csv';

% Fixed values / Top-HTC
t_env_top = 293;                  % Ambient free stream temperature [K]
h_top = 10;                       % HTC at the top edge [W/(m^2*K)]

%% === Initial guess of the model parameter ===
theta0.Lambda = 50;  % Heat conduction [W/(m*K)]
theta0.rho    = 7800;% Density [kg/m^3]
theta0.cp     = 460; % Heat capacity [J/(kg*K)]
if modellTyp=="2D"
    h_bot0 = 20*ones(nx,1);
else
    h_bot0 = 20*ones(nx*nz,1);
end

%% === Read experimental data ===
if modellTyp=="2D"
    Tdata = readtable(csv2D);
    Texp = Tdata{:,3:end};
    Nstates = nx*ny;
else
    Tdata = readtable(csv3D);
    Texp = Tdata{:,4:end};
    Nstates = nx*ny*nz;
end
Nmeas = size(Texp,2);

%% === Calculate side indices ===
sideNodes = false(Nstates,1);
for idx=1:Nstates
    tmp=idx-1; ix=mod(tmp,nx)+1; iy=floor(tmp/nx)+1;
    if strcmp(modellTyp,"2D")
        sideNodes(idx) = ix==1 || ix==nx;
    else
        iz=floor(tmp/(nx*ny))+1;
        sideNodes(idx) = ix==1||ix==nx||iy==1||iy==ny;
    end
end
sideIdx = find(sideNodes);

%% === Build input-matrix u ===
T_sideBC = Texp(sideIdx,:);         % Side nodes - BC (experimental value)
T_topBC  = t_env_top * ones(1,Nmeas);
u = [T_sideBC' T_topBC' zeros(Nmeas,length(h_bot0))];

%% === idgrey Setup ===
param0 = [theta0.Lambda theta0.rho theta0.cp h_bot0'];
order  = [Nstates Nstates Nstates 0];

[A,B,C,D] = fHeat2D3D(param0, dt);


sysi = idgrey('fHeat2D3D', param0, 'c', order, ...
    'InputName',[repmat({'T_sideBC'},1,length(sideIdx)) {'T_top'} repmat({'T_bot'},1,length(h_bot0))], ...
    'OutputName', repmat({'Tnode'},1,Nstates));

% Release parameters (all)
sysi.Structure.Parameters.Free = true;

%% === Parameter-estimation ===
data_id = iddata(Texp', u, dt);
opt = greyestOptions('Display','on','InitialCondition','zero');
sys_est = greyest(data_id, sysi, opt);

%% === Extract estimated values ===
th = sys_est.Structure.Parameters;
Lambda_est = th(1).Value; rho_est = th(2).Value; cp_est = th(3).Value;
hBot_est = [th(4:end).Value]';

%% === Plots ===
Tsim = sim(sys_est, data_id);

% Plot h_bot
figure;
if strcmp(modellTyp,"2D")
    plot(unique(Tdata.x), hBot_est, 'o-', 'LineWidth', 2);
    xlabel('x [m]'); ylabel('h_{bot} [W/m^2K]');
    title('HTC at the bottom edge (2D)'); grid on;
else
    % 3D-Surface plot of h_bot at bottom surface
    [Xq, Zq] = meshgrid(unique(Tdata.x), unique(Tdata.z));
    hq = griddata(Tdata.x, Tdata.z, hBot_est, Xq, Zq, 'natural');

    surf(Xq, Zq, hq);
    xlabel('x [m]'); ylabel('z [m]'); zlabel('h_{bot} [W/m^2K]');
    title('HTC at bottom surface (3D)'); colorbar;
    shading interp; view(30, 30);
end

% Comparison temperature-simulation vs. temperature-measurement
figure;
Nrow = ceil(sqrt(Nstates));
for i = 1:Nstates
    subplot(Nrow, Nrow, i);
    plot(times, Texp(i,:), '-', 'LineWidth', 1.5); hold on;
    plot(times, Tsim.y(i,:), ':', 'LineWidth', 1.5);
    title(sprintf('Nodes %d', i));
    xlabel('Time [s]'); ylabel('T [°C]');
    legend('Experiment', 'Simulation');
    grid on;
end

% 3D-Point: final temperature distribution
if strcmp(modellTyp,"3D")
    figure;
    scatter3(Tdata.x, Tdata.y, Tdata.z, 50, Tsim.y(:,end), 'filled');
    xlabel('x'); ylabel('y'); zlabel('z');
    title('Simulated temperature distribution at t=tend');
    colorbar; view(30,30);
end

%% End of the script

The call of the script "greybox_heat_diffusion.m" is not working properly and the function "fHeat2D3D" is not called correct.

For the averaging of the corne nodes (2D) respectively the corner & edge nodes (3D). I made this script:

function [A,B,C,D] = fHeat2D3D(theta, Ts)
% fHeat2D3D Continuous-time grey-box state-space for 2D/3D heat diffusion
% Signature for idgrey: fHeat2D3D(theta, Ts)
% Inputs:
%   theta = [Lambda, rho, cp, hBot1, ..., hBotN]
%   Ts    = sample time (ignored in continuous-time model)
% Globals: nx, ny, nz, dx, dy, dz, modellTyp, h_top

global nx ny nz dx dy dz modellTyp h_top

% Extract physical parameters
Lambda = theta(1);
rho    = theta(2);
cp     = theta(3);

% Determine expected number of bottom convection coefficients
if strcmp(modellTyp,'2D')
    expectedNbot = nx;
else
    expectedNbot = nx * nz;
end
% Validate theta length
if length(theta) ~= 3 + expectedNbot
    error('Theta must have length 3+%d (got %d)', expectedNbot, length(theta));
end
% Extract bottom convection coefficients
hBot = theta(4:end);

% Number of states
if strcmp(modellTyp,'2D')
    N = nx * ny;
else
    N = nx * ny * nz;
end

% Build A matrix (2D/3D Laplacian)
A = zeros(N);
dists = [dx, dx, dy, dy];
for idx = 1:N
    tmp = idx - 1;
    ix  = mod(tmp, nx) + 1;
    iy  = floor(tmp / nx) + 1;
    if strcmp(modellTyp,'2D')
        iz = 1;
    else
        iz = floor(tmp / (nx * ny)) + 1;
    end
    diagSum = 0;
    nbr = [-1,0; 1,0; 0,-1; 0,1];
    for k = 1:4
        x2 = ix + nbr(k,1);
        y2 = iy + nbr(k,2);
        if x2>=1 && x2<=nx && y2>=1 && y2<=ny
            z2 = iz;
            idx2 = x2 + (y2-1)*nx + (z2-1)*nx*ny;
            K = Lambda / (rho * cp * dists(k)^2);
            A(idx, idx2) = K;
            diagSum = diagSum + K;
        end
    end
    A(idx, idx) = -diagSum;
end

% Build B matrix
% Identify side (Dirichlet) nodes
side = false(N,1);
for idx = 1:N
    tmp = idx - 1;
    x   = mod(tmp, nx) + 1;
    y   = floor(tmp / nx) + 1;
    if strcmp(modellTyp,'2D')
        side(idx) = (x == 1) || (x == nx);
    else
        z = floor(tmp / (nx * ny)) + 1;
        side(idx) = (x==1)||(x==nx)||(y==1)||(y==ny);
    end
end
Nside = sum(side);
B = zeros(N, Nside + 1 + expectedNbot);

% Side BC inputs
sidx = find(side);
for i = 1:Nside
    B(sidx(i), i) = 1;
end
% Top convection BC
if strcmp(modellTyp,'2D')
    topIdx = find(~side & floor((0:N-1)/nx)==0);
else
    topIdx = find(~side & floor((0:N-1)/(nx*ny))==nz-1);
end
for i = 1:numel(topIdx)
    B(topIdx(i), Nside+1) = h_top;
end
% Bottom convection BC
if strcmp(modellTyp,'2D')
    botIdx = find(~side & floor((0:N-1)/nx)==ny-1);
else
    botIdx = find(~side & floor((0:N-1)/(nx*ny))==0);
end
for i = 1:numel(botIdx)
    B(botIdx(i), Nside+1+i) = hBot(i);
end

% Outputs and feedthrough
C = eye(N);
% === Corner and edge nodes averaging ===
for idx = 1:N
    bcCount = nnz(B(idx,:));  % Count of BC-Input
    if bcCount > 1
        % Balanced averaging: scale A- and B-Line
        A(idx,:) = A(idx,:) / bcCount;
        B(idx,:) = B(idx,:) / bcCount;
    end
end

D = zeros(N, size(B,2));(N, size(B,2));
end

I'm not sure if this is the rigth approach.

Maybe someone has a hint how can I get rid of the errors.

With best regards

Steffen

Answer 1

[Vanita]

Hi Steffen B.

Based on your data generation code, it appears that 16 states are defined, but your model only uses 9 states. Below are the corrections to ensure compatibility with 9 states:

Correction1: in Top and Bottom Convection BC snippet replace variable side with side transpose as shown in the code below:

In fHeat2D3D(theta, Ts, varargin) :

% Top convection BC()
if strcmp(modellTyp,'2D')
    topIdx = find(~side' & floor((0:N-1)/nx)==0);
else
    topIdx = find(~side' & floor((0:N-1)/(nx*ny))==nz-1);
end
for i = 1:numel(topIdx)
    B(topIdx(i), Nside+1) = h_top;
end


% Bottom convection BC
if strcmp(modellTyp,'2D')
    botIdx = find(~side' & floor((0:N-1)/nx)==ny-1);
else
    botIdx = find(~side' & floor((0:N-1)/(nx*ny))==0);
end
for i = 1:numel(botIdx)
    B(botIdx(i), Nside+1+i) = hBot(i);
end

Correction 2: using only 9 states out of 16.

In main file:

%% === Read experimental data ===
if modellTyp=="2D"
    Tdata = readtable(csv2D);
    Texp = Tdata{1:9,3:end};
    Nstates = nx*ny;

Correction 3: add OutputNames.

OutputNames = arrayfun(@(k) sprintf('Tnode%d', k), 1:Nstates, 'UniformOutput', false);
sysi = idgrey('fHeat2D3D', param0, 'c', order, ...
    'InputName',[repmat({'T_sideBC'},1,length(sideIdx)) {'T_top'} repmat({'T_bot'},1,length(h_bot0))], ...
    'OutputName', OutputNames); %repmat({'Tnode'},1,Nstates));

Correction 4: replace InitialCondition with InitialState

%% === Parameter-estimation ===
data_id = iddata(Texp', u, dt);
opt = greyestOptions('Display','on','InitialState','zero');
sys_est = greyest(data_id, sysi, opt);

Correction 5: replace th(x).Value with th.Value(x);

%% === Extract estimated values ===
th = sys_est.Structure.Parameters;
Lambda_est = th.Value(1); rho_est = th.Value(2); cp_est = th.Value(3);
hBot_est = [th.Value((4:end))]';

Hope this helps!

[Torsten]

@Vanita

I tried to make the changes to the code, but didn't succeed in producing a running code.

Could you make the necessary changes or include the complete corrected code you used ?

% greybox_heat_diffusion.m
% Parametric continuous-time grey-box model for 2D/3D heat conduction
clear; clc; close all;
global nx ny nz dx dy dz modellTyp h_top t_env_top
%% === User setup ===
modellTyp     = "2D";            % "2D" or "3D"
nx = 3; ny = 3; nz = 3;          % Count of nodes in x, y, z direction(ignore nz for 2D)
dx = 0.02; dy = 0.03; dz = 0.04; % Distance between nodes [m]
dt = 1;                           % Time step [s]
times = 0:dt:60;                  % Time range [s]
csv2D = 'example_tempdata2D.csv';
csv3D = 'example_tempdata3D.csv';
% Fixed values / Top-HTC
t_env_top = 293;                  % Ambient free stream temperature [K]
h_top = 10;                       % HTC at the top edge [W/(m^2*K)]


% generate_dummy_tempdata_and_model_fixed.m
% Generate Dummy-data and prepare (2D/3D)

T0 = 300;                         % Initial temperature [K]
%% === 1) Generate dummy-data and save as CSV ===
data = [];
switch modellTyp
    case "2D"
        for iy = 0:ny-1
            for ix = 0:nx-1
                x = ix * dx;
                y = iy * dy;
                Tcurve = T0 + 10 * sin(2*pi*times/60 + ix*0.3 + iy*0.5);
                data = [data; x, y, Tcurve];
            end
        end
        % Header ass string array
        header = ["x","y", "t" + string(times)];
        csvFile = 'example_tempdata2D.csv';
    case "3D"
        for iz = 0:nz-1
            for iy = 0:ny-1
                for ix = 0:nx-1
                    x = ix * dx;
                    y = iy * dy;
                    z = iz * dz;
                    Tcurve = T0 + 10 * sin(2*pi*times/60 + ix*0.2 + iy*0.3 + iz*0.4);
                    data = [data; x, y, z, Tcurve];
                end
            end
        end
        % Header as string array
        header = ["x","y","z", "t" + string(times)];
        csvFile = 'example_tempdata3D.csv';
end
% Bulid table and save
Ttbl_out = array2table(data, 'VariableNames', header);
writetable(Ttbl_out, csvFile);
disp(['Dummy-data saved: ', csvFile]);

Dummy-data saved: example_tempdata2D.csv

%% === Initial guess of the model parameter ===
theta0.Lambda = 50;  % Heat conduction [W/(m*K)]
theta0.rho    = 7800;% Density [kg/m^3]
theta0.cp     = 460; % Heat capacity [J/(kg*K)]
if modellTyp=="2D"
    h_bot0 = 20*ones(nx,1);
else
    h_bot0 = 20*ones(nx*nz,1);
end
%% === Read experimental data ===
if modellTyp=="2D"
    Tdata = readtable(csv2D);
    Texp = Tdata{1:9,3:end};
    Nstates = nx*ny;
else
    Tdata = readtable(csv3D);
    Texp = Tdata{:,4:end};
    Nstates = nx*ny*nz;
end
Nmeas = size(Texp,2);
%% === Calculate side indices ===
sideNodes = false(Nstates,1);
for idx=1:Nstates
    tmp=idx-1; ix=mod(tmp,nx)+1; iy=floor(tmp/nx)+1;
    if strcmp(modellTyp,"2D")
        sideNodes(idx) = ix==1 || ix==nx;
    else
        iz=floor(tmp/(nx*ny))+1;
        sideNodes(idx) = ix==1||ix==nx||iy==1||iy==ny;
    end
end
sideIdx = find(sideNodes);
%% === Build input-matrix u ===
T_sideBC = Texp(sideIdx,:);         % Side nodes - BC (experimental value)
T_topBC  = t_env_top * ones(1,Nmeas);
u = [T_sideBC' T_topBC' zeros(Nmeas,length(h_bot0))];
%% === idgrey Setup ===
param0 = [theta0.Lambda theta0.rho theta0.cp h_bot0'];
order  = [Nstates Nstates Nstates 0];
[A,B,C,D] = fHeat2D3D(param0, dt);

A = 9×9

-0.0503 0.0348 0 0.0155 0 0 0 0 0 0.0348 -0.0852 0.0348 0 0.0155 0 0 0 0 0 0.0348 -0.0503 0 0 0.0155 0 0 0 0.0155 0 0 -0.0658 0.0348 0 0.0155 0 0 0 0.0155 0 0.0348 -0.1006 0.0348 0 0.0155 0 0 0 0.0155 0 0.0348 -0.0658 0 0 0.0155 0 0 0 0.0155 0 0 -0.0503 0.0348 0 0 0 0 0 0.0155 0 0.0348 -0.0852 0.0348 0 0 0 0 0 0.0155 0 0.0348 -0.0503

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

B = 9×10

1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 10 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 20 0 0 0 0 0 0 0 1 0 0 0 0

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

C = 9×9

1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

D = 9×10

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

OutputNames = arrayfun(@(k) sprintf('Tnode%d', k), 1:Nstates, 'UniformOutput', false);
sysi = idgrey('fHeat2D3D', param0, 'c', order, ...
    'InputName',[repmat({'T_sideBC'},1,length(sideIdx)) {'T_top'} repmat({'T_bot'},1,length(h_bot0))], ...
    'OutputName', OutputNames); %repmat({'Tnode'},1,Nstates));

Error using idgrey (line 375)
The ODE function "fHeat2D3D" could not be evaluated successfully using the given set of parameters, sample time and optional arguments. The error message generated during the evaluation was:
Undefined function 'fHeat2D3D' for input arguments of type 'double'.

% Release parameters (all)
sysi.Structure.Parameters.Free = true;
%% === Parameter-estimation ===
data_id = iddata(Texp', u, dt);
opt = greyestOptions('Display','on','InitialState','zero');
sys_est = greyest(data_id, sysi, opt);
%% === Extract estimated values ===
th = sys_est.Structure.Parameters;
Lambda_est = th.Value(1); rho_est = th.Value(2); cp_est = th.Value(3);
hBot_est = [th.Value((4:end))]';
%% === Plots ===
Tsim = sim(sys_est, data_id);
% Plot h_bot
figure;
if strcmp(modellTyp,"2D")
    plot(unique(Tdata.x), hBot_est, 'o-', 'LineWidth', 2);
    xlabel('x [m]'); ylabel('h_{bot} [W/m^2K]');
    title('HTC at the bottom edge (2D)'); grid on;
else
    % 3D-Surface plot of h_bot at bottom surface
    [Xq, Zq] = meshgrid(unique(Tdata.x), unique(Tdata.z));
    hq = griddata(Tdata.x, Tdata.z, hBot_est, Xq, Zq, 'natural');
    surf(Xq, Zq, hq);
    xlabel('x [m]'); ylabel('z [m]'); zlabel('h_{bot} [W/m^2K]');
    title('HTC at bottom surface (3D)'); colorbar;
    shading interp; view(30, 30);
end
% Comparison temperature-simulation vs. temperature-measurement
figure;
Nrow = ceil(sqrt(Nstates));
for i = 1:Nstates
    subplot(Nrow, Nrow, i);
    plot(times, Texp(i,:), '-', 'LineWidth', 1.5); hold on;
    plot(times, Tsim.y(i,:), ':', 'LineWidth', 1.5);
    title(sprintf('Nodes %d', i));
    xlabel('Time [s]'); ylabel('T [°C]');
    legend('Experiment', 'Simulation');
    grid on;
end
% 3D-Point: final temperature distribution
if strcmp(modellTyp,"3D")
    figure;
    scatter3(Tdata.x, Tdata.y, Tdata.z, 50, Tsim.y(:,end), 'filled');
    xlabel('x'); ylabel('y'); zlabel('z');
    title('Simulated temperature distribution at t=tend');
    colorbar; view(30,30);
end
%% End of the script

function [A,B,C,D] = fHeat2D3D(theta, Ts)
% fHeat2D3D Continuous-time grey-box state-space for 2D/3D heat diffusion
% Signature for idgrey: fHeat2D3D(theta, Ts)
% Inputs:
%   theta = [Lambda, rho, cp, hBot1, ..., hBotN]
%   Ts    = sample time (ignored in continuous-time model)
% Globals: nx, ny, nz, dx, dy, dz, modellTyp, h_top
global nx ny nz dx dy dz modellTyp h_top
% Extract physical parameters
Lambda = theta(1);
rho    = theta(2);
cp     = theta(3);
% Determine expected number of bottom convection coefficients
if strcmp(modellTyp,'2D')
    expectedNbot = nx;
else
    expectedNbot = nx * nz;
end
% Validate theta length
if length(theta) ~= 3 + expectedNbot
    error('Theta must have length 3+%d (got %d)', expectedNbot, length(theta));
end
% Extract bottom convection coefficients
hBot = theta(4:end);
% Number of states
if strcmp(modellTyp,'2D')
    N = nx * ny;
else
    N = nx * ny * nz;
end
% Build A matrix (2D/3D Laplacian)
A = zeros(N);
dists = [dx, dx, dy, dy];
for idx = 1:N
    tmp = idx - 1;
    ix  = mod(tmp, nx) + 1;
    iy  = floor(tmp / nx) + 1;
    if strcmp(modellTyp,'2D')
        iz = 1;
    else
        iz = floor(tmp / (nx * ny)) + 1;
    end
    diagSum = 0;
    nbr = [-1,0; 1,0; 0,-1; 0,1];
    for k = 1:4
        x2 = ix + nbr(k,1);
        y2 = iy + nbr(k,2);
        if x2>=1 && x2<=nx && y2>=1 && y2<=ny
            z2 = iz;
            idx2 = x2 + (y2-1)*nx + (z2-1)*nx*ny;
            K = Lambda / (rho * cp * dists(k)^2);
            A(idx, idx2) = K;
            diagSum = diagSum + K;
        end
    end
    A(idx, idx) = -diagSum;
end
% Build B matrix
% Identify side (Dirichlet) nodes
side = false(N,1);
for idx = 1:N
    tmp = idx - 1;
    x   = mod(tmp, nx) + 1;
    y   = floor(tmp / nx) + 1;
    if strcmp(modellTyp,'2D')
        side(idx) = (x == 1) || (x == nx);
    else
        z = floor(tmp / (nx * ny)) + 1;
        side(idx) = (x==1)||(x==nx)||(y==1)||(y==ny);
    end
end
Nside = sum(side);
B = zeros(N, Nside + 1 + expectedNbot);
% Side BC inputs
sidx = find(side);
for i = 1:Nside
    B(sidx(i), i) = 1;
end
% Top convection BC()
if strcmp(modellTyp,'2D')
    topIdx = find(~side' & floor((0:N-1)/nx)==0);
else
    topIdx = find(~side' & floor((0:N-1)/(nx*ny))==nz-1);
end
for i = 1:numel(topIdx)
    B(topIdx(i), Nside+1) = h_top;
end
% Bottom convection BC
if strcmp(modellTyp,'2D')
    botIdx = find(~side' & floor((0:N-1)/nx)==ny-1);
else
    botIdx = find(~side' & floor((0:N-1)/(nx*ny))==0);
end
for i = 1:numel(botIdx)
    B(botIdx(i), Nside+1+i) = hBot(i);
end
% Outputs and feedthrough
C = eye(N);
% === Corner and edge nodes averaging ===
for idx = 1:N
    bcCount = nnz(B(idx,:));  % Count of BC-Input
    if bcCount > 1
        % Balanced averaging: scale A- and B-Line
        A(idx,:) = A(idx,:) / bcCount;
        B(idx,:) = B(idx,:) / bcCount;
    end
end
D = zeros(N, size(B,2));
end

function [A,B,C,D] = fHeat2D3D_2(theta, Ts)
% fHeat2D3D Continuous-time grey-box state-space for 2D/3D heat diffusion
% Signature for idgrey: fHeat2D3D(theta, Ts)
% Inputs:
%   theta = [Lambda, rho, cp, hBot1, ..., hBotN]
%   Ts    = sample time (ignored in continuous-time model)
% Globals: nx, ny, nz, dx, dy, dz, modellTyp, h_top
global nx ny nz dx dy dz modellTyp h_top
% Extract physical parameters
Lambda = theta(1);
rho    = theta(2);
cp     = theta(3);
% Determine expected number of bottom convection coefficients
if strcmp(modellTyp,'2D')
    expectedNbot = nx;
else
    expectedNbot = nx * nz;
end
% Validate theta length
if length(theta) ~= 3 + expectedNbot
    error('Theta must have length 3+%d (got %d)', expectedNbot, length(theta));
end
% Extract bottom convection coefficients
hBot = theta(4:end);
% Number of states
if strcmp(modellTyp,'2D')
    N = nx * ny;
else
    N = nx * ny * nz;
end
% Build A matrix (2D/3D Laplacian)
A = zeros(N);
dists = [dx, dx, dy, dy];
for idx = 1:N
    tmp = idx - 1;
    ix  = mod(tmp, nx) + 1;
    iy  = floor(tmp / nx) + 1;
    if strcmp(modellTyp,'2D')
        iz = 1;
    else
        iz = floor(tmp / (nx * ny)) + 1;
    end
    diagSum = 0;
    nbr = [-1,0; 1,0; 0,-1; 0,1];
    for k = 1:4
        x2 = ix + nbr(k,1);
        y2 = iy + nbr(k,2);
        if x2>=1 && x2<=nx && y2>=1 && y2<=ny
            z2 = iz;
            idx2 = x2 + (y2-1)*nx + (z2-1)*nx*ny;
            K = Lambda / (rho * cp * dists(k)^2);
            A(idx, idx2) = K;
            diagSum = diagSum + K;
        end
    end
    A(idx, idx) = -diagSum;
end
% Build B matrix
% Identify side (Dirichlet) nodes
side = false(N,1);
for idx = 1:N
    tmp = idx - 1;
    x   = mod(tmp, nx) + 1;
    y   = floor(tmp / nx) + 1;
    if strcmp(modellTyp,'2D')
        side(idx) = (x == 1) || (x == nx);
    else
        z = floor(tmp / (nx * ny)) + 1;
        side(idx) = (x==1)||(x==nx)||(y==1)||(y==ny);
    end
end
Nside = sum(side);
B = zeros(N, Nside + 1 + expectedNbot);
% Side BC inputs
sidx = find(side);
for i = 1:Nside
    B(sidx(i), i) = 1;
end
% Top convection BC()
if strcmp(modellTyp,'2D')
    topIdx = find(~side' & floor((0:N-1)/nx)==0);
else
    topIdx = find(~side' & floor((0:N-1)/(nx*ny))==nz-1);
end
for i = 1:numel(topIdx)
    B(topIdx(i), Nside+1) = h_top;
end
% Bottom convection BC
if strcmp(modellTyp,'2D')
    botIdx = find(~side' & floor((0:N-1)/nx)==ny-1);
else
    botIdx = find(~side' & floor((0:N-1)/(nx*ny))==0);
end
for i = 1:numel(botIdx)
    B(botIdx(i), Nside+1+i) = hBot(i);
end
% Outputs and feedthrough
C = eye(N);
D = zeros(N, size(B,2));
end

[Vanita]

Hi Torsten

I am not encountering the error as you mentioned. Please try using the below code.

I changed the function signature to: function [A,B,C,D] = fHeat2D3D(theta, Ts,varargin).

% greybox_heat_diffusion.m
% Parametric continuous-time grey-box model for 2D/3D heat conduction
clear; clc; close all;
global nx ny nz dx dy dz modellTyp h_top t_env_top
%% === User setup ===
modellTyp     = "2D";            % "2D" or "3D"
nx = 3; ny = 3; nz = 3;          % Count of nodes in x, y, z direction(ignore nz for 2D)
dx = 0.02; dy = 0.03; dz = 0.04; % Distance between nodes [m]
dt = 1;                           % Time step [s]
times = 0:dt:60;                  % Time range [s]
csv2D = 'example_tempdata2D.csv';
csv3D = 'example_tempdata3D.csv';
% Fixed values / Top-HTC
t_env_top = 293;                  % Ambient free stream temperature [K]
h_top = 10;                       % HTC at the top edge [W/(m^2*K)]
% generate_dummy_tempdata_and_model_fixed.m
% Generate Dummy-data and prepare (2D/3D)
T0 = 300;                         % Initial temperature [K]
%% === 1) Generate dummy-data and save as CSV ===
data = [];
switch modellTyp
    case "2D"
        for iy = 0:ny-1
            for ix = 0:nx-1
                x = ix * dx;
                y = iy * dy;
                Tcurve = T0 + 10 * sin(2*pi*times/60 + ix*0.3 + iy*0.5);
                data = [data; x, y, Tcurve];
            end
        end
        % Header ass string array
        header = ["x","y", "t" + string(times)];
        csvFile = 'example_tempdata2D.csv';
    case "3D"
        for iz = 0:nz-1
            for iy = 0:ny-1
                for ix = 0:nx-1
                    x = ix * dx;
                    y = iy * dy;
                    z = iz * dz;
                    Tcurve = T0 + 10 * sin(2*pi*times/60 + ix*0.2 + iy*0.3 + iz*0.4);
                    data = [data; x, y, z, Tcurve];
                end
            end
        end
        % Header as string array
        header = ["x","y","z", "t" + string(times)];
        csvFile = 'example_tempdata3D.csv';
end
% Bulid table and save
Ttbl_out = array2table(data, 'VariableNames', header);
writetable(Ttbl_out, csvFile);
disp(['Dummy-data saved: ', csvFile]);

Dummy-data saved: example_tempdata2D.csv

%% === Initial guess of the model parameter ===
theta0.Lambda = 50;  % Heat conduction [W/(m*K)]
theta0.rho    = 7800;% Density [kg/m^3]
theta0.cp     = 460; % Heat capacity [J/(kg*K)]
if modellTyp=="2D"
    h_bot0 = 20*ones(nx,1);
else
    h_bot0 = 20*ones(nx*nz,1);
end
%% === Read experimental data ===
if modellTyp=="2D"
    Tdata = readtable(csv2D);
    Texp = Tdata{1:9,3:end};
    Nstates = nx*ny;
else
    Tdata = readtable(csv3D);
    Texp = Tdata{:,4:end};
    Nstates = nx*ny*nz;
end
Nmeas = size(Texp,2);
%% === Calculate side indices ===
sideNodes = false(Nstates,1);
for idx=1:Nstates
    tmp=idx-1; ix=mod(tmp,nx)+1; iy=floor(tmp/nx)+1;
    if strcmp(modellTyp,"2D")
        sideNodes(idx) = ix==1 || ix==nx;
    else
        iz=floor(tmp/(nx*ny))+1;
        sideNodes(idx) = ix==1||ix==nx||iy==1||iy==ny;
    end
end
sideIdx = find(sideNodes);
%% === Build input-matrix u ===
T_sideBC = Texp(sideIdx,:);         % Side nodes - BC (experimental value)
T_topBC  = t_env_top * ones(1,Nmeas);
u = [T_sideBC' T_topBC' zeros(Nmeas,length(h_bot0))];
%% === idgrey Setup ===
param0 = [theta0.Lambda theta0.rho theta0.cp h_bot0'];
order  = [Nstates Nstates Nstates 0];
[A,B,C,D] = fHeat2D3D(param0, dt);
% OutputNames = arrayfun(@(k) sprintf('Tnode%d', k), 1:Nstates, 'UniformOutput', false);
% sysi = idgrey('fHeat2D3D', param0, 'c', order, ...
%     'InputName',[repmat({'T_sideBC'},1,length(sideIdx)) {'T_top'} repmat({'T_bot'},1,length(h_bot0))], ...
%     'OutputName', OutputNames); %repmat({'Tnode'},1,Nstates));
%adding
OutputNames = arrayfun(@(k) sprintf('Tnode%d', k), 1:Nstates, 'UniformOutput', false);
sysi = idgrey('fHeat2D3D', param0, 'c', order, ...
    'InputName',[repmat({'T_sideBC'},1,length(sideIdx)) {'T_top'} repmat({'T_bot'},1,length(h_bot0))], ...
    'OutputName', OutputNames); %repmat({'Tnode'},1,Nstates));

Error using idgrey (line 375)
The ODE function "fHeat2D3D" could not be evaluated successfully using the given set of parameters, sample time and optional arguments. The error message generated during the evaluation was:
Undefined function 'fHeat2D3D' for input arguments of type 'double'.

% Release parameters (all)
sysi.Structure.Parameters.Free = true;
%% === Parameter-estimation ===
data_id = iddata(Texp', u, dt);
opt = greyestOptions('Display','on','InitialState','zero');
sys_est = greyest(data_id, sysi, opt);
%% === Extract estimated values ===
th = sys_est.Structure.Parameters;
Lambda_est = th.Value(1); rho_est = th.Value(2); cp_est = th.Value(3);
hBot_est = [th.Value((4:end))]';

function [A,B,C,D] = fHeat2D3D(theta, Ts,varargin)
% fHeat2D3D Continuous-time grey-box state-space for 2D/3D heat diffusion
% Signature for idgrey: fHeat2D3D(theta, Ts)
% Inputs:
%   theta = [Lambda, rho, cp, hBot1, ..., hBotN]
%   Ts    = sample time (ignored in continuous-time model)
% Globals: nx, ny, nz, dx, dy, dz, modellTyp, h_top
global nx ny nz dx dy dz modellTyp h_top
% Extract physical parameters
Lambda = theta(1);
rho    = theta(2);
cp     = theta(3);
% Determine expected number of bottom convection coefficients
if strcmp(modellTyp,'2D')
    expectedNbot = nx;
else
    expectedNbot = nx * nz;
end
% Validate theta length
if length(theta) ~= 3 + expectedNbot
    error('Theta must have length 3+%d (got %d)', expectedNbot, length(theta));
end
% Extract bottom convection coefficients
hBot = theta(4:end);
% Number of states
if strcmp(modellTyp,'2D')
    N = nx * ny;
else
    N = nx * ny * nz;
end
% Build A matrix (2D/3D Laplacian)
A = zeros(N);
dists = [dx, dx, dy, dy];
for idx = 1:N
    tmp = idx - 1;
    ix  = mod(tmp, nx) + 1;
    iy  = floor(tmp / nx) + 1;
    if strcmp(modellTyp,'2D')
        iz = 1;
    else
        iz = floor(tmp / (nx * ny)) + 1;
    end
    diagSum = 0;
    nbr = [-1,0; 1,0; 0,-1; 0,1];
    for k = 1:4
        x2 = ix + nbr(k,1);
        y2 = iy + nbr(k,2);
        if x2>=1 && x2<=nx && y2>=1 && y2<=ny
            z2 = iz;
            idx2 = x2 + (y2-1)*nx + (z2-1)*nx*ny;
            K = Lambda / (rho * cp * dists(k)^2);
            A(idx, idx2) = K;
            diagSum = diagSum + K;
        end
    end
    A(idx, idx) = -diagSum;
end
% Build B matrix
% Identify side (Dirichlet) nodes
side = false(N,1);
for idx = 1:N
    tmp = idx - 1;
    x   = mod(tmp, nx) + 1;
    y   = floor(tmp / nx) + 1;
    if strcmp(modellTyp,'2D')
        side(idx) = (x == 1) || (x == nx);
    else
        z = floor(tmp / (nx * ny)) + 1;
        side(idx) = (x==1)||(x==nx)||(y==1)||(y==ny);
    end
end
Nside = sum(side);
B = zeros(N, Nside + 1 + expectedNbot);
% Side BC inputs
sidx = find(side);
for i = 1:Nside
    B(sidx(i), i) = 1;
end
% Top convection BC()
if strcmp(modellTyp,'2D')
    topIdx = find(~side' & floor((0:N-1)/nx)==0);
else
    topIdx = find(~side' & floor((0:N-1)/(nx*ny))==nz-1);
end
for i = 1:numel(topIdx)
    B(topIdx(i), Nside+1) = h_top;
end
% Bottom convection BC
if strcmp(modellTyp,'2D')
    botIdx = find(~side' & floor((0:N-1)/nx)==ny-1);
else
    botIdx = find(~side' & floor((0:N-1)/(nx*ny))==0);
end
for i = 1:numel(botIdx)
    B(botIdx(i), Nside+1+i) = hBot(i);
end
% Outputs and feedthrough
C = eye(N);
% === Corner and edge nodes averaging ===
for idx = 1:N
    bcCount = nnz(B(idx,:));  % Count of BC-Input
    if bcCount > 1
        % Balanced averaging: scale A- and B-Line
        A(idx,:) = A(idx,:) / bcCount;
        B(idx,:) = B(idx,:) / bcCount;
    end
end
D = zeros(N, size(B,2));
end
function [A,B,C,D] = fHeat2D3D_2(theta, Ts)
% fHeat2D3D Continuous-time grey-box state-space for 2D/3D heat diffusion
% Signature for idgrey: fHeat2D3D(theta, Ts)
% Inputs:
%   theta = [Lambda, rho, cp, hBot1, ..., hBotN]
%   Ts    = sample time (ignored in continuous-time model)
% Globals: nx, ny, nz, dx, dy, dz, modellTyp, h_top
global nx ny nz dx dy dz modellTyp h_top
% Extract physical parameters
Lambda = theta(1);
rho    = theta(2);
cp     = theta(3);
% Determine expected number of bottom convection coefficients
if strcmp(modellTyp,'2D')
    expectedNbot = nx;
else
    expectedNbot = nx * nz;
end
% Validate theta length
if length(theta) ~= 3 + expectedNbot
    error('Theta must have length 3+%d (got %d)', expectedNbot, length(theta));
end
% Extract bottom convection coefficients
hBot = theta(4:end);
% Number of states
if strcmp(modellTyp,'2D')
    N = nx * ny;
else
    N = nx * ny * nz;
end
% Build A matrix (2D/3D Laplacian)
A = zeros(N);
dists = [dx, dx, dy, dy];
for idx = 1:N
    tmp = idx - 1;
    ix  = mod(tmp, nx) + 1;
    iy  = floor(tmp / nx) + 1;
    if strcmp(modellTyp,'2D')
        iz = 1;
    else
        iz = floor(tmp / (nx * ny)) + 1;
    end
    diagSum = 0;
    nbr = [-1,0; 1,0; 0,-1; 0,1];
    for k = 1:4
        x2 = ix + nbr(k,1);
        y2 = iy + nbr(k,2);
        if x2>=1 && x2<=nx && y2>=1 && y2<=ny
            z2 = iz;
            idx2 = x2 + (y2-1)*nx + (z2-1)*nx*ny;
            K = Lambda / (rho * cp * dists(k)^2);
            A(idx, idx2) = K;
            diagSum = diagSum + K;
        end
    end
    A(idx, idx) = -diagSum;
end
% Build B matrix
% Identify side (Dirichlet) nodes
side = false(N,1);
for idx = 1:N
    tmp = idx - 1;
    x   = mod(tmp, nx) + 1;
    y   = floor(tmp / nx) + 1;
    if strcmp(modellTyp,'2D')
        side(idx) = (x == 1) || (x == nx);
    else
        z = floor(tmp / (nx * ny)) + 1;
        side(idx) = (x==1)||(x==nx)||(y==1)||(y==ny);
    end
end
Nside = sum(side);
B = zeros(N, Nside + 1 + expectedNbot);
% Side BC inputs
sidx = find(side);
for i = 1:Nside
    B(sidx(i), i) = 1;
end
% Top convection BC()
if strcmp(modellTyp,'2D')
    topIdx = find(~side' & floor((0:N-1)/nx)==0);
else
    topIdx = find(~side' & floor((0:N-1)/(nx*ny))==nz-1);
end
for i = 1:numel(topIdx)
    B(topIdx(i), Nside+1) = h_top;
end
% Bottom convection BC
if strcmp(modellTyp,'2D')
    botIdx = find(~side' & floor((0:N-1)/nx)==ny-1);
else
    botIdx = find(~side' & floor((0:N-1)/(nx*ny))==0);
end
for i = 1:numel(botIdx)
    B(botIdx(i), Nside+1+i) = hBot(i);
end
% Outputs and feedthrough
C = eye(N);
D = zeros(N, size(B,2));
end

.

[Torsten]

Hello Vanita,

unfortunately the same problem as before (see above). With which MATLAB version did you make it work on your PC ?

I found that simple MATLAB examples don't work here under "Answers" - I don't know why.

par = [-1; 0.25; 0];
aux = {};
T = 0;
m = idgrey('myfunc',par,'c',aux,T);

Error using idgrey (line 375)
The ODE function "myfunc" could not be evaluated successfully using the given set of parameters, sample time and optional arguments. The error message generated during the evaluation was:
Undefined function 'myfunc' for input arguments of type 'double'.

load dcmotordata
data = iddata(y,u,0.1);
m_est = greyest(data,m);

function [A,B,C,D,K,x0] = myfunc(par,T)
  A = [0 1; 0 par(1)]; 
  B = [0;par(2)];
  C = eye(2);
  D = zeros(2,1);
  K = zeros(2,2);
  x0 = [par(3);0];
end

Example is taken from

https://uk.mathworks.com/help/ident/ug/estimating-linear-grey-box-models.html

[Vanita]

Hi Torsten,

It seems that the code runs without any issues when it is copied and pasted directly into MATLAB Desktop. The error might be due to a formatting issue introduced when opening the code in MATLAB online directly.

[Torsten]

@Vanita

Thank you for your help.