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

}

文章转自:http://blog.csdn.net/superwen_go/article/details/7669429

基于遗传算法的新安江模型参数优化率定(四)相关推荐

  1. 新安江遗传算法c语言,基于遗传算法的新安江模型参数优化率定(四)

    4.3.1新安江三水源模型 //新安江三水源模型.hios #include 算法 #include 函数 #include 优化 #include spa const intVariableNum ...

  2. 基于遗传算法的新安江模型参数优化率定(一)

    1参数率定选用的流域 1.1流域概况 竹溪坡流域慨况及暴雨洪水一般特性: 竹溪坡流域位于湖南资水支流伊溪的上游,集水面积686km2,属山区,河长49km,干流坡降2.42‰.流域内有三个雨量站,即竹 ...

  3. 基于遗传算法的新安江模型参数优化率定(三)

  4. 【PEST++】02 新安江模型参数自动率定

    文章目录 PEST++系列文章 一.背景 1.1 模型简介 1.2 PEST++简介 1.3 所用程序 二.原理 2.1 目标函数 2.2 参数范围 三.过程 3.1 准备文件 3.1.1 实测值文件 ...

  5. 基于遗传算法的风电储能蓄电池容量优化配置matlab优化程序

    基于遗传算法的风电储能蓄电池容量优化配置 风电+储能蓄电池微电网配置(基于matlab的遗传算法微电网配置优化程序) 参考文献:基于遗传算法的风电储能蓄电池容量优化配置 摘要:为了降低独立风力发电系统 ...

  6. 基于遗传算法的微电网经济运行优化matlab程序

    基于遗传算法的微电网经济运行优化matlab程序 摘 要: 微电网作为智能电网的一部分,是分布式电源接入电网的一种有效手段,微电网经济运行是其中一个重要研究方面.考察微电网经济性,通常是从最小运行成本 ...

  7. 基于遗传算法的电动汽车有序充电优化调度

    基于遗传算法的电动汽车有序充电优化调度 软件:Matlab 利用遗传算法对电动汽车有序充电进行优化:优化目标包括充电费用最低,电动汽车充到足够的电,负荷峰谷差最小. 分别利用传统.精英和变异遗传算法进 ...

  8. 基于遗传算法的电动汽车有序充电优化调度 利用遗传算法对电动汽车有序充电进行优化;优化目标包括充电费用最低,电动汽车充到足够的电,负荷峰谷差最小

    基于遗传算法的电动汽车有序充电优化调度 软件:Matlab 利用遗传算法对电动汽车有序充电进行优化:优化目标包括充电费用最低,电动汽车充到足够的电,负荷峰谷差最小. 分别利用传统.精英和变异遗传算法进 ...

  9. 模型参数优化(一):遗传算法

    参数是指算法中的未知数,有的需要人为指定,比如神经网络算法中的学习效率,有的是从数据中拟合而来,比如线性回归中的系数,如此等等.在使用选定算法进行建模时,设定或得到的参数很可能不是最优或接近最优的,这 ...

最新文章

  1. 监控员工离职倾向系统已被下架,网友:劝你善良
  2. 嵌入式Linux操作系统学习规划 (转)
  3. 【百度】大型网站的HTTPS实践(一)——HTTPS协议和原理
  4. 前端学习(2769):发送网络请求
  5. mac linux 蓝牙键盘,还在纠结Mac版键盘?试试KeyRemap4MacBook吧!
  6. Win-MASM64汇编语言-CMPXCHG指令
  7. vue3中这几个变化你要注意了
  8. WPF 使用自定义的TTF字体
  9. 快速入门Maxwell基本操作流程(2D部分)
  10. Openssl CA证书生成以及双向认证,及windows系统证书批量导出,android cer转bks
  11. C语言入门递归算法——汉诺塔(简单易懂,最后还有汉诺塔游戏)
  12. ubuntu php mysql 乱码,ubuntu 服务器字符乱码问题
  13. 社群团购到底有哪些优势?为什么社群团购很容易打造爆品?
  14. 2022.07.10 第九小组 高小涵 学习笔记
  15. The Cook and the Chef: Musk’s Secret Sauce
  16. AMD三核、六核安装SQL2000
  17. 计算机课程设计参考文献,近几年课程设计参考文献 课程设计参考文献有哪些...
  18. 记录使用pytest测试UI自动化遇到的self = <script.test01_user_login.TestUserLogin object at 0x000001A8BE16E430>问题解决
  19. Qt基于定时器实现简单动图展示(2例)
  20. window自带的计算机应用程序,如何设置Window 10系统中默认的启动程序

热门文章

  1. Mac电脑下载的google chrome无法使用
  2. Security Best Practices+Klocwork
  3. 段码液晶屏学习应用笔谈
  4. document.documentElement对象
  5. SQL查询语句注入实战(手注,显注)
  6. Android libyuv应用系列(二)libyuv在Android中的使用
  7. Matlab学习日记(2)输入与输出
  8. 免费超大量邮件发送服务Amazon SES和Mailgun提供SMTP和API支持
  9. matlab画梯形并平移,matlab 批量处理梯形变形
  10. [置顶]乔布斯的斯坦福演讲(双语)