Worse result when using "lsqnonlin"

11 views (last 30 days)
I have to estimate some parameters in my work. I have used "lsqnonlin" to estimate the parameters. Unfortunately, the answer of "lsqnonlin" is worse than the first simulation's answer (before parameter estimation). Actually, the estimated parameters aren't good. Would you please let me know why I face this problem?
========================================== Main Program
%Global Parameters
% Lower bound (Reaction rate constant should not be negative)
lb = 0.0001*ones(5,1);
% Initial guess
HH0 = [8 5 17 12 2]';
options = optimoptions('lsqnonlin', 'Algorithm','trust-region-reflective', 'TolFun',1e-03, 'TolX',1e-07);
options = optimoptions(options, 'Display','iter', 'Diagnostics','on', 'DiffMaxChange',1e+15, 'DiffMinChange',1e-03, 'MaxFunEvals',1000);
% Parameter estimation
[k,resnorm,res,exitflag,output] = lsqnonlin(@objfun,HH0,lb,[],options);
======================================== Objective Function
function final_f = objfun(HH)
% Global Parameters
% Molecular weight, g/mol (so far no overall mass balance implemented)
mC = 12.011;
mH = 1.0079;
mO = 15.999;
mCo = 58.933;
% Import data from excel
ndata = xlsread('dimension.xls');
% Mean velocity of RH, m/s
URH = MF/(rhoRH*A);
% Initial concentrations in the liquid phase, mol/m³
c0 = a column vector with size 14*1
% Interval of integration vector
L = 100; % Number of grid points
zspan = linspace(0,LR,L);
% Options for the ODE solver
options = odeset('RelTol',1.0e-10, 'AbsTol',1.0e-12, 'InitialStep',0.001);
% Integration
[z,f] = ode23tb(@(z,g)PFR_DGL(g,n,URH,HH),zspan,f0,options);
L_E = ndata(:,1);
t_E = L_E./(URH*1000);
% Matrix of experimental data
final_E = [O2_E, ROOH_E, ROH_E, RO_E];
% Standard deviation of experimental data
std_dev = std(final_E,1,1);
% Finding the simulation points which are correspond to measured point
[hh,aa] = size(ndata);
new_f = zeros(hh,n+3);
for i = 1:hh
for ii = 1:L
if t(ii,:) - t_E(i,:) <= 0.2
new_f(i,:) = f(ii,:);
end
end
end
% Matrix of simulation results which are corespond to experimental data
final_new_f = [new_f(:,2), new_f(:,6), new_f(:,8), new_f(:,11)];
% Difference between simulation results and experimental data
final = final_new_f - final_E;
% It comes parameters in the same range
ddd = zeros(hh,4);
eee = zeros(hh,4);
for j = 1:4
for i = 1:hh
ddd(i,j) = final(i,j)/std_dev(1,j);
eee(i,j) = ddd(i,j).^2;
end
end
% Transfer from matrix to a column vector
final_f = eee(:);
========================================================= Function for ODE
function f = PFR_DGL(g,n,URH,HH)
% Global Parameters
%%Value transfer
c(1:n) = g(1:n); % Concentrations in the liquid phase, mol/m3
T = g(n+1); % Temperature, K
p0 = g(n+2); % Pressure, Pa
VR = g(n+3); % Reactor volume, m3
% The stoichiometric coefficient matrix
n_u = nu_ij_Po_u(0);
% Reaction rates [mol/(m3 s)]
ru = rate_Po2_u(c,T,HH);
%%Change of Molar rates in the liquid phase, mol/(m3 s)
Ri = n_u'*ru;
%%Output
f(1:n-1) = Ri/URH; % Concentrations of the components 1 to n-1, mol/m3
f(14) = 0; % Concentration of the Cat., mol/m3
f(n+1) = 0; % Temperature change, K/s
f(n+2) = dpdz; % Pressure change, Pa/s
f(n+3) = dVdz; % Volume change, m3/s
% Transfer of f as a column vector
f = f';
============================================ Needed function for ODE function
function ru = rate_Po2_u(c,T,HH)
% Global Paranmeter
% Collision factors,
k0u = a row vector with size 1*22
% Activation Energies, J/mol
Eu = a row vector with size 1*22
% Reaction rate constant at temperature T
kTu = k0u.*exp(-Eu*RTi);
% Reaction rates
ru = [ HH(1)*c(1)*c(2)
kTu(2)*c(3)*c(2)
HH(2)*c(1)*c(5)
kTu(4)*c(1)*c(7)
HH(3)*c(6)
HH(4)*c(6)*c(8)
kTu(7)*c(6)*c(11)
kTu(8)*c(6)*c(1)
kTu(9)*c(6)*c(7)
kTu(10)*c(6)*c(5)
kTu(11)*c(8)*c(5)
kTu(12)*c(11)*c(5)
kTu(13)*c(3)*c(6)
kTu(14)*c(3)*c(6)
HH(5)*c(5)*c(5)
kTu(16)*c(12)*c(6)
kTu(17)*c(8)*c(12)
kTu(18)*c(5)*c(5)
kTu(19)*c(1)*c(9)
kTu(20)*c(1)*c(7)
kTu(21)*c(1)*c(7)
kTu(22)*c(12)];
  4 Comments
Matt Kindig
Matt Kindig on 14 Aug 2013
Edited: Matt Kindig on 14 Aug 2013
Can you post the full contents of output.message ?
Mehr Markazi
Mehr Markazi on 14 Aug 2013
Hey Matt Kinding, output.message says: Solver stopped prematurely. lsqnonlin stopped because it exceeded the function evaluation limit, options.MaxFunEvals = 500 (the default value).
I know I should increase the MaxFunEvals. I increased it until 1000, but again it stops because of the same massage.

Sign in to comment.

Accepted Answer

Matt J
Matt J on 14 Aug 2013
Edited: Matt J on 14 Aug 2013
The algorithm is not converging. It is running until the limits on MaxFunEvals and MaxIter are reached and then stopping.
One reason this might occur is if your function is not scaled well, e.g., its derivatives even very close to the solution are orders of magnitude larger than the ratio TolFun/TolX.

More Answers (0)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!