Data structure for estimation of continuous-Time Grey-Box model (continuous-time state-space model)

Hello together,

I am trying to understand the grey-box modelling in matlab with continuous-time state-space model structure using the example of 1d heat conduction in the rod (once with insulated end and convection on the other side and once with convection on both sides).

I want to estimate the parameters of the model with greyest, by applying it to a known temperature field that I define as a reference.

The script is shown below:

clc,close all,clear all

Tdata=xlsread('TempData.xlsx');
save TempData.mat Tdata;

time=Tdata(:,1);
T1=Tdata(:,3); % at x1=0;         P1
T2=Tdata(:,5); % at x2=0.1667;    P2
T3=Tdata(:,7); % at x3=2*0.1667;  P3
T4=Tdata(:,9); % at x3=3*0.1667;  P4
T5=Tdata(:,11); % at x3=4*0.1667; P4
T6=Tdata(:,13); % at x3=5*0.1667; P5

rho=7850;     % density
cp=420;       % heat capacity
lambda=45;    % conductivity
kappa=1;      % Initial guess for kappa=lambda/(rho*cp);

htf=1;        % Initial guess of htc at the right side
m_o_g_s=6;    % model order count = number of measuring points
r_l=1;        % rod length
Tini=22;      % Initial temperature

dL=r_l/m_o_g_s;

m=idgrey('heatd',{kappa,htf},'c',{m_o_g_s,r_l,Tini})
m2=idgrey('heatd2h',{kappa,htf},'c',{m_o_g_s,r_l,Tini})
me=greyest(Tdata,m)

The state-space system is defined below for insulation at the left side of the rod (x=0m)

function [A,B,C,D,K,x0] = heatd(kappa,htf,T,Ngrid,L,temp)
% ODE file parameterizing the heat diffusion model

% kappa (first parameter) - heat diffusion coefficient
% htf (second parameter) - heat transfer coefficient 
%                          at the far end of rod

% Auxiliary variables for computing state-space matrices:
% Ngrid: Number of points in the space-discretization
% L: Length of the rod
% temp: Initial room temperature (uniform)

% Compute space interval
deltaL = L/Ngrid;

% A matrix
A = zeros(Ngrid,Ngrid);
for kk = 2:Ngrid-1
  A(kk,kk-1) = 1;
  A(kk,kk) = -2;
  A(kk,kk+1) = 1;
end

% A matrix
A(1,1) = -1; A(1,2) = 1;
A(Ngrid,Ngrid-1) = 1;
A(Ngrid,Ngrid) = -1;
A = A*kappa/deltaL/deltaL;

% B matrix = Boundary condition on left end (insulated at x=0) and right end (convection at x=L)
B = zeros(Ngrid,1);
B(Ngrid,1) = htf/deltaL;

% C matrix
C = zeros(1,Ngrid);
C(1,1) = 1;

% D matrix (fixed to zero)
D = 0;

% K matrix: fixed to zero
K = zeros(Ngrid,1);

% Initial states: fixed to room temperature
x0 = temp*ones(Ngrid,1);

And slightly changed for free convection to air at the left side of the rod (x=0m):

function [A,B,C,D,K,x0] = heatd2h(kappa,htf,T,Ngrid,L,temp)
% ODE file parameterizing the heat diffusion model

% kappa (first parameter) - heat diffusion coefficient
% htf (second parameter) - heat transfer coefficient 
%                          at the far end of rod

% Auxiliary variables for computing state-space matrices:
% Ngrid: Number of points in the space-discretization
% L: Length of the rod
% temp: Initial room temperature (uniform)
hta = 20;
% Compute space interval
deltaL = L/Ngrid;

% A matrix
A = zeros(Ngrid,Ngrid);
for kk = 2:Ngrid-1
  A(kk,kk-1) = 1;
  A(kk,kk) = -2;
  A(kk,kk+1) = 1;
end

% A matrix 
A(1,1) = -1; A(1,2) = 1;
A(Ngrid,Ngrid-1) = 1;
A(Ngrid,Ngrid) = -1;
A = A*kappa/deltaL/deltaL;

% B matrix = Boundary condition on left end (free convection at x=0) and right end (convection at x=L)
B = zeros(Ngrid,1);
B(1,1) = hta/deltaL;
B(Ngrid,1) = htf/deltaL;

% C matrix 
C = zeros(1,Ngrid);
C(1,1) = 1;

% D matrix (fixed to zero)
D = 0;

% K matrix: fixed to zero
K = zeros(Ngrid,1);

% Initial states: fixed to room temperature
x0 = temp*ones(Ngrid,1);

The models seems okay and the output is:

m =
  Continuous-time linear grey box model defined by "heatd" function:
      dx/dt = A x(t) + B u(t) + K e(t)
       y(t) = C x(t) + D u(t) + e(t)
 
  A = 
        x1   x2   x3   x4   x5   x6
   x1  -36   36    0    0    0    0
   x2   36  -72   36    0    0    0
   x3    0   36  -72   36    0    0
   x4    0    0   36  -72   36    0
   x5    0    0    0   36  -72   36
   x6    0    0    0    0   36  -36
 
  B = 
       u1
   x1   0
   x2   0
   x3   0
   x4   0
   x5   0
   x6   6
 
  C = 
       x1  x2  x3  x4  x5  x6
   y1   1   0   0   0   0   0
 
  D = 
       u1
   y1   0
 
  K = 
       y1
   x1   0
   x2   0
   x3   0
   x4   0
   x5   0
   x6   0
 
  Model parameters:
   Par1 = 1
   Par2 = 1
 
Parameterization:
   ODE Function: heatd
   Disturbance component: parameterized by the ODE function
   Initial state: parameterized by the ODE function
   Number of free coefficients: 2
   Use "getpvec", "getcov" for parameters and their uncertainties.

Status:                                                         
Created by direct construction or transformation. Not estimated.

and for the second variant (with free convection at the left side):

m2 =
  Continuous-time linear grey box model defined by "heatd2h" function:
      dx/dt = A x(t) + B u(t) + K e(t)
       y(t) = C x(t) + D u(t) + e(t)
 
  A = 
        x1   x2   x3   x4   x5   x6
   x1  -36   36    0    0    0    0
   x2   36  -72   36    0    0    0
   x3    0   36  -72   36    0    0
   x4    0    0   36  -72   36    0
   x5    0    0    0   36  -72   36
   x6    0    0    0    0   36  -36
 
  B = 
        u1
   x1  120
   x2    0
   x3    0
   x4    0
   x5    0
   x6    6
 
  C = 
       x1  x2  x3  x4  x5  x6
   y1   1   0   0   0   0   0
 
  D = 
       u1
   y1   0
 
  K = 
       y1
   x1   0
   x2   0
   x3   0
   x4   0
   x5   0
   x6   0
 
  Model parameters:
   Par1 = 1
   Par2 = 1
 
Parameterization:
   ODE Function: heatd2h
   Disturbance component: parameterized by the ODE function
   Initial state: parameterized by the ODE function
   Number of free coefficients: 2
   Use "getpvec", "getcov" for parameters and their uncertainties.

Status:                                                         
Created by direct construction or transformation. Not estimated.

When i run the script i get this message:

Error using idgrey/greyest (line 46)
The model size is not consistent with the number of columns in the data matrix.

Error in untitled2 (line 48)
me=greyest(Tdata,m)
 

My plan was to estmate the two parameter (kappa and htf) with the data set.

Unfortunately the data structure is not right.

I don't know how the data input for estimation should look like.

Maybe someone has a clue how the structure should be.