用fmincon和lsqnonlin优化问题

identify1100

实际上做的就是合成气制甲醇动力学参数拟合问题，我编了下面的程序，得出的结果实在是不理想，哪位高手能指导我一下，在此先谢过了，我的QQ：258738020，邮箱chenpeng1100@hotmail.com。

function Kinetics_global
clear all
clc
format long
k0 = [214.7 -41835 12940 -60968 -36.81 -1000 -0.145 10560 -1.720 -1220];  %参数初值 C302
lb = [10 -100000 1000 -100000 -100 -10000 -10 1000 -10 -10000];            %参数下限
ub = [10000 -1000 100000 -1000 0 -100 0 100000 0 -100];                    %参数上限
%k0 = [9110 -5842 10480 -77872 -36.81 -1000 -0.145 10560 -1.720 -1220];
%--------------------------------------------------------------
W = 3.0901;                                                                %实验数据
%W = 2.0337;
P = 5;
R = 8.3145;
AT = 508.9;                                                                %反应平均温度
%--------------------------------------------------------------
KineticsDate_ying;
%Data09_03_31;
%Data09_03;
%Data09_31;
T = ExpData(:,1)+273.15 ;                                                  %实验温度                                            
yCO2_in = ExpData(:,5);
yH2 = ExpData(:,6);                                                        %进口组分
yCO = ExpData(:,7);
yCO2 = ExpData(:,8);
ym = ExpData(:,9);
yH2O = ExpData(:,10);
Nin = ExpData(:,1);                                                        %进口流量
Nout = Nin./(1+2.*ExpData(:,9));                                           %出口流量 Nout=Nin/(1+2*ym,out)
                                                            
% 使用函数fmincon()进行参数估计
[k,fval,flag] = fmincon(@Gk_fmincon,k0,[],[],[],[],lb,ub,[],[],ExpData,T,yCO,P,yCO2,yH2,ym,yH2O,R,AT,Nin,Nout,W,yCO2_in);
fprintf('\n使用函数fmincon()估计得到的参数值为:\n')
fprintf('\tk1 = %.4f\n',k(1))
fprintf('\tk2 = %.4f\n',k(2))
fprintf('\tk3 = %.4f\n',k(3))
fprintf('\tk4 = %.4f\n',k(4))
fprintf('\tk5 = %.4f\n',k(5))
fprintf('\tk6 = %.4f\n',k(6))
fprintf('\tk7 = %.4f\n',k(7))
fprintf('\tk8 = %.4f\n',k(8))
fprintf('\tk9 = %.4f\n',k(9))
fprintf('\tk10 = %.4f\n',k(10))
fprintf('  The sum of the squares is: %.1e\n\n',fval)
k_fmincon = k;

% 使用函数lsqnonlin()进行参数估计
[k,resnorm,residual,exitflag,output,lambda,jacobian] = ...
    lsqnonlin(@Gk_lsqnonlin,k0,lb,ub,[],ExpData,T,yCO,P,yCO2,yH2,ym,yH2O,R,AT,Nin,Nout,W,yCO2_in);      
ci = nlparci(k,residual,jacobian);
fprintf('\n\n使用函数lsqnonlin()估计得到的参数值为:\n')
fprintf('\tk1 = %.4f ± %.4f\n',k(1),ci(1,2)-k(1))
fprintf('\tk2 = %.4f ± %.4f\n',k(2),ci(2,2)-k(2))
fprintf('\tk3 = %.4f ± %.4f\n',k(3),ci(3,2)-k(3))
fprintf('\tk4 = %.4f ± %.4f\n',k(4),ci(4,2)-k(4))
fprintf('\tk5 = %.4f ± %.4f\n',k(5),ci(5,2)-k(5))
fprintf('\tk6 = %.4f ± %.4f\n',k(6),ci(6,2)-k(6))
fprintf('\tk7 = %.4f ± %.4f\n',k(7),ci(7,2)-k(7))
fprintf('\tk8 = %.4f ± %.4f\n',k(8),ci(8,2)-k(8))
fprintf('\tk9 = %.4f ± %.4f\n',k(9),ci(9,2)-k(9))
fprintf('\tk10 = %.4f ± %.4f\n',k(10),ci(10,2)-k(10))
fprintf('  The sum of the squares is: %.1e\n\n',resnorm)

% 以函数fmincon()估计得到的结果为初值，使用函数lsqnonlin()进行参数估计
k0 = k_fmincon;
[k,resnorm,residual,exitflag,output,lambda,jacobian] = ...
    lsqnonlin(@Gk_lsqnonlin,k0,lb,ub,[],ExpData,T,yCO,P,yCO2,yH2,ym,yH2O,R,AT,Nin,Nout,W,yCO2_in);      
ci = nlparci(k,residual,jacobian);
fprintf('\n\n以fmincon()的结果为初值，使用函数lsqnonlin()估计得到的参数值为:\n')
fprintf('\tk1 = %.4f ± %.4f\n',k(1),ci(1,2)-k(1))
fprintf('\tk2 = %.4f ± %.4f\n',k(2),ci(2,2)-k(2))
fprintf('\tk3 = %.4f ± %.4f\n',k(3),ci(3,2)-k(3))
fprintf('\tk4 = %.4f ± %.4f\n',k(4),ci(4,2)-k(4))
fprintf('\tk5 = %.4f ± %.4f\n',k(5),ci(5,2)-k(5))
fprintf('\tk6 = %.4f ± %.4f\n',k(6),ci(6,2)-k(6))
fprintf('\tk7 = %.4f ± %.4f\n',k(7),ci(7,2)-k(7))
fprintf('\tk8 = %.4f ± %.4f\n',k(8),ci(8,2)-k(8))
fprintf('\tk9 = %.4f ± %.4f\n',k(9),ci(9,2)-k(9))
fprintf('\tk10 = %.4f ± %.4f\n',k(10),ci(10,2)-k(10))
fprintf('  The sum of the squares is: %.1e\n\n',resnorm)

%--------------------------------------------------------------------------
function f=Gk_lsqnonlin(k,ExpData,T,yCO,P,yCO2,yH2,ym,yH2O,R,AT,Nin,Nout,W,yCO2_in)
%----------------------计算组成逸度------------------
jCO = exp((-0.09326+189.156./T-399940./(T.^3)-181.527.*yCO./T+140.001.*(yCO.^2)./T).*P./(0.101325.*T));
jCO2 = exp((-0.343605+428.452./T-6.92177*(10^7)./(T.^3)-327.402.*yCO2./T-374.954.*(yCO2.^2)./T).*P./(0.101325.*T));
jH2 = exp((0.110785+35.3324./T-5005.74./(T.^3)-19.6109.*yH2./T-20.9799.*(yH2.^2)./T).*P./(0.101325.*T));
jm = exp((-1.49696+997.85./T-(10^8)./(T.^3)-792.109.*ym./T-803.4.*(ym.^2)./T).*P./(0.101325.*T));
jH2O = exp((-1.78527+1408.49./T-1.83959*(10^8)./(T.^3)-3648.32.*yH2O./T-3116.5.*(yH2O.^2)./T).*P./(0.101325.*T));
fCO = yCO.*P.*jCO;
fCO2 = yCO2.*P.*jCO2;
fH2 = yH2.*P.*jH2;
fm = ym.*P.*jm;
fH2O = yH2O.*P.*jH2O;
Kf1 = exp(13.1652+9203.26./T-5.92839.*log(T)-0.352404*0.01.*T+0.102264*(10^(-4)).*(T.^2)-0.769446*(10^(-8)).*(T.^3)+0.238583*(10^(-11)).*(T.^4)).*(0.101325^(-2));
Kf2 = exp(1.6654+4553.34./T-2.72613.*log(T)-1.422914*0.01.*T+0.172060*(10^(-4)).*(T.^2)-1.106294*(10^(-8)).*(T.^3)+0.319698*(10^(-11)).*(T.^4)).*(0.101325^(-2));;
B1 = fm./(fCO.*(fH2.^2).*Kf1);
B2 = fm.*fH2O./(fCO2.*(fH2.^3).*Kf2);
%----------------------计算组成平衡常数--------------
K1 = k(1).*exp(k(2)./(R.*T));
K2 = k(3).*exp(k(4)./(R.*T));
Kco = exp(k(5)+k(6).*(1./T-1/AT));
Kco2 = exp(k(7)+k(8).*(1./T-1/AT));
KH2 = exp(k(9)+k(10).*(1./T-1/AT));
%--------------------出口浓度计算值----------------------
yCO2_c=(Nin.*yCO2_in-W.*K2.*fCO2.*(fH2.^3).*(1-B2)./((1+Kco.*fCO+Kco2.*fCO2+KH2.*fH2).^4))./Nout;
ym_c=W.*(K2.*fCO2.*(fH2.^3).*(1-B2)./((1+Kco.*fCO+Kco2.*fCO2+KH2.*fH2).^4)+K1.*fCO.*(fH2.^2).*(1-B1)./((1+Kco.*fCO+Kco2.*fCO2+KH2.*fH2).^3))./Nout;
f1 = ExpData(:,8)-yCO2_c;
f2 = ExpData(:,9)-ym_c;
f = [f1;f2];
%-----------------------------------------------------------------------
function f=Gk_fmincon(k,ExpData,T,yCO,P,yCO2,yH2,ym,yH2O,R,AT,Nin,Nout,W,yCO2_in)
%----------------------计算组成逸度------------------
jCO = exp((-0.09326+189.156./T-399940./(T.^3)-181.527.*yCO./T+140.001.*(yCO.^2)./T).*P./(0.101325.*T));
jCO2 = exp((-0.343605+428.452./T-6.92177*(10^7)./(T.^3)-327.402.*yCO2./T-374.954.*(yCO2.^2)./T).*P./(0.101325.*T));
jH2 = exp((0.110785+35.3324./T-5005.74./(T.^3)-19.6109.*yH2./T-20.9799.*(yH2.^2)./T).*P./(0.101325.*T));
jm = exp((-1.49696+997.85./T-(10^8)./(T.^3)-792.109.*ym./T-803.4.*(ym.^2)./T).*P./(0.101325.*T));
jH2O = exp((-1.78527+1408.49./T-1.83959*(10^8)./(T.^3)-3648.32.*yH2O./T-3116.5.*(yH2O.^2)./T).*P./(0.101325.*T));
fCO = yCO.*P.*jCO;
fCO2 = yCO2.*P.*jCO2;
fH2 = yH2.*P.*jH2;
fm = ym.*P.*jm;
fH2O = yH2O.*P.*jH2O;
Kf1 = exp(13.1652+9203.26./T-5.92839.*log(T)-0.352404*0.01.*T+0.102264*(10^(-4)).*(T.^2)-0.769446*(10^(-8)).*(T.^3)+0.238583*(10^(-11)).*(T.^4)).*(0.101325^(-2));
Kf2 = exp(1.6654+4553.34./T-2.72613.*log(T)-1.422914*0.01.*T+0.172060*(10^(-4)).*(T.^2)-1.106294*(10^(-8)).*(T.^3)+0.319698*(10^(-11)).*(T.^4)).*(0.101325^(-2));;
B1 = fm./(fCO.*(fH2.^2).*Kf1);
B2 = fm.*fH2O./(fCO2.*(fH2.^3).*Kf2);
%----------------------计算组成平衡常数--------------
K1 = k(1).*exp(k(2)./(R.*T));
K2 = k(3).*exp(k(4)./(R.*T));
Kco = exp(k(5)+k(6).*(1./T-1/AT));
Kco2 = exp(k(7)+k(8).*(1./T-1/AT));
KH2 = exp(k(9)+k(10).*(1./T-1/AT));
%--------------------出口浓度计算值----------------------
yCO2_c =( Nin.*yCO2_in-W.*K2.*fCO2.*(fH2.^3).*(1-B2)./((1+Kco.*fCO+Kco2.*fCO2+KH2.*fH2).^4))./Nout;
ym_c = W.*(K2.*fCO2.*(fH2.^3).*(1-B2)./((1+Kco.*fCO+Kco2.*fCO2+KH2.*fH2).^4)+K1.*fCO.*(fH2.^2).*(1-B1)./((1+Kco.*fCO+Kco2.*fCO2+KH2.*fH2).^3))./Nout;
f1 = ExpData(:,8)-yCO2_c;
f2 = ExpData(:,9)-ym_c;
f = sum(f1.^2)+sum(f2.^2);

identify1100

哎。。。期待好心人出现。。

messenger

弄这么大段程序，也没时间看呀。只贴程序，不说明，谁也不知道你要做什么。结果哪不理想呀？

identify1100

3# messenger
首先谢谢回复。
恩 我的失误。。。程序里调用函数中的式子也就是一些计算公式，应该是没做么问题的。
我的问题就是说要求得参数K使这个目标函数：
f = sum((ExpData(:,8)-yCO2_c).^2)+sum((ExpData(:,9)-ym_c).^2);值最小，结果用fmincon和lsqnonlin得出的结果不同，差别很大。

shamohu

如果可能，用1stOpt对比计算一下。