Too many output arguments error when using lsqcurvefit

조회 수: 3 (최근 30일)
Dursman Mchabe
Dursman Mchabe 2018년 10월 1일
댓글: Dursman Mchabe 2018년 10월 1일

I have a set of 7 ordinary differential equations and 1 algebraic which I use to calculate the concentrations of 8 species with time. I have also conducted experiments to measure the concentrations of such 8 species. The challenge I have now that the concentrations predicted by the equations are not the same as the concentrations measured during experiments. As a result, I am looking at ways the regress the 4 parameters and 8 initial conditions, such that the calculated concentrations can match the experimental concentrations. I have found an example on the internet, from this link:

https://www.mathworks.com/matlabcentral/answers/329901-solving-coupled-differential-equations

in which the approach of using matlab's lsqcurvefit was followed. I Have tried to follow the step by step, however, I get the error message:

Error using EstimationOfParametersAndInitialConditionsUsingLSQCURVEFIT3/ode15ifun
Too many output arguments.
Error in lsqcurvefit (line 213)
            initVals.F = feval(funfcn_x_xdata{3},xCurrent,XDATA,varargin{:});
Error in EstimationOfParametersAndInitialConditionsUsingLSQCURVEFIT3 (line 155)
B = lsqcurvefit(@ode15ifun, eqns0, ExperimentalTime, ExperimentalConcentrations)
Caused by:
    Failure in initial objective function evaluation. LSQCURVEFIT cannot continue.

My current code is:

function EstimationOfParametersAndInitialConditionsUsingLSQCURVEFIT3
global ExperimentalTime ExperimentalConcentrations
function ode15ifun
ExperimentalTime = [600
                    1200
                    1800
                    2400
                    3000
                    10200
                    17400
                    24600
                    31800
                    39000
                    46200
                    53400
                    60600
                    67800
                    75000
                    82200
                    89400
                    96600
                    103800
                    111000
                    118200
                    125400
                    132600
                    139800
                    147000
                    154200
                    161400
                    168600
                    175800
                    183000];
ExperimentalConcentrations = [  0.00E+00    6.97E-08    3.34E-01    4.88E-06    4.88E-06    4.86E+01    0.00E+00    7.41E-06
                                1.72E-02    8.76E-08    6.79E-01    4.88E-06    4.88E-06    4.92E+01    0.00E+00    9.33E-06
                                1.71E-02    1.18E-07    1.04E+00    5.11E-06    5.11E-06    4.94E+01    0.00E+00    1.20E-05
                                1.70E-02    1.03E-07    1.41E+00    5.11E-06    5.11E-06    4.96E+01    0.00E+00    1.05E-05
                                1.68E-02    1.77E-07    6.00E+00    5.35E-06    5.35E-06    3.78E+01    6.30E+00    1.74E-05
                                1.74E-02    3.65E-07    1.07E+01    5.35E-06    5.35E-06    3.09E+01    9.33E+00    3.72E-05
                                1.65E-02    3.50E-07    1.56E+01    5.36E-06    5.36E-06    2.67E+01    1.13E+01    3.55E-05
                                1.67E-02    8.83E-07    2.07E+01    5.36E-06    5.36E-06    1.98E+01    1.43E+01    1.00E-04
                                1.67E-02    5.01E-06    2.59E+01    5.07E-06    5.07E-06    1.14E+01    1.79E+01    4.07E-02
                                1.70E-02    5.06E-06    3.14E+01    5.07E-06    5.07E-06    9.58E+00    1.91E+01    2.45E-01
                                1.80E-02    4.93E-06    3.67E+01    4.93E-06    4.93E-06    5.65E+00    2.09E+01    6.17E-01
                                2.15E-02    5.29E-06    4.18E+01    5.30E-06    5.30E-06    3.32E+00    2.21E+01    1.32E+00
                                2.68E-02    5.30E-06    4.66E+01    5.30E-06    5.30E-06    0.00E+00    2.37E+01    2.29E+00
                                3.07E-02    5.30E-06    5.09E+01    5.30E-06    5.30E-06    0.00E+00    2.40E+01    2.34E+00
                                3.33E-02    9.45E-06    5.51E+01    9.46E-06    9.46E-06    0.00E+00    2.43E+01    2.40E+00
                                3.68E-02    9.45E-06    5.90E+01    9.46E-06    9.46E-06    0.00E+00    2.47E+01    1.82E+00
                                4.14E-02    1.13E-05    6.23E+01    1.13E-05    1.13E-05    0.00E+00    2.50E+01    1.38E+00
                                4.48E-02    1.15E-05    6.56E+01    1.15E-05    1.15E-05    0.00E+00    2.54E+01    2.09E+00
                                4.87E-02    1.10E-05    6.87E+01    1.10E-05    1.10E-05    0.00E+00    2.58E+01    1.82E+00
                                5.19E-02    1.10E-05    7.12E+01    1.10E-05    1.10E-05    0.00E+00    2.62E+01    1.58E+00
                                5.47E-02    1.04E-05    7.36E+01    1.04E-05    1.04E-05    0.00E+00    2.67E+01    2.29E+00
                                5.75E-02    1.04E-05    7.59E+01    1.04E-05    1.04E-05    0.00E+00    2.71E+01    1.62E+00
                                5.97E-02    1.04E-05    7.78E+01    1.04E-05    1.04E-05    0.00E+00    2.76E+01    1.12E+00
                                6.17E-02    1.03E-05    7.95E+01    1.03E-05    1.03E-05    0.00E+00    2.81E+01    8.91E-01
                                6.33E-02    1.04E-05    8.11E+01    1.04E-05    1.04E-05    0.00E+00    2.86E+01    8.51E-01
                                6.44E-02    1.05E-05    8.24E+01    1.06E-05    1.06E-05    0.00E+00    2.90E+01    7.59E-01
                                6.48E-02    1.13E-05    8.24E+01    1.13E-05    1.13E-05    0.00E+00    2.90E+01    8.71E-01
                                6.53E-02    1.13E-05    8.23E+01    1.13E-05    1.13E-05    0.00E+00    2.90E+01    1.12E+00
                                6.57E-02    1.13E-05    8.21E+01    1.13E-05    1.13E-05    0.00E+00    2.90E+01    1.00E+00
                                6.59E-02    1.03E-05    8.20E+01    1.03E-05    1.03E-05    0.00E+00    2.90E+01    8.51E-01];
%-------------------------------------------------------------------------------------------------------------------------------------
   % b(1:4) = Parameters,  b(5:12) = Initial Conditions
   % MAPPING: b(1) = kga,  b(2) = kLa_SO2,  b(3) = ktot, b(4) = Kad, 
   %          b(5)= CSO2_gas(t0), b(6)= CCO2_gas(t0), b(7)= S_total(t0), b(8)= C_total(t0), b(9)= Ca2_total(t0),b(10)= CCaCO3(t0),b(11)= CCaSO3(t0),b(12)= CH(t0)   
syms  b V_Headspace F_rate CSO2_in CCO2_in R T HSO2 DCa2 DSO2 DHSO3 DSO32 kLa_CO2 HCO2 DCO2 DHCO3 DCO32 KSPCaSO3 BETCaSO3 BETCaCO3 MWCaCO3 KCO2 KHCO3 KSO2 KHSO3 Kw CCa2
eqn1 = diff(b(5),t)== F_rate/ V_Headspace * ( CSO2_in -  b(5)) - (b(5) * R * T - HSO2 * ((b(7) * b(12)^2)/(b(12)^2 + KSO2 * b(12) + KSO2 * KHSO3))) / ((1/b(1)) + 1/(b(2) * (((DSO2 * (HSO2 * b(5) * R * T - ((b(7) * b(12)^2)/(b(12)^2 + KSO2 * b(12) + KSO2 * KHSO3))) + DHSO3 * (((KSO2 * HSO2 *  b(5) * R * T / b(12)) - ((b(7) * KSO2 * b(12))/(b(12)^2 + KSO2 * b(12) + KSO2 * KHSO3)))) + DSO32 * ((((KSO2 * KHSO3 * HSO2 * b(5) * R * T/b(12)^2) - ((b(7) * KSO2 * KHSO3)/(b(12)^2 + KSO2 * b(12) + KSO2 * KHSO3)))))) / (DSO2 * (HSO2 * b(5) * R * T - ((b(7) * b(12)^2)/(b(12)^2 + KSO2 * b(12) + KSO2 * KHSO3))))))));
eqn2 = diff(b(6) ,t)== F_rate/ V_Headspace * ( CCO2_in -  b(6)) - (kLa_CO2 * ((((DCO2*(HCO2 * b(6) * R * T - ((b(8) * b(12)^2)/(b(12)^2 + KCO2 * b(12) + KCO2 * KHCO3))) + DHCO3*(((KCO2 * HCO2 *  b(6) * R * T/b(12)) - ((C_total (t) * KCO2 * b(12))/(b(12)^2 + KCO2 * b(12) + KCO2 * KHCO3)))) + DCO32 * ((((KCO2 * KHCO3 * HCO2 * b(6) * R * T/b(12)^2) - ((b(8) * KCO2 * KHCO3)/(b(12)^2 + KCO2 * b(12) + KCO2 * KHCO3)))))) / (DCO2 * (HCO2 * b(6) * R * T - ((C_total (t) * b(12)^2)/(b(12)^2 + KCO2 * b(12) + KCO2 * KHCO3)))))))*(b(6) * R * T/HCO2 - ((b(8) * b(12)^2)/(b(12)^2 + KCO2 * b(12) + KCO2 * KHCO3))));
eqn3 = diff(b(7),t)== (b(5) * R * T - HSO2*((b(7) * b(12)^2)/(b(12)^2 + KSO2 * b(12) + KSO2 * KHSO3))) / ((1/b(1)) + 1/(b(2) * (((DSO2*(HSO2* b(5) * R * T - ((b(7) * b(12)^2)/(b(12)^2 + KSO2 * b(12) + KSO2 * KHSO3))) + DHSO3*(((KSO2 * HSO2 *  b(5) * R * T/b(12)) - ((b(7) * KSO2 * b(12))/(b(12)^2 + KSO2 * b(12) + KSO2 * KHSO3)))) + DSO32 * ((((KSO2 * KHSO3 *HSO2* b(5) * R * T/b(12)^2) - ((b(7) * KSO2 * KHSO3)/(b(12)^2 + KSO2 * b(12) + KSO2 * KHSO3)))))) / (DSO2 * (HSO2 * b(5) * R * T - ((b(7) * b(12)^2)/(b(12)^2 + KSO2 * b(12) + KSO2 * KHSO3))))))) - 0.162 * exp(-5153/T) * (((CCa2 * ((b(7) * KSO2 * KHSO3)/(b(12)^2 + KSO2 * b(12) + KSO2 * KHSO3))) / KSPCaSO3) - 1)^2 * (BETCaSO3/((CCa2 * (((((KSO2 * KHSO3 *HSO2* b(5) * R * T/b(12)^2) - ((b(7) * KSO2 * KHSO3)/(b(12)^2 + KSO2 * b(12) + KSO2 * KHSO3)))))) / (DSO2 * (HSO2 * b(5) * R * T - ((b(7) * b(12)^2)/(b(12)^2 + KSO2 * b(12) + KSO2 * KHSO3))))))))/ KSPCaSO3;
eqn4 = diff(b(8),t)== (kLa_CO2 * ((((DCO2*(HCO2 * b(6) * R * T - ((b(8) * b(12)^2)/(b(12)^2 + KCO2 * b(12) + KCO2 * KHCO3))) + DHCO3*(((KCO2 * HCO2 *  b(6) * R * T/b(12)) - ((C_total (t) * KCO2 * b(12))/(b(12)^2 + KCO2 * b(12) + KCO2 * KHCO3)))) + DCO32 * ((((KCO2 * KHCO3 * HCO2 * b(6) * R * T/b(12)^2) - ((b(8) * KCO2 * KHCO3)/(b(12)^2 + KCO2 * b(12) + KCO2 * KHCO3)))))) / (DCO2 * (HCO2 * b(6) * R * T - ((C_total (t) * b(12)^2)/(b(12)^2 + KCO2 * b(12) + KCO2 * KHCO3)))))))*(b(6) * R * T /HCO2 - ((b(8) * b(12)^2)/(b(12)^2 + KCO2 * b(12) + KCO2 * KHCO3))))  + (b(3) * BETCaCO3 * MWCaCO3 * b(10) * b(12) *(1 - (b(4) * b(12))/(1 + b(4) * b(12))));
eqn5 = diff(b(9),t)== (-1) * (b(3) * BETCaCO3 * MWCaCO3 * b(10) * b(12) *(1 - (b(4) * b(12))/(1 + b(4) * b(12)))) - (0.162 * exp(-5153/T) * (((CCa2 * ((b(7) * KSO2 * KHSO3)/(b(12)^2 + KSO2 * b(12) + KSO2 * KHSO3))) / KSPCaSO3) - 1)^2 * (BETCaSO3/((CCa2 * (((((KSO2 * KHSO3 *HSO2* b(5) * R * T/b(12)^2) - ((b(7) * KSO2 * KHSO3)/(b(12)^2 + KSO2 * b(12) + KSO2 * KHSO3)))))) / (DSO2 * (HSO2 * b(5) * R * T - ((b(7) * b(12)^2)/(b(12)^2 + KSO2 * b(12) + KSO2 * KHSO3))))))))/ KSPCaSO3;
eqn6 = diff(b(10),t)== (-1) * (b(3) * BETCaCO3 * MWCaCO3 * b(10) * b(12) * (1 - (b(4) * b(12))/(1 + b(4) * b(12))));
eqn7 = diff(b(11),t)== 0.162 * exp(-5153/T) * ((((CCa2 * ((b(7) * KSO2 * KHSO3)/(b(12)^2 + KSO2 * b(12) + KSO2 * KHSO3))) / KSPCaSO3) - 1)^2 * (BETCaSO3/((CCa2 * (((((KSO2 * KHSO3 *HSO2* b(5) * R * T/b(12)^2) - ((b(7) * KSO2 * KHSO3)/(b(12)^2 + KSO2 * b(12) + KSO2 * KHSO3)))))) / (DSO2 * (HSO2 * b(5) * R * T - ((b(7) * b(12)^2)/(b(12)^2 + KSO2 * b(12) + KSO2 * KHSO3))))))))/ KSPCaSO3;
eqn8 = diff(b(12),t)== b(12) + 2 * CCa2 - ((b(7) *KSO2 * b(12))/(b(12)^2 + KSO2 * b(12) + KSO2 * KHSO3))-  2*((b(7) * KSO2 * KHSO3)/(b(12)^2 + KSO2 * b(12) + KSO2 * KHSO3)) - ((b(8) * KCO2 * b(12))/(b(12)^2 + KCO2 * b(12) + KCO2 * KHCO3)) - 2 * ((b(8) * KCO2 * KHCO3)/(b(12)^2 + KCO2 * b(12) + KCO2 * KHCO3))- Kw / b(12);
eqns = [eqn1 eqn2 eqn3 eqn4 eqn5 eqn6 eqn7 eqn8];
vars = [b(5); b(6); b(7); b(8); b(9); b(10); b(11); b(12)];
origVars = length(vars);
M = incidenceMatrix(eqns, vars);
[eqns, vars] = reduceDifferentialOrder(eqns, vars);
isLowIndexDAE(eqns,vars);
f = daeFunction(eqns,vars, V_Headspace, F_rate, CSO2_in, CCO2_in, CCa2, R, T, HSO2, b(1), DCa2, DSO2, DHSO3, DSO32, b(2), kLa_CO2, HCO2, DCO2, DHCO3, DCO32, KSPCaSO3, BETCaSO3, b(3), BETCaCO3, MWCaCO3, b(4), KCO2, KHCO3, KSO2, KHSO3, Kw);
V_Headspace = 1.5e-6;       
F_rate      = 1.66667e-5;    
CSO2_in     = 6.51332e-2;   
CCO2_in     = 0;            
CCa2        = 10;           
R           = 8.314;        
T           = 323.15;       
HSO2        = 149;                                              
b(1)        = 4.14e-6;                                          
DCa2        = 1.39e-9;                                          
DSO2        = 2.89e-9;                                    
DHSO3       = 2.89e-9;                                     
DSO32       = 2.89e-9;                                    
b(2)        = 8.4e-4;                                           
kLa_CO2     = 9.598e-4;                                         
HCO2        = 5.15e3;                                           
DCO2        = 3.53e-9;                                    
DHCO3       = 2.89e-9;                                    
DCO32       = 2.89e-9;                                    
KSPCaSO3    = 1.07e-7;                                          
BETCaSO3    = 10;                                              
b(3)        = 8.825e-3;                                         
BETCaCO3    = 12.54;                                            
MWCaCO3     = 100.0869;                                         
b(4)         = 0.84;                                             
KCO2        = 1.7e-3;                                           
KHCO3       = 6.55e-8;                                          
KSO2        = 6.24;                                             
KHSO3       = 5.68e-5;                                          
Kw          = 5.3e-8;                                          
F = @(t, Y, YP) f(t, Y, YP, V_Headspace, F_rate, CSO2_in, CCO2_in, CCa2, R, T, HSO2, b(1), DCa2, DSO2, DHSO3, DSO32, b(2), kLa_CO2, HCO2, DCO2, DHSO3, DSO32, KSPCaSO3, BETCaSO3, b(3), BETCaCO3, MWCaCO3, b(4), KCO2, KHCO3, KSO2, KHSO3, Kw);
vars;                                   
%b0est = [0; 0; 3.3351e-1; 4.879e-6; 4.879e-6; 4.862e+1; 0; 7.413e-6];
 b0est = [1;1;1;1;1;1;1;1;];                                   
  bp0est = zeros(8,1);                                  
 opt = odeset('RelTol', 10.0^(-7), 'AbsTol' , 10.0^(-7));                                   
 [b0, bp0] = decic(F, 0, b(5:12), [], bp0est, [], opt);                                   
 [tSol,eqns] = ode15i(F, ExperimentalTime, b0, bp0, opt);                                   
 for k = 1:origVars                                   
 S{k} = char(vars(k));                                  
 end                                      
end
eqns0 = [4.14e-6;  8.4e-4;  8.825e-3;  0.84; 0; 0; 3.3351e-1; 4.879e-6; 4.879e-6; 4.862e+1; 0; 7.413e-6];
%B = lsqcurvefit(@ode15ifun, eqns0, ExperimentalTime, ExperimentalConcentrations, zeros(size(eqns0)), inf(size(eqns0)))
B = lsqcurvefit(@ode15ifun, eqns0, ExperimentalTime, ExperimentalConcentrations)
fprintf(1, '\n\tParameters:\n\t\tkga = %11.3E\n\t\tkLa = %11.3E\n\t\tktot = %11.3E\n\t\tKad = %11.3E\n\tInitial Conditions:\n\t\tCSO2gas = %11.3E\n\t\tCCO2gas = %11.3E\n\t\tStotal = %11.3E\n\t\tCtotal = %11.3E\n\t\tCatotal = %11.3E\n\t\tCCaCO3 = %11.3E\n\t\tCCaSO3 = %11.3E\n\t\tCH = %11.3E\n', B);
F = @(t, B, BP) f(t, B, BP, V_Headspace, F_rate, CSO2_in, CCO2_in, CCa2, R, T, HSO2, DCa2, DSO2, DHSO3, DSO32, kLa_CO2, HCO2, DCO2, DHSO3, DSO32, KSPCaSO3, BETCaSO3, BETCaCO3, MWCaCO3, KCO2, KHCO3, KSO2, KHSO3, Kw);
[tSol,eqns] = ode15i(F, ExperimentalTime, B(5:12), bp0, opt);
figure(2)
plot(ExperimentalTime, ExperimentalConcentrations, 'p')
hold on
plot(tSol, ySol, '--')
hold off
grid
legend('Exp-SO_{2-gas}','Exp-CO_{2-gas}','Exp-S_{total}','Exp-C_{total}','Exp-Ca^{2+}_{total}','Exp-CaCO_{3}','Exp-CaSO_{3}.^{1}/_{2}H_{2}O','Exp-H^{+}', 'Mod-SO_{2-gas}','Mod-CO_{2-gas}','Mod-S_{total}','Mod-C_{total}','Mod-Ca^{2+}_{total}','Mod-CaCO_{3}','Mod-CaSO_{3}.^{1}/_{2}H_{2}O','Mod-H^{+}', 'Location','best')
end

I can not spot where I am going wrong. I do not know where to fix. What can I try?

When I run the code without including the lsqcurvefit, the code does not give errors, however, the experimental results does not match the calculated values. I think that when I regress 4 input values, the calculated values will fit the experimental values. Below is the code on which I have commented-out the lsqcurvefit:

%function EstimationOfParametersAndInitialConditionsUsingLSQCURVEFIT2
function ode15ifun
ExperimentalTime = [600
                    1200
                    1800
                    2400
                    3000
                    10200
                    17400
                    24600
                    31800
                    39000
                    46200
                    53400
                    60600
                    67800
                    75000
                    82200
                    89400
                    96600
                    103800
                    111000
                    118200
                    125400
                    132600
                    139800
                    147000
                    154200
                    161400
                    168600
                    175800
                    183000];
ExperimentalConcentrations = [  0.00E+00    6.97E-08    3.34E-01    4.88E-06    4.88E-06    4.86E+01    0.00E+00    7.41E-06
                                1.72E-02    8.76E-08    6.79E-01    4.88E-06    4.88E-06    4.92E+01    0.00E+00    9.33E-06
                                1.71E-02    1.18E-07    1.04E+00    5.11E-06    5.11E-06    4.94E+01    0.00E+00    1.20E-05
                                1.70E-02    1.03E-07    1.41E+00    5.11E-06    5.11E-06    4.96E+01    0.00E+00    1.05E-05
                                1.68E-02    1.77E-07    6.00E+00    5.35E-06    5.35E-06    3.78E+01    6.30E+00    1.74E-05
                                1.74E-02    3.65E-07    1.07E+01    5.35E-06    5.35E-06    3.09E+01    9.33E+00    3.72E-05
                                1.65E-02    3.50E-07    1.56E+01    5.36E-06    5.36E-06    2.67E+01    1.13E+01    3.55E-05
                                1.67E-02    8.83E-07    2.07E+01    5.36E-06    5.36E-06    1.98E+01    1.43E+01    1.00E-04
                                1.67E-02    5.01E-06    2.59E+01    5.07E-06    5.07E-06    1.14E+01    1.79E+01    4.07E-02
                                1.70E-02    5.06E-06    3.14E+01    5.07E-06    5.07E-06    9.58E+00    1.91E+01    2.45E-01
                                1.80E-02    4.93E-06    3.67E+01    4.93E-06    4.93E-06    5.65E+00    2.09E+01    6.17E-01
                                2.15E-02    5.29E-06    4.18E+01    5.30E-06    5.30E-06    3.32E+00    2.21E+01    1.32E+00
                                2.68E-02    5.30E-06    4.66E+01    5.30E-06    5.30E-06    0.00E+00    2.37E+01    2.29E+00
                                3.07E-02    5.30E-06    5.09E+01    5.30E-06    5.30E-06    0.00E+00    2.40E+01    2.34E+00
                                3.33E-02    9.45E-06    5.51E+01    9.46E-06    9.46E-06    0.00E+00    2.43E+01    2.40E+00
                                3.68E-02    9.45E-06    5.90E+01    9.46E-06    9.46E-06    0.00E+00    2.47E+01    1.82E+00
                                4.14E-02    1.13E-05    6.23E+01    1.13E-05    1.13E-05    0.00E+00    2.50E+01    1.38E+00
                                4.48E-02    1.15E-05    6.56E+01    1.15E-05    1.15E-05    0.00E+00    2.54E+01    2.09E+00
                                4.87E-02    1.10E-05    6.87E+01    1.10E-05    1.10E-05    0.00E+00    2.58E+01    1.82E+00
                                5.19E-02    1.10E-05    7.12E+01    1.10E-05    1.10E-05    0.00E+00    2.62E+01    1.58E+00
                                5.47E-02    1.04E-05    7.36E+01    1.04E-05    1.04E-05    0.00E+00    2.67E+01    2.29E+00
                                5.75E-02    1.04E-05    7.59E+01    1.04E-05    1.04E-05    0.00E+00    2.71E+01    1.62E+00
                                5.97E-02    1.04E-05    7.78E+01    1.04E-05    1.04E-05    0.00E+00    2.76E+01    1.12E+00
                                6.17E-02    1.03E-05    7.95E+01    1.03E-05    1.03E-05    0.00E+00    2.81E+01    8.91E-01
                                6.33E-02    1.04E-05    8.11E+01    1.04E-05    1.04E-05    0.00E+00    2.86E+01    8.51E-01
                                6.44E-02    1.05E-05    8.24E+01    1.06E-05    1.06E-05    0.00E+00    2.90E+01    7.59E-01
                                6.48E-02    1.13E-05    8.24E+01    1.13E-05    1.13E-05    0.00E+00    2.90E+01    8.71E-01
                                6.53E-02    1.13E-05    8.23E+01    1.13E-05    1.13E-05    0.00E+00    2.90E+01    1.12E+00
                                6.57E-02    1.13E-05    8.21E+01    1.13E-05    1.13E-05    0.00E+00    2.90E+01    1.00E+00
                                6.59E-02    1.03E-05    8.20E+01    1.03E-05    1.03E-05    0.00E+00    2.90E+01    8.51E-01];
%-------------------------------------------------------------------------------------------------------------------------------------
syms  CSO2_gas(t) CCO2_gas(t) S_total(t) C_total(t) Ca2_total(t) CCaCO3(t) CCaSO3(t) CH(t) V_Headspace F_rate CSO2_in CCO2_in R T HSO2 kga DCa2 DSO2 DHSO3 DSO32 kLa_SO2 kLa_CO2 HCO2 DCO2 DHCO3 DCO32 KSPCaSO3 BETCaSO3 ktot BETCaCO3 MWCaCO3 Kad KCO2 KHCO3 KSO2 KHSO3 Kw CCa2
eqn1 = diff(CSO2_gas(t),t)== (F_rate/ V_Headspace) * ( CSO2_in -  CSO2_gas(t)) - (( CSO2_gas(t) * R * T - HSO2 * ((S_total(t) * CH(t)^2)/(CH(t)^2 + KSO2 * CH(t) + KSO2 * KHSO3))) / ((1/kga) + 1/( kLa_SO2 * ((( DSO2 * (HSO2 * CSO2_gas(t) * R * T - ((S_total(t) * CH(t)^2)/(CH(t)^2 + KSO2 * CH(t) + KSO2 * KHSO3))) + DHSO3 * (((KSO2 * HSO2 *  CSO2_gas(t) * R * T/CH(t)) - ((S_total(t) * KSO2 * CH(t))/(CH(t)^2 + KSO2 * CH(t) + KSO2 * KHSO3)))) + DSO32 * (((( KSO2 * KHSO3 * HSO2 * CSO2_gas(t) * R * T/CH(t)^2) - ((S_total(t) * KSO2 * KHSO3)/(CH(t)^2 + KSO2 * CH(t) + KSO2 * KHSO3)))))) / ( DSO2 * ( HSO2 * CSO2_gas(t) * R * T - ((S_total(t) * CH(t)^2)/(CH(t)^2 + KSO2 * CH(t) + KSO2 * KHSO3)))))))));
eqn2 = diff(CCO2_gas(t) ,t)== (F_rate/ V_Headspace) * ( CCO2_in -  CCO2_gas(t)) - ( kLa_CO2 * ((((DCO2 * (HCO2 * CCO2_gas(t) * R * T - ((C_total(t) * CH(t)^2)/(CH(t)^2 + KCO2 * CH(t) + KCO2 * KHCO3))) + DHCO3 * (((KCO2 * HCO2 *  CCO2_gas(t) * R * T/CH(t)) - ((C_total (t) * KCO2 * CH(t))/(CH(t)^2 + KCO2 * CH(t) + KCO2*KHCO3)))) + DCO32 * ((((KCO2 * KHCO3 * HCO2 * CCO2_gas(t) * R * T/CH(t)^2) - ((C_total(t) * KCO2 * KHCO3)/(CH(t)^2 + KCO2 * CH(t) + KCO2 * KHCO3)))))) / (DCO2 * (HCO2 * CCO2_gas(t) * R * T - ((C_total (t) * CH(t)^2)/(CH(t)^2 + KCO2 * CH(t) + KCO2 * KHCO3))))))) * (CCO2_gas(t) * R * T / HCO2 - (( C_total(t) * CH(t)^2)/(CH(t)^2 + KCO2 * CH(t) + KCO2 * KHCO3))));
eqn3 = diff(S_total(t),t)== (( CSO2_gas(t) * R * T - HSO2 * ((S_total(t) * CH(t)^2)/(CH(t)^2 + KSO2 * CH(t) + KSO2 * KHSO3))) / ((1/kga) + 1/(kLa_SO2 * ((( DSO2 * (HSO2 * CSO2_gas(t) * R * T - ((S_total(t) * CH(t)^2)/(CH(t)^2 + KSO2 * CH(t) + KSO2 * KHSO3))) + DHSO3 * ((( KSO2 * HSO2 *  CSO2_gas(t) * R * T/CH(t)) - ((S_total(t) * KSO2 * CH(t))/(CH(t)^2 + KSO2 * CH(t) + KSO2 * KHSO3)))) + DSO32 * (((( KSO2 * KHSO3 * HSO2 * CSO2_gas(t) * R * T/CH(t)^2) - ((S_total(t) * KSO2 * KHSO3)/(CH(t)^2 + KSO2 * CH(t) + KSO2 * KHSO3)))))) / ( DSO2 * ( HSO2 * CSO2_gas(t) * R * T - (( S_total(t) * CH(t)^2)/(CH(t)^2 + KSO2 * CH(t) + KSO2 * KHSO3)))))))) - ( 0.162 * exp(-5153 / T) * ((( CCa2 * ((S_total(t) * KSO2 * KHSO3)/(CH(t)^2 + KSO2 * CH(t) + KSO2 * KHSO3))) / KSPCaSO3) - 1)^2 * (BETCaSO3/((CCa2 * ((((( KSO2 * KHSO3 * HSO2 * CSO2_gas(t) * R * T/CH(t)^2) - ((S_total(t) * KSO2 * KHSO3)/(CH(t)^2 + KSO2 * CH(t) + KSO2 * KHSO3)))))) / (DSO2 * (HSO2 * CSO2_gas(t) * R * T - ((S_total(t) * CH(t)^2)/(CH(t)^2 + KSO2 * CH(t) + KSO2 * KHSO3))))))))/ KSPCaSO3);
eqn4 = diff(C_total(t),t)== ( kLa_CO2 * ((((DCO2 * ( HCO2 * CCO2_gas(t) * R * T - ((C_total(t) * CH(t)^2) / (CH(t)^2 + KCO2 * CH(t) + KCO2 * KHCO3))) + DHCO3 * (((KCO2 * HCO2 *  CCO2_gas(t) * R * T / CH(t)) - ((C_total (t) * KCO2 * CH(t))/(CH(t)^2 + KCO2 * CH(t) + KCO2 * KHCO3)))) + DCO32 * ((((KCO2 * KHCO3 * HCO2 * CCO2_gas(t) * R * T / CH(t)^2) - (( C_total(t) * KCO2 * KHCO3)/(CH(t)^2 + KCO2 * CH(t) + KCO2 * KHCO3)))))) / ( DCO2 * (HCO2 * CCO2_gas(t) * R * T - ((C_total (t) * CH(t)^2)/(CH(t)^2 + KCO2 * CH(t) + KCO2 * KHCO3))))))) * (CCO2_gas(t) * R * T / HCO2 - ((C_total(t) * CH(t)^2)/(CH(t)^2 + KCO2 * CH(t) + KCO2 * KHCO3))))  + (ktot * BETCaCO3 * MWCaCO3 * CCaCO3(t) * CH(t) *(1 - (Kad * CH(t))/(1 + Kad * CH(t))));
eqn5 = diff(Ca2_total(t),t)== (-1) * ( ktot * BETCaCO3 * MWCaCO3 * CCaCO3(t) * CH(t) * (1 - (Kad * CH(t))/(1 + Kad * CH(t)))) - (0.162 * exp(-5153 / T) * ((( CCa2 * ((S_total(t) * KSO2 * KHSO3)/(CH(t)^2 + KSO2 * CH(t) + KSO2 * KHSO3))) / KSPCaSO3) - 1)^2 * (BETCaSO3/((CCa2 * (((((KSO2 * KHSO3 * HSO2 * CSO2_gas(t) * R * T/CH(t)^2) - ((S_total(t) * KSO2 * KHSO3)/(CH(t)^2 + KSO2 * CH(t) + KSO2*KHSO3)))))) / (DSO2 * (HSO2 * CSO2_gas(t) * R * T - ((S_total(t) * CH(t)^2)/(CH(t)^2 + KSO2 * CH(t) + KSO2 * KHSO3))))))))/ KSPCaSO3;
eqn6 = diff(CCaCO3(t),t)== (-1) * (ktot * BETCaCO3 * MWCaCO3 * CCaCO3(t) * CH(t) * (1 - (Kad * CH(t)) / (1 + Kad * CH(t))));
eqn7 = diff(CCaSO3(t),t)== 0.162 * exp(-5153 / T) * ((((CCa2 * ((S_total(t) * KSO2 * KHSO3)/(CH(t)^2 + KSO2 * CH(t) + KSO2 * KHSO3))) / KSPCaSO3) - 1)^2 * (BETCaSO3 / ((CCa2 * (((((KSO2 * KHSO3 * HSO2 * CSO2_gas(t) * R * T/CH(t)^2) - ((S_total(t) * KSO2 * KHSO3) / (CH(t)^2 + KSO2 * CH(t) + KSO2 * KHSO3)))))) / (DSO2 * (HSO2 * CSO2_gas(t) * R * T - ((S_total(t) * CH(t)^2)/(CH(t)^2 + KSO2 * CH(t) + KSO2 * KHSO3))))))))/ KSPCaSO3;
eqn8 = diff(CH(t),t)== CH(t) + 2 * CCa2 - ((S_total(t) * KSO2 * CH(t)) / (CH(t)^2 + KSO2 * CH(t) + KSO2 * KHSO3)) -  2 * ((S_total(t) * KSO2 * KHSO3)/(CH(t)^2 + KSO2 * CH(t) + KSO2 * KHSO3)) - ((C_total(t) * KCO2 * CH(t))/(CH(t)^2 + KCO2 * CH(t) + KCO2 * KHCO3)) - 2 * ((C_total(t) * KCO2 * KHCO3)/(CH(t)^2 + KCO2 * CH(t) + KCO2 * KHCO3))- Kw / CH(t);
eqns = [eqn1 eqn2 eqn3 eqn4 eqn5 eqn6 eqn7 eqn8];
vars = [CSO2_gas(t); CCO2_gas(t); S_total(t); C_total(t); Ca2_total(t); CCaCO3(t); CCaSO3(t); CH(t)];
origVars = length(vars);
M = incidenceMatrix(eqns, vars);
[eqns, vars] = reduceDifferentialOrder(eqns, vars);
isLowIndexDAE(eqns,vars);
f = daeFunction(eqns,vars, V_Headspace, F_rate, CSO2_in, CCO2_in, CCa2, R, T, HSO2, kga, DCa2, DSO2, DHSO3, DSO32, kLa_SO2, kLa_CO2, HCO2, DCO2, DHCO3, DCO32, KSPCaSO3, BETCaSO3, ktot, BETCaCO3, MWCaCO3, Kad, KCO2, KHCO3, KSO2, KHSO3, Kw);
V_Headspace = 1.5e-6;       
F_rate      = 1.66667e-5;    
CSO2_in     = 6.51332e-2;   
CCO2_in     = 0;            
CCa2        = 10;           
R           = 8.314;        
T           = 323.15;       
HSO2        = 149;          
kga         = 4.14e-6;      
DCa2        = 1.39e-9;      
DSO2        = 2.89e-9;
DHSO3       = 2.89e-9; 
DSO32       = 2.89e-9;
kLa_SO2     = 8.4e-4;       
kLa_CO2     = 9.598e-4;     
HCO2        = 5.15e3;       
DCO2        = 3.53e-9;
DHCO3       = 2.89e-9;
DCO32       = 2.89e-9;
KSPCaSO3    = 1.07e-7;      
BETCaSO3    = 10;           
ktot        = 8.825e-3;     
BETCaCO3    = 12.54;        
MWCaCO3     = 100.0869;     
Kad         = 0.84;         
KCO2        = 1.7e-3;       
KHCO3       = 6.55e-8;      
KSO2        = 6.24;         
KHSO3       = 5.68e-5;      
Kw          = 5.3e-8;      
F = @(t, Y, YP) f(t, Y, YP, V_Headspace, F_rate, CSO2_in, CCO2_in, CCa2, R, T, HSO2, kga, DCa2, DSO2, DHSO3, DSO32, kLa_SO2, kLa_CO2, HCO2, DCO2, DHSO3, DSO32, KSPCaSO3, BETCaSO3, ktot, BETCaCO3, MWCaCO3, Kad, KCO2, KHCO3, KSO2, KHSO3, Kw);
vars;
y0est = [0; 0; 3.3351e-1; 4.879e-6; 4.879e-6; 4.862e+1; 0; 7.413e-6];
yp0est = zeros(8,1);
opt = odeset('RelTol', 10.0^(-7), 'AbsTol' , 10.0^(-7));
[y0, yp0] = decic(F, 0, y0est, [], yp0est, [], opt);
[tSol,ySol] = ode15i(F, [600, 183000], y0, yp0, opt);
for k = 1:origVars
  S{k} = char(vars(k));
end
figure(1)
plot(tSol,ySol)
hold on
plot(ExperimentalTime, ExperimentalConcentrations, 'p')
legend('Exp-SO_{2-gas}','Exp-CO_{2-gas}','Exp-S_{total}','Exp-C_{total}','Exp-Ca^{2+}_{total}','Exp-CaCO_{3}','Exp-CaSO_{3}.^{1}/_{2}H_{2}O','Exp-H^{+}', 'Mod-SO_{2-gas}','Mod-CO_{2-gas}','Mod-S_{total}','Mod-C_{total}','Mod-Ca^{2+}_{total}','Mod-CaCO_{3}','Mod-CaSO_{3}.^{1}/_{2}H_{2}O','Mod-H^{+}', 'Location','best')
hold off
end
%eqns0 = [4.14e-6;  8.4e-4;  8.825e-3;  0.84; 0; 0; 3.3351e-1; 4.879e-6; 4.879e-6; 4.862e+1; 0; 7.413e-6];
%B = lsqcurvefit(@ode15ifun, eqns0, ExperimentalTime, ExperimentalConcentrations, zeros(size(eqns0)), inf(size(eqns0)))
%fprintf(1, '\n\tParameters:\n\t\tkga = %11.3E\n\t\tkLa = %11.3E\n\t\tktot = %11.3E\n\t\tKad = %11.3E\n\tInitial Conditions:\n\t\tCSO2gas = %11.3E\n\t\tCCO2gas = %11.3E\n\t\tStotal = %11.3E\n\t\tCtotal = %11.3E\n\t\tCatotal = %11.3E\n\t\tCCaCO3 = %11.3E\n\t\tCCaSO3 = %11.3E\n\t\tCH = %11.3E\n', B);
%F = @(t, Y, YP) f(t, B, BP, V_Headspace, F_rate, CSO2_in, CCO2_in, CCa2, R, T, HSO2, B(1), DCa2, DSO2, DHSO3, DSO32, B(2), kLa_CO2, HCO2, DCO2, DHSO3, DSO32, KSPCaSO3, BETCaSO3, B(3), BETCaCO3, MWCaCO3, B(4), KCO2, KHCO3, KSO2, KHSO3, Kw);
%[tSol,eqns] = ode15i(F, ExperimentalTime, B(5:12), bp0, opt);
%figure(2)
%plot(ExperimentalTime, ExperimentalConcentrations, 'p')
%hold on
%plot(tSol, ySol, '--')
%hold off
%grid
%legend('Exp-SO_{2-gas}','Exp-CO_{2-gas}','Exp-S_{total}','Exp-C_{total}','Exp-Ca^{2+}_{total}','Exp-CaCO_{3}','Exp-CaSO_{3}.^{1}/_{2}H_{2}O','Exp-H^{+}', 'Mod-SO_{2-gas}','Mod-CO_{2-gas}','Mod-S_{total}','Mod-C_{total}','Mod-Ca^{2+}_{total}','Mod-CaCO_{3}','Mod-CaSO_{3}.^{1}/_{2}H_{2}O','Mod-H^{+}', 'Location','best')
%end
  댓글 수: 2
Stephen23
Stephen23 2018년 10월 1일
@ Dursman Mchabe: your question is very difficult to read with so much code. Please edit your question, delete the code, upload the code by clicking the paperclip button.
Dursman Mchabe
Dursman Mchabe 2018년 10월 1일
Hi Stephen Cobeldick, Thanks for the comment. In a nutshell, I have a set of 7 ordinary differential equations and 1 algebraic which I use to calculate the concentrations of 8 species with time. I have also conducted experiments to measure the concentrations of such 8 species. The challenge I have now that the concentrations predicted by the equations are not the same as the concentrations measured during experiments. As a result, I am looking at ways the regress the 4 parameters and 8 initial conditions, such that the calculated concentrations can match the experimental concentrations. I attach a code where I attempt to use lsqcurvefitting and another code where I comment out the lsqcurvefit. The later does not give errors.

댓글을 달려면 로그인하십시오.

답변 (1개)

Steven Lord
Steven Lord 2018년 10월 1일
The relevant line of your code included in the error message:
Error in EstimationOfParametersAndInitialConditionsUsingLSQCURVEFIT3 (line 155)
B = lsqcurvefit(@ode15ifun, eqns0, ExperimentalTime, ExperimentalConcentrations)
What requirements does lsqcurvefit place on the function handle you pass into it as the first input?
"Function you want to fit, specified as a function handle or the name of a function. fun is a function that takes two inputs: a vector or matrix x, and a vector or matrix xdata. fun returns a vector or matrix F, the objective function evaluated at x and xdata."
Does your function ode15ifun satisfy those requirements? Here's how it was defined:
function ode15ifun
How many input arguments does ode15ifun accept? Does it accept two?
How many output arguments does ode15ifun return? Does it return one?
Even if you corrected the signature of your ode15ifun, if you're calling lsqcurvefit inside ode15ifun with ode15ifun as input you'll get into an infinite loop unless you're very careful.
ode15ifun calls lsqcurvefit which calls ode15ifun which calls lsqcurvefit which calls ode15ifun ...
  댓글 수: 2
Dursman Mchabe
Dursman Mchabe 2018년 10월 1일
Hi Steven Lord, Thanks a lot for pointing me in the right direction, I took the syntax of lsqcurvefit from:
https://www.mathworks.com/help/optim/ug/lsqcurvefit.html
as
x = lsqcurvefit(fun,x0,xdata,ydata);
However, the function I want to fit, ode15ifun, does not takes two inputs: a vector or matrix x, and a vector or matrix xdata. Also, it does not returns a vector or matrix F, the objective function evaluated at x and xdata.
I now realise my first mistake. My first correction will be to input fun as a function handle for an anonymous function, F. And remove line 3, function ode15ifun.
Dursman Mchabe
Dursman Mchabe 2018년 10월 1일
The above changes triggers this error:
Index exceeds array bounds.
Error in sym/subsref (line 859)
R_tilde = builtin('subsref',L_tilde,Idx);
Error in EstimationOfParametersAndInitialConditionsUsingLSQCURVEFIT3 (line 73)
eqn1 = diff(b(5),t)== F_rate/ V_Headspace * ( CSO2_in - b(5)) - (b(5) * R * T - HSO2 * ((b(7) * b(12)^2)/(b(12)^2 + KSO2 * b(12) + KSO2 * KHSO3)))
/ ((1/b(1)) + 1/(b(2) * (((DSO2 * (HSO2 * b(5) * R * T - ((b(7) * b(12)^2)/(b(12)^2 + KSO2 * b(12) + KSO2 * KHSO3))) + DHSO3 * (((KSO2 * HSO2 *
b(5) * R * T / b(12)) - ((b(7) * KSO2 * b(12))/(b(12)^2 + KSO2 * b(12) + KSO2 * KHSO3)))) + DSO32 * ((((KSO2 * KHSO3 * HSO2 * b(5) * R * T/b(12)^2)
- ((b(7) * KSO2 * KHSO3)/(b(12)^2 + KSO2 * b(12) + KSO2 * KHSO3)))))) / (DSO2 * (HSO2 * b(5) * R * T - ((b(7) * b(12)^2)/(b(12)^2 + KSO2 * b(12) +
KSO2 * KHSO3))))))));

댓글을 달려면 로그인하십시오.

카테고리

Help CenterFile Exchange에서 Conversion Between Symbolic and Numeric에 대해 자세히 알아보기

Community Treasure Hunt

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

Start Hunting!

Translated by