所得出的结果等于初始值,这个问题折磨我很久了 不知道怎么解决
求助大虾帮忙
愿意奉上BB
QQ403188374
Optimization terminated: norm of the current step is less
than OPTIONS.TolX.
global theta yexp S HI
k0 = [0.4 2.5 3.6 1000 1000 50 200 2500 5000 200 5000 10 90 ]; % 参数初值
lb = [0 0 0 0 0 0 0 0 0 0 0 0 0 ]; % 参数下限
ub = [ +inf +inf +inf +inf +inf +inf +inf +inf +inf +inf +inf +inf +inf ]; % 参数上限
x0=[0.000434037 0.000696735 0.001125753 0.00034319 0.000297428 0.000414565 6.96812E-06 2.79345E-07 0];
tspan=[0 1];
Kinetics =...
[1 0.00047 0.000722807 0.001132813 0.000320087 0.000290179 0.000399206 4.86E-06 9.52E-08 8.63E-06
2 0.000437 0.000721053 0.001149219 0.000331224 0.000289286 0.000407937 4.51E-06 1.56E-07 9.24E-06
3 0.000478 0.000712281 0.001141406 0.000322478 0.000283036 0.000400794 5.21E-06 1.84E-07 7.41E-06
4 0.000463 0.000714912 0.001182031 0.00032137 0.000286607 0.000410635 4.98E-06 1.05E-07 8.28E-06
5 0.000483 0.000705263 0.001135625 0.00033344 0.000296429 0.000400794 3.52E-06 1.18E-07 1.17E-05
6 0.000462 0.000776316 0.001135156 0.000308688 0.000283929 0.00036746 4.03E-06 1.63E-07 1.03E-05
7 0.000491 0.000716667 0.001105469 0.000336181 0.000284821 0.000361111 3.65E-06 1.44E-07 1.13E-05
8 0.000499 0.000819298 0.001078125 0.000325277 0.000283036 0.000352381 3.36E-06 1.49E-07 1.20E-05
9 0.000495 0.000767544 0.001126563 0.000336589 0.000290179 0.000331746 4.00E-06 1.28E-07 1.05E-05
10 0.000547 0.000798246 0.001145313 0.000315831 0.000289286 0.000365873 4.00E-06 1.12E-07 1.05E-05
11 0.000509 0.000729825 0.001130313 0.000323761 0.000277679 0.000387302 3.56E-06 9.84E-08 1.17E-05
12 0.000473 0.000722807 0.001157031 0.000324257 0.000288393 0.00036619 4.91E-06 8.74E-08 8.36E-06
13 0.0005 0.000813158 0.001153125 0.000334461 0.000283036 0.000373016 1.67E-06 5.48E-08 1.64E-05
14 0.000571 0.000823684 0.001151563 0.000192187 0.000253571 0.000343651 2.77E-06 7.74E-08 1.36E-05
15 0.000546 0.000718421 0.001185938 0.000333032 0.000282143 0.000394444 3.62E-06 9.64E-08 1.14E-05
];
yexp=Kinetics(:,2:10);
theta=[259.7663376 389.6495064 519.5326752 649.415844 259.7663376 389.6495064 519.5326752 649.415844 259.7663376 389.6495064 519.5326752 649.415844 259.7663376 389.6495064 519.5326752]';
HI=[0.008232711 0.008232711 0.008232711 0.008232711 0.013572204 0.013572204 0.013572204 0.013572204 0.016198704 0.016198704 0.016198704 0.016198704 0.021367521 0.021367521 0.021367521]';
S=[1.525111607 2.033482143 2.541852679 3.55859375 2.055803571 1.541852679 3.59765625 2.569754464 2.583705357 3.6171875 1.550223214 2.066964286 3.65625 2.611607143 2.089285714]';
[k,resnorm,residual,exitflag,output,lambda,jacobian] = ...
lsqnonlin(@ZObjFunc,k0,lb,ub,[],tspan,x0,theta,S,HI,yexp)
function f = ZObjFunc(k,tspan,x0,theta,S,HI,yexp) % 目标函数
for i=1:length(theta)
[t Xsim] = ode45(@Zhydrosulf,tspan,x0,[],k,theta(i),S(i),HI(i));
ysim(i,1) = Xsim(end,1);
ysim(i,2) = Xsim(end,2);
ysim(i,3) = Xsim(end,3);
ysim(i,4) = Xsim(end,4);
ysim(i,5) = Xsim(end,5);
ysim(i,6) = Xsim(end,6);
ysim(i,7) = Xsim(end,7);
ysim(i,8) = Xsim(end,8);
ysim(i,9) = Xsim(end,9);
end
f = [ysim(:,1)-yexp(:,1); ysim(:,2)-yexp(:,2); ysim(:,3)-yexp(:,3);ysim(:,4)-yexp(:,4); ysim(:,5)-yexp(:,5); ysim(:,6)-yexp(:,6);ysim(:,7)-yexp(:,7); ysim(:,8)-yexp(:,8); ysim(:,9)-yexp(:,9)];
function dCdt = Zhydrosulf (t,C,k,theta,S,HI) % ODE模型方程
dens = 1/(C(1)+C(2)+C(3)+C(4)+C(5)+C(6)+0.000641572+0.004526409+HI);
denom1 = (1+theta*dens*(k(6)*(C(4)+C(5)+C(6))+k(7)*0.004526409+k(8)*(C(7)+C(8))+k(9)*C(9)))*S;
denom2 =S*(1+theta*dens*(k(10)*(C(4)+C(5)+C(6))+k(11)*(C(7)+C(8))+(k(12)*HI)^2+k(13)*C(9)))^3
;
dCdt = ...
[theta*dens*theta*dens*k(1)*C(4)*HI / denom1
theta*dens*theta*dens*k(2)*C(5)*HI / denom1
theta*dens*theta*dens*k(3)*C(6)*HI / denom1
-theta*dens*theta*dens*k(1)*C(4)*HI / denom1
-theta*dens*theta*dens*k(2)*C(5)*HI / denom1
-theta*dens*theta*dens*k(3)*C(6)*HI / denom1
-theta*dens*theta*dens*k(4)*C(7)*HI / denom2
-theta*dens*theta*dens*k(5)*C(8)*HI / denom2
theta*dens*theta*dens*k(4)*C(7)*HI / denom2+theta*dens*theta*dens*k(5)*C(8)*HI / denom2];
[Last edited by billduke on 2011-4-10 at 17:11]