24小时热门版块排行榜

返回列表

飞鸿印雪jay

银虫 (小有名气)

应助: 1 (幼儿园)
金币: 340.5
散金: 60
帖子: 69
在线: 60.3小时
虫号: 2338094
注册: 2013-03-11
性别: GG
专业: 生物化工与食品化工

[求助] 遗传算法和最小二乘法结合拟合非线性方程组遇到问题求大神帮忙看看时间紧急了

function GA_LS2
clear all;close all;clc
data=[
0 0.590890215 0.431553748 0.111114174 0.164307696
5 0.477300105 0.495979271 0.139568658 0.181064681
10 0.396894324 0.536060613 0.167735717 0.247594176
15 0.335733264 0.543725871 0.182513065 0.235030697
20 0.268980912 0.543446351 0.199082168 0.263018724
25 0.218558133 0.543414959 0.218928771 0.301457754
30 0.16267805 0.540071789 0.249769619 0.363496193
35 0.128636265 0.523054713 0.264234746 0.415444552
40 0.097241151 0.483990339 0.266649889 0.440978169
45 0.073957065 0.449848122 0.270023805 0.468755311
50 0.062685596 0.431687277 0.278595408 0.504694097
60 0.040334681 0.386251269 0.292986268 0.596506937
70 0.022374728 0.320132412 0.28173439 0.638552746
80 0.015131643 0.27930472 0.279228981 0.695497244
90 0.008490777 0.227685253 0.266323912 0.768221053
100 0.003855323 0.182959369 0.254308566 0.820826006
110 0.002596683 0.147439246 0.234401738 0.875385242
120 0.000949164 0.121096727 0.224282059 0.931174132
];
p=data(:,1);    %pi X轴
Lexp=data(:,2:5); %Li Y轴

%-------------遗传算法-----------------------------------------------------

options = gaoptimset('Generations',1000,'StallGenLimit',300,...
'StallTimeLimit',50,'TolFun',1e-12,'TolCon',1e-12);
[k1,fva,reason,output,final_pop]=ga(@objfun,4,options);
options = gaoptimset('InitialPopulation',final_pop,'Generations',1000,'StallGenLimit',300,...
'StallTimeLimit',50,'TolFun',1e-12,'TolCon',1e-12);
[k2,fva,reason,output,final_pop2]=ga(@objfun,4,options);
fprintf('\n\n遗传算法的初始估计数值:\n');
fprintf('\n\t参数 A0 = %.9f',k2(1));
fprintf('\n\t参数 B0 = %.9f',k2(2));
fprintf('\n\t参数 C0 = %.9f',k2(3));
fprintf('\n\t参数 D0 = %.9f',k2(4));
fprintf('\n\t参数 E0 = %.9f',k2(5));
fprintf('\n\t参数 F0 = %.9f',k2(6));

%------------------------------最小二乘法 ----------------------------------
k0=k2;
lb=[0 0 0 0 0 0];
ub=[1 1 1 1 1 1]*1e6;
OPTIONS=optimset('MaxFunEvals',1000,'TolFun',1e-12,'Algorithm','trust-region-reflective','Display','Off');
[k3,resnorm,residual,exitflag,output,lambda,jacobian] = ...
lsqnonlin(@ObjFunc,k0,lb,ub,OPTIONS,p,Lexp);
%-----------------------------------------------------------------

%------------------结果输出与图形化----------------------------------------------
k=k3;
y=Lorentz(p,k,C);
fprintf('\n\n最终拟合结果:\n')
fprintf(' \t残差平方和= %.6e\n\n',resnorm);
fprintf('\n\t参数 a = %.9f',k(1))
fprintf('\n\t参数 b = %.9f',k(2))
fprintf('\n\t参数 c = %.9f',k(3))
fprintf('\n\t参数 d = %.9f',k(4))
fprintf('\n\t参数 c = %.9f',k(5))
fprintf('\n\t参数 d = %.9f',k(6))
n=length(p);
R2A=1-sum((Lexp(:,1)-y).^2)./sum((Lexp(:,1)-mean(Lexp(:,1))).^2);
R2B=1-sum((Lexp(:,2)-y).^2)./sum((Lexp(:,2)-mean(Lexp(:,2))).^2);
R2C=1-sum((Lexp(:,3)-y).^2)./sum((Lexp(:,3)-mean(Lexp(:,3))).^2);
R2D=1-sum((Lexp(:,4)-y).^2)./sum((Lexp(:,4)-mean(Lexp(:,4))).^2);

MSEA=1/n*sum((Lexp(:,1)-y).^2);
MSEB=1/n*sum((Lexp(:,2)-y).^2);
MSEC=1/n*sum((Lexp(:,3)-y).^2);
MSED=1/n*sum((Lexp(:,4)-y).^2);

% MAE=1/n*sum(abs(Lexp-y));
%
%
%
%
% MAS=max(abs(Lexp-y));

fprintf('\n\t决定系数 R-Square = %.9f',R2A);
fprintf('\n\t决定系数 R-Square = %.9f',R2B);
fprintf('\n\t决定系数 R-Square = %.9f',R2C);
fprintf('\n\t决定系数 R-Square = %.9f',R2D);

fprintf('\n\t均方误差 MSE = %.6f',MSEA);
fprintf('\n\t均方误差 MSE = %.6f',MSEB);
fprintf('\n\t均方误差 MSE = %.6f',MSEC);
fprintf('\n\t均方误差 MSE = %.6f',MSED);

% fprintf('\n\t平均绝对误差 MAE = %.6f',MAE);
%
%
%
%
% fprintf('\n\t最大绝对误差 MAS = %.6f',MAS);

figure
plot(p,y,'b',p,Lexp,'or'),axis([0 1.05 -0.05 1.05]),...
text(0.05,0.95,['决定系数 R-Square =' num2str(R2)]),...
text(0.05,0.85,['均方误差 MSE =' num2str(MSE)]),...
%    text(0.05,0.75,['平均绝对误差 MAE =' num2str(MAE)]),...
%    text(0.05,0.65,['最大绝对误差 MAS =' num2str(MAS)]),...
xlabel('p'),ylabel('L'),
legend('拟合结果','原始数据','Location','Best')
%-------------------------------------------------------------------------

function f = ObjFunc(k,t,x0,Lexp)                         % 目标函数
[t, Xsim] = ode45(@Lorentz,t,x0,[],k); %四阶,五级Runge-Kutta单步算法
Xsim1=Xsim(:,1);
Xsim2=Xsim(:,2);
Xsim3=Xsim(:,3);
Xsim4=Xsim(:,4);
Xsim5=Xsim(:,5);
Xsim6=Xsim(:,6);
Xsim7=Xsim(:,7);
Xsim8=Xsim(:,8);

ysim(:,1) = Xsim1(2:end);
ysim(:,2) = Xsim2(2:end);
ysim(:,3) = Xsim3(2:end);
ysim(:,4) = Xsim4(2:end);
ysim(:,5) = Xsim5(2:end);
ysim(:,6) = Xsim6(2:end);
ysim(:,7) = Xsim7(2:end);
ysim(:,8) = Xsim8(2:end);

f=Lorentz(p,k,C)-Lexp(i);

%----------------构造拟合目标函数-------------------------------------------
function fun=objfun(k)
data=[
0 0.590890215 0.431553748 0.111114174 0.164307696
5 0.477300105 0.495979271 0.139568658 0.181064681
10 0.396894324 0.536060613 0.167735717 0.247594176
15 0.335733264 0.543725871 0.182513065 0.235030697
20 0.268980912 0.543446351 0.199082168 0.263018724
25 0.218558133 0.543414959 0.218928771 0.301457754
30 0.16267805 0.540071789 0.249769619 0.363496193
35 0.128636265 0.523054713 0.264234746 0.415444552
40 0.097241151 0.483990339 0.266649889 0.440978169
45 0.073957065 0.449848122 0.270023805 0.468755311
50 0.062685596 0.431687277 0.278595408 0.504694097
60 0.040334681 0.386251269 0.292986268 0.596506937
70 0.022374728 0.320132412 0.28173439 0.638552746
80 0.015131643 0.27930472 0.279228981 0.695497244
90 0.008490777 0.227685253 0.266323912 0.768221053
100 0.003855323 0.182959369 0.254308566 0.820826006
110 0.002596683 0.147439246 0.234401738 0.875385242
120 0.000949164 0.121096727 0.224282059 0.931174132];
p=data(:,1);    %pi
C=data(:,2:5); %Li
n=length(t);
g1=k(1)+k(2)+k(3)+k(4);
if (k(1)<0||k(2)<0||k(3)<0||k(4)<0||g1<1)
fun=inf;
else
for i=1:n
FF(i)=(Lexp(i)-Lorentz(p(i),k(i),C(i)))^2;
end
fun=sum(FF);
end
%------------------------------构造待拟合函数------------------------------
function f = Lorentz(p,k,C)
dCAdt =-(k(1)+k(2)+k(3))*C(1);
dCBdt =k(1)*C(1)-(k(4)+k(5))*C(2);
dCCdt =k(2)*C(1)+k(4)*C(2)-k(6)*C(3);
dCDdt =k(3)*C(1)+k(5)*C(2)+k(6)*C(3);
f = [dCAdt;dCBdt;dCCdt;dCDdt];

回复此楼