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

墨蓝 2022-08-24 11:18 168阅读 0赞

## 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

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

发表评论取消回复

还没有评论，来说两句吧...

相关阅读

相关基于遗传算法优化SVM参数的热负荷预测

相关基于遗传算法的装配线调度优化

相关基于遗传算法优化的BP神经网络

相关 matlab 遗传算法参数优化,Matlab基于BP神经网络模型的遗传算法多参数优化程序

相关【优化求解】基于遗传算法的PID参数优化设计matlab代码

相关【优化求解】基于遗传算法的PID参数优化设计matlab代码

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

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

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

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