xaj.asv

这个同样是利用matlab编程

function [fit,dc,result]=XAJ(XX)
% XAJ是新安江的运行程序,用于单纯形和遗传算法调用,也用于新安江模型的预报
% XX是调用的优化参数
% fit 返回目标函数的适值
% dc返回有效性系数.
% result是一个数组,返回格式为[时间,雨量,实测流量,计算流量];
% 输入起始值 W,WU,WL,WD,QG
WU=20;WL=50;WD=10;FR=0.5; S=2; AREA=1640;U=AREA/3.6;
W=WU+WL+WD;

%输入雨量P,蒸散发能力E,实测流量QS
global DATA;
TIME=DATA(:,1);
P=DATA(:,2);
EM=DATA(:,3);
QS=DATA(:,4);

QI0=0.3.*QS(1);    
QG0=0.4.*QS(1);
%  参数处理
[num,numvars]=size(XX);
% 优化参数
A_KC=XX(:,1);
A_SM=XX(:,2);
A_KG=XX(:,3);
A_KI=XX(:,4);
A_CG=XX(:,5);
A_CI=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     %%%% %%% 对每组数计算

    KC=A_KC(I);
    SM=A_SM(I);
    KG=A_KG(I);
    KI=A_KI(I);
    CG=A_CG(I);
    CI=A_CI(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-KC.*EM;
    for T=1:M               %%  T以时段为单位计算
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %以下为产流计算
        if PE(T)<0
            R=0;
        else
            if  W>=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
            if WU+PE(T)>0
                EU=K*EM(T);
                ED=0;
                EL=0;
                WU=WU+PE(T);
            else
                EU=WU+P(T);
                WU=0;
                if WL>C*WLM
                    EL=(K.*EM(T)-EU).*WL/WLM;
                    WL=WL-EL;
                    ED=0;
                else
                    if WL>C.*(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)-R<WUM
                WU=WU+PE(T)-R;
            else
                if WU+WL+PE(T)-WUM>WLM
                    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<=0)
            RS=0;
            RG=S.*KG.*FR;
            RSS=RG.*KI./KG;
        else
            X=FR;
            FR=(R-PE(T).*IMP)./PE(T);
            S=X.*S./FR;
            SS=S;
            Q=R./FR;
            G=fix(Q./5)+1;
            Q=Q./G;
            KID=KI.^(1/G);
            KGD=KID.*KG./KI;
            RS=0;
            RG=0;
            RI=0;
            for J=1:G
                if S>=SM
                    AU=SMM;
                else
                    AU=SMM.*(1-(1-S./SM).^(1./(1+EX)));
                end

                if AU+Q<SMM
                    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;
                RI=S.*KID.*FR+RI;
                S=J.*Q+SS-(RS+RI+RG)./FR;
            end
        end
        OUT(T,:)=[RS,RSS,RG];
    end     % 一次数据演算完
%以下为汇流计算%%%%%%%%%%%%%%%%%%%%%%%%%%
    RS=OUT(:,1); RI=OUT(:,2);RG=OUT(:,3);
    QS(1)=RS(1).*U;
    QI(1)=QI0 ;
    QG(1)=QG0 ;
    Q(1)=QS(1)+QI(1)+QG(1);
    for T=2:M
        QS(T)=RS(T).*U;
        QI(T)=QI(T-1).*CI+RI(T).*(1-CI).*U;
        QG(T)=QG(T-1).*CG+RG(T).*(1-CG).*U;
        Q(T)=QS(T)+QI(T)+QG(T);
    end
    QJ=Q;
    if L<0  L=0;end
    for T=L+2:M
        QJ(T)=CS.*QJ(T-1)+(1-CS).*Q(T-L);
    end
%以下为目标函数计算%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    alf=0.6;
    y1=0;y2=0;
    n1=1;n2=1;
    for T=1:M
        if QJ(T)>800
            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);
    f2=sum((QS-mean(QS).*ones(M,1)).^2);
    dq=1-f1/f2;
    dc(I)=dq;
    result =[TIME,P,QS,QJ'];
end  %一组参数计算结束I

fit=-fit'; %遗传算法为了求最大值,在此加负号.
dc=dc';