MATLAB实现基于遗传算法/引力搜索算法优化新安江水文模型

%% 根据求得最优参数计算流量 进行绘图 % 蒸散发参数: KC、WUM、WLM、C； % 产流量参数： WM、B； % 水源划分参数： SM、EX、 KSS、KG % 坡地汇流参数：KKSS、KKG KC=g_best(max_it,1); WUM=g_best(max_it,2); WLM=g_best(max_it,3); C=g_best(max_it,4); WM=g_best(max_it,5); B=g_best(max_it,6); SM=g_best(max_it,7); EX=g_best(max_it,8); KSS=g_best(max_it,9); KG=g_best(max_it,10); KKSS=g_best(max_it,11); KKG=g_best(max_it,12); % 率定期 % 将原始的流域统计资料进行赋值存储 Pfile=fopen('P87-89.DAT','rt'); %降雨 Pcell=textscan(Pfile,'%s %f'); fclose(Pfile); Qfile=fopen('Q87-89.DAT','rt'); %流量 Qcell=textscan(Qfile,'%s %f'); fclose(Qfile); Efile=fopen('E87-89.DAT','rt'); %蒸发 Ecell=textscan(Efile,'%s %f'); fclose(Efile); DATE=Pcell{1}'; P=cell2mat(Pcell(:,2)); Q=cell2mat(Qcell(:,2)); E0=cell2mat(Ecell(:,2)); %蒸发皿测得多年平均蒸发量 FDATE_FILE=fopen('F87-89.DAT','rt'); FDATE_cell=textscan(FDATE_FILE,'%d %s %s'); fclose(FDATE_FILE); FNUM=FDATE_cell{1}; %FNUM 洪号 FDATES=FDATE_cell{2}; %FDATES 实测流量起涨时间 FDATEE=FDATE_cell{3}; %FDATEE 实测流量退水时间 %将给定的单位线进行赋值存储 load('q.mat') %24h10min单位线 UH=q; %单位线 K=48; % 蓄泄系数为48h(小时) FE=0.8; % 流域初始含水量可在0 ~ 0.2WM 之间选择，即初始土壤含水容量折算系数FE=0.8~1 F=6448; % 流域面积单位为km^2 IMP=0.001; % 流域不透水面积比重 % 初始数据处理 EP=KC.*E0; %水面蒸发乘系数KC转换成流域蒸发能力值 WDM=WM-WUM-WLM; S_0=SM*FE; %初始自由水蓄量 WU_0=FE*WUM; WL_0=FE*WLM; WD_0=FE*WDM; %初始土壤含水量 %模型运算 [PO,RCE,DATE_CE,QCE,RCE1]=Datastatistics(DATE,P,Q,F,K,FNUM,FDATES,FDATEE); [R,E,PE,W,WMM,WU,WL,WD,a]=runoffgenerate(P , EP , WM , WUM , WLM , WDM , WU_0 , WL_0 , WD_0 , C , B , IMP); [RS,RI,RG,FR,AU,S]=three_water_sources_division(P,EP,PE,W,WM,B,EX,KG,KSS,S_0,SM,R,IMP); [QRS,QRI,QRG,Q_all,R_all]=confluence(UH,F,RS,RI,RG,KKG,KKSS); % 预报评价 [DC,DC_all ]=EVALUATE_DC(DATE,FDATES,FDATEE,Q_all,Q); % 计算洪峰确定性系数和流量相对误差RE [RE,DC_qmax]=lastcheck(Q ,Q_all , DATE , FDATES , FDATEE); RE=RE/(3*365); fprintf('率定期RE值为 %5.3f\n',RE);%目标函数值 fprintf('率定期DC值为 %5.3f\n',DC_all);%目标函数值 %% 1987年 Q1=Q(1:365);P1=P(1:365);Q_all1=Q_all(1:365); figure(12) n1=length(Q1); hold on;box on;grid off; [AX,H1,H2]=plotyy(1:n1,Q1,1:n1,P1,'plot','bar'); %画双轴，AX(1)左轴，AX（2）右轴，H为曲线本身 set(AX(2),'YDir','reverse','Ylim',[0,max(P1)*3],'YTick',[0:30:max(P1)*2.5],'FontSize',10); %设置右边轴为倒立 set(gca,'box','off','Ytick',[]) set(AX(1),'YLim',[0,1.5*max(Q1)],'YTick',[0:500:1.5*max(Q1)],'Fontsize',10,'YColor','k'); H3=plot(1:n1,Q_all1,'r-','LineWidth',0.8); %设置坐标轴的标题 set(get(AX(1),'Xlabel'),'String','日期(d)'); set(get(AX(1),'Ylabel'),'String','流量Q（m^3/s)'); set(get(AX(2),'Ylabel'),'string','降雨量P/mm'); %设置figuer中线，柱的属性 axis([1 365 0 1.5*max(Q1)]); set(H2,'BarWidth',0.5 ,'facecolor',[0.2 0.4 1]); set(H1,'LineWidth',0.8,'Color',[1 0.25 0]); %设置H1的曲线宽为4，颜色为g绿色 set(H3,'LineWidth',0.8,'Color',[0.28235 0.87843 0.8]); set(gca,'xtick',[1 93 210 300],'xticklabel',{'1987/1/1' ,'1987/4/1' ,'1987/7/1' ,'1987/10/1'}); set(gcf,'color','w'); %设置当前figure的背景颜色 legend('模拟值','实测值','降水量');%设置图例 title('1987年降雨径流图（率定期）') %% 1988年 Q2=Q(365+1:365+365);P2=P(365+1:365+365);Q_all2=Q_all(365+1:365+365); figure(13) n2=length(Q2); hold on;box on;grid off; [AX,H1,H2]=plotyy(1:n2,Q2,1:n2,P2,'plot','bar'); %画双轴，AX(1)左轴，AX（2）右轴，H为曲线本身 set(AX(2),'YDir','reverse','Ylim',[0,max(P2)*3],'YTick',[0:30:max(P2)*2.5],'FontSize',10); %设置右边轴为倒立 set(gca,'box','off','Ytick',[]) set(AX(1),'YLim',[0,1.5*max(Q2)],'YTick',[0:200:1.5*max(Q2)],'Fontsize',10,'YColor','k'); H3=plot(1:n1,Q_all2,'r-','LineWidth',0.8); %设置坐标轴的标题 set(get(AX(1),'Xlabel'),'String','日期(d)'); set(get(AX(1),'Ylabel'),'String','流量Q（m^3/s)'); set(get(AX(2),'Ylabel'),'string','降雨量P/mm'); %设置figuer中线，柱的属性 axis([1 365 0 1.5*max(Q2)]); set(H2,'BarWidth',0.5 ,'facecolor',[0.2 0.4 1]); set(H1,'LineWidth',0.8,'Color',[1 0.25 0]); %设置H1的曲线宽为4，颜色为g绿色, set(H3,'LineWidth',0.8,'Color',[0.28235 0.87843 0.8]); set(gca,'xtick',[1 93 210 300],'xticklabel',{'1988/1/1' ,'1988/4/1' ,'1988/7/1' ,'1988/10/1'}); set(gcf,'color','w'); %设置当前figure的背景颜色 legend('模拟值','实测值','降水量');%设置图例 title('1988年降雨径流图（率定期）') %% 1989年 Q3=Q(365*2+1:365*2+365);P3=P(365*2+1:365*2+365);Q_all3=Q_all(365*2+1:365*2+365); figure(14) n3=length(Q3); hold on;box on;grid off; [AX,H1,H2]=plotyy(1:n3,Q3,1:n3,P3,'plot','bar'); %画双轴，AX(1)左轴，AX（2）右轴，H为曲线本身 set(AX(2),'YDir','reverse','Ylim',[0,max(P3)*3],'YTick',[0:30:max(P3)*2.5],'FontSize',10); %设置右边轴为倒立 set(gca,'box','off','Ytick',[]) set(AX(1),'YLim',[0,1.5*max(Q3)],'YTick',[0:200:1.5*max(Q3)],'Fontsize',10,'YColor','k'); H3=plot(1:n1,Q_all3,'r-','LineWidth',0.8); %设置坐标轴的标题 set(get(AX(1),'Xlabel'),'String','日期(d)'); set(get(AX(1),'Ylabel'),'String','流量Q（m^3/s)'); set(get(AX(2),'Ylabel'),'string','降雨量P/mm'); %设置figuer中线，柱的属性 axis([1 365 0 1.5*max(Q3)]); set(H2,'BarWidth',0.5 ,'facecolor',[0.2 0.4 1]); set(H1,'LineWidth',0.8,'Color',[1 0.25 0]); %设置H1的曲线宽为4，颜色为g绿色, set(H3,'LineWidth',0.8,'Color',[0.28235 0.87843 0.8]); set(gca,'xtick',[1 93 210 300],'xticklabel',{'1989/1/1' ,'1989/4/1' ,'1989/7/1' ,'1989/10/1'}); set(gcf,'color','w'); %设置当前figure的背景颜色 legend('模拟值','实测值','降水量');%设置图例 title('1989年降雨径流图（率定期）') %% 检验期 % inputs 输入 % 将原始的流域统计资料进行赋值存储 Pfile=fopen('P90.DAT','rt'); %降雨 Pcell=textscan(Pfile,'%s %f'); fclose(Pfile); Qfile=fopen('Q90.DAT','rt'); %流量 Qcell=textscan(Qfile,'%s %f'); fclose(Qfile); Efile=fopen('E90.DAT','rt'); %蒸发 Ecell=textscan(Efile,'%s %f'); fclose(Efile); DATE_check=Pcell{1}'; P_check=cell2mat(Pcell(:,2)); Q_check=cell2mat(Qcell(:,2)); E0_check=cell2mat(Ecell(:,2)); %蒸发皿测得多年平均蒸发量 FDATE_FILE=fopen('F90.DAT','rt'); FDATE_cell=textscan(FDATE_FILE,'%d %s %s'); fclose(FDATE_FILE); FNUM_check=FDATE_cell{1}; %FNUM 洪号 FDATES_check=FDATE_cell{2}; %FDATES 实测流量起涨时间 FDATEE_check=FDATE_cell{3}; %FDATEE 实测流量退水时间 % 初始数据处理 EP=KC.*E0_check; %水面蒸发乘系数KC转换成流域蒸发能力值 WDM=WM-WUM-WLM; S_0=SM*FE; %初始自由水蓄量 WU_0=FE*WUM; WL_0=FE*WLM; WD_0=FE*WDM; %初始土壤含水量 %场次洪水统计 [PO,RCE,DATE_CE,QCE,RCE1]=Datastatistics(DATE_check,P_check,Q_check,F,K,FNUM_check,FDATES_check,FDATEE_check); % 模型运算 [R,E,PE,W,WMM,WU,WL,WD,a]=runoffgenerate(P_check , EP , WM , WUM , WLM , WDM , WU_0 , WL_0 , WD_0 , C , B , IMP); [RS,RI,RG,FR,AU,S]=three_water_sources_division(P_check,EP,PE,W,WM,B,EX,KG,KSS,S_0,SM,R,IMP); [QRS,QRI,QRG,Q_all_check,R_all_check]=confluence(UH,F,RS,RI,RG,KKG,KKSS); % 预报评价 [DC_check,DC_all_check ]=EVALUATE_DC(DATE_check,FDATES_check,FDATEE_check,Q_all_check,Q_check); % 计算洪峰确定性系数和流量相对误差RE [RE_check,DC_qmax_check]=lastcheck(Q_check ,Q_all_check , DATE_check , FDATES_check , FDATEE_check); RE_check=RE_check/365; fprintf('检验期RE值为 %5.3f\n',RE_check); fprintf('检验期DC值为 %5.3f\n',DC_all_check); %% 1990年 figure(15) hold on;box on;grid off; [AX,H1,H2]=plotyy(1:n3,Q_check,1:n3,P_check,'plot','bar'); %画双轴，AX(1)左轴，AX（2）右轴，H为曲线本身 set(AX(2),'YDir','reverse','Ylim',[0,max(P_check)*4],'YTick',[0:20:max(P_check)*3],'FontSize',10); %设置右边轴为倒立 set(gca,'box','off','Ytick',[]) set(AX(1),'YLim',[0,2000],'YTick',[0:200:2000],'Fontsize',