应《网络安全法》要求，自2017年10月1日起，未进行实名认证将不得使用互联网跟帖服务。为保障您的帐号能够正常使用，请尽快对帐号进行手机号验证，感谢您的理解与支持！

24小时热门版块排行榜    

小木虫论坛-学术科研互动平台 » 计算模拟区 » 计算模拟 » 遗传算法和最小二乘法结合 拟合非线性方程组遇到问题 求大神帮忙看看 时间紧急了
11/1返回列表
查看: 903  |  回复: 0
只看楼主 @他人 存档 新回复提醒 (忽略) 收藏    在APP中查看

飞鸿印雪jay

银虫 (小有名气)

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

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];
回复此楼

» 猜你喜欢
1楼 2014-11-20 20:49:48
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

智能机器人

Robot (super robot)

我们都爱小木虫

找到一些相关的精华帖子，希望有用哦~

科研从小木虫开始，人人为我，我为人人
相关版块跳转 我要订阅楼主 飞鸿印雪jay 的主题更新
11/1返回列表
普通表情 龙 兔 虎 猫
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[基金申请] 要骂人了，新模版改版就是要淡化问题凝练这种虚的东西，结果有个评委还在说凝练得不够 +5 瞬息宇宙 2026-08-31 8/400 2026-09-02 09:37 by 昵称告罄
[基金申请] 科研人应该花精力去思考如何解决问题，而不是去凝练问题 +3 瞬息宇宙 2026-09-01 5/250 2026-09-02 08:01 by anjeeshine
[文学芳草园] 梦想 +5 myrtle 2026-08-26 8/400 2026-09-02 02:42 by Leogzhya
[基金申请] 为什么国自然不能直接公布 +5 bjdxyxy 2026-08-26 5/250 2026-09-01 17:54 by 小伟大博士
[基金申请] 学科评审组评审是指会评吗？ +4 瞬息宇宙 2026-08-31 4/200 2026-09-01 14:58 by jiaoxg
[基金申请] 为什么资助数各大高校都创新高，自己申请怎么就这么难 +13 Kittylucky 2026-08-27 14/700 2026-09-01 11:06 by feng6531
[论文投稿] 小白求助 投论文要求的highlights应该如何写 5+3 l1963982152 2026-08-29 4/200 2026-09-01 09:04 by 北京莱茵编辑
[基金申请] 麻烦专家们看看评委们的意见（F口面上） +7 gdd2018 2026-08-28 12/600 2026-09-01 08:32 by 尼古拉斯小虫
[基金申请] 基金未中，这种答复是模板吗？ +6 zhaosm1982 2026-08-27 7/350 2026-08-31 21:18 by qdxxmc
[基金申请] 哪位高人中了，把查询到的截图贴出来让我看看，让我长长见识 +6 yuleib84 2026-08-26 7/350 2026-08-31 19:46 by 鱼翔浅底1
[基金申请] 面上意见出来了 +12 黄鸟于飞Chao 2026-08-29 23/1150 2026-08-31 18:57 by 黄鸟于飞Chao
[基金申请] 中青基了要发朋友圈吗？ +7 349506619 2026-08-28 7/350 2026-08-31 13:39 by 冼亮淀粉酶
[基金申请] 29号明天会评吗 +4 笨笨唐 2026-08-28 4/200 2026-08-31 09:30 by huixian257
[基金申请] 为什么到现在没收到通知？ +5 tannykie 2026-08-29 5/250 2026-08-30 21:05 by purplejack
[基金申请] 国自然评审意见 +13 wangmingqi 2026-08-28 19/950 2026-08-29 10:22 by Poppy1104
[基金申请] 2026年叶企孙基金 +4 bud_bud 2026-08-27 7/350 2026-08-29 07:23 by foolishmani
[基金申请] 基金系统什么内容也没有 30+4 winsaint 2026-08-27 9/450 2026-08-28 11:06 by maolC
[基金申请] 怎么查啊 +6 huang1991js 2026-08-26 6/300 2026-08-28 08:42 by winsaint
[基金申请] 看板上这么多中的，有点像50人群里49个人都是骗子的那种感觉…… +5 a089 2026-08-26 6/300 2026-08-27 14:05 by jonewore
[基金申请] 我不理解！ +15 Edward_pc 2026-08-26 23/1150 2026-08-26 20:34 by zzuzxg
信息提示
关闭
请填处理意见
关闭