基于遗传算法的新安江模型参数优化率定（四）

4.3  C++程序代码

4.3.1  新安江三水源模型

//新安江三水源模型.h

#include <fstream>

#include <iostream>

#include <iomanip>

#include <cmath>

 

const intVariableNum = 11, //待率定参数个数

              M = 485,//其为降雨起止间时段数(应该为)

              Mq = 14;//单位线中q[i]的个数，此课设中为

const double  F = 686, dt = 6;//面积、系列数据的时间间隔

const double FE =0.8;

const int  Days[4] = {30, 31, 30, 31}; //四至七月份每月天数

const doubleAveE[2][4] = {1.57, 2.29, 2.65, 3.41,0.97, 1.49, 1.71,     2.34};//四至七月份每月的多年平均日蒸发量

double K, WUM, WLM, C, //蒸发(蒸散发能力折算系数、分层蓄水容量、深层蒸散发系数)

        WM, b,    //产流(蓄水容量、蓄水容量曲线指数)

       SM,  EX, KI,   //水源划分(自由水蓄水容量、自由水蓄水容量曲线指数、自由水蓄水库出流系数)

        KKI, KKG;             //汇流(壤中流、地下径流消退系数)

 doubleP[M], Ep[M], EU[M], EL[M], ED[M], E[M], PE[M], WU[M + 1], WL[M + 1], WD[M + 1],W[M], a[M], R[M];

 doubleaf[M];//af[i]指产流面积比率(%),三水源划分中需要的数据

 doubleP1[M], P2[M], P3[M], Qo[ M ];     //三个雨量站实测雨量，实测流域流量

 doubleS[M], AU[M], RS[M], RI[M], RG[M];

 doubleq[Mq], QS[Mq + M - 1],Qs[Mq + M - 1][M];

 doubleQ[Mq + M - 1], QI[Mq + M - 1], QG[Mq + M - 1];

 doubleSumQo, AveQo, S1, S2;

 doubleU = F/3.6/dt;//折算系数

 doubleDC;//确定性系数

 

void inputQ();//读入原始数据

double FuntionDC(doubleCanshu[]);//计算确定性系数

void outputQ(doubleCanshu[]);//输出模拟流量

 

//**********读入原始数据(函数)************   

void inputQ()

{

       usingnamespace std;

        ifstream infile;//读人三个雨量站实测流域流量

        infile.open("infile_3P_Qo.txt");

        for(int i = 0; i < M; i++)

               infile>>P1[i]>>P2[i]>>P3[i]>>Qo[i];

 

       SumQo = 0, AveQo;

        for (int i = 0; i < M; i++)  

              SumQo += Qo[i];

       AveQo = SumQo/M;

        S2 = 0;

        for (int i = 0; i < M; i++)  

              S2 += pow(Qo[i] -AveQo, 2);

       infile.close();

 

        infile.open("infile_q.txt");

        for (int i = 0; i < Mq; i++)  

        {//读人时段单位线数据

               infile>>q[i];

        }

        infile.close();

}

 

//**********计算确定性系数(函数)************      

 doubleFuntionDC(double Canshu[])

 {

       usingnamespace std;

       K = Canshu[10]; WM =Canshu[0]; WUM = Canshu[1]; WLM = Canshu[2]; C = Canshu[3];

       b = Canshu[4];

       SM = Canshu[5]; EX = Canshu[6];KI = Canshu[7];             

       KKI = Canshu[8]; KKG  = Canshu[9];     

 

 

   //******三层蒸发模式下的蓄满产流模型(开始)

        double WDM =WM - WUM - WLM, KG = 0.7 - KI;

        double WMM =WM*(1 + b);

        WU[0] = FE*WUM;

        WL[0] = FE*WLM;

        WD[0] = FE*WDM;

 

       //********计算蒸发能力

       intSumTime1 = 0, SumTime2 = 0;

       for(int j = 0; j < 4; j++)

       {

              SumTime1 =  SumTime2;

              SumTime2+=4*Days[j];

              if(SumTime2 > M) SumTime2 = M;

               for(int i = SumTime1; i < SumTime2; i++)

               {

                      P[i] = (P1[i]+P2[i]+P3[i])/3;

                     if(P[i] < 3)  Ep[i] = AveE[0][j] * K;

                     else   Ep[i] =AveE[1][j] * K;

               }

       }

        for (int i = 0; i < M; i++)

        {

               W[i] = WU[i] + WL[i]+ WD[i];

               if((1 -W[i]/WM)<0)

                      a[i] = WMM;

               else a[i] =WMM*(1 - pow(1 - W[i]/WM,1/(b + 1)));

 

               if ( P[i] ==0)

               {

                      if (WU[i] <Ep[i])

                      {

                             EU[i] = WU[i];

                             if (WL[i]>= C*WLM)

                             {

                                    EL[i] = (Ep[i] - EU[i])*WL[i]/WLM;

                                    ED[i] = 0;

                             }

                             else

                                    if (WL[i] <(Ep[i] - EU[i])*C)

                                    {

                                           EL[i] = WL[i];

                                           ED[i] = (Ep[i] - EU[i])*C - EL[i];

                                    }

                                    else

                                    {

                                           EL[i] = (Ep[i] - EU[i])*C;

                                           ED[i] = 0;

                                    }

                             E[i] = EU[i] + EL[i] + ED[i];

                             WU[i + 1] = 0;

                             WL[i + 1] = WL[i] - EL[i];

                             WD[i + 1] = WD[i] - ED[i];

                      }

                      else

                      {

                             EL[i] = ED[i] = 0;

                             E[i] = EU[i] = Ep[i];

                             WU[i + 1] = WU[i] - Ep[i];

                             WL[i + 1] = WL[i];

                             WD[i + 1] = WD[i];

                      }

                      PE[i] = -E[i];   //PE为负值

                      R[i] = 0;

               }

               else

               {//3

                      if (P[i] +WU[i] < Ep[i])

                      {

                             EU[i] =P[i] + WU[i];

                             if (WL[i] >=C*WLM)

                             {

                                    EL[i] = (Ep[i] - EU[i])*WL[i]/WLM;

                                    ED[i] = 0;

                             }

                             else

                                    if (WL[i] <(Ep[i] - EU[i])*C)

                                    {

                                           EL[i] = WL[i];

                                           ED[i] = (Ep[i] - EU[i])*C - EL[i];

                                    }

                                    else

                                    {

                                           EL[i] = (Ep[i] - EU[i])*C;

                                           ED[i] = 0;

                                    }

                             E[i] = EU[i] + EL[i] + ED[i];

                             PE[i] = P[i] - E[i];    //PE为负值

                             R[i] = 0;

                             WU[i + 1] = 0;

                             WL[i + 1] = WL[i] - EL[i];

                             WD[i + 1] = WD[i] - ED[i];

                      }

                      else

                      {

                             EL[i] = ED[i] = 0;

                             E[i] = EU[i] = Ep[i];

                             PE[i] = P[i] - E[i];

                             if (P[i] <Ep[i])    //PE为负值

                             {

                                    R[i] = 0;

                                WU[i + 1] = WU[i] +  P[i] - E[i];

                                WL[i + 1] = WL[i];

                                   WD[i+ 1] = WD[i];

                             }//到此,因PE为负而无产流的情况全部讨论完毕

                             else

                             {//以下情况出现产流，注意此时蒸发只发生在WU，即EL=ED=0

                                    if (a[i] +PE[i] <= WMM)

                                    {

                                           R[i] = PE[i] + W[i]- WM  + WM*(pow(1 - (a[i] + PE[i])/WMM,b + 1));

                                           WU[i + 1] = PE[i] + WU[i]- R[i];

                                           WL[i + 1] = WL[i ];

                                           WD[i + 1] = WD[i];

                                           if (WU[i + 1]> WUM)

                                           {

                                                 WL[i+ 1] = WU[i + 1] - WUM + WL[i];

                                                 WU[i+ 1] =  WUM;

                                                 if (WL[i + 1] > WLM)

                                                 {

                                                        WD[i+ 1] = WL[i + 1] - WLM + WD[i];

                                                        WL[i+ 1] = WLM;

                                                        if (WD[i + 1] > WDM) WD[i + 1] = WDM;

                                                   }

                                             }

                                    }

                                    else

                                    {

                                           R[i] = PE[i] + W[i]- WM;

                                           WU[i + 1] = WUM;

                                           WL[i + 1] = WLM;

                                           WD[i + 1] = WDM;

                                    }

                             }

                      }

               }

               if((a[i]+PE[i])>WMM)

                      af[i] = 1;

               else

                     af[i] = 1 -pow(1 - (a[i]+PE[i])/WMM,b);

       }

     //**********三水源划分(开始)

        double SMM =SM*(1 + EX);

        double S0 =FE*SM, af0 = 1 - pow(1 - a[0]/WMM,b);

                                                        //af0指W[0]对应的面积比率(%)

        for (int i = 0; i < M; i++)

        {

               if(i == 0)     S[i] = S0*af0/af[i];

               else S[i] =S[i - 1]*af[i - 1]/af[i]+( R[i - 1]- RS[i - 1] - RI[i - 1] - RG[i - 1])/af[i];

               if(S[i]>SM)S[i] = SM;

               AU[i] = SMM*(1 - pow(1 - S[i]/SM,1/(EX + 1)));

               if(R[i] == 0)

                      RS[i]  =0;

               else

               {

                      if(AU[i] +PE[i] > SMM)

                            RS[i] =(S[i] + PE[i] - SM)*af[i];

                  else

                            RS[i] =(S[i] + PE[i] - SM + SM*(pow(1 - (AU[i] + PE[i])/SMM,EX + 1)))*af[i];

               }                

               RI[i] = KI*(S[i] + (R[i] -RS[i])/af[i])*af[i];

               RG[i]= KG/KI*RI[i];

        }

 

     //**********汇流(开始)

        for(int  j = 0; j< M; j++)

               for(int i = 0; i < Mq + M - 1; i++)

                      if(i <j)  Qs[i][j] = 0;

                      else

                      {

                             if(i < j +Mq)  Qs[i][j] = RS[j]/10*q[i - j];

                             else Qs[i][j]= 0;

                      }

 

        for(int  i = 0; i< Mq + M - 1; i++)

        {

              QS[i] = 0;//一定要初始化

               for(int j = 0; j < M ; j++)

                       QS[i] += Qs[i][j];

        }

        QI[0] = RI[0]*(1 - KKI)*U;

        QG[0] = RG[0]*(1 - KKG)*U;

        for(int i = 1; i < Mq + M - 1; i++)

               if(i < M )

               {

                      QI[i] = RI[i]*(1 - KKI)*U + QI[i - 1]*KKI;

                      QG[i] = RG[i]*(1 - KKG)*U + QG[i - 1]*KKG;

               }

               else

               {

                      QI[i] = QI[i - 1]*KKI;

                      QG[i] = QG[i - 1]*KKG;

               }  

 

        for(int  i = 0; i< Mq + M - 1; i++)

               Q[i] = QS[i] + QI[i] + QG[i];

       //*****确定性系数计算(开始)

        S1 = 0;

        for (int i = 0; i < M; i++)  

              S1 += pow(Q[i] - Qo[i],2);

       DC = 1 - S1/S2;

       returnDC;

 }

//**********输出模拟流量过程(函数)************

 voidoutputQ(double Canshu[])

 {

       usingnamespace std;

        ofstream outfile;

        outfile.open("outfile_Q.txt");

        outfile<<"模拟效率系数"<<FuntionDC(Canshu)<<endl

               <<setw(10)<<"P雨量"<<setw(10)<<"Q实测"<<setw(10)<<"Q模拟"<<endl;

        for (int i = 0; i < M; i++)  

              outfile<<setw(10)<<P[i]<<setw(10)<<Qo[i]<<setw(10)<<Q[i]<<endl;     

        outfile.close();

 }

 

4.3.2  基因遗传算法

//遗传算法.h

#include <ctime>

#include "新安江三水源模型.h"

const intGenerationNum = 200,//最大演算世代次数

            SumNum = 60,//当最优适应度重复次数超过此值时停止演算

                     IndividualNum =21,//该种群的个体数目

                     ChromosomeNum =11;//每个个体的基因（待率定参数）数目

                           

const double ChrTop[ChromosomeNum]//基因（待率定参数）的阈值上限

                                                 ={200,20, 90, 0.20, 0.4, 25, 1.5, 0.7, 1.0, 1.0, 1.5},

                     ChrBottom[ChromosomeNum]//基因（待率定参数）的阈值下限

                                   ={120,10, 60, 0.09, 0.1,  5,  1.0, 0.0, 0.0, 0.0, 0.5},

                     Pc = 0.5,//个体的交叉率（crossoverrate）

                     PcChr = 0.7,//交叉对交叉的基因比率

                     PmInd = 0.7,//个体变异率（mutationrate）

                     PmChr = 0.5,//变异个体的基因变异率（mutationrate)

                  Bm = 4;//变异运算的形状系数

 intnc =  (int)((IndividualNum - 1)*Pc/2),//nc对个体的基因参与交叉

       ncChr = (int) (ChromosomeNum*PcChr),//两个体交叉的基因数

       nmInd = (int) ((IndividualNum - 1)*PmInd),//nmInd个个体发生变异

       nmChr = (int) (ChromosomeNum*PmChr),//个体的nmChr个基因发生变异

       x, y,tx1,tx2, Best,Worst,//挑出最优及最差的个体

      CountNum= 1;//计数最优适应度重复次数

 doubleIndChr[IndividualNum][ChromosomeNum],//每代种群的所有个体的基因

           Fitness[IndividualNum],//个体的适应度

              BestFitness = 0,//备份最优个体的适应度

              BestIndChr[ChromosomeNum],//备份最优个体的基因

              SumFitness = 0,//累积适应度

              SumPs[IndividualNum] ={0}, //累积选择概率

              dSumPs = 0,//用来求累积选择概率的

               r = 0,//伪随机数,交叉变异时使用

               rs[IndividualNum] ={ 0},//伪随机数,选择时使用

               temp;//中间变量

 

void YiChuanSuanFa()

{

       usingnamespace std;

        ofstream outfile,outtext;

        outfile.open("outfile_BestIndividual.txt");//写出文件，用于绘制遗传算法的进化过程

        outfile<<setw(10)<<"WM"<<setw(10)<<"WUM"<<setw(10)<<"WLM"<<setw(10)<<"C"<<setw(10)<<"b"

            <<setw(10)<<"SM"<<setw(10)<<"EX"<<setw(10)<<"KI"

               <<setw(10)<<"KKI"<<setw(10)<<"KKG"

               <<setw(10)<<"K"

               <<setw(10)<<"BestFit" <<setw(10)<<"WorstFit" <<setw(10)<<"AverageFit"<<endl;

 

        //**********初始化

        srand( (unsigned)time(NULL ) ); 

        for(int i=0; i<IndividualNum; i++)

              for (int j=0; j<ChromosomeNum ;j++)

                     IndChr[i][j] = (double)rand()/RAND_MAX*(ChrTop[j]-ChrBottom[j])+ChrBottom[j];

       //**********世代更替演算(开始)

       for(int g = 1; g< GenerationNum; g++)

       {

           //**********适应度(开始)

              for(int i = 0; i <IndividualNum - 1; i++)

                     Fitness[i]=  FuntionDC(IndChr[i]);

              if(g == 1) //计算初始化的最后一个个体的适应度

                     Fitness[IndividualNum- 1] =  FuntionDC(IndChr[IndividualNum -1]);

 

              for(tx1 = 0; tx1 < IndividualNum; tx1++)

              {//找出最优个体

                     Best = tx1;

                     for( tx2 = tx1+1; tx2< IndividualNum; tx2++)

                             if(Fitness[tx1]< Fitness[tx2])

                            {

                                   Best= tx2;

                                   break;

                            }

                     if(tx2 ==IndividualNum) break;

              }

              for(tx1 = 0; tx1 < IndividualNum; tx1++)

              {//找出最差个体

                     Worst = tx1;

                     for( tx2 = tx1+1; tx2< IndividualNum; tx2++)

                             if(Fitness[tx1]> Fitness[tx2])

                            {

                                   Worst= tx2;

                                   break;

                            }

                     if(tx2 ==IndividualNum) break;

              }

              for (int k=0; k<ChromosomeNum;k++)

              {//将最优个体排至最后,

                     outfile<<setw(10)<<IndChr[Best][k];

                     temp =IndChr[IndividualNum - 1][k] ;

                     IndChr[IndividualNum- 1][k] = IndChr[Best][k];

                     IndChr[Best][k]= temp;

              }//最优个体不参加选择、交叉

              outfile<<setw(10)<<Fitness[Best]<<setw(10)<<Fitness[Worst];

              temp =  Fitness[IndividualNum - 1];

              Fitness[IndividualNum -1] = Fitness[Best];

              Fitness[Best] = temp;

              SumFitness = 0;//初始化

               for(int i = 0; i < IndividualNum - 1; i++)

              {

                     rs[i] = (double)rand()/RAND_MAX;

                     SumFitness +=Fitness[i];

              }

              outfile<<setw(10)<<(SumFitness+ Fitness[IndividualNum - 1])/IndividualNum<< endl;

              if(BestFitness == Fitness[Best])//进化停滞

              {

                     if((++CountNum) >= SumNum)break;

              }

              else //进化成功

              {

                     BestFitness =Fitness[Best];

                     CountNum = 1;

              }

                     //**********选择(轮盘开始)

              dSumPs = 0;

              for(int i = 0; i <IndividualNum - 1; i++)

              {

                     SumPs[i] =Fitness[i]/SumFitness;

                     dSumPs +=SumPs[i];

                     SumPs[i] =dSumPs;

                     for(int j = i+1; j< IndividualNum - 1; j++)

                     {//按升序排列随机数

                            if(rs[j]<rs[i])

                            {

                                   temp= rs[i];

                                   rs[i]= rs[j];

                                   rs[j]= temp;

                            }

                     }

              }

           for(int i = 0; i < IndividualNum - 1; i++)

              {//最优个体不参加选择

                     for(int j = 0 ; j< IndividualNum - 1; j++)

                            if (SumPs[j] > rs[i])

                            {

                                   for (int k=0;k<ChromosomeNum; k++)

                                          IndChr[i][k]= IndChr[j][k];

                                   break;

                            }

              }

       //**********交叉(开始)************

              for(int i = 0; i < nc; i++)

              {//随机交叉个体对

                     x = y = 0;//最优个体不参加交叉

                     while(x == y)

                     {//+0.5四舍五入

                            x = (int)((double)rand()/RAND_MAX*(IndividualNum - 2) + 0.5);

                            y = (int)((double)rand()/RAND_MAX*(IndividualNum- 2) + 0.5);

                     }

                     for(int j = 0; j <ncChr; j++)

                     {//随机交叉基因对

                            r = (double)rand()/RAND_MAX;

                            int k = (int)((double)rand()/RAND_MAX*(ChromosomeNum - 1) + 0.5);

                            temp =IndChr[x][k] ;

                            IndChr[x][k]= r*temp+(1-r)*IndChr[y][k];

                            IndChr[y][k]= r*IndChr[y][k]+(1-r)*temp;

                     }

              }

                     //**********变异(开始)************

              double t = g/GenerationNum;//变异运算的进化标记

              for(int i = 0; i < nmInd; i++)//随机变异个体

              {//最优个体只能进行优化变异

                     x = (int)((double)rand()/RAND_MAX*(IndividualNum- 1) + 0.5);

                     if(x== (IndividualNum - 1))

                            for (int k=0;k<ChromosomeNum ;k++)

                                   BestIndChr[k]= IndChr[x][k];//备份最优基因

                     for(int j = 0; j <nmChr; j++)//随机变异个体上的随机变异基因

                     {

                            int k = (int)((double)rand()/RAND_MAX*(ChromosomeNum - 1) + 0.5);

                            r = (double)rand()/RAND_MAX;

                            if((rand()%2) ==0)

                                   IndChr[x][k]+= (ChrTop[k] - IndChr[x][k])*pow(r*(1-t),Bm);

                            else

                                   IndChr[x][k]-= (IndChr[x][k]  -ChrBottom[k])*pow(r*(1-t),Bm);

                     }

                     if(x== (IndividualNum - 1)) 

                     {//判断最优基因变异后是否优化了

                            Fitness[x]=  FuntionDC(IndChr[x]);

                            if(Fitness[x] < BestFitness)//变异退化了

                            {

                                   for (int k=0;k<ChromosomeNum ;k++)

                                     IndChr[x][k] = BestIndChr[k];//换回最优基因

                                   Fitness[x]= BestFitness;

                            }

                            else

                            {

                                   BestFitness= Fitness[x];//变异优化了

                                   CountNum= 0;

                            }

                     }

              }    

       }

        outfile.close();

        outputQ(IndChr[Best]);//输出模拟流量

}

4.3.3  主函数

//新安江模型参数率定.cpp

#include "stdafx.h"

#include "遗传算法.h"

void main()

{

       inputQ();

      YiChuanSuanFa();

}