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

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();

}