4.3.1新安江三水源模型
//新安江三水源模型.hios
#include 算法
#include 函数
#include 优化
#include spa
const intVariableNum = 11, //待率定参数个数.net
M = 485,//其为降雨起止间时段数(应该为)blog
Mq = 14;//单位线中q[i]的个数，此课设中为ip
const double F = 686, dt = 6;//面积、系列数据的时间间隔get
const double FE =0.8;input
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]
{
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]
{
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]
{
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]
{
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
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<
<
for (int i = 0; i < M; i++)
outfile<
outfile.close();
}
4.3.2基因遗传算法
//遗传算法.h
#include
#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<
<
<
<
<
//**********初始化
srand( (unsigned)time(NULL ) );
for(int i=0; i
for (int j=0; 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
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
{//将最优个体排至最后,
outfile<
temp =IndChr[IndividualNum - 1][k] ;
IndChr[IndividualNum- 1][k] = IndChr[Best][k];
IndChr[Best][k]= temp;
}//最优个体不参加选择、交叉
outfile<
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<
if(BestFitness == Fitness[Best])//进化停滞
{
if((++CountNum) >= SumNum)break;
}
else //进化成功
{
BestFitness =Fitness[Best];
CountNum = 1;
}
//**********选择(轮盘开始)
dSumPs = 0;
for(int i = 0; i
{
SumPs[i] =Fitness[i]/SumFitness;
dSumPs +=SumPs[i];
SumPs[i] =dSumPs;
for(int j = i+1; j< IndividualNum - 1; j++)
{//按升序排列随机数
if(rs[j]
{
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
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
{//随机交叉基因对
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
BestIndChr[k]= IndChr[x][k];//备份最优基因
for(int j = 0; 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
IndChr[x][k] = BestIndChr[k];//换回最优基因
Fitness[x]= BestFitness;
}
else
{
BestFitness= Fitness[x];//变异优化了
CountNum= 0;
}
}
}
}
outfile.close();
outputQ(IndChr[Best]);//输出模拟流量
}
//新安江模型参数率定.cpp
#include "stdafx.h"
#include "遗传算法.h"
void main()
{
inputQ();
YiChuanSuanFa();
}
文章转自：http://blog.csdn.net/superwen_go/article/details/7669429
扫一扫
6905
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
