gould.cpp

古尔德法

 /*
 x(i,j)--水库入流（年，月径流量）
 A(i)--蓄水状态变量                            
 S(i,j)--月初蓄水量 
 S1(i)--年末蓄水量
 Q(i,j)--年转移概率矩阵
 IND(j)--年末出现小于等于某一状态的总次数
 PI(i)--水库稳定蓄水概率分布
 FY(i)--条件供水破坏概率（相对破坏年数）
 FM(i)--条件供水破坏概率（相对破坏月数）
 QP--调节流量
 Vn--有效库容
 M--状态总数
 K1--计算总年教
 PY--年供水保证率
 PM--月供水保证率
*/
#include<math.h>
#include<fstream.h>
#define  M    60
#define  K1   49
#define  m    18
void main()
{
/*状态离散*/
	int i,j,k,l,e;
	float Z[M],V[M],A[m];
	float Z0=162.0,Zn=207.5;
	float V0,Vn,DV;
	ifstream ZV("D:\\水位库容关系.txt");
		for(i=0;i<M;i++)
			ZV>>Z[i]>>V[i];
    ZV.close();
 	for(i=0;i<M-1;i++)
	{
        if(Z0>=Z[i]&&Z0<Z[i+1])
			V0=V[i]+(V[i+1]-V[i])*(Z0-Z[i])/(Z[i+1]-Z[i]);
		if(Zn>=Z[i]&&Zn<Z[i+1])
			Vn=V[i]+(V[i+1]-V[i])*(Zn-Z[i])/(Z[i+1]-Z[i]);
	}
	DV=Vn/(m-2);
	A[0]=0;	
	A[m-1]=Vn;
	for(i=1;i<m-1;i++)
		A[i]=DV*i-DV/2;
/*年末蓄水状态及推求条件破坏概率*/
	float QP=7.8;
	float FY[m],FM[m],IND[m];
	float S[K1][13],x[K1][12],S1[K1],Q[m][m];
	int   IY,IM,INF;
    ifstream Qin("D:\\径流资料.txt");
	for(i=0;i<K1;i++)
		for(j=0;j<12;j++)
			Qin>>x[i][j];
	Qin.close();
		for(i=0;i<m;i++)
	{
		IY=0;
		IM=0;
		for(k=0;k<K1;k++)
		{
				S[k][0]=A[i];
			    INF=0;
				for(l=0;l<12;l++)
				{
					S[k][l+1]=S[k][l]+x[k][l]*0.0263-QP;
			    	if(S[k][l+1]<=0)
					{
				     	S[k][l+1]=0;
					    INF=INF+1;
					}
			    	else if(S[k][l+1]>=Vn)
				    	S[k][l+1]=Vn;
				}
			    if(INF>0)
				{
				   IY=IY+1;
				   IM=IM+INF;
				}
			S1[k]=S[k][12];
		}
			FY[i]=IY/float(K1);
			FM[i]=IM/float(K1*12);
/*推求水库蓄水状态的年转移概率*/
		for(j=0;j<m;j++)
			IND[j]=0;
		for(k=0;k<K1;k++)
		{
			if(S1[k]==0.0)
				IND[0]+=1;
			else if(S1[k]==Vn)
				IND[m-1]+=1;
			else
				{
					e=int(S1[k]/DV)+1;
                    IND[e]+=1;
				}
		}
		for(j=0;j<m;j++)
				Q[i][j]=IND[j]/float(K1);
	}
/*推求水库稳定的蓄水概率分布*/
		float PAI[m],PI[m];
		for(j=0;j<m;j++)
		{
			if(j==0)
				PAI[j]=1.0;
			else
				PAI[j]=0.0;
		}
		for(j=0;j<m;j++)
		{
			PI[j]=0.0;
			for(i=0;i<m;i++)
				PI[j]=PI[j]+PAI[i]*Q[i][j];
    		for(e=0;e<m;e++)
			if(fabs(PI[e]-PAI[e])>=0.001)
				PAI[j]=PI[j];
		}
/*计算水库正常供水保证率*/
		float HY=0.0,HM=0.0;
		float PY,PM;
		for(i=0;i<m;i++)
		{
			HY+=PI[i]*FY[i];
			HM+=PI[i]*FM[i];
		}
		PY=1-HY;
		PM=1-HM;
/*输出文件*/
		ofstream file1("D:\\状态转移矩阵.txt");
			for(i=0;i<m;i++)
			{
				for(j=0;j<m;j++)
					file1<<Q[i][j]<<"\t";
				    file1<<"\n";
			}
			file1.close();
		ofstream file2("D:\\条件破坏概率.txt");
			for(i=0;i<m;i++)
				file2<<FY[i]<<"\t"<<FM[i]<<endl;
			file2.close();
		ofstream file3("D:\\水库稳定蓄水概率分布.txt");
			for(i=0;i<m;i++)
				file3<<PI[i]<<endl;
			file3.close();
		ofstream file4("D:\\状态点各值.txt");
		    for(i=0;i<m;i++)
			    file4<<A[i]<<endl;
			file4.close();
		ofstream file5("D:\\正常供水保证率.txt");
		    file5<<PY<<"\t"<<PM<<endl;
			file5.close();
}