水库优化调度动态规划程序

#include<iostream> #include<cmath> #include<iomanip> #include<fstream> #include<cstdlib> #include<string> using namespace std; const int N = 100; const int M = 372; const int M1= 6; const int M2= 6; float chazhi(float x0,float y1[],float y2[],int num) { int m1; for(m1=1;m1<=num-1;m1++) if(x0<=y1[m1]&&x0>y1[m1-1]) return(y2[m1 - 1] + (x0 - y1[m1 - 1]) * (y2[m1] - y2[m1 - 1]) / (y1[m1] -y1[m1 - 1])); if(x0<=y1[0]) return(y2[0]); if(x0>y1[num-1]) return(y2[num-2] + (x0- y1[num-2]) * (y2[num-1] - y2[num-2]) / (y1[num-1] - y1[num-2])); } int main() { //-----------变量的定义--------------------------------------------------- int cfxs,i,j,k,js; int num1=M,num2=M1,num3=M2; float vmax,vmin,sume,xxz1,xxz2,zcz,vjun,bv,zdq,zxq,zcv,xxv1,xxv2,sv,Nzh,Np,max[N]={0},Q,clxs; float qt[M],sz,abc; //来流量 float yv[M][N],yz[M][N],yzx[M][N],fz[M1],fv[M1]; //水库库容水位，水位库容曲线 float fqx[M2],fzx[M2],yqx[M][N],yh[M][N],yr[M][N],yw[M][N]={0},ye[M][N];//下泄流量曲线，各阶段出力 int ml[M][N],xx[M][N],xxx[N][N]; float c0=400.0; //惩罚系数 float le[N][N],lq[N][N],lzx[N][N],lh[N][N],lr[N][N],lw[N][N],lyz[N][N];//临时阶段变量 //-------输入输出数据--------------- char *fname="result.txt"; //physical file name ofstream fout(fname,ios::out); //opening file for output if( !fout ) //handles errors in file opening { cerr<<"Error in opening file for output!"; exit(1); } string ifname; ifstream fin; cout<<"Enter the filename of in file"<<endl; cin>>ifname; fin.open(ifname.c_str()); if( fin.fail() ) { cout<<"ERROR open the in file "<<endl; exit(1); } else { for( i = 0; i < num1; i++ ) fin>>qt[i]; //来流量 for(i=0;i<num2;i++) fin>>fz[i]; for(i=0;i<num2;i++) fin>>fv[i]; //上游库容水位 for(i=0;i<num3;i++) fin>>fzx[i]; for(i=0;i<num3;i++) fin>>fqx[i]; //下游水位和下泄流量 fin>>zdq>>zxq>>Nzh>>Np>>xxz1>>xxz2>>zcz>>sz>>clxs; }//end open infile and intialize tha array x[] fin.close(); //--------------------------------------------------------------------------------------- sv=chazhi(sz,fz,fv,num2); xxv1=chazhi(xxz1,fz,fv,num2); xxv2=chazhi(xxz2,fz,fv,num2); zcv=chazhi(zcz,fz,fv,num2); for(i=0;i<num1;i++) { if (((i+1)%12)<3) { if(xxv1>sv) { vmax = xxv1; vmin=sv; } else { vmax=xxv1; vmin=xxv1; } } if (((i+1)%12)==3) { if(xxv2>sv) { vmax = xxv2; vmin=sv; } else { vmax=xxv2; vmin=xxv2; } } if (((i+1)%12)>=4) { vmax = zcv; vmin=sv; } if(i==num1-1) { if(xxv1>sv) vmax=vmin=sv; else vmax=vmin=xxv1; } vjun=(vmax-vmin)/(N-1); for(j=0;j<N;j++) { yv[i][j]=vmin+j*vjun; //各阶段各状态库容 } } //--------------------------------------------第一阶段 从死水位开始 for(i=0,j=0;j<N;j++) { if(xxv1>sv) vmin=sv; else vmin=xxv1; cfxs=0; bv=yv[i][j]/2.0+vmin/2.0; //第一阶段平均库容 yz[i][j]=chazhi(bv,fv,fz,num2); //第一阶段平均上游水位 yqx[i][j]=qt[i]+(vmin-yv[i][j])*100.0/2.63; //第一阶段平均下泄流量 if(yqx[i][j]>zxq&&yqx[i][j]<=zdq) { xx[i][j]=1; yzx[i][j]=chazhi(yqx[i][j],fqx,fzx,num3);//下游尾水位 yh[i][j]=yz[i][j]-yzx[i][j]; //毛水头 yr[i][j]=clxs*yqx[i][j]*yh[i][j]/1000; //第一阶段出力 MW if(yr[i][j]>Nzh) { yr[i][j]=Nzh; yw[i][j]=yqx[i][j]-yr[i][j]/clxs/yh[i][j]*1000; yqx[i][j]=yr[i][j]/clxs/yh[i][j]; } } if(yqx[i][j]>zdq) { xx[i][j]=1; yzx[i][j]=chazhi(yqx[i][j],fqx,fzx,num3);//下游尾水位 yh[i][j]=yz[i][j]-yzx[i][j]; //毛水头 Q=yqx[i][j]; yqx[i][j]=zdq; yr[i][j]=clxs*yqx[i][j]*yh[i][j]/1000; //各状态出力 if(yr[i][j]>Nzh) { yr[i][j]=Nzh; yw[i][j]=Q-yr[i][j]/clxs/yh[i][j]*1000; yqx[i][j]=yr[i][j]/clxs/yh[i][j]; } else yw[i][j]=Q-yqx[i][j]; } if(yqx[i][j]<=zxq) { yr[i][j]=-50000.0; xx[i][j]=0; } if(yr[i][j]<Np) cfxs=1; ye[i][j]=yr[i][j]-cfxs*c0; //加上惩罚项的目标函数值 // ml[i][j]=0;//记录上一阶段的状态点 } //第二阶段到第M个阶段－－－－－－－－－－－－－－－－－－－－－－－－－－ for(i=1;i<num1;i++) { for(j=0;j<N;j++) { max[j]=-5000; for(k=0;k<N;k++) { cfxs=0; bv=yv[i][j]/2.0+yv[i-1][k]/2.0; //平均蓄水量 lyz[j][k]=chazhi(bv,fv,fz,num2); lq[j][k]=qt[i]+(yv[i-1][k]-yv[i][j])*100.0/2.63; if(lq[j][k]>zxq) { if(lq[j][k]<=zdq) { xxx[j][k]=1; lzx[j][k]=chazhi(lq[j][k],fqx,fzx,num3);//下游尾水位 lh[j][k]=lyz[j][k]-lzx[j][k]; //毛水头 lr[j][k]=clxs*lq[j][k]*lh[j][k]/1000; //出力 if(lr[j][k]>Nzh) { lr[j][k]=Nzh; lw[j][k]=lq[j][k]-lr[j][k]/clxs/lh[j][k]*1000; lq[j][k]=lr[j][k]/clxs/lh[j][k]*1000; } } if(lq[j][k]>zdq) { xxx[j][k]=1; lzx[j][k]=chazhi(lq[j][k],fqx,fzx,num3);//下游尾水位 lh[j][k]=lyz[j][k]-lzx[j][k]; //毛水头 Q=lq[j][k]; lq[j][k]=zdq; lr[j][k]=clxs*lq[j][k]*lh[j][k]/1000; //出力 if(lr[j][k]>Nzh) { lr[j][k]=Nzh; lw[j][k]=Q-lr[j][k]/clxs/lh[j][k]*1000; lq[j][k]=lr[j][k]/clxs/lh[j][k]*1000; } else yw[i][j]=Q-lq[j][k]; } } if(lq[j][k]<=zxq) { xxx[j][k]=0; lr[j][k]=-50000.0; } if(lr[j][k]<Np) cfxs=1; le[j][k]=lr[j][k]+ye[i-1][k]-cfxs*c0; if(max[j]<=le[j][k]&&xx[i-1][k]==1&&xxx[j][k]==1) //此处有错 { max[j]=le[j][k]; ye[i][j]=le[j][k]; //总出力 yr[i][j]=lr[j][k]; //阶段出力 yqx[i][j]=lq[j][k]; yw[i][j]=lw[j][k]; yzx[i][j]=lzx[j][k]; yz[i][j]=lyz[j][k]; yh[i][j]=lh[j][k]; ml[i][j]=k; //该阶段的状态点 xx[i][j]=1; } } } } //------------------------------------------------------------ bv=0; js=0; sume=0; float e[M],r[M]; float qx[M],w[M],zx[M],zs[M],z[M+1],h[M],v[M+1]; //z[M]是阶段水位 for(i=num1-1;i>=0;i--) { if (i==num1-1 ) { e[i]=ye[i][0]; r[i]=yr[i][0]; qx[i]=yqx[i][0]; w[i]=yw[i][0]; zx[i]=yzx[i][0]; zs[i]=yz[i][0]; h[i]=yh[i][0]; v[i+1]=yv[i][0]; z[i+1]=chazhi(v[i+1],fv,fz,num2); } if( i==num1-2 ) { j = ml[i + 1][0]; e[i] = e[i + 1] - r[i + 1]; r[i] = yr[i][j]; qx[i] = yqx[i][j]; w[i] = yw[i][j]; zx[i] = yzx[i][j]; zs[i] = yz[i][j]; h[i]=yh[i][j]; v[i+1 ] = yv[i][j]; z[i+1]=chazhi(v[i+1],fv,fz,num2); } if (i < num1-2) { j = ml[i + 1][j]; e[i] = e[i + 1] - r[i + 1]; r[i] = yr[i][j]; qx[i] = yqx[i][j]; w[i] = yw[i][j]; zx[i] = yzx[i][j]; zs[i] = yz[i][j]; h[i]=yh[i][j]; v[i+1 ] = yv[i][j]; z[i+1]=chazhi(v[i+1],fv,fz,num2); } sume=sume+r[i]; v[0]=xxv1; z[0]=xxz1; if( r[i]>=Np) js=js+1; } float phy,Pj,phc,Epj; if(sume>0) { Pj=sume*1.0/num1; phc=js*1.0/num1; Epj=sume*365*24/num1; } else cout<<"无可行解,请重新输入"; for(i=1;i<num1;i++) { if(qt[i]==qx[i]) w[i]=0; if(qt[i]!=qx[i]) { w[i]=qt[i]-qx[i]-(v[i+1]-v[i])*100/2.63; if(w[i]<0.1) w[i]=0; } } fout<<"平均出力值P="<<Pj<<"\n"<<"出力保证率p="<<phc<<"\n"<<"多年平均发电量E="<<Epj<<"\n"; fout<<"qt"<<"\t"<<"qx"<<"\t"<<"w"<<"\t"<<"v"<<"\t"<<"z"<<"\t"<<"zs"<<"\t"<<"zx"<<"\t"<<"h"<<"\t"<<"r"<<"\n"; for(i=0;i<num1;i++) { if(i%12==0) fout<<"\n"; fout<<qt[i]<<"\t"<<qx[i]<<"\t"<<w[i]<<"\t"<<v[i+1]<<"\t"<<z[i+1]<<"\t"<<zs[i]<<"\t"<<zx[i]<<"\t"<<h[i]<<"\t"<<r[i]<<"\n"; } // for(i=0;i<M;i++) // { // fout<<"\n"; // for(j=0;j<N;j++) // { // fout<<ye[i][j]<<"\t"; // } // }