matlab曲线优化,matlab，lsqnonlin拟合曲线，数优化，的修改 - 计算模拟 - 小木虫 - 学术 科研 互动社区...

function parafitN2

clear all

clc

%        t/s       Fe(2)     H2O2       4-Cp    dilos  / mol/L

Kinetics=[0       0.125    2.5        1        0

100      0.1      1.2       0.5       0.1000

200      0.1      0.5       0.2      0.22

400      0.1      0.15      0.1      0.24

900      0.1      0         0.1       0.25   ]*1e-3;  %%用于拟合的实验数据

B0=[4e9 0.6e9 (2.99)*1e11  (1.00)*1e9 (2.00)*1e7 (4.90)*1e9 (9.58)*1e5  400  40];

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

x0=(1e-3)*[0.125  0  2.5  0 0 0 1 0 0 0  0 0];

yexp = Kinetics  ;               % yexp: 实验数据[x1        x4        x5        x6]

% 使用函数lsqnonlin()进行参数估计

options = optimset('largescale','off','display','iter');

options=optimset(options,'tolx',1e-100);

options=optimset(options,'tolfun',1e-100);

[B,resnorm,residual,exitflag,output,lambda,jacobian] = ...

lsqnonlin(@ObjFunc7LNL,B0,lb,[],options,x0,yexp);

ci = nlparci(B,residual,jacobian);

fprintf('\n\n使用函数lsqnonlin()估计得到的参数值为:\n')

fprintf('\tk16 = %.11f\n',B(1))

fprintf('\tk17 = %.11f\n',B(2))

fprintf('\tk18 = %.11f\n',B(3))

fprintf('\tk21 = %.11f\n',B(4))

fprintf('\tk22 = %.11f\n',B(5))

fprintf('\tk23 = %.11f\n',B(6))

fprintf('\tk24 = %.11f\n',B(7))

fprintf('\tk25 = %.11f\n',B(8))

fprintf('\tk26 = %.11f\n',B(9))

fprintf('  所求残差为: %.1e\n\n',resnorm)

% ------------------------------------------------------------------

function f = ObjFunc7LNL(B,x0,yexp) %目标函数=实验值和和理论值之差

tspan = [0 100 200  400 900];

[t x] = ode23s(@KineticEqs,tspan,x0,[],B);

f1 = x(:,1) - yexp(:,2);

f2 = x(:,3) - yexp(:,3);

f3 = x(:,7) - yexp(:,4);

f4 = x(:,8) - yexp(:,5);

f = [f1 f2 f3 f4];

% ------------------------------------------------------------------

function dxdt = KineticEqs(t,x,B)

% 反应模型方程

k(1)=76; k(2)=0.01; k(3)=(1.0e+7)*2.7;k(4)=(1.0e+5)*1.58; k(5)=1.0e+10;k(6)=3.2*(1.0e+8); k(7)=1.2*(1.0e+6); k(8)=3.1*(1.0e+5);k(9)=(1.0e+7) ;k(10)=5*(1.0e+7); k(11)=4.2*(1.0e+9);  k(12)=8.3*(1.0e+5); k(13)=(1.0e+10); k(14)=(1.0e+10); k(15)=9.7*(1.0e+7);

r1=k(1)*x(1)*x(3);r2 = k(2)*x(2)*x(3); r3 = k(3)*x(4)*x(3);   r4 = k(4)*x(5);

r5 = k(5)*(1e-3)*x(6);r6 = k(6)*x(4)*x(1);r7 = k(7)*x(5)*x(1)*1e-3;  r8 = k(8)*x(5)*x(2)*1e-3;

r9 = k(9)*x(6)*x(1)*1e-6; r10 = k(10)*x(6)*x(2);r11 = k(11)*x(4)*x(4);r12 = k(12)*x(5)*x(5);

r13 = k(13)*x(4)*x(5);r14 = k(14)*x(4)*x(6);r15 = k(15)*x(5)*x(6);

r16 = B(1)*x(7)*x(4);

r17 = B(2)*x(7)*x(4);

r18 = B(3)*x(9)*x(4);

r21 = B(4)*x(8)*x(4);

r22 = B(5)*x(8)*x(2)*x(2);

r23 = B(6)*x(11)*x(4);

r24 = B(7)*x(12)*x(4);

r25 = B(8)*x(12)*x(2);

r26 = B(9)*x(2)*x(12)*0.3;

dx(1)=-r1+r2-r6-r7+r8-r9+r10+2*r22+r25;

%检测指标[Fe2 +]:

dx(2)= -(-r1+r2-r6-r7+r8-r9+r10+2*r22+r25)-r26;

%[Fe3+ ]

dx(3)=-r1-r2-r3+r7+r9+r11+r12+r15;% 原文“-r15”

%检测指标[h2o2]

dx(4)=r1-r3-r6-r11-r13-r14-r16-r21-r23-r24-r18;

%[·OH  ]:

dx(5)=r2-r4+r5-r7-r8-r12-r13-r15;

%[HO2·]

dx(6)=r4-r5-r9-r10-r14-r15;

%x(o2)

dx(7)=-r16-r17;

%检测指标[4v-xP]

dx(8)=r17+r18-r21-r22;

%检测指标[diols ]

dx(9)=r16-r18;

%[xlDHxD ]

dx(10)=0;

%[xlDHxDP ]

dx(11)=r21+r22-r23;

%[B Q]

dx(12)=3*r23-r24-r25-r26;

%aSS

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