matlap新安江三水源模型程序.doc

.新安江模型程序核心源代码function fit,dc,result=XAJ(XX)% XAJ是新安江的运行程序,用于单纯形和遗传算法调用,也用于新安江模型的预报% XX是调用的优化参数% fit 返回目标函数的适值% dc返回有效性系数.% result是一个数组,返回格式为时间,雨量,实测流量,计算流量;.% $Date: 2005/5/25 $%email:damao_% 输入起始值 W,WU,WL,WD,QGWU=20;WL=50;WD=10;FR=0.89; S=2; AREA=7547;U=AREA/3.6;W=WU+WL+WD;%输入雨量E,蒸散发能力P,实测流量QSglobal DATA;TIME=DATA(:,1);P=DATA(:,2);EM=DATA(:,3);QS=DATA(:,4);TRSS0=0.3.*QS(1);TRG0=0.4.*QS(1);% 参数处理num,numvars=size(XX);% 优化参数A_K=XX(:,1);A_SM=XX(:,2);A_KG=XX(:,3);A_KSS=XX(:,4);A_KKG=XX(:,5);A_KKSS=XX(:,6);A_CS=XX(:,7);A_WUM=XX(:,8);A_WLM=XX(:,9);A_WDM=XX(:,10);A_IMP=XX(:,11);A_B=XX(:,12);A_C=XX(:,13);A_EX=XX(:,14);A_L=XX(:,15);A_WM=A_WUM+A_WLM+A_WDM;for I=1:num % % 对每组数计算 K=A_K(I); SM=A_SM(I); KG=A_KG(I); KSS=A_KSS(I); KKG=A_KKG(I); KKSS=A_KKSS(I); CS=A_CS(I); WUM=A_WUM(I); WLM=A_WLM(I); WDM=A_WDM(I); WM=WUM+WLM+WDM; IMP=A_IMP(I); B=A_B(I); C=A_C(I); EX=A_EX(I); L=A_L(I); L=round(L); WMM=(1+B).*WM/(1-IMP); M=size(P,1); PE=P-K.*EM; for T=1:M % T以时段为单位计算% %以下为产流计算 if PE(T)=WM A=WMM; else A=WMM*(1-(1-W/WM).(1/(1+B); end if A+PE(T)0 if A+PE(T)WMM R=PE(T)-WM+W+WM.*(1-(PE(T)+A)./WMM).(1+B); else R=PE(T)+W-WM; end else R=0; end end %5 % 以下为蒸发计算zhengfa if PE(T)0 EU=K*EM(T); ED=0; EL=0; WU=WU+PE(T); else EU=WU+P(T); WU=0; if WLC*WLM EL=(K.*EM(T)-EU).*WL/WLM; WL=WL-EL; ED=0; else if WLC.*(K.*EM(T)-EU) EL=C.*(K.*EM(T)-EU); WL=WL-EL; ED=0; else EL=WL; WL=0; ED=C.*(K*EM(T)-EU)-EL; WD=WD-ED; end end end else EU=K.*EM(T); ED=0; EL=0; if WU+PE(T)-RWLM WU=WUM; WL=WLM; WD=W+PE(T)-R-WU-WL; else WU=WUM; WL=WU+WL+PE(T)-R-WUM; end end end E=EU+EL+ED; W=WU+WL+WD;% 以下为分水计算% SMM=(1+EX).*SM; if (PE(T)=0)|(R=SM AU=SMM; else AU=SMM.*(1-(1-S./SM).(1./(1+EX); end if AU+QSMM RS=(Q-SM+S+SM.*(1-(Q+AU)./SMM).(1+EX).*FR+RS; else RS=(Q+S-SM).*FR+RS; end S=J.*Q-RS./FR+S; RG=S.*KGD.*FR+RG; RSS=S.*KSSD.*FR+RSS; S=J.*Q+SS-(RS+RSS+RG)./FR; end end OUT(T,:)=RS,RSS,RG; end % 一次数据演算完%以下为汇流计算% RS=OUT(:,1); RSS=OUT(:,2);RG=OUT(:,3); TRS(1)=RS(1).*U; TRSS(1)=TRSS0 ; TRG(1)=TRG0 ; TR(1)=TRS(1)+TRSS(1)+TRG(1); for T=2:M TRS(T)=RS(T).*U; TRSS(T)=TRSS(T-1).*KKSS+RSS(T).*(1-KKSS).*U; TRG(T)=TRG(T-1).*KKG+RG(T).*(1-KKG).*U; TR(T)=TRS(T)+TRSS(T)+TRG(T); end QJ=TR; if L800 y1=(QJ(T)-QS(T).2+y1; n1=n1+1; else y2=(QJ(T)-QS(T).2+y2; n2=n2+1; end end q0=mean(QS); q1=mean(QJ); y=(y1*alf/n1+y2*(1-alf)/n2)*(1+abs(q0-q1)/q0);fit(I)=y;%以下为(有效性系数)确定性系数计算% f1=sum( (QS-QJ).2);