基于遗传算法的新安江模型参数优化率定（四）
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();
}
