动态规划进行单库优化调度计算C++程序代码

1、#ifndef SRDP_H#define SRDP_H class SRDPprivate:/1.阶段变量 stage variablelong m_StageNum;double m_Delta_t;double m_DeltaHour;/2.决策变量 decision variable/double *m_OptDeciVar_q;/3.状态变量 state variablelong *m_StateNum;/double *m_StateV;/4.指标函数 target functiondouble *m_TarFun;/5.目标函数object functiondouble m_Ob

2、jFun;最优轨迹记录数组long *m_OptArray;原始资料记录数组double *m_RunoffArray;double *m_EvapArray;时间（年月）记录数组long *m_YearArray;long *m_MonthArray;计算结果double *m_OptVolumnArray;double *m_OptNArray;double *m_OptAbandonQ;阶段数存储量时间间隔T=2626560s一个月的小时数：h库容离散点数库容离散阶段指标函数/目标函数长系列径流蒸发量/最优轨迹对应库容数据数组序列/最优轨迹对应出力数据数组序列/最优轨迹对应的弃水量序列d

3、ouble *m_OptUpWaterLevelArray; 最优轨迹对应的水库上游水位序列 double *m_OptDownWaterLevelArray;/最优轨迹对应的水库下游水位序列 double *m_OptVarialbe_q;最优轨迹对应的发电流量序列最优轨迹对应的预想出力最优轨迹中低于保证出力时的破坏深度/最优轨迹对应的保证率double *m_OptMeanWaterHeadArray; /最优轨迹对应的平均发电水头 double *m_OptNmaxArray;double *m_OptDamageDepth;double m_OptReliability;double

4、m_OptGenerationProduction; 发电量double m_OptAnnualMaxPower;/Parameters of Hydropower Station/死水位m/死库容.03亿m3正常蓄水位835m汛期防洪限制水位m汛期最大库容.03亿m3/double m_DeadWaterLevel;double m_DeadVolumn;double m_NormalWaterLevel;double m_FloodLimitLevel;double m_FloodVmax;非汛期最大库容.37亿m3机组台数装机容量*2MW/保证出力.2MWdouble m_NFloodV

5、max;int m_UnitsNum;double m_InstallCapacity;double m_FirmCapacity_N;预想出力double m_PotentialWaterLevel1;double m_PotentialCapacity_N1;double m_PotentialWaterLevel2;double m_PotentialCapacity_N2;double m_DesignAnnualGeneration;设计年发电量.99xkWe V);上游水位计算double DownWaterLevelCalculation(double q); 下游水位计算voi

6、d ReverseSearchCalculation();反向寻优计算void TerminalOptionalResultCalculation();/最终结果计算;#endif类的实现/ SRDP.cpp: implementation of the SRDP class./#include stdafx.h#include SRDP.h#include #include assert.h#include fstream#include iomanip#include iostream#include using namespace std;/ Construction/Destructi

7、on/SRDP:SRDP() 原始资料长系列径流数据/蒸发数据this - m_RunoffArray = NULL;this-m_EvapArray = NULL;/动态规划模型数据阶段数时间间隔一个月的小时数/状态数阶段指标/目标函数时间（年）数组时间（月）数组/最优轨迹记录/最优轨迹对应的预想出力/最优轨迹中低于保证出力时的破坏深度/this-m_StageNum=100;/this-m_Delta_t=2626560;/this-m_DeltaHour=730;this-m_StateNum = NULL;this-m_TarFun = NULL;/this-m_ObjFun=0.00

8、;this-m_YearArray = NULL;this - m_MonthArray = NULL;this-m_OptArray = NULL;this - m_OptNmaxArray = NULL;this-m_OptDamageDepth = NULL;SRDP:SRDP()原始初始资料数据析构delete this- m_RunoffArray;delete this-m_EvapArray;中间变量析构delete this-m_StateNum;delete this-m_YearArray;delete this-m_MonthArray;计算结果数组析构delete th

9、is-m_OptVolumnArray;delete this-m_OptNArray;delete this-m_OptAbandonQ;delete this-m_OptUpWaterLevelArray;/最优轨迹对应的水库上游水位序列delete this-m_OptDownWaterLevelArray; 最优轨迹对应的水库上游水位序列 delete this-m_OptVarialbe_q;/最优轨迹对应的发电流量序列/最优轨迹对应的预想出力/最优轨迹中低于保证出力时的破坏深度delete this-m_OptMeanWaterHeadArray;delete this-m_Opt

10、NmaxArray;delete this-m_OptDamageDepth;for(int i=0;im_StageNum;i+)/delete this-m_OptDeciVar_qi;/delete this-m_StateVi;delete this-m_OptArrayi;delete this-m_TarFuni;/delete this-m_OptDeciVar_q;/delete this-m_StateV;delete this-m_OptArray;delete this-m_TarFun;/数据初始化内存分配/void SRDP:InitionalizingSRDP(ch

11、ar *DataFile1, char *DataFile2,char *DataFile3)this-m_StageNum = 516;this-m_BeginMonth=6;this-m_BeginYear=1957;this-m_OptGenerationProduction=0.00;this-m_OptAnnualMaxPower = 0.00;/this-m_DeltaHour=730;/获取时间序列this-m_YearArray = new longthis-m_StageNum+1;this- m_MonthArray = new longthis-m_StageNum+1;

12、this-m_StateNum = new longthis-m_StageNum+1;for(int i=0;im_StageNum;i+)this-m_YearArrayi = 0;this-m_MonthArrayi = 0;this-m_StateNumi = 10;this-m_YearArray0=this-m_BeginYear;/水文年起始年份this-m_MonthArray0=this-m_BeginMonth; /水文年起始月份for(int i = 1;im_StageNum;i+ )if(this-m_MonthArrayi-1+1)m_MonthArrayi = t

13、his-m_MonthArrayi-1+1;this-m_YearArrayi = this-m_YearArrayi-1;else if(this-m_MonthArrayi-1+1)12)this-m_MonthArrayi = this-m_MonthArrayi-1-11;this-m_YearArrayi = this-m_YearArrayi-1+1;/coutvvivv m_YearArrayi-m_MonthArrayiendl;/获取状态离散点数/for(int i=0;im_StageNum;i+)if(this-m_MonthArrayi=6&this-m_MonthAr

14、rayim_StateNumi=1;elsethis-m_StateNumi=51;/coutvvi+1vvtvvthis-m_YearArrayivvVvvthis-m_MonthArrayivvtvvthis-m_StateNumivve ndl;/动态分配内存并初始化径流资料、蒸发数据、状态数据等/一维数组this-m_RunoffArray = new doublethis-m_StageNum+1;this-m_EvapArray = new doublethis-m_StageNum+1;/this-m_OptVolumnArray = new doublethis-m_Stage

15、Num+1;this-m_OptNArray = new doublethis-m_StageNum+1;this-m_OptAbandonQ = new doublethis-m_StageNum+1;this-m_OptUpWaterLevelArray = new doublethis-m_StageNum+1;/最优轨迹对应的水库上游水位序列this-m_OptDownWaterLevelArray = new doublethis-m_StageNum+1; /最优轨迹对应的水库上游水位 序列this-m_OptVarialbe_q = new doublethis-m_StageN

16、um+1;/最优轨迹对应的发电流量序列this-m_OptMeanWaterHeadArray= new doublethis-m_StageNum+1;this-m_OptNmaxArray= new doublethis-m_StageNum+1;/最优轨迹对应的预想出力this-m_OptDamageDepth= new doublethis-m_StageNum+1;/最优轨迹中低于保证出力时的破坏深度/二维数组/this-m_OptDeciVar_q = new double*this-m_StageNum+1;/this-m_StateV = new double*this-m_S

17、tageNum+1;this-m_OptArray = new long*this-m_StageNum+1;this-m_TarFun = new double*this-m_StageNum+1;for(int i=0;im_StageNum;i+)/this-m_OptDeciVar_qi = new doublethis-m_StateNumi;/this-m_StateVi = new doublethis-m_StateNumi;this-m_OptArrayi = new longthis-m_StateNumi;this-m_TarFuni = new doublethis-m

18、_StateNumi;for(int i=0;im_StageNum;i+)this-m_RunoffArrayi = 0.00;this-m_EvapArrayi = 0.00;this-m_OptVolumnArrayi = 0.00;this-m_OptNArray i = 0.00;this-m_OptAbandonQi=0.00;this-m_OptUpWaterLevelArrayi = 0.00;this-m_OptDownWaterLevelArrayi = 0.00;this-m_OptVarialbe_qi = 0.00;this-m_OptMeanWaterHeadArr

19、ayi = -1;this-m_OptNmaxArrayi= -1;最优轨迹对应的预想出力this-m_OptDamageDepthi= 0;/最优轨迹中低于保证出力时的破坏深度for(long j=0;jm_StateNumi;j+)/this-m_OptDeciVar_qij = 0.00;/this-m_StateVij =0.00;this-m_OptArrayij = 0;this-m_TarFunij = -1e+6;/Input System Parameteters from file/ifstream myFileIn1,myFileIn2,myFileIn3;myFileI

20、n1.open(DataFile1);系统特性参数文件if(!myFileIn1)coutCannot Open System Parameters File this-m_DeadWaterLevelthis-m_DeadVolumnthis-m_NormalWaterLevelthis-m_FloodLimitLevel this-m_FloodVmaxthis-m_NFloodVmax this-m_UnitsNumthis-m_InstallCapacitythis-m_FirmCapacity_Nthis-m_PotentialWaterLevel1this-m_PotentialC

21、apacity_N1this-m_PotentialWaterLevel2this-m_PotentialCapacity_N2this-m_DesignAnnualGenerationthis-m_CCthis-m_LHC this-m_LHmaxthis-m_Delta_t;myFileIn1.close();/*cout”系统特性参数endl;cout|时段长度cout|死水位cout|死库容cout|正常蓄水位cout |汛期防洪限制水位cout|汛期最大库容cout|非汛期最大库容 cout|系统装机容量cout|=|endl;| Delta_T = m_Delta_tsendl;|

22、 DeadWaterLevel = m_DeadWaterLevelmendl;| DeadVolumn=m_DeadVolumn亿立方米endl;| NormalWaterLevel = m_NormalWaterLevelmendl;| FloodLimitLevel = m_FloodLimitLevelmendl;| FloodPeriodVmax = m_FloodVmax亿立方米endl;| NonFloodPeriodVmax = m_NFloodVmax 亿立方米endl;| InstallCapacity = m_UnitsNumxm_InstallCapacity/2MWe

23、ndl;cout| 系统保证出力| FirmCapacity = m_FirmCapacity_NMWendl;cout|水电站预想出力|;cout| 水头|setw(6)m_PotentialWaterLevel1m|setw(4)m_PotentialWaterLevel2m|endl;cout|出力 |m_PotentialCapacity_N1MW|m_PotentialCapacity_N2MW|endl;cout| 设计年发电量| DesignAnnualGeneration = m_DesignAnnualGenerationkW.hendl;cout|水电站出力系数cout|水

24、头损失系数cout|最大水头损失cout|径流系列起始年月| OutputCoefficient=m_CCendl;| HeadLossCoefficient=m_LHCendl;| HeadLossMax=m_LHmaxendl;| BeginTime=m_BeginYear-0m_BeginMonthendl;cout|=|endl; */Input runoff data from file/从文件中读取长系列径流数据 myFileIn2.open(DataFile2);/径 流文件 if(!myFileIn2)coutCannot Open Runoff_data File ; exi

25、t(0);elsewhile(myFileIn2.good()for(int i=0;im_StageNum;i+)myFileIn2m_RunoffArrayi;/*coutsetiosflags(ios:left)setw(6)m_RunoffArrayi;if(i+1)%12=0) coutendl; */ break;myFileIn2.close();/Input evaporation data from file/从文件中读取蒸发数据myFileIn3.open(DataFile3);/ 蒸发数据文件if(!myFileIn3)coutCannot Open Evaporatio

26、n_data File ;exit(0);elsewhile(myFileIn3.good()for(int i=0;im_StageNum;i+)myFileIn3this-m_EvapArrayi;if(0=(i+1)%12)myFileIn3.close();myFileIn3.open(DataFile3);break;myFileIn3.close();/运行计算/void SRDP:RunningSRDP()反向逆推计算寻优this-ReverseSearchCalculation();正向顺推计算this-TerminalOptionalResultCalculation();/

27、保存计算结果/void SRDP:SaveCaluResults(char *FileName)ofstream myFileOut;myFileOut.open(FileName);/FileNam e 是保存数据的文件名myFileOutsetiosflags(ios:fixed)setprecision(3);myFileOut=endl;myFileOutsetw(35)”系统特性参数endl;myFileOut=endl;myFileOut时段长度m_Delta_tsendl;myFileOut死水位m_DeadWaterLevelmendl;myFileOut死库容m_DeadVo

28、lumn亿立方米endl;myFileOut正常蓄水位m_NormalWaterLevelmendl;myFileOut汛期防洪限制水位m_FloodLimitLevelmendl;myFileOut汛期最大库容”m_FloodVmax”亿立方米endl;myFileOut非汛期最大库容m_NFloodVmax亿立方米endl;myFileOut系统装机容量m_UnitsNumx m_InstallCapacity/2MW endl;myFileOut系统保证出力m_FirmCapacity_NMWendl;myFileOut水电站预想出力水头(m)setw(6)m_PotentialWate

29、rLevel1setw(8)m_PotentialWaterLevel2endl;myFileOut出力(MW)setw(7)m_PotentialCapacity_N1setw(6)m_PotentialCapacity_N2endl;myFileOut设计年发电量m_DesignAnnualGenerationkW.hendl;myFileOut水电站出力系数m_CCendl;myFileOut水头损失系数m_LHCendl;myFileOut最大水头损失m_LHmaxendl;myFileOut径流系列起始年月m_BeginYear/0m_BeginMonthendl;myFileOut

30、=endl;myFileOutendl;myFileOut= =endl;myFileOutsetw(35)最优化计算结果endl;myFileOut=endl;myFileOutendl;myFileOutsetiosflags(ios:fixed)setprecision(2);myFileOut总发电量：m_OptGenerationProduction/100000000亿KW.hendl;myFileOut多年平均发电量：m_OptGenerationProduction/43.0/100000000”亿KW.hendl;myFileOut可靠保证率：m_OptReliability

31、*100%endl;myFileOut年最大发电量：m_OptAnnualMaxPower/100000000亿KW.hendl;myFileOutendl;myFileOutsetiosflags(ios:fixed)setiosflags(ios:right)setw(10)时间(年/月)setw(11)来流量”setw(8)蒸发量”setw(10)初库容”setw(10)末库容”setw(10)初水位”setw(10)末水位”setw(10)”下游水位”setw(10)发电水头”setw(10)发电流量setw(10)弃水量”setw(10)出力值”setw(13)预想出力值”setw(

32、11)”破坏深度”endl;for(int i=0;im_StageNum-1;i+)时间(年/月)输出myFileOutsetw(5)m_YearArrayim_MonthArrayi10)myFileOut0m_MonthArrayit;elsemyFileOutm_MonthArrayit;数据输出myFileOutsetiosflags(ios:fixed)setprecision(2)setw(6)m_RunoffArrayisetw(7)m_EvapArrayi;myFileOutsetiosflags(ios:fixed)setprecision(3)setiosflags(io

33、s:right)setw(11)m_OptVolumnArrayisetw(10)m_OptVolumnArrayi+1setw(10)m_OptUpWaterLevelArrayisetw(10)m_OptUpWaterLevelArrayi+1setw(10)m_OptDownWaterLevelArrayisetw(10)m_OptMeanWaterHeadArrayisetw(10)m_OptVarialbe_qisetw(10)m_OptAbandonQisetw(10)m_OptNArrayisetw(12)m_OptNmaxArrayisetw(10)m_OptDamageDep

34、thi;myFileOutsetw(5)m_YearArrayim_MonthArrayi10)myFileOut0 m_MonthArrayi;elsemyFileOutm_MonthArrayi;myFileOutendl;myFileOut.close();ofstream myfileout;myfileout.open(WaterLevel.txt);for(int i=0;im_StageNum-1;i+)myfileoutsetw(5)m_YearArrayim_MonthArrayi10)myfileout0m_MonthArrayit;elsemyfileoutm_Month

35、Arrayit; myfileoutm_OptUpWaterLevelArrayim_TarFunthis-m_StageNum0=0.00; for(int i=this-m_StageNum-1;i=0;i-) 对库容进行离散 if(1=this-m_StateNumi) Delta_V1 = this-m_FloodVmax - this-m_DeadVolumn; elseDelta_V1 = (this-m_NFloodVmax - this-m_DeadVolumn)/(this-m_StateNumi-1); if(1=this-m_StateNumi+1)Delta_V2 =

36、this-m_FloodVmax - this-m_DeadVolumn; elseDelta_V2 = (this-m_NFloodVmax - this-m_DeadVolumn )/(this-m_StateNumi+1-1);/*coutvvivvsetw(4)vvthis-m_MonthArrayivvsetw(10)vvDelta_V1vvsetw(4)vvthis-m_StateNumivv setw(10)vvthis-m_MonthArrayi+1vvsetw(10)vvDelta_V2vvsetw(4)vvthis-m_StateNumi+1vvendl;*/ 进行逆推迭代

37、计算for(int j=0;jm_StateNumi;j+) V1 = this-m_DeadVolumn+j*Delta_V1;单位：亿立方米for(int k=0;km_StateNumi+1;k+)V2 = this-m_DeadVolumn+k*Delta_V2;单位：亿立方米水量平衡验算求出发电流量q=this-m_RunoffArrayi-(V2 -V1)*100000000/this-m_Delta_t-this-m_EvapArrayi;/ 单 位：立方米/秒if(qUpWaterLevelCalculation(V1)+this-UpWaterLevelCalculation

38、(V2)/2.0; 水位流量关系曲线求下游水位值Z_down = this-DownWaterLevelCalculation(q); 水头损失计算和发电水头计算 HeadLoss= this-m_LHC*q*q; if(HeadLossthis-m_LHmax)HeadLoss = this-m_LHmax; GWH = Z_up-Z_down-HeadLoss; if(GWH=0)N = 0.00; else 最大出力值计算if(GWH=57) Nmax = 120.0;elseNmax=2.94*GWH-47.58;/coutsetw(10)Nmaxm_CC*q*GWH/1000;if(

39、Nm_FirmCapacity_N)/低于保证出力 N=N-2000*(this-m_FirmCapacity_N-N); else if(N=Nmax)/大于最大出力限制时按最大出力计算N = Nmax; /couti”t”j”t”k”t”Nm_TarFuni+1k;if(tempTarFun=this-m_TarFunij)this-m_TarFunij=tempTarFun;this-m_OptArrayij = k;this-m_ObjFun = this-m_TarFun00;/coutm_ObjFunm_OptVolumnArray0=this-m_DeadVolumn;for(i

40、nt i=0;im_StageNum;i+)this-m_OptUpWaterLevelArrayi=this-UpWaterLevelCalculation(this-m_OptVolumnArrayi);对库容进行离散if(1=this-m_StateNumi+1)Delta_V = this-m_FloodVmax - this-m_DeadVolumn;elseDelta_V = (this-m_NFloodVmax - this-m_DeadVolumn)/(this-m_StateNumi+1-1);m=this-m_OptArrayim;this-m_OptVolumnArray

41、i+1=this-m_DeadVolumn+m*Delta_V;时段初末水位值计算Z1=this-UpWaterLevelCalculation(this-m_OptVolumnArrayi);Z2=this-UpWaterLevelCalculation(this-m_OptVolumnArrayi+1);水量平衡计算单位：立方米/秒q=this-m_RunoffArrayi-this-m_EvapArrayi-(this-m_OptVolumnArrayi+1-this-m_OptVolumnArrayi)*100000000/this-m_Delta_t;if(qm_OptNArrayi

42、=0.00;this-m_OptVarialbe_qi=0.00;this-m_OptAbandonQi=0.00;this-m_OptDamageDepthi=-10000;count+; else /平均发电水头计算Z_up = (Z1+Z2)/2;Z_down = this-DownWaterLevelCalculation(q);Delta_H = this-m_LHC*q*q;if(Delta_H5.59)Delta_H = 5.59;GWH = Z_up-Z_down-Delta_H;预想出力值计算if(GWH=57)Nmax = 120.0;elseNmax=2.94*GWH-4

43、7.58;N=this-m_CC*q*GWH/1000; 出力计算if(NNmax) 出力大于最大出力double q0,q1,q2;q0=q; q1=0.00; q2=q;while(abs(N-Nmax)0.0001)q=(q1+q2)/2;Z_down = this-DownWaterLevelCalculation(q);Delta_H = this - m_LHC*q*q;if(Delta_H5.59)Delta_H = 5.59;GWH = Z_up-Z_down-Delta_H;/预想出力值计算if(GWH=57)Nmax = 120.0;elseNmax=2.94*GWH-47

44、.58;N=this-m_CC*q*GWH/1000; 出力计算if(NNmax)q2=q;elseq1=q;计算结果this-m_OptNArrayi=N;this-m_OptAbandonQi=q0-q;this-m_OptVarialbe_qi=q;/最优轨迹对应出力数据数组序列/最优轨迹对应的弃水量序列/最优轨迹对应的发电流量序列this-m_OptMeanWaterHeadArrayi=GWH; /最优轨迹对应的平均发电水头 this-m_OptNmaxArrayi=Nmax;/最优轨迹对应的预想出力this-m_OptDownWaterLevelArrayi=Z_down最优轨迹对应的水库下游水位序列 else if(Nm_FirmCapa