Parameter identification with nlgreyest

I have a non-linear model I want to fit into an available to identify system parameters. I have used nlgreyest but the result is very different from the behaviour of the avaialable data.

My ODE model . I want A,B and D to be the free parameters.

I have tried tightening the tolerance but the no improvement.

What I am looking for is how to set upper and lower bounds for these parameters since I have an idea of what they could be.

Answer 1

[Alex Sha]

Hi, would you please attach your data.

Answer 2

[A Afonja]

clc
clear all
% i want to use this script to estimate damping
fprintf('ensure all PRT parameters are correct\n')

xt=0.1; %% ART width [m]
wr=0.17; %% reservoir length [m]
wd=0.514; %% duct length [m]
hd=0.17; %% Duct height [m]
w=wd+wr; %% reservoir distance [m]
hr= 0.235; %% Datum level [m]
ht=0.610; %% total height of tank [m]
rd= 0.058; %% COG distance from duct c/l []

% r1=0.065; %% inner fllet radius
% r2= 0.115; %% external fillet raduis 
rho=997.561; %% mass density of fluid in U ART [kg/m^3]
g=9.81; %% acceleration due to gravity[m/s^2]
% q=0.022; %% this is the coefficient of resistance derived from experiement or estimated
% WWeight=(((hr+hd/2)*wr)*2 + (wd*hd))*xt*rho; %% weight of water in tank
%% tank properties
Vt=wr*xt*(2*hr+w*wr/hd);
alfa3=g*rho*xt*wr*w;
alfa4=g*rho*xt*wr;
dn=0.02; %% this is an initial guess for nonlinear damping value
dl=0.02; %% this is an inital guess for linear damping value
%% Simulation result upload
TTRange=[1995,12.611455481821041]; %% in case you want to remove initial part part of the simulation use this to specify the index and time from where you want result to be used, if not use [1 0]
% TTRange=[1095,7.019779593687162];
%%% note the first item in TTRange is the index from where will begin and
%%% second item is time at that point
% TTRange=[1 0]; %if you dont want to cut result at anypart use this
Sim_R=readtable('leftdH_009.csv');
dH=Sim_R.leftdHMonitor_LeftdHMonitor(TTRange(1):end); %1742
Times=Sim_R.Time(TTRange(1):end)-TTRange(2);% this operation shift the simulation time to zero if you have chosen to remove the early part of result. 
Sim_Vel=readtable('leftdHVel_009.csv');
Vel=Sim_Vel.VelocityLevelLeftMonitor_VelocityLevelLeftMonitor_m_s_(TTRange(1):end);
% calculate acceleartion from velocity
% Accel=zeros(length(Vel),1);
Accel=diff(Vel)./(diff(Times));
Accel(end+1)=Accel(end);


%%% CFD result not in fixed time step, use cubline spline to get result at
%%% equally spaced time step
ddt=0.01;%% sampel interval for spline curve
Inter_time=(0:ddt:floor(Times(end)))';
Inter_dH=interp1(Times,dH,Inter_time,'spline');
Inter_Vel=interp1(Times,Vel,Inter_time,'spline');
Inter_Accel=diff(Inter_Vel)./(diff(Inter_time));
Inter_Accel(end+1)=Inter_Accel(end);

%% concatenate the vector of accel, velocity and product of velocity
%%%  a_tt*dh" + dl*dh' + dn*|dh'|*dh' + c_tt*dh=0
%%% lets keep only c_tt fixed
% A=[Accel*rho,Vel,Vel.*abs(Vel)];
%% Matrix MEthod 1 use to get some initial estimate of the parameters
%%
A=[Accel*rho,Vel,Vel.*abs(Vel)];
C=-2*alfa4*dH;
Coeef=pinv(A)*C;
VtFit=Coeef(1);
dlFit=Coeef(2);
dnFit=Coeef(3);
tspan=0:0.001:Times(end);
y0=[dH(1) Vel(1)];
% %% reconstruct using ODE45 for Matrix method 1
% [tfit,dHfit]=ode45(@(t,y) nonlinearODE(t,y,rho*VtFit,dlFit,dnFit,2*alfa4), tspan, y0);
% figNo =  findobj('type','figure');figure(length(figNo)+1)
% plot(tfit,dHfit(:,1)); hold on
% plot(Times,dH)
% legend('data','CFD')
% hold off


%% Set up programe for optimization

data = iddata(Inter_dH,[],ddt,'Name','CFD_Model3D_009');%% iddata requires that the Output of experiement be a column vector 
data.OutputName = 'Water Level';
data.OutputUnit = 'mm';
data.Tstart = 0;
data.TimeUnit = 's'; 
order         = [1 0 2];
initial_states = [dH(1); Vel(1)];
Ts=0;
parameters    = {rho,alfa4,VtFit,dlFit,dnFit};
% create the nonlinear model for nlgreyest
nonlinear_model = idnlgrey('Tank_non_linearLagrangian',order,parameters,initial_states,Ts);

% set free and fixed parameterts
setpar(nonlinear_model,'Fixed',{true true false false false});
%
setpar(nonlinear_model,'Minimum',{rho,alfa4,0.001, 0.001 dnFit});
setpar(nonlinear_model,'Maximum',{rho,alfa4,0.1, 5 0.1});

nonlinear_model = nlgreyest(data,nonlinear_model,'Display','Full');
% extract solution
Vt_est = nonlinear_model.Parameters(3).Value; %% rename the free parameter to what you want
dl_est = nonlinear_model.Parameters(4).Value; %% rename the free parameter to what you want
dn_est = nonlinear_model.Parameters(5).Value; %% rename the free parameter to what you want
% extract solution and standard deviation
[nonlinear_b_est, S_deviation] = getpvec(nonlinear_model,'free')
compare(data,nonlinear_model)

Answer 3

[A Afonja]

I was able to ge the system identification that matches the CFD response but I would like to set some comtraints on the maximum and mininum valus for the free parameters