Python 逆調整池運用解析修正版（Flood Routine）
逆調整池運用解析プログラムの修正版
- 関数FLOODの中の収束計算で、前回貯水池水位のみを既知としていたが、本来貯水池水位と放流量はペアなので、前回放流量も既知とした。プログラム作成時になぜ収束計算の中で前回流量を既知として確定していなかったのかは思い出せない。
- 関連出力は、データフレームにして一気に書き出すことにした。
以下最新版プログラム（2020.08.08）
# ============================================== # Flood Routine Analysis # ============================================== import sys import numpy as np import pandas as pd from scipy import interpolate def FLOOD(iOFC,fnameR1,fnameR2,fnameR3,ram,qrel): nr,rh,rv=INP_HV(fnameR1) # reservoir H-V nd,ti,q_in,ELini,Qoini=INP_INFLOW(fnameR2) # inflow time history no,elof,qof=INP_OUTFLOW(fnameR3) # outflow capacity pEL=np.zeros(nd,dtype=np.float64) # storage of reervoir water level pQo=np.zeros(nd,dtype=np.float64) # storage of outflow pDE=np.zeros(nd,dtype=np.float64) # storage of error pVO=np.zeros(nd,dtype=np.float64) # storage of reservoir volume pRA=np.zeros(nd,dtype=np.float64) # storage of ramp up rate pUD=np.zeros(nd,dtype=int) # storage of indicator od water level Vin=np.sum(q_in)*(ti[2]-ti[1])/2 # total volume of inflow Qimax=np.max(q_in) # Initial setting elv=ELini # elv=WLINI(q_in[0],elof,qof,no) # Initial reservoir water level EL=elv VOL=RET_V(nr,rh,rv,elv) # Initial reservoir volume q_out=Qoini qout0=Qoini print('q_in0={0:10.3f}'.format(q_in[0])) print('elv_ini={0:10.3f}'.format(EL)) print('q_ini={0:10.3f}'.format(q_out)) i=0 iUD=0 pUD[i]=iUD pEL[i]=EL pQo[i]=q_out pDE[i]=EL-elv pVO[i]=VOL pRA[i]=0 # Iterative calculation iUD=1 # Assume reservoir water level increase dh=0.0005 # Water level increment for searching equilibrium point eps=0.001 # Convergence criteria itmax=int(1.0/dh)*100 Vf=RET_V(nr,rh,rv,60.0) Vz=0 Vin2=0 Vot2=0 q1=q_out for i in range(0,nd-1): dt=ti[i+1]-ti[i] # Time interval Qin=0.5*(q_in[i]+q_in[i+1]) # Average inflow rate during dt hh=0.0 # Zero set of water level increment k=0 j=0 while 1: k=k+1 j=j+1; hh=hh+float(iUD)*dh # Update waler level increment elv1=EL elv2=EL+hh if iOFC==8: q2=OFC8(q1,qout0,elv2,elof,qof,no,ti[i],ram,qrel) Qout=0.5*(q1+q2) # Average outflow rate during dt dS=(Qin-Qout)*dt*3600.0 # Storage accumulated during dt R=VOL+dS # Total reservoir volume elv=RET_H(nr,rh,rv,R) # Reservoir water level at reservoir volume R if iUD==1 and j==10: # 10 times trial assuming reservoir water level increase if EL+hh>elv: # if reservoir water level is less than assumed water level, iUD=-1 # assume water level decrease hh=0.0 j=0 if iUD==-1 and j==10: # 10 times trial assuming reservoir water level decrease if EL+hh<elv: # if reservoir water level is greater than assumed water level, iUD=1 # assume water level decrease hh=0.0 j=0 if np.abs(EL+hh-elv)<eps: break # Judgement of convergence if itmax<k: break VOL=R # Cumulative volume EL=EL+hh # Elevation q_out=q2 # Outflow q1=q2 pUD[i+1]=iUD pEL[i+1]=EL pQo[i+1]=q_out pDE[i+1]=EL-elv pVO[i+1]=VOL pRA[i+1]=(pQo[i+1]-pQo[i])/(ti[i+1]-ti[i]) hmax=max(pEL)-ELini sys.stdout.write('Time: %10.3f\n'% np.max(ti)) sys.stdout.write('h : %10.3f\n'% hmax) sys.stdout.write('Qin : %10.3f %10.3f\n'% (np.min(q_in),np.max(q_in))) sys.stdout.write('Qout: %10.3f %10.3f\n'% (np.min(pQo),np.max(pQo))) sys.stdout.write('EL : %10.3f %10.3f\n'% (np.min(pEL),np.max(pEL))) sys.stdout.write('\n') # TWC of 80m downstream of Re-Regulating Dam fnameR='inp_twc.txt' data=np.loadtxt(fnameR,skiprows=0,usecols=(0,1)) qq=data[:,0] ee=data[:,1] ff=interpolate.interp1d(qq,ee) el_ds=ff(pQo) df = pd.DataFrame({ 'i' : np.arange(0,nd), # numper 'iUD': pUD, # indicator for water level 'ti' : ti, # time 'EL' : pEL, # water level of reservoir 'dEL': pDE, # error 'VOL': pVO, # reservoir volume 'Qi' : q_in, # inflow to reservoir 'Qo' : pQo, # outflow from reservoir 'RA' : pRA, # ramp up rate 'ELD': el_ds # downstream water level }) df = df.set_index('i') return df def WLINI(q,elof,qof,no): # To obtain initial water level at oitflow discharge q for i in range(0,no-1): if qof[i]<=q and q<=qof[i+1]: break if qof[no-1]<q: i=no-2 x1=qof[i] y1=elof[i] x2=qof[i+1] y2=elof[i+1] a=(y2-y1)/(x2-x1) b=(x2*y1-x1*y2)/(x2-x1) elv=a*q+b if q<0.0: elv=elof[0] return elv def OFC8(q1,qout0,elv2,elof,qof,no,tii,ram,qrel): Qlim=300.0 if 5.0<tii and tii<= 29.0 :qout0=qrel if 29.0<tii and tii<= 53.0 :qout0=qrel if 53.0<tii and tii<= 77.0 :qout0=qrel if 77.0<tii and tii<=101.0 :qout0=qrel if 101.0<tii: qout0=39.0 q2=OFC(elv2,elof,qof,no) rq2=q2 if q2<qout0: qout0=q2 q2=qout0 st=0.01 ts1=0.0 ts2=24.0 ts3=48.0 ts4=72.0 ts5=96.0 if ts1+st<=tii and tii<=ts1+5.0+st: if elv2<68.0: q2=q1+ram*st if rq2<q2: q2=rq2 if Qlim<=q2: q2=Qlim else: q2=rq2 if ts2+st<=tii and tii<=ts2+5.0+st: if elv2<68.0: q2=q1+ram*st if rq2<q2: q2=rq2 if Qlim<=q2: q2=Qlim else: q2=rq2 if ts3+st<=tii and tii<=ts3+5.0+st: if elv2<68.0: q2=q1+ram*st if rq2<q2: q2=rq2 if Qlim<=q2: q2=Qlim else: q2=rq2 if ts4+st<=tii and tii<=ts4+5.0+st: if elv2<68.0: q2=q1+ram*st if rq2<q2: q2=rq2 if Qlim<=q2: q2=Qlim else: q2=rq2 if ts5+st<=tii and tii<=ts5+5.0+st: if elv2<68.0: q2=q1+ram*st if rq2<q2: q2=rq2 if Qlim<=q2: q2=Qlim else: q2=rq2 return q2 def OFC(elv,elof,qof,no): # Free overflow discharge for i in range(0,no-1): if elof[i]<=elv and elv<=elof[i+1]: break if elof[no-1]<elv: i=no-2 x1=elof[i] y1=qof[i] x2=elof[i+1] y2=qof[i+1] a=(y2-y1)/(x2-x1) b=(x2*y1-x1*y2)/(x2-x1) q=a*elv+b if elv<=elof[0]: q=0.0 return q def RET_V(nr,rh,rv,elv): # To obtain reservoir cumulative volume from the water level for i in range(0,nr-1): if rh[i]<=elv and elv<=rh[i+1]: break if rh[nr-1]<elv: i=nr-2 x1=rv[i] y1=rh[i] x2=rv[i+1] y2=rh[i+1] a=(y2-y1)/(x2-x1) b=(x2*y1-x1*y2)/(x2-x1) v=(elv-b)/a return v def RET_H(nr,rh,rv,v): # To obtain reservoir water level from cumulative volume for i in range(0,nr-1): if rv[i]<=v and v<=rv[i+1]: break if rv[nr-1]<v: i=nr-2 x1=rv[i] y1=rh[i] x2=rv[i+1] y2=rh[i+1] a=(y2-y1)/(x2-x1) b=(x2*y1-x1*y2)/(x2-x1) elv=a*v+b return elv def INP_HV(fnameR1): # Input reservoir H-V data fin=open(fnameR1,'r') dat=fin.readlines() fin.close() n=len(dat) text=dat[0] text=text.strip() text=text.split() vcoef=float(text[0]) rh=np.zeros(n-1,dtype=np.float64) rv=np.zeros(n-1,dtype=np.float64) for i in range(0,n-1): text=dat[i+1] text=text.strip() text=text.split() rh[i]=float(text[0]) rv[i]=float(text[1])*vcoef nr=len(rh) return nr,rh,rv def INP_INFLOW(fnameR2): # Input time sequence of inflow fin=open(fnameR2,'r') dat=fin.readlines() fin.close() n=len(dat) text=dat[0] text=text.strip() text=text.split() ELini=float(text[0]) Qoini=float(text[1]) ti =np.zeros(n-1,dtype=np.float64) q_in=np.zeros(n-1,dtype=np.float64) for i in range(0,n-1): text=dat[i+1] text=text.strip() text=text.split() ti[i] =float(text[0]) q_in[i]=float(text[1]) nd=len(ti) return nd,ti,q_in,ELini,Qoini def INP_OUTFLOW(fnameR3): # Input outflow capacity (water level-discharge curve) fin=open(fnameR3,'r') dat=fin.readlines() fin.close() n=len(dat) elof=np.zeros(n,dtype=np.float64) qof =np.zeros(n,dtype=np.float64) for i in range(0,n): text=dat[i] text=text.strip() text=text.split() elof[i]=float(text[0]) qof[i] =float(text[1]) no=len(elof) return no,elof,qof def main(): lout=['out_150.csv','out_175.csv','out_200.csv','out_225.csv','out_250.csv'] lram=[150,175,200,225,250] lqrel=[89.0,87.1,82.9,81.5,80.7] for out,ram,qrel in zip(lout,lram,lqrel): iOFC=8 fnameR1='inp_RHV.txt' # H-V data of reservoir fnameR2='inp_c00.txt' # Inflow time history (Hydrograph) fnameR3='inp_OF8.txt' # Outflow capacity fnameW = out # Output file name df=FLOOD(iOFC,fnameR1,fnameR2,fnameR3,ram,qrel) df.to_csv(fnameW, sep=",") #============================ if __name__ == '__main__': main()
以 上
