Python和Julia河流湖泊沿海水域特征数值算法模型

🎯要点

  1. 一维水流场景计算和绘图:
    1. 🎯恒定透射率水头和流量计算:🖊两条完全穿透畜水层理想河流之间 | 🖊无承压畜水层两侧及两条完全穿透畜水层的补给 | 🖊分水岭或渗透性非常低的岩体的不渗透边界补给 | 🖊两个不同且恒定透射率的畜水层区域补给。
    2. 🎯半承压畜水层水头和流量矢量计算:🖊运河到圩田区域 | 🖊湖泊和排水区域 | 🖊流向一条又长又直且宽的河流 | 🖊一条河流两种不同的畜水层 | 🖊两河流两种不同畜水层补给区​。
    3. 🎯可变饱和厚度无侧限水流:🖊不可渗透的边界和河流之间补给区域 | 🖊畜水层底部水流 | 🖊承压流变为无承压流的畜水层​。
    4. 🎯沿海承压畜水层:界面位置和流量:🖊稳定界面流 | 🖊无侧限界面流 | 🖊渗漏层与畜水层隔开的海洋界面流​。
    5. 🎯瞬态流:🖊地表水位突变响应 | 🖊地下水对地表水位周期性变化响应 | 🖊两条河流之间的瞬时补给 | 🖊拉普拉斯变换解地表水流微分方程 | 🖊恒定透射率和存储期间的线性化饱和厚度变化的响应。
  2. 二维水流场景计算和绘图:
    1. 🎯抽水孔、注入孔以及靠近河流和不渗透边界的孔:🖊无约束流向半径的孔 | 🖊抽水孔和注入孔 | 🖊不均匀边界抽水孔 | 🖊半承压畜水层中抽水孔 | 🖊两畜水层底部抽水的孔。
    2. 🎯均匀流情况下的孔:🖊单孔 | 🖊涡流孔 | 🖊沿河流域的单孔 | 🖊具有渗流层的沿河流域单孔 | 🖊沿海流域的单孔。
    3. 🎯解析元建模 | 🎯瞬态流 | 🎯垂直水流​。
  3. 🎯Python统计模型处理流体数据 | 🎯Julia模拟逼近非饱和区水域有限差分传输模型。

🍇Python计算畜水层补给成本

import time
start = time.time()
from SAb.sample import latin
from SAb.util import read_param_file
from SAb.analyze import delta
import numpy as np
import numpy_financial as npf
import pandas as pd
import math
from matplotlib import pyplot as plt

p = read_param_file('inputs_asr.txt')
p1 = dict(p)
p2 = dict(p)
p3 = dict(p)
p4 = dict(p)

bnds = p.get("bounds")
bnds1 = list(bnds)
bnds2 = list(bnds)
bnds3 = list(bnds)
bnds4 = list(bnds)

del bnds, p

f = 1.3 
bnds1[1:3] = [25,25*f],[.95,.95*1.05] 
bnds2[1:3] = [20,20*f],[.90,.90*1.05]
bnds3[1:3] = [15,15*f],[.85,.85*1.05]
bnds4[1:3] = [10,10*f],[.80,.80*1.05] 

n = 10000

rechMLy = [250,500,1000,2000,3000,4000,5000]

p1.update({"bounds":bnds1})
x1 = latin.sample(p1, n)
p2.update({"bounds":bnds2})
x2 = latin.sample(p2, n)
p3.update({"bounds":bnds3})
x3 = latin.sample(p3, n)
p4.update({"bounds":bnds4})
x4 = latin.sample(p4, n)

x = np.stack((x1,x2,x3,x4))
xn = x.shape[0] 


th = x[:,:,0].astype(int) 
boreyieldLs = x[:,:,1].round() 
re = x[:,:,2] 
treatcap = x[:,:,3] 
treatop = x[:,:,4] 
rate = x[:,:,5]
pumpinstrcap = x[:,:,6].astype(int) 
alloc_cost = x[:,:,7].astype(int) 
pipecostkm = x[:,:,8] 
wellredevpct = x[:,:,9] t 
rivpumplift = x[:,:,10].astype(int) 
transpump = x[:,:,11] 
pumpop = x[:,:,12] 
monitoring = x[:,:,13].astype(int) 
obsbore = x[:,:,14].astype(int)  
injwell = x[:,:,15].astype(int) 

Feas = x[:,:,16].astype(int) 
dy = x[:,:,17].astype(int) 
Dummy = x[:,:,18] 

bmin = np.amin(boreyieldLs,axis=1).astype(int) 
bmax = np.amax(boreyieldLs,axis=1).astype(int) 
remax = np.amax(re,axis=1).round(2)
remin = np.amin(re,axis=1).round(2)

trig = pd.read_csv('lock_1_avg_ann_dish_72_21.txt', header=None)
trig = np.asarray(trig).reshape(50,)
t25 = np.asarray([1] + [0] * 23 + [1]) 

def timef (th, rechMLy, re):
    rech = np.zeros(50,)
    recov = np.zeros(50,)
    thres = np.array([th]*50)
    rech[trig > thres] = 1
    rech[0] = 0 
    recov[trig < thres] = 1
    aquifML = rechMLy * 1e6 
    
    a = np.array([rechMLy * rech]).reshape(50,)
    b = np.array([rechMLy * recov * 1. * -1]).reshape(50,)
    d = np.zeros_like(a)

    if a[0]+b[0]>aquifML:
        d[0] = aquifML
    else:
        d[0] = a[0]+b[0]
    for j in range(1,len(a)):
        c = (d[j-1] * re) + a[j] + b[j]
        if c>aquifML:
            d[j] = aquifML
        elif c<0:
            d[j] = 0
        elif c<aquifML:
            d[j] = (d[j-1] * re) + a[j] + b[j]
    # calc rech recov vols from d
    e = np.zeros_like(d)
    e[0] = d[0]
    for k in range(1,len(d)):
            e[k] = d[k] - (d[k-1] * re)
        
    vi = e * rech # vols in
    vo = e * recov # vols out
    vi_tot = np.sum(vi.round(1))
    vo_tot = np.sum(vo.round(1)*-1)
    vi_vo = vi + vo
    vi_vo_0 = 49 - np.count_nonzero(vi + vo) 
    return vi, vo, vi_tot, vo_tot, rech, vi_vo, vi_vo_0

#%% run vols times series (vol, scen, yrs, iter)
vol_in = np.zeros([len(rechMLy),xn,50,n])
vol_out = np.zeros([len(rechMLy),xn,50,n])
vi_vo = np.zeros([len(rechMLy),xn,50,n])
vi_vo_0 = np.zeros([len(rechMLy),xn,50,n])
Vol_in_tot_ML = np.zeros([len(rechMLy),xn,n])
Vol_out_tot_ML = np.zeros([len(rechMLy),xn,n])
rech_yrs = np.zeros([len(rechMLy),xn,50,n])

for r in range(len(rechMLy)):
    for l in range(xn):
        for i in range(n):
            vol_in[r,l,:,i], vol_out[r,l,:,i], Vol_in_tot_ML[r,l,i], Vol_out_tot_ML[r,l,i], rech_yrs[r,l,:,i], vi_vo[r,l,:,i], vi_vo_0[r,l,:,i] = timef(th[l,i], rechMLy[r], re[l,i])     

Vol_in_tot_m3 = Vol_in_tot_ML * 1000
Vol_out_tot_m3 = Vol_out_tot_ML * 1000
Pct_vol_recovered = Vol_out_tot_ML / Vol_in_tot_ML *100

def rech_capex (rechMLy, injwell, pumpinstrcap, boreyieldLs, dy, rate):
        Nbores = math.ceil(rechMLy / (boreyieldLs * 86400 * dy * 0.000001))
        Wellconstr = Nbores * injwell
        Pumpinstrcap = npf.npv(rate, (pumpinstrcap * Nbores) * t25) 
        Rech_cap = Pumpinstrcap + Wellconstr
        Rech_capML = Rech_cap / rechMLy
        return Rech_cap, Rech_capML, Nbores, Wellconstr, Pumpinstrcap

def rech_opex (vol_in, pumpop, Wellconstr, wellredevpct, redev5y, rate):
        rech_pump = vol_in * pumpop * 50 
        Rech_pump = npf.npv(rate, rech_pump)
        wellredev = (Wellconstr * wellredevpct) * redev5y 
        Rech_maint = npf.npv(rate, wellredev) 
        return Rech_pump, Rech_maint

#%% recharge capex
Rech_cap = np.zeros([len(rechMLy),x.shape[0],x.shape[1]])
Rech_capML = np.zeros([len(rechMLy),x.shape[0],x.shape[1]])
Nbores = np.zeros([len(rechMLy),x.shape[0],x.shape[1]])
Wellconstr = np.zeros([len(rechMLy),x.shape[0],x.shape[1]])
Pumpinstrcap = np.zeros([len(rechMLy),x.shape[0],x.shape[1]])
redev5y = np.zeros_like(trig)
redev5y = redev5y + ([0,0,0,0,1] * 10) 

for r in range(len(rechMLy)):
    for l in range(xn):
        for i in range (n):
            Rech_cap[r,l,i], Rech_capML[r,l,i], Nbores[r,l,i], Wellconstr[r,l,i], Pumpinstrcap[r,l,i] = rech_capex(rechMLy[r], injwell[l,i], pumpinstrcap[l,i], boreyieldLs[l,i], dy[l,i], rate[l,i])

#%% recharge opex
Rech_pump = np.zeros([len(rechMLy),x.shape[0],x.shape[1]])
Rech_maint = np.zeros([len(rechMLy),x.shape[0],x.shape[1]])
for r in range(len(rechMLy)):
    for l in range (xn):
        for i in range (n):
            Rech_pump[r,l,i], Rech_maint[r,l,i] = rech_opex(vol_in[r,l,:,i], pumpop[l,i], Wellconstr[r,l,i], redev5y, rate[l,i], wellredevpct[l,i])

def monitor(monitoring, rate, obsbore, vol_in, rechMLy):
    Monitor_op = npf.npv(rate, ([monitoring] * (vol_in / 50)))
    Obsbore_cap = math.ceil(rechMLy / 1000) * (obsbore * 50)
    return Monitor_op, Obsbore_cap

Monitor_op = np.zeros([len(rechMLy),x.shape[0],x.shape[1]])
Obsbore_cap = np.zeros([len(rechMLy),x.shape[0],x.shape[1]])
for r in range(len(rechMLy)):
    for l in range(xn):
        for i in range (n):
            Monitor_op[r,l,i], Obsbore_cap[r,l,i] = monitor(monitoring[l,i], rate[l,i], obsbore[l,i], vol_in[r,l,:,i], rechMLy[r])

def pipe_cost(pipecostkm, rechMLy, rate):
    Pipe_cap = pipecostkm * rechMLy * 1e-03
    pipeop = Pipe_cap * 0.02 
    Pipe_op = npf.npv(rate, ([pipeop] * 50)) - pipeop # no pipe maint in 1st yr
    Pipe_tot = Pipe_cap + Pipe_op
    return Pipe_cap, Pipe_op, Pipe_tot

#%% run pipe costs
Pipe_cap = np.zeros([len(rechMLy),x.shape[0],x.shape[1]]) 
Pipe_op = np.zeros([len(rechMLy),x.shape[0],x.shape[1]])  
Pipe_tot = np.zeros([len(rechMLy),x.shape[0],x.shape[1]]) 
for r in range(len(rechMLy)):
    for l in range(xn):
        for i in range (n):
            Pipe_cap[r,l,i], Pipe_op[r,l,i], Pipe_tot[r,l,i] = pipe_cost(pipecostkm[l,i], rechMLy[r], rate[l,i])

def water_alloc_cost (rate, alloc_cost, vol_in):
    water_alloc = vol_in * alloc_cost
    Water_alloc_op = npf.npv(rate, water_alloc) - alloc_cost 
    return Water_alloc_op

Water_alloc_op = np.zeros([len(rechMLy),x.shape[0],x.shape[1]]) 
for r in range(len(rechMLy)):
    for l in range (xn):
        for i in range (n):
            Water_alloc_op[r,l,i] = water_alloc_cost(rate[l,i], alloc_cost[l,i], vol_in[r,l,:,i])
    
def recov_cost (vol_out, pumpop, rate):
    recov = vol_out * pumpop * 50 * -1
    Recov_op = npf.npv(rate, recov)
    return Recov_op

Recov_op = np.zeros([len(rechMLy),x.shape[0],x.shape[1]]) 
for r in range(len(rechMLy)):
    for l in range (xn):
        for i in range (n):
            Recov_op[r,l,i] = recov_cost(vol_out[r,l,:,i], pumpop[l,i], rate[l,i])

def river_pump (pumpop, rivpumplift, vol_in, rate, rechMLy, dy, transpump):
    riverpump = pumpop * rivpumplift * vol_in
    River_pump_op = npf.npv(rate, riverpump)
    river_pump_cap = math.ceil(rechMLy / (dy * 1.296)) * transpump
    River_pump_cap = npf.npv(rate, river_pump_cap * t25)
    return River_pump_op, River_pump_cap

River_pump_op = np.zeros([len(rechMLy),x.shape[0],x.shape[1]])
River_pump_cap =  np.zeros([len(rechMLy),x.shape[0],x.shape[1]])
for r in range(len(rechMLy)):
    for l in range (xn):
        for i in range (n):
            River_pump_op[r,l,i], River_pump_cap[r,l,i] = river_pump(pumpop[l,i], rivpumplift[l,i], vol_in[r,l,:,i], rate[l,i], rechMLy[r], dy[l,i], transpump[l,i]) 
       
def pretreatment (rechMLy, rate, rech_yrs):
    rf = np.random.uniform(low=0.8, high=1.0) 
    Pretreat_cap = (-588.8 * np.log(rechMLy) + 6528.6) * rf * rechMLy
    treat =  (-465.3 * np.log(rechMLy) + 4302.2) * rechMLy * rf * rech_yrs
    Pretreat_op = npf.npv(rate, treat)
    return Pretreat_cap, Pretreat_op

Pretreat_cap = np.zeros([len(rechMLy),x.shape[0],x.shape[1]])
Pretreat_op = np.zeros([len(rechMLy),x.shape[0],x.shape[1]])
for r in range(len(rechMLy)):
    for l in range (xn):
        for i in range (n):
            Pretreat_cap[r,l,i], Pretreat_op[r,l,i] = pretreatment(rechMLy[r], rate[l,i], rech_yrs[r,l,:,i])
        

dir()
Tot_cost = Feas+Monitor_op+Obsbore_cap+Pipe_cap+Pipe_op+Pretreat_cap+Pretreat_op+Rech_cap+Rech_maint+Rech_pump+Recov_op+River_pump_op+River_pump_cap+Water_alloc_op
Rech_cost = Feas+Monitor_op+Obsbore_cap+Pipe_cap+Pipe_op+Pretreat_cap+Pretreat_op+Rech_cap+Rech_maint+Rech_pump+River_pump_op+River_pump_cap+Water_alloc_op
LC_rech_m3 = (Rech_cost / Vol_in_tot_m3).round(2)
LC_recov_m3 = (Tot_cost / Vol_out_tot_m3).round(2)
Tot_capex = Feas+Obsbore_cap+Pipe_cap+Pretreat_cap+Rech_cap+River_pump_cap
Tot_opex = Monitor_op+Pipe_op+Pretreat_op+Rech_maint+Rech_pump+Recov_op+River_pump_op+Water_alloc_op
Ann_opex = (Monitor_op+Pipe_op+Pretreat_op+Rech_maint+Rech_pump+Recov_op+River_pump_op+Water_alloc_op) / 49
Tot_opex_capex = Tot_opex / Tot_capex
Ann_opex_capex = Ann_opex / Tot_capex
Pretreat_op_cap = Pretreat_op / Pretreat_cap
Prop_rech_lost = (Vol_in_tot_ML - Vol_out_tot_ML) / Vol_in_tot_ML

Pretreat_prop_cap = (Pretreat_cap / Vol_out_tot_m3) / LC_recov_m3
Injwell_prop_cap = (Rech_cap / Vol_out_tot_m3) / LC_recov_m3
Pump_pipe_prop_cap = ((River_pump_cap+Pipe_cap) / Vol_out_tot_m3) / LC_recov_m3
Obsbore_prop_cap = (Obsbore_cap / Vol_out_tot_m3) / LC_recov_m3
Feas_prop_cap = (Feas / Vol_out_tot_m3) / LC_recov_m3

Pretreat_prop_op = (Pretreat_op / Vol_out_tot_m3) / LC_recov_m3
Pump_prop_op = ((River_pump_op + Recov_op + Rech_pump) / Vol_out_tot_m3) / LC_recov_m3
Wateralloc_prop_op = (Water_alloc_op / Vol_out_tot_m3) / LC_recov_m3
Monitor_prop_op = (Monitor_op / Vol_out_tot_m3) / LC_recov_m3
Maint_prop_op = ((Pipe_op+Rech_maint) / Vol_out_tot_m3) / LC_recov_m3

dp = dict([(key, []) for key in ['ind','MLy','RE_Yield',
           'Pre-treatment plant','Injection well construction','Pumps & pipes','Observation bores','Feasibility studies',
           'Pre-treatment','Pumping cost','Opportunity cost water','Monitoring','Maintenance']])
for r in range(len(rechMLy)): # vol
    for m in range (int(xn)): # scen
        dp['ind'].append(str(rechMLy[r])+", "+str(remin[m])+"-"+str(remax[m])+", "+str(bmin[m])+"-"+str(bmax[m]))
        dp['MLy'].append(rechMLy[r])
        dp['RE_Yield'].append(str(remin[m])+"-"+str(remax[m])+", "+str(bmin[m])+"-"+str(bmax[m])+"L/s")
        dp['Pre-treatment plant'].append(np.mean(Pretreat_prop_cap[r,m]).round(2))
        dp['Injection well construction'].append(np.mean(Injwell_prop_cap[r,m]).round(2))
        dp['Pumps & pipes'].append(np.mean(Pump_pipe_prop_cap[r,m]).round(2))
        dp['Observation bores'].append(np.mean(Obsbore_prop_cap[r,m]).round(2))
        dp['Feasibility studies'].append(np.mean(Feas_prop_cap[r,m]).round(2))
        dp['Pre-treatment'].append(np.mean(Pretreat_prop_op[r,m]).round(2))
        dp['Pumping cost'].append(np.mean(Pump_prop_op[r,m]).round(2))
        dp['Opportunity cost water'].append(np.mean(Wateralloc_prop_op[r,m]).round(2))
        dp['Monitoring'].append(np.mean(Monitor_prop_op[r,m]).round(2))
        dp['Maintenance'].append(np.mean(Maint_prop_op[r,m]).round(2))
dfp = pd.DataFrame.from_dict(dp)        
dfp.to_csv("asr_prop_costs.csv")

dfp1 = dfp.groupby('MLy').mean()
labels = rechMLy * np.repeat([0.001], len(rechMLy))
plt.rcParams['font.size'] = '7'
plt.pcolormesh(dfp1,cmap='Reds')
plt.yticks(np.arange(0.5, len(dfp1), 1),labels=labels)
plt.ylabel('Scheme capacity (MCM/y)')
plt.xticks(np.arange(0.5, len(dfp1.columns), 1), dfp1.columns, rotation=45)
plt.colorbar()
plt.text(0,1,'  <-------------------- Capital costs -------------------->')
plt.text(5,1,'  <---------------- Operational costs ------------------>')
plt.savefig('asr_heatmap.png',dpi=300)
plt.show()
plt.close()

d = dict([(key, []) for key in ['mar_type','MLy','re','yield',
          'RE_Yield','recharged_vol','recovered_vol','prop_lost', 'asr_well_capML', 'capex$M_min','capex$M_max','capex$M_mean','capex$M_std',
          'ann_opex$K','ann_opex$K_std','tot_opex$M','tot_opex:capex','ann_opex:capex','treat_cap$M','ann_treat_op$K','treat_op:cap','lc_rech',
          'lc_rech_std','lc_recov','lc_recov_std','lc_recov_cov']])
for r in range(len(rechMLy)): # vol
    for m in range (int(xn)): # scen
        d['mar_type'].append(str("ASR"))
        d['MLy'].append(rechMLy[r])
        d['re'].append(remin[m])
        d['yield'].append(bmin[m])
        d['RE_Yield'].append(str(remin[m])+"-"+str(remax[m])+", "+str(bmin[m])+"-"+str(bmax[m])+" L/s")
        d['recharged_vol'].append(np.mean(Vol_in_tot_ML[r,m]))
        d['recovered_vol'].append(np.mean(Vol_out_tot_ML[r,m]))
        d['prop_lost'].append(np.mean(Prop_rech_lost[r,m]))
        d['asr_well_capML'].append(np.mean(Rech_capML[r,m]))
        d['capex$M_min'].append((np.min(Tot_capex[r,m])*0.000001))
        d['capex$M_max'].append((np.max(Tot_capex[r,m])*0.000001))
        d['capex$M_mean'].append((np.mean(Tot_capex[r,m])*0.000001))
        d['capex$M_std'].append((np.std(Tot_capex[r,m])*0.000001))
        d['ann_opex$K'].append((np.mean(Ann_opex[r,m])*0.001))
        d['ann_opex$K_std'].append((np.std(Ann_opex[r,m])*.001))
        d['ann_opex:capex'].append((np.mean(Ann_opex_capex[r,m])*0.001))
        d['tot_opex:capex'].append(np.mean(Tot_opex_capex[r,m]))
        d['tot_opex$M'].append((np.mean(Tot_opex[r,m])*0.000001))
        d['treat_cap$M'].append((np.mean(Pretreat_cap[r,m])*0.000001))
        d['ann_treat_op$K'].append((np.mean(Pretreat_op[r,m])*0.001/49))
        d['treat_op:cap'].append(np.mean(Pretreat_op_cap[r,m]))
        d['lc_rech'].append(np.mean(LC_rech_m3[r,m]))
        d['lc_rech_std'].append(np.std(LC_rech_m3[r,m]))
        d['lc_recov'].append(np.mean(LC_recov_m3[r,m]))
        d['lc_recov_std'].append(np.std(LC_recov_m3[r,m]))
        d['lc_recov_cov'].append(np.std(LC_recov_m3[r,m]) / np.mean(LC_recov_m3[r,m]))
        
df = pd.DataFrame.from_dict(d)

#%%pivots for plots
df2 = df.pivot(index='MLy',columns='RE_Yield', values='lc_rech')
df3 = df.pivot(index='MLy',columns='RE_Yield', values='lc_rech_std') #yerr std
df4 = df.pivot(index='MLy',columns='RE_Yield', values='capex$M_mean')
df5 = df.pivot(index='MLy',columns='RE_Yield', values='capex$M_std') #yerr std
df6 = df.pivot(index='MLy',columns='RE_Yield', values='ann_opex$K')
df7 = df.pivot(index='MLy',columns='RE_Yield', values='ann_opex$K_std') #yerr std
df8 = df.pivot(index='MLy',columns='RE_Yield', values='tot_opex:capex')
df9 = df.pivot(index='MLy',columns='RE_Yield', values='ann_opex:capex')
df10 = df.pivot(index='MLy',columns='RE_Yield', values='lc_recov')
df11 = df.pivot(index='MLy',columns='RE_Yield', values='lc_recov_std') #yerr std
df.to_csv("asr_results.csv")

townsx = [250,750]
townsy = [1.50,1.10]
ticks = rechMLy
labels = rechMLy * np.repeat([0.001], len(rechMLy))
plt.rcParams['font.size'] = '7'
fig, axes = plt.subplots(nrows=1, ncols=2, figsize=(17/2.54,9/2.54))

df2.plot(ax=axes[0], fontsize=8,yerr=df3,linewidth = 1, 
              alpha = 0.6, capsize =2, xticks=ticks, xlim=(0,np.max(rechMLy)+100),
              legend=False, rot=90)
axes[0].set_xlabel('Scheme capacity (MCM/y)')
axes[0].set_ylabel('Levelised cost recharge $/m3')
axes[0].set_xticklabels(labels)
axes[0].text(-500,0,'a')
df10.plot(ax=axes[1], fontsize=8,yerr=df11,linewidth = 1, 
              alpha = 0.9, capsize =2, xticks=ticks, xlim=(0,np.max(rechMLy)+100),
              rot=90)
axes[1].set_xlabel('Scheme capacity (MCM/y)')
axes[1].set_xticklabels(labels)
axes[1].set_ylabel('Levelised cost recovery $/m3')
axes[1].legend(title='Recovery efficiency, bore yield')
axes[1].text(-500,-0.2,'b')
axes[1].scatter(townsx, townsy, alpha=0.3, c='k')
axes[1].text(350,1.55,'Tamworth')
axes[1].text(850,1.15,'Bathurst')
plt.tight_layout()
plt.savefig('asr_lev_cost_rech_recov.png',dpi=300)
plt.show()
plt.close()

rot=90
df4.plot.line(figsize=(9/2.54,9/2.54),fontsize=8,yerr=df5,linewidth = 1., 
              alpha = 0.9, capsize =2)
plt.xlabel('Scheme capacity (MCM/y)')
plt.ylabel('Total capital cost AUD$M')
plt.xlim(0,np.max(rechMLy)+100)
plt.xticks(ticks, labels, rotation=rot)
plt.legend(title="Recovery efficiency, bore yield", title_fontsize=7, loc=0, fontsize=7)
plt.text(-500,-0.05,'a')
plt.tight_layout()
plt.savefig('asr_capex_cost_v07.png',dpi=300)
plt.show()
plt.close()

x1[:,3] = Pretreat_cap[0,0,:] # 1st vol, 1st scen
x1[:,4] = Pretreat_op[0,0,:]
bnds1[3:5] = [np.min(x1[:,3]),np.max(x1[:,3])],[np.min(x1[:,4]),np.max(x1[:,4])]
p1l = p1
p1l.update({"bounds":bnds1})

x4[:,3] = Pretreat_cap[0,3,:] 
x4[:,4] = Pretreat_op[0,3,:]
bnds4[3:5] = [np.min(x4[:,3]),np.max(x4[:,3])],[np.min(x4[:,4]),np.max(x4[:,4])]
p4l = p4
p4l.update({"bounds":bnds4})

x1h = x1
x1h[:,3] = Pretreat_cap[-1,0,:] 
x1h[:,4] = Pretreat_op[-1,0,:]
bnds1[3:5] = [np.min(x1h[:,3]),np.max(x1h[:,3])],[np.min(x1h[:,4]),np.max(x1h[:,4])]
p1h = p1
p1h.update({"bounds":bnds1})

x4h = x4
x4h[:,3] = Pretreat_cap[-1,3,:] 
x4h[:,4] = Pretreat_op[-1,3,:]
bnds4[3:5] = [np.min(x4h[:,3]),np.max(x4h[:,3])],[np.min(x4h[:,4]),np.max(x4h[:,4])]
p4h = p4
p4h.update({"bounds":bnds4})

Si1l = delta.analyze(p1l, x1, LC_recov_m3[0,0,:], print_to_console=False)
Si1l_df = pd.DataFrame.from_dict(Si1l)
Var1l_df = pd.DataFrame.from_dict(p1l)
Si1l_df['variable'] = Var1l_df['names']

Si1h = delta.analyze(p1h, x1h, LC_recov_m3[-1,0,:], print_to_console=False)
Si1h_df = pd.DataFrame.from_dict(Si1h)
Var1h_df = pd.DataFrame.from_dict(p1l)
Si1h_df['variable'] = Var1h_df['names']

Si4l = delta.analyze(p4l, x4, LC_recov_m3[0,3,:], print_to_console=False)
Si4l_df = pd.DataFrame.from_dict(Si4l)
Var4l_df = pd.DataFrame.from_dict(p4l)
Si4l_df['variable'] = Var4l_df['names']

Si4h = delta.analyze(p4h, x4h, LC_recov_m3[-1,3,:], print_to_console=False)
Si4h_df = pd.DataFrame.from_dict(Si4h)
Var4h_df = pd.DataFrame.from_dict(p4l)
Si4h_df['variable'] = Var4h_df['names']

s = dict([(key, []) for key in []])
s["var"] = list(Si1l_df['variable'])
s["s1l"] = list(Si1l_df['S1'])
s["s1lc"] = list(Si1l_df['S1_conf'])

s["s1h"] = list(Si1h_df['S1'])
s["s1hc"] = list(Si1h_df['S1_conf'])

s["s4l"] = list(Si4l_df['S1'])
s["s4lc"] = list(Si4l_df['S1_conf'])

s["s4h"] = list(Si4h_df['S1'])
s["s4hc"] = list(Si4h_df['S1_conf'])

sdf = pd.DataFrame.from_dict(s)
sdf.to_excel('asr_sens.xlsx')

#%% plt sens
plt.rcParams['font.size'] = '7'
plt.rcParams['figure.constrained_layout.use'] = True
fig, ax = plt.subplots(2,2,figsize=(19/2.54 , 12/2.54))
xlim = [0,0.9]

#s1 best case
ax[0,0].barh(sdf['var'], sdf['s1l'], align='center', color='grey', xerr=sdf['s1lc'], label='s1l')
ax[0,0].set_title('Best case low vol', fontsize=7)
ax[0,0].set_xlim(xlim)

ax[0,1].barh(sdf['var'], sdf['s1h'], align='center', color='grey', xerr=sdf['s1hc'], label='s1h')
ax[0,1].get_yaxis().set_visible(False)
ax[0,1].set_title('Best case high vol', fontsize=7)
ax[0,1].set_xlim(xlim)

#s4 worst case
ax[1,0].barh(sdf['var'], sdf['s4l'], align='center', color='grey', xerr=sdf['s4lc'], label='s4l')
ax[1,0].set_title('Worst case low vol', fontsize=7)
ax[1,0].set_xlim(xlim)
ax[1,0].set_xlabel('S1')

ax[1,1].barh(sdf['var'], sdf['s4h'], align='center', color='grey', xerr=sdf['s4hc'], label='s4h')
ax[1,1].set_title('Worst case high vol', fontsize=7)
ax[1,1].get_yaxis().set_visible(False)
ax[1,1].set_xlim(xlim)
ax[1,1].set_xlabel('S1')

fig.savefig("asr_sens.png", dpi=300)
fig.show()

#%%
end = time.time()
print("run time", end - start, "s")

参阅:亚图跨际

最近更新

  1. TCP协议是安全的吗?

    2024-04-24 18:50:03       16 阅读
  2. 阿里云服务器执行yum,一直下载docker-ce-stable失败

    2024-04-24 18:50:03       16 阅读
  3. 【Python教程】压缩PDF文件大小

    2024-04-24 18:50:03       15 阅读
  4. 通过文章id递归查询所有评论(xml)

    2024-04-24 18:50:03       18 阅读

热门阅读

  1. 笨蛋学C++【C++基础第二弹】

    2024-04-24 18:50:03       12 阅读
  2. OneFlow 概念清单

    2024-04-24 18:50:03       11 阅读
  3. MySQL 自建数据库慢日志分析

    2024-04-24 18:50:03       12 阅读
  4. 栈的简单应用:括号匹配

    2024-04-24 18:50:03       11 阅读
  5. Vue.js(过滤器(Filter))

    2024-04-24 18:50:03       12 阅读
  6. class094 贪心经典题目专题6【左程云算法】

    2024-04-24 18:50:03       10 阅读
  7. c# 连接数据库、excel数据批量导入到数据库

    2024-04-24 18:50:03       11 阅读
  8. Semaphore

    Semaphore

    2024-04-24 18:50:03      9 阅读
  9. Dubbo

    Dubbo

    2024-04-24 18:50:03      11 阅读
  10. jvm学习笔记

    2024-04-24 18:50:03       9 阅读