#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";
// }
// }