基于简化运行策略的水电站中长期常规调度

//毕设之常规调度程序---简化运行
#include<stdio.h>
#include<stdlib.h>
#include<math.h>

#define z_ss  328.5    //死水位
#define z_qfh  346.30   //防洪限制水位-前汛期
#define z_hfh  346.80   //防洪限制水位-后汛期
#define z_zc  347      //正常蓄水位
#define Y   48        //流量资料年数
#define M   12        //一年月数
#define T   576       //入流流量数据个数
#define NUM1   24        //水位库容关系点数
#define NUM2   12        //下泄流量与水位关系点数
#define NUM3   3         //水头限制出力关系点数

//插值函数部分
float chazhi(float x,float y1[],float y2[],int n)  //插值函数
 {
     int i;
     for(i=0;i<n-1;i++)
        {
            if(x>=y1[i]&&x<=y1[i+1])
                return (y2[i]+(y2[i+1]-y2[i])/(y1[i+1]-y1[i])*(x-y1[i]));
        }
        if(x<y1[0])
            return y1[0];
        if(x>y1[n-1])
            return (y2[n-2]+(y2[n-1]-y2[n-2])/(y1[n-1]-y1[n-2])*(x-y1[n-2]));
 }

int main(void)
{
    float Q[T+1],vcz[NUM1],zsy_cz[NUM1],qcz[NUM2],zxy_cz[NUM2],z_mao[NUM3],n_yx[NUM3];
    int i,j;
    int n=1,m=1;
    Q[0]=0;

//数据输入部分
    FILE *fQ1i;      //流量数据输入
    if((fQ1i=fopen("流量.txt","r"))==NULL)
 {
  printf("can't open 流量.txt file\n");
  exit(0);
 }
 for(i=0;i<Y;i++)
 {
  for(j=0;j<M;j++)
  {
   fscanf(fQ1i,"%f",&Q[n++]);
  }
 }

 FILE *fQsw;
 fQsw=fopen("水文年排列的流量.txt","w");//将原始流量数据按照水文年排列
 for(i=5;i<T+1;i++)
    {
        fprintf(fQsw,"%6.2f\t",Q[i]);
        if((i-4)%12==0)
            fprintf(fQsw,"\n");
    }
    for(i=1;i<5;i++)
    {
        fprintf(fQsw,"%6.2f\t",Q[i]);
    }

    FILE *fQi;      //水文年排列的流量数据输入
    if((fQi=fopen("水文年排列的流量.txt","r"))==NULL)
 {
  printf("can't open 水文年排列的流量.txt file\n");
  exit(0);
 }
 for(i=0;i<Y;i++)
 {
  for(j=0;j<M;j++)
  {
   fscanf(fQi,"%f",&Q[m++]);
  }
 }

 FILE *fv_z;    //水位-库容数据输入
 if((fv_z=fopen("水位库容.txt","r"))==NULL)
 {
  printf("can't open 水位库容.txt file\n");
  exit(0);
 }
 for(i=0;i<NUM1;i++)
 {
   fscanf(fv_z,"%f%f",&zsy_cz[i],&vcz[i]);
   //printf("zsy_cz[%d]=%f\tvcz[%d]=%f\n",i,zsy_cz[i],i,vcz[i]);
 }

 FILE *fq_z;   //水位-尾流量数据输入
 if((fq_z=fopen("水位流量.txt","r"))==NULL)
 {
  printf("can't open 水位流量.txt file\n");
  exit(0);
 }
 for(i=0;i<NUM2;i++)
 {
        fscanf(fq_z,"%f%f",&zxy_cz[i],&qcz[i]);
        //printf("zxy_cz[%d]=%f\tqcz[%d]=%f\n",i,zxy_cz[i],i,qcz[i]);
 }

 FILE *fn_z;//水头出力数据导入
 if((fn_z=fopen("水头出力.txt","r"))==NULL)
 {
     printf("can't open 水头出力.txt file\n");
     exit(0);
 }
 for(i=0;i<NUM3;i++)
        {
            fscanf(fn_z,"%f""%f",&z_mao[i],&n_yx[i]);
            //printf("z_mao[%d]=%f\tn_yx[%d]=%f\n",i,z_mao[i],i,n_yx[i]);
        }
/*//库容系数计算
 float v_zc,v_sk;    //正常蓄水位对应的库容
 float v_xl;    //兴利库容
 float Q_ave;   //多年平均流量
    float w_yong;  //多年平均水量
 float Q_total=0;//多年平均流量之和
 float b;  //库容系数
 for(i=1;i<T+1;i++)
    {
        Q_total+=Q[i];
    }
    Q_ave=Q_total/T*12;
    w_yong=Q_ave*2.63*pow(10,-2);
    v_sk=chazhi(z_ss,zsy_cz,vcz,NUM1);
    v_zc=chazhi(z_zc,zsy_cz,vcz,NUM1);
    v_xl=v_zc-v_sk;
    b=v_xl/w_yong;
    printf("多年平均流量Q_ave=%f\n",Q_ave);
    printf("多年平均年径流量Q_ave=%f\n",w_yong);
    printf("多兴利库容v_xl=%f\n",v_xl);
    printf("库容系数b=%f\n",b);
    if(b>0.3)
        printf("此水库属于多年调节水库\n");
    else if(b<=0.3&&b>=0.08)
        printf("此水库属于年调节水库\n");
    else if(b<0.08&&b>=0.03)
        printf("此水库属于不完全年调节水库\n");
    else
        printf("此水库属于日调节水库\n");*/

//常规调度部分
    float Qfd[T+1],Qfd2[T+1],Qfd3;
    float Qfd1[T+1]={0};//试算发电流量部分
    float zsy[T+1],zxy[T+1];//上游水位
    float dzs[T+1],dz[T+1],dzs1[T+1];//水头损失
    float q_qs[T+1]={0};//弃水
    float v[T+2];//库容
    float r=0.00001;//检验值
    float s=0.3;//加权数
    float k=7.74;//出力系数
    float NP=1.51;//保证出力
    float N[T+1];//出力序列
    float Vmax,Vmin;//库容上下限
    float Nyu;//预想出力
    int n_yt[T+1];//每月天数
    float  yn=0;//破坏的年份数
    float b_y;//保证率
    float E[T+1];//发电量

    v[1]=chazhi(328.5,zsy_cz,vcz,NUM1);//定义库容起始值
    Vmin=chazhi(z_ss,zsy_cz,vcz,NUM1);//计算库容下限
    for(i=1;i<T+1;i++)
    {
        if(i%12<=3&&i%12>=1)//汛期非汛期库容的限制值
            Vmax=chazhi(z_qfh,zsy_cz,vcz,NUM1);
        else if(i%12<=6&&i%12>=4)
            Vmax=chazhi(z_hfh,zsy_cz,vcz,NUM1);
        else
            Vmax=chazhi(z_zc,zsy_cz,vcz,NUM1);
       //确定每月小时数
       if(i%12==1||i%12==3||i%12==4||i%12==6||i%12==8||i%12==9||i%12==11)//每月31天的月份
           n_yt[i]=744;
       else if(i%12==10)//2月
           n_yt[i]=696;
       else//平月
           n_yt[i]=720;
        for(;;)
        {
            zsy[i]=chazhi(v[i],vcz,zsy_cz,NUM1);//插值计算上游水位
            zxy[i]=chazhi(Qfd1[i],qcz,zxy_cz,NUM2);//插值计算下游水位
            dzs1[i]=8.057*pow(10,-4)*Qfd1[i]*Qfd1[i]*4;//水头损失值

            if(dzs1[i]<7) //水头损失确认
                 dzs[i]=dzs1[i];
            else
                 dzs[i]=7;

            dz[i]=zsy[i]-zxy[i]-dzs[i];//实际水头
            Qfd[i]=NP*pow(10,4)/(k*dz[i]);//由保证出力及水头计算发电流量

            if(fabs(Qfd[i]-Qfd1[i])>=r)//判断假设的发电流量是否合理
                 Qfd1[i]=s*Qfd[i]+(1-s)*Qfd1[i];//重新定义发电流量的试算值
            else
                break;
        }
        Qfd2[i]=(Qfd[i]+Qfd1[i])/2;
        v[i+1]=v[i]+(Q[i]-Qfd2[i])*2.63/100;//由水量平衡计算末时段库容
        Nyu=chazhi(dz[i],z_mao,n_yx,NUM3);//预想出力限制部分

        if(v[i+1]>=Vmin&&v[i+1]<=Vmax)//库容为限制条件，计算不同情况下的实际发电流量
        {
            N[i]=NP;
            Qfd[i]=Qfd2[i];
        }
        else if (v[i+1]<Vmin)  //降低出力
        {
            Qfd[i]=Q[i]+(v[i]-Vmin)*100/2.63;
            N[i]=Qfd[i]*k*dz[i]*pow(10,-4);
            v[i+1]=v[i]+(Q[i]-Qfd[i])*2.63/100;
        }
        else
        {
            Qfd3=Q[i]+(v[i]-Vmax)*100/2.63;   //加大出力
            N[i]=Qfd3*k*dz[i]*pow(10,-4);
            if(N[i]>Nyu)
            {
                N[i]=Nyu;
                Qfd[i]=Nyu*pow(10,4)/(k*dz[i]);
                q_qs[i]=Qfd3-Qfd[i];
            }
            else
                Qfd[i]=Qfd3;
            v[i+1]=v[i]+(Q[i]-Qfd3)*2.63/100;
        }
        E[i]=N[i]*n_yt[i];
    }

    for(i=0;i<Y;i++)
    {
        for(j=1;j<M+1;j++)
        {
           if(N[i*12+j]<NP)
           {
               yn++;
               break;
           }
        }
    }
    b_y=(Y-yn)/Y;
    printf("%3.0f\n",yn);
    printf("%f",b_y);

    FILE *fZSY,*fVY,*fNS,*fQS,*fQQ,*fdz,*fFD;
 fZSY=fopen("常规调度下的上游水位.txt","w");    //上游水位
 fVY=fopen("常规调度下的库容.txt","w");    //库容
 fQS=fopen("常规调度下的发电流量.txt","w");    //发电流量
    fNS=fopen("常规调度下的出力.txt","w");       //出力
    fQQ=fopen("常规调度下的弃水.txt","w");    //弃水
    fdz=fopen("常规调度下的水头.txt","w");     //水头
    fFD=fopen("常规调度下的发电量.txt","w");   //发电量

 for(i=1;i<T+1;i++)
 {
       fprintf(fZSY,"%6.2f\t",zsy[i]);
       fprintf(fVY,"%6.2f\t",v[i]);
       fprintf(fQS,"%6.2f\t",Qfd[i]);
       fprintf(fNS,"%6.2f\t",N[i]);
       fprintf(fQQ,"%6.2f\t",q_qs[i]);
       fprintf(fdz,"%6.2f\t",dz[i]);
       fprintf(fFD,"%6.2f\t",E[i]);
       if(i%12==0)
       {
          fprintf(fZSY,"\n");
          fprintf(fVY,"\n");
          fprintf(fQS,"\n");
          fprintf(fNS,"\n");
          fprintf(fQQ,"\n");
          fprintf(fdz,"\n");
          fprintf(fFD,"\n");
       }
 }
 for(i=1;i<T+1;i++)
 {
       fprintf(fZSY,"%6.2f\n",zsy[i]);
       fprintf(fVY,"%6.2f\n",v[i]);
       fprintf(fQS,"%6.2f\n",Qfd[i]);
       fprintf(fNS,"%6.2f\n",N[i]);
       fprintf(fQQ,"%6.2f\n",q_qs[i]);
       fprintf(fdz,"%6.2f\n",dz[i]);
       fprintf(fFD,"%6.2f\n",E[i]);
 }
   return 0;
}