How to get variable value in nonlinear model fitting

Hi everyone,

I'm currently facing a seemingly simple but confusing problem in running nonlinear model fitting with LSQNONLIN. The entire code is attached as ZIP file, which consists of multiple .m files. Parallel computing was used to speed up the calculation, as seen in parfor in compute_y_hat_NTA_H2.m. The code itself is working (presumably as intended) and executable through optimize_NTA_H2.m.

The particularly problematic .m file in question is calc_r_r_NTA_H.m, which is shown below. I want to obtain/print the values of r_r = [r_r_NH3_1, r_r_NH3_2, r_r_N2D1_1, r_r_N2D1_2] (shown in Line 154 of the code), which is supposed to be in form of matrix with four rows, each representing _r_NH3_1, r_r_NH3_2, r_r_N2D1_1, r_r_N2D1_2, respectively. However, no matter how I tried, so far it only yields error messages (i.e. variable not found etc.).

Could someone help me to find the solution? Your efforts are much appreciated.

function r_r = calc_r_r_NTA_H(T, x_NO_1, x_NO_2)

%  Rate equation is Langmuir-Hinshelwood mechanism
%  Reverse reaction not assumed in the model
%  
%  A(i) and E(i) to be fitted

global R
global searchParams 
global Kf_NH3_1 Kf_NH3_2  Kf_N2_D1_1 Kf_N2_D1_2 
global Ka_H2 Ka_NO Ka_CO Ka_OH_1 Ka_OH_2 Ka_COOH Ka_NHx_1 Ka_NHx_2 %Ka_H1 
global P_t
global F_NO_0 F_H2_0 F_NH3_0 F_N2_0 F_CO_0 F_CO2_0 F_H2O_0 F_Ar_0
global p_NO p_H2 p_NH3 p_N2 p_CO p_CO2 p_H2O p_Ar
global F_total
global x_NH3_1 x_NH3_2 x_N2_1 x_N2_2

global ADS 
%int ratio_r_NH3 = 0.5;
%int ratio_r_N2 = 0.5;
global ratio_r_NH3 ratio_r_N2 %ratio_r

   if isempty(ratio_r_NH3) || isempty(ratio_r_N2)
       ratio_r_NH3 = 0.5;
       ratio_r_N2  = 0.3;
       
       % conversion to NH3 with CO
        x_NH3_1 = x_NO_1*ratio_r_NH3;
        % conversion to NH3 withOUT CO
        x_NH3_2 = x_NO_1*(1-ratio_r_NH3);
        % conversion to N2 with CO
        x_N2_1 = x_NO_2*ratio_r_N2;
        % conversion to N2 withOUT CO
        x_N2_2 = x_NO_2*(1-ratio_r_N2);
   end
    
  ratio_r = [ratio_r_NH3, ratio_r_N2];
   
R = 8.314;


%x_N2 = y(2);
% reaction to form NH3, with CO
% 4NO + 5CO + H2O + 5H2 <-> 4NH3 + 5CO2
% 2NO + 5CO + 3H2O <-> 2NH3 + 5CO2     x_NH3_1
% or without CO
% 2NO + 5H2 <-> 2NH3 + 2H2O            x_NH3_2

% reaction to form N2, with CO
% 2NO + 2CO <-> N2 + 2CO2              x_N2_1
% without CO
% 2NO + 2H2 <-> N2 + 2H2O              x_N2_2

F_NO = F_NO_0 - F_NO_0 * (x_NO_1 + x_NO_2);
F_H2 = F_H2_0 - F_NO_0*(2.5* x_NH3_2 + x_N2_2); % with H2 feed
%F_H2 = F_H2_0 - F_NO_0*(1.25 * x_NH3_1 + 2.5* x_NH3_2 + x_N2_2); % with H2 feed
%F_H2 = F_H2_0;  % without H2 feed
F_CO = F_CO_0 - F_NO_0*(2.5 * x_NH3_1 + x_N2_1);
%F_CO = F_CO_0 - F_NO_0*(1.25 * x_NH3_1 + 0.5*x_N2_1);
F_H2O = F_H2O_0 - F_NO_0*(1.5 * x_NH3_1 - x_NH3_2 - x_N2_2);
%F_H2O = F_H2O_0 - F_NO_0*(0.25 * x_NH3_1 - x_N2_2);
F_NH3 = F_NH3_0 + F_NO_0 * (x_NH3_1 + x_NH3_2); 
%F_NH3 = F_NH3_0 + F_NO_0 * (x_NH3_1 + 0.5*x_NH3_2); 
F_N2 = F_N2_0 + F_NO_0 * (0.5*x_N2_1 + 0.5*x_N2_2);
%F_N2 = F_N2_0 + F_NO_0 * (x_N2_1 + 0.5*x_N2_2);
%F_CO2 = F_CO2_0 + 1.75* F_NO_0 * x_NO_1 + 0.5* F_NO_0 * x_NO_2;
F_CO2 = F_CO2_0 + F_NO_0 * (2.5 * x_NH3_1 + x_N2_1);
F_Ar = F_Ar_0;

F_total = F_NO + F_CO + F_H2O + F_H2 + F_NH3 + F_N2 + F_CO2 + F_Ar;

p_NO = F_NO / F_total * P_t;
%p_NO = F_NO ./ F_total .* P_t;
p_CO = F_CO / F_total * P_t;
p_H2O = F_H2O / F_total * P_t;
p_H2 = F_H2 / F_total * P_t;
p_NH3 = F_NH3 / F_total * P_t;
p_N2 = F_N2 / F_total * P_t;
p_CO2 = F_CO2 / F_total * P_t;
p_Ar = F_Ar / F_total * P_t;
%p_Ar = F_Ar ./ F_total .* P_t;



%  reaction coefficients (starting from SP3)
%  reaction coefficients
Kf_NH3_1 = (searchParams(1)*10000)*exp(-searchParams(2)*10000/(R*T));
Kf_N2_D1_1 =(searchParams(3)*10000)*exp(-searchParams(4)*10000/(R*T));
%Ka_H1 = searchParams(5)*exp(searchParams(6)*10000/(R*T));   %NOT USED, REACTION WITH DISS H2 ALREADY IN NO+H2
Ka_H2 = searchParams(5)*exp(searchParams(6)*10000/(R*T));
Ka_NO = searchParams(7)*exp(searchParams(8)*10000/(R*T));
Ka_CO = searchParams(9)*exp(searchParams(10)*10000/(R*T));
Ka_COOH = searchParams(11)*exp(searchParams(12)*10000/(R*T));
Ka_OH_1 = searchParams(13)*exp(searchParams(14)*10000/(R*T));
%Ka_N_1 = searchParams(17)*exp(searchParams(18)*10000/(R*T));
Ka_NHx_1 = searchParams(15)*exp(searchParams(16)*10000/(R*T));

Kf_NH3_2 = (searchParams(17)*10000)*exp(-searchParams(18)*10000/(R*T));
Kf_N2_D1_2 =(searchParams(19)*10000)*exp(-searchParams(20)*10000/(R*T));
Ka_OH_2 = searchParams(21)*exp(searchParams(22)*10000/(R*T));
%Ka_N_2 = searchParams(25)*exp(searchParams(26)*10000/(R*T));
Ka_NHx_2 = searchParams(23)*exp(searchParams(24)*10000/(R*T));

%A_H = Ka_H1*(p_H2^0.5)+Ka_H2*(p_H2O^0.5)*(p_CO^0.5);
A_H = Ka_H2*(p_H2O^0.5)*(p_CO^0.5);
%A_H = Ka_H1*(p_H2^0.5)+Ka_H2*(p_H2O^0.5)*(p_CO^0.5)*((p_CO2+0.00001)^-0.5);

% REDUCTION WITH CO AND H2O
Bf_NH3_1 = Kf_NH3_1 * p_NO*(p_CO^0.5)*((p_H2O^-0.5))*A_H^4;        %forward numerator NH3 formation
Bf_N2_1 = Kf_N2_D1_1 * (p_NO^2)*(p_CO^searchParams(25));        %forward numerator N* -> N2

%Ba_H = Ka_H1*(p_H2^0.5)+Ka_H2*(p_H2O^0.5)*(p_CO^0.5)*((p_CO2+0.00001)^-0.5);
Ba_H = Ka_H2*(p_H2O^0.5)*(p_CO^0.5)*((p_CO2+0.00001)^-0.5);
Ba_NO = Ka_NO*p_NO;
Ba_CO = Ka_CO*p_CO;
Ba_COOH = Ka_COOH*(p_CO^0.5)*(p_H2O^0.5)*((p_CO2+0.00001)^0.5);
Ba_OH_1 = Ka_OH_1*((p_CO+0.00001)^-0.5)*(p_H2O^0.5)*((p_CO2+0.00001)^0.5);
%Ba_N_1 = Ka_N_1*p_NO*(p_H2O^-0.5)*(p_CO^0.5)*(p_CO2^-0.5)*A_H;
Ba_NHx_1 = Ka_NHx_1*p_NO*((p_H2O+0.00001)^-0.5)*(p_CO^0.5)*((p_CO2+0.00001)^-0.5)*A_H^searchParams(26);

% REDUCTION WITH ONLY H2
Bf_NH3_2 = Kf_NH3_2 * p_NO*(p_H2^2.5);                             %forward numerator NH3 formation
Bf_N2_2 = Kf_N2_D1_2 * (p_NO^2)*(p_H2^searchParams(27));        %forward numerator N* -> N2

Ba_OH_2 = Ka_OH_2*(p_H2O^1)*((p_H2+0.00001)^-1.5);
%Ba_N_2 = Ka_N_2*p_NO*(p_H2O^-0.5)*(p_CO^0.5)*(p_CO2^-0.5)*A_H;
Ba_NHx_2 = Ka_NHx_2*p_NO*((p_H2O+0.00001)^-1)*(p_H2^searchParams(28));


% reaction rate equation

% reaction rate total of the following reactions
% NH2* + H* <-> NH3(g) + 2*                Rx C3
% 2N* <-> N2(g) + 2*                       Rx D1

ADS = 1+Ba_H+Ba_NO+Ba_CO+Ba_COOH+Ba_OH_1+Ba_NHx_1+Ba_OH_2+Ba_NHx_2;

%r_r_X_1 is rate or X formation reaction, CO-dependent; ..._2 is
%non-CO-dependent (reduction with H2)

% NH3 formation rate
r_r_NH3_1 = Bf_NH3_1/ADS^4;
r_r_NH3_2 = Bf_NH3_2/ADS^4;

% N2 formation rate
r_r_N2D1_1 = Bf_N2_1/ADS^2;
r_r_N2D1_2 = Bf_N2_2/ADS^2;

%ratio_r_X is ratio of X formation rate CO-dependent per overall 
ratio_r_NH3 = r_r_NH3_1/(r_r_NH3_1 + r_r_NH3_2);
ratio_r_N2 = r_r_N2D1_1/(r_r_N2D1_1 + r_r_N2D1_2);

% Rx C3 rate
r_r = [r_r_NH3_1, r_r_NH3_2, r_r_N2D1_1, r_r_N2D1_2];

end