🥦Python和Julia河流湖泊沿海水域特征数值算法模型
关键词
Python | Julia | 一维 | 二维 | 透射率 | 水头 | 流量 | 畜水层 | 河流 | 海域 | 承压层 | 补给 | 分水岭 | 渗透 | 矢量 | 算法 | 模型 | 建模 | 湖泊 | 排水 | 稳定 | 渗漏层 | 瞬态流 | 阶跃响应 | 周期性 | 拉普拉斯变换 | 微分方程 | 饱和 | 线性 | 厚度 | 孔 | 不均匀 | 垂直 | 有限差分
🎯要点
一维水流场景计算和绘图：
🎯恒定透射率水头和流量计算：🖊两条完全穿透畜水层理想河流之间 | 🖊无承压畜水层两侧及两条完全穿透畜水层的补给 | 🖊分水岭或渗透性非常低的岩体的不渗透边界补给 | 🖊两个不同且恒定透射率的畜水层区域补给。
🎯半承压畜水层水头和流量矢量计算：🖊运河到圩田区域 | 🖊湖泊和排水区域 | 🖊流向一条又长又直且宽的河流 | 🖊一条河流两种不同的畜水层 | 🖊两河流两种不同畜水层补给区。
🎯可变饱和厚度无侧限水流：🖊不可渗透的边界和河流之间补给区域 | 🖊畜水层底部水流 | 🖊承压流变为无承压流的畜水层。
🎯沿海承压畜水层：界面位置和流量：🖊稳定界面流 | 🖊无侧限界面流 | 🖊渗漏层与畜水层隔开的海洋界面流。
🎯瞬态流：🖊地表水位突变响应 | 🖊地下水对地表水位周期性变化响应 | 🖊两条河流之间的瞬时补给 | 🖊拉普拉斯变换解地表水流微分方程 | 🖊恒定透射率和存储期间的线性化饱和厚度变化的响应。
二维水流场景计算和绘图：
🎯抽水井、注入井以及靠近河流和不渗透边界的井：🖊无约束流向半径的井 | 🖊抽水井和注入井 | 🖊不均匀边界抽水井 | 🖊半承压畜水层中抽水井 | 🖊两畜水层底部抽水的井。
🎯均匀流情况下的孔：🖊单孔 | 🖊涡流孔 | 🖊沿河流域的单孔 | 🖊具有渗流层的沿河流域单孔 | 🖊沿海流域的单孔。
🎯解析元建模 | 🎯瞬态流 | 🎯垂直水流。
🎯Python统计模型处理流体数据 | 🎯Julia模拟逼近非饱和区水域有限差分传输模型。
🍇Python计算畜水层补给成本
畜水层管理补给是有意、有管理性地对畜水层进行补给，以供日后使用或环境效益。
可以根据不同的目标进行操作，例如平衡季节性需求和供应或促进水循环利用。
还可用于增强长期水安全和抗旱能力。
含水层储存实际上是一个不需要建造的水库，主要的基础设施成本仅限于进入含水层进行补给和恢复。
根据基础设施项目经济评估指南，通常需要测试输入参数和假设的敏感性，以确保该过程对潜在变化和情景范围具有鲁棒性。通过概率建模和单变量敏感性分析来量化畜水层不确定性因素。
分析框架和评估模型：
分析框架如下图所示。使用针对一系列计划规模和运行条件的水平衡和财务评估函数，将输入变量传递到迭代模型中。 产出包括平准化成本分布、敏感性分析以及分类资本和运营成本。
用于计算不同规模和运营条件下的畜水层计划成本分布的子模型和可变范围场景：
子模型 井注入 渗透盆地 情况 12341234 补给 (MCM/y)∗0.25,0.5,1,2,3,4,50.25,0.5,1,2,3,4,50.25,0.5,1,2,3,4,50.25,0.5,1,2,3,4,50.25,0.5,1,2,3,4,50.25,0.5,1,2,3,4,50.25,0.5,1,2,3,4,50.25,0.5,1,2,3,4,5 井注入或渗透率 25−32.5L/s20−26L/s15−19.5L/s10−13L/s0.8−1m/d0.6−0.78m/d0.4−0.52m/d0.2−0.26m/d 回收效率 (%)95−10090−94.585−8980−8495−10090−94.585−8980−84
计算畜水层补给成本
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")Last updated