Worse result when using "lsqnonlin"

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)];

[Matt J]

What is the rationale behind TolFun=1e-3 instead of its default 1e-6 and similary why DiffMaxChange=1e+15 instead of Inf?

Also, what were the results of the lsqnonlin output arguments "exitflag" and "output"?

[Mehr Markazi]

Thanks for your answer, Matt.

I tried different value for them. I also used the default value for them. But I couldn't get the answer that I want.

The result for "exitflag" is "0". Here is the results for output:

output.firstorderopt: Value(0.439204516011054) Min(0.439204516011054) Max(0.439204516011054)

output.iterations: Value(83) Min(83) Max(83)

output.funcCount: Value(504) Min(504) Max(504)

output.cgiterations: Value(0) Min(0) Max(0)

output.algorithm: 'trust-region-reflective'

output.message: char

[Matt Kindig]

Can you post the full contents of output.message ?

[Mehr Markazi]

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.

Answer 1

[Matt J]

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.