运行 一直提示 Not enough input arguments 以下是我的代码

1 Call lsqnonlin

tic % time check

clear all;

%Experimental data

tspan = [0,600,900,1200,1500,1800,2100,2400,2700,3000]';

% Concentration of Mg2+  mol/L

y_1 = [0,0.186,0.271,0.36,0.384,0.311,0.262,0.257,0.252,0.245]';

% PH

y_2 = [10.35,8.64,8.45,8.35,8.21,8.09,8,7.83,7.65,7.42]';

data = [tspan y_1 y_2];

yexp = [data(:,2) data(:,3)];

% initial parameter identification

kn = 10^2.6; %Crystallization kinetics

kg = 8.06e-2;

n = 6.5 ;

g = 6.5 ;

X_1=1; % mass transfer coefficient enhancement factor

X_2=1;

X_3=0.05;

teta_0 = [kn kg n g X_1 X_2 X_3]; %just collect in 1 vector sampling instants t and measured B

%initial values for ODE

x0 = [0,0,3.3e-4,3.03e-11,0,0,1.65e-4,100,0,0,0,0,0];

lb = [0 0 0 0 0 0 0 0 0 0 0 0 0];

ub = [1 1 1 1 1 1 1 101 1 1 1 1 1];

% Call the optimizer:

[teta_opt,resnorm,residual,exitflag,output,lambda,jacobian] ...

= lsqnonlin(@Objfun_ffMgOH2_CO2_H2O_with_precipitation,teta_0,lb,ub,[],tspan,x0,yexp);

ci = nlparci(teta_opt,residual,jacobian);

kn = teta_opt(1);

kg = teta_opt(2);

n  = teta_opt(3);

g  = teta_opt(4);

X_1 = teta_opt(5);

X_2 = teta_opt(6);

X_3 = teta_opt(7);

%Optimized solution

tsim = [0,3000];

[t,x] = ode15s(@ode_ffMgOH2_CO2_H2O_with_precipitation,tsim,x0,[],kn,kg,n,g,X_1,X_2,X_3);

figure(1)

plot(tspan,y_1,'o',tsim,x(:,7)+x(:,9)) %plot the data vs solution

title('Parameter fitting of Mg2+ in the bulk,mol/L');

figure(2)

PH=-log10(x(:,4));

plot(tspan,y_2,'o',tsim,PH) %plot the data vs solution

title('Parameter fitting of PH');

toc %time check

2 Call Objfun

%Call 'Objfun' function:

function lsq = Objfun_ffMgOH2_CO2_H2O_with_precipitation(teta,tspan,x0,yexp)

[~,x] = ode15s(@ode_ffMgOH2_CO2_H2O_with_precipitation,tspan,x0,[],teta);

%Concentration of Mg2+

y_cal(:,1) = x(:,7)+x(:,9);

%PH

y_cal(:,2) = -log10(x(:,4));

lsq = [y_cal(:,1)-(yexp(:,1)) (y_cal(:,2)-yexp(:,2))];

end

3 Call ode_function

% the reaction of Calc hydroxide[Mg(OH)2] and Carbon dioxide

function dx=ode_ffMgOH2_CO2_H2O_with_precipitation(t,x,kn,kg,n,g,X_1,X_2,X_3)

dx = zeros(13,1);

% x(1) is the concentration of liquid CO2(l)

% initial value = 0 mol/L

% x(2) is the concentration of gas CO2(g)

% initial value = 0 mol/L

% x(3) is the concentration of Hydroxy [OH-] in bulk

% initial value = 3.3e-4 mol/L

% x(4) is the concentration of Hydrogen [H+] in bulk

% initial value = 3.03e-11 mol/L

% x(5) is the concentration of Bicarbonate [HCO3-]

% initial value = 0 mol/L

% x(6) is the concentration of carbonate [CO3 2-]

% initial value = 0 mol/L

% x(7) is the concentration of Magnesium [Mg 2+]

% initial value = 1.65e-4  mol/L

% x(8) is the mass of Mg(OH)2

% initial value = 100  g

% x(9) is the concentraiton of liquid MgCO3

% initial value = 0  mol/L

% x(10) is o monent which is number density of crystal n(L)

% initial value = 0

% x(11) is 1 monent which is number density of crystal n(L)

% initial value = 0

% x(12) is 2 monent which is number density of crystal n(L)

% initial value = 0

% x(13) is 3 monent which is number density of crystal n(L)

% initial value = 0

Qg = 1.67e-2;        % L/s          1L/min

N  = 650;

vl = 2.826;          % vl is the dispersion volum                  L

vg = 2.047*10^(-5)*N^(0.579)*Qg^0.63*vl;

% vg is the gas volum in the gas-liquid mixture m3

%agitating speed rpm

P = 1;               % P is the atmospheric pressure                atm

pi = 3.14;

kla = 0.0304;        % kla is the mass transfer coefficient of CO2  s^(-1)

%(9.6218e-5)*N^0.4316*Qg^0.3840*c^0.1117*E;

kH = 0.035;          % kh is the Henry's constant               mol/(L*atm)

M_CO2 = 44;          % M_CO2 is the co2 molar mass              g/mol

M_MgOH = 58;         % molar mass of Mg(OH)2               g/mol

M_MgCO3 = 84;        % molar mass of M_MgCO3    g/mol

dd_CO2 = 1.8;        % dd_CO2 is the density of CO2

dd_MgOH = 2.34e3;    % dd_Mg(OH)2 is the density of Mg(OH)2          g/L

dd_MgCO3 = 2.96e3;   % dd_Mg(OH)2 is the density of MgCO3          g/L

dd_mol_MgCO3 = dd_MgCO3/M_MgCO3; % mol/L

ddslurry = 1e3;      % density of slurry g/L                     % ddCO2 = 1960; ddCO2 is the density of CO2

c = 1;               % volume fraction of CO2 in the inlet gas

CO2in = dd_CO2/M_CO2*c;% CO2in is the concentration of CO2 entering the reactor mol/L

% CO2in = ddCO2/MCO2*0.15 =1800/44*0.05 =6.75  mol/m3

%Ksp1 = 1.8e-11;      % Ksp[Mg(OH)2]          mol2/L2

Ksp2 = 3.5e-8;       % Ksp[MgCO3]            mol2/L2

k11 = 8.4e3;         % L/(mol*s) k11 = 8.4*10^3             L/mol*s

k12 = 2e-4;          % s^(-1)     k12 = 2.0*10^(-4)           s^(-1)

k21 = 6e9;           % L/(mol*s) k21 = 6*10^(9)             L/mol*s

k22 = 1.2e6;         % s^(-1)     k22 = 1.2*10^(6)            s^(-1)

k31 = 1.4e11;        % L/(mol*s) k31 = 1.4*10^(11)          L/mol*s

k32 = 1.3e-3;        % mol/(L*s) k32 = 1.3*10^(-3)          mol/L*s

k41 = 2.4e-2;        % s^(-1)     k41 = 2.4*10^(-2)          s^(-1)

k42 = 5.7e2;         % L/(mol*s) k42 = 5.7*10^4             L/mol*s

k52 = 9e-3;          % s^(-1) k52 = 9e-3                     s^(-1)

k51 = 1.88e6;         %k52/Ksp2;      % L/(mol*s)

%D_CO2_l  = 1.91e-7;  %diffusion coefficeint of CO2(l) dm2/s

D_OH     = 5.27e-7;  %diffusion coefficeint of OH- dm2/s

D_H      = 9.31e-7;  %diffusion coefficeint of H+ dm2/s

%D_HCO3   = 1.19e-7;  %diffusion coefficeint of HCO3 dm2/s

%D_CO3    = 9.23e-8;  %diffusion coefficeint of CO3 dm2/s

D_Mg     =1.06e-9; %7.06e-8;  %diffusion coefficeint of Mg2+ dm2/s

%zk_CO2_l  = 0;

zk_OH     = -1;

zk_H      = 1;

%zk_HCO3   = -1;

%zk_CO3    = -2;

zk_Mg     = 2;

%C_CO2_l  = 0;       % concentration on the interface i

C_OH     = 3.3e-4;

C_H      = 3.03e-11;

%C_HCO3   = 0;

%C_CO3    = 0;

C_Mg     = 1.65e-4; % concentration on the interface i

%d_CO2_l   = C_CO2_l-x(1); % Ci-Cb

d_OH      = C_OH-x(3);

d_H       = C_H-x(4);

%d_HCO3    = C_HCO3-x(5);

%d_CO3     = C_CO3-x(6);

d_Mg      = C_Mg-x(7); % Ci-Cb

% d_CO2_l   = x(1)-C_CO2_l; % Cb-Ci

% d_OH      = x(3)-C_OH;

% d_H       = x(4)-C_H;

% d_HCO3    = x(5)-C_HCO3;

% d_CO3     = x(6)-C_CO3;

% d_Mg      = x(7)-C_Mg; % Cb-Ci

%A_CO2_l   = (C_CO2_l+x(1))/2; % (Ci+Cb)/2

A_OH      = (C_OH+x(3))/2;

A_H       = (C_H+x(4))/2;

%A_HCO3    = (C_HCO3+x(5))/2;

%A_CO3     = (C_CO3+x(6))/2;

A_Mg      = (C_Mg+x(7))/2;  % (Ci+Cb)/2

Ntot1 = 8.17e10;      %total number of Mg(OH)2 100g = 8.17e10

%Ntot1 = 100/(3.14*d^3*2340/6)

UL = 1.01e-1; % dynamic viscosity 1.01e-3 Pas = kg*s/m = 1.01e-1 g*s/dm

VL = 1.01e-4; % kinematic viscosity 1.01e-6 m2/s = 1.01e-4 dm2/s

eT = 14028; % total energy dissipation w/kg = m2/s3*100 = dm2/s3

dp = (6*x(8)/(Ntot1*dd_MgOH*3.14))^(1/3); %diameter of solid dm

a = (3.3*x(8)^(2/3)*(3.14*Ntot1)^(1/3)/dd_MgOH^(2/3)/vl);%specific area

if a > 0

% mass transfer %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Re = eT^(1/3)*dp^(4/3)/VL;

%Sc_CO2_l  = UL/(ddslurry*D_CO2_l);

Sc_OH     = UL/(ddslurry*D_OH);

Sc_H      = UL/(ddslurry*D_H);

%Sc_HCO3   = UL/(ddslurry*D_HCO3);

%Sc_CO3    = UL/(ddslurry*D_CO3);

Sc_Mg     = UL/(ddslurry*D_Mg);

%Sh_CO2_l  = (2^5.8+(0.61*(Re^0.58)*(Sc_CO2_l^(1/3)))^5.8)^(1/5.8);

Sh_OH     = (2^5.8+(0.61*(Re^0.58)*(Sc_OH^(1/3)))^5.8)^(1/5.8);

Sh_H      = (2^5.8+(0.61*(Re^0.58)*(Sc_H^(1/3)))^5.8)^(1/5.8);

%Sh_HCO3   = (2^5.8+(0.61*(Re^0.58)*(Sc_HCO3^(1/3)))^5.8)^(1/5.8);

%Sh_CO3    = (2^5.8+(0.61*(Re^0.58)*(Sc_CO3^(1/3)))^5.8)^(1/5.8);

Sh_Mg     = (2^5.8+(0.61*(Re^0.58)*(Sc_Mg^(1/3)))^5.8)^(1/5.8);

% Parameter fitting!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

%Ks_CO2_l  = Sh_CO2_l*D_CO2_l/dp;

Ks_OH     = Sh_OH*D_OH/dp*X_1;

Ks_H      = Sh_H*D_H/dp*X_2;

%Ks_HCO3   = Sh_HCO3*D_HCO3/dp;

%Ks_CO3    = Sh_CO3*D_CO3/dp;

Ks_Mg     = Sh_Mg*D_Mg/dp*X_3;

% K = (zk_CO2_l*Ks_CO2_l*d_CO2_l+zk_OH*Ks_OH*d_OH+zk_H*Ks_H*d_H+zk_HCO3*Ks_HCO3*d_HCO3+zk_CO3*Ks_CO3*d_CO3+zk_Mg*Ks_Mg*d_Mg)/...

%     (zk_CO2_l^2*Ks_CO2_l*A_CO2_l+zk_OH^2*Ks_OH*A_OH+zk_H^2*Ks_H*A_H+zk_HCO3^2*Ks_HCO3*A_HCO3+zk_CO3^2*Ks_CO3*A_CO3+zk_Mg^2*Ks_Mg*A_Mg);

K = (zk_OH*Ks_OH*d_OH+zk_H*Ks_H*d_H+zk_Mg*Ks_Mg*d_Mg)/...

(zk_OH^2*Ks_OH*A_OH+zk_H^2*Ks_H*A_H+zk_Mg^2*Ks_Mg*A_Mg);

%NCO2_l  = Ks_CO2_l*(d_CO2_l-A_CO2_l*zk_CO2_l*K);

NOH     = Ks_OH*(d_OH-A_OH*zk_OH*K);

NH      = Ks_H*(d_H-A_H*zk_H*K);

%NHCO3   = Ks_HCO3*(d_HCO3-A_HCO3*zk_HCO3*K);

%NCO3    = Ks_CO3*(d_CO3-A_CO3*zk_CO3*K);

NMg     = Ks_Mg*(d_Mg-A_Mg*zk_Mg*K);

% K = (zk_CO2_l*D_CO2_l*d_CO2_l+zk_OH*D_OH*d_OH+zk_H*D_H*d_H+zk_HCO3*D_HCO3*d_HCO3+zk_CO3*D_CO3*d_CO3+zk_Mg*D_Mg*d_Mg)/...

%     (zk_CO2_l^2*D_CO2_l*A_CO2_l+zk_OH^2*D_OH*A_OH+zk_H^2*D_H*A_H+zk_HCO3^2*D_HCO3*A_HCO3+zk_CO3^2*D_CO3*A_CO3+zk_Mg^2*D_Mg*A_Mg);

%

% NCO2_l  = Ks_CO2_l*d_CO2_l-Ks_CO2_l*A_CO2_l*zk_CO2_l*K;

% NOH     = Ks_OH*d_OH-Ks_OH*A_OH*zk_OH*K;

% NH      = Ks_H*d_H-Ks_H*A_H*zk_H*K;

% NHCO3   = Ks_HCO3*d_HCO3-Ks_HCO3*A_HCO3*zk_HCO3*K;

% NCO3    = Ks_CO3*d_CO3-Ks_CO3*A_CO3*zk_CO3*K;

% NMg     = Ks_Mg*d_Mg-Ks_Mg*A_Mg*zk_Mg*K;

% reaction kinetics%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

r1 = k11*x(1)*x(3)-k12*x(5);

r2 = k21*x(5)*x(3)-k22*x(6);

r3 = k31*x(3)*x(4)-k32;

r4 = k41*x(1)-k42*x(5)*x(4);

r5 = k51*x(7)*x(6)-k52;

% Crystallization kinetics%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Parameter fitting!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

% kn = 10^2.6;  % kn = 10^2.6;%mol^(-n)*dm^(3n-3)*s^-1

% kg = 8.06e-2; % kg = 8.06e-2;%mol^(-g)*dm^(3g+1)*s^-1

% n = 6.5;      % n = 4.2;

% g = 6.5;        % g = 2;

S = sqrt((x(6)+x(9))*(x(7)+x(9)))-sqrt(Ksp2);

B = kn*(S^n);

G = kg*(S^g);

%current = 2*NMg+NH-NOH

dx(1) = kla*(kH*M_CO2*P*x(2)/dd_CO2-x(1))*(vl+vg)/vl-r1-r4;%+NCO2_l*a;

dx(2) = (Qg*(CO2in-x(2))-kla*(kH*M_CO2*P*x(2)/dd_CO2-x(1))*(vl+vg))/vg;

dx(3) = -r1-r2-r3+NOH*a;

dx(4) = -r3+r4+NH*a;

dx(5) = r1-r2+r4;%+NHCO3*a;

dx(6) = r2-r5;%+NCO3*a;

dx(7) = NMg*a-r5;

dx(8) = -NMg*a*vl*M_MgOH;

dx(9) = r5-pi/6*(1e-21)*dd_mol_MgCO3*B-G/2*pi*dd_mol_MgCO3*x(12); %MgCO3(l) mol/L

dx(10) = B; % m0

dx(11) = G*x(10);%m1

dx(12) = 2*G*x(11);%m2

dx(13) = 3*G*x(12);%m3

dx = [dx(1);dx(2);dx(3);dx(4);dx(5);dx(6);dx(7);dx(8);dx(9);dx(10);dx(11);dx(12);dx(13)];

end

end