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

关键词

Python | Julia | 一维 | 二维 | 透射率 | 水头 | 流量 | 畜水层 | 河流 | 海域 | 承压层 | 补给 | 分水岭 | 渗透 | 矢量 | 算法 | 模型 | 建模 | 湖泊 | 排水 | 稳定 | 渗漏层 | 瞬态流 | 阶跃响应 | 周期性 | 拉普拉斯变换 | 微分方程 | 饱和 | 线性 | 厚度 | 孔 | 不均匀 | 垂直 | 有限差分

Python和Julia河流湖泊沿海水域特征数值算法模型 | Notion亚图跨际 on Notion

🏈指点迷津 | Brief

🎯要点

1. 一维水流场景计算和绘图：

   1. 🎯恒定透射率水头和流量计算：🖊两条完全穿透畜水层理想河流之间 | 🖊无承压畜水层两侧及两条完全穿透畜水层的补给 | 🖊分水岭或渗透性非常低的岩体的不渗透边界补给 | 🖊两个不同且恒定透射率的畜水层区域补给。

   2. 🎯半承压畜水层水头和流量矢量计算：🖊运河到圩田区域 | 🖊湖泊和排水区域 | 🖊流向一条又长又直且宽的河流 | 🖊一条河流两种不同的畜水层 | 🖊两河流两种不同畜水层补给区​。

   3. 🎯可变饱和厚度无侧限水流：🖊不可渗透的边界和河流之间补给区域 | 🖊畜水层底部水流 | 🖊承压流变为无承压流的畜水层​。

   4. 🎯沿海承压畜水层：界面位置和流量：🖊稳定界面流 | 🖊无侧限界面流 | 🖊渗漏层与畜水层隔开的海洋界面流​。

   5. 🎯瞬态流：🖊地表水位突变响应 | 🖊地下水对地表水位周期性变化响应 | 🖊两条河流之间的瞬时补给 | 🖊拉普拉斯变换解地表水流微分方程 | 🖊恒定透射率和存储期间的线性化饱和厚度变化的响应。

2. 二维水流场景计算和绘图：

   1. 🎯抽水井、注入井以及靠近河流和不渗透边界的井：🖊无约束流向半径的井 | 🖊抽水井和注入井 | 🖊不均匀边界抽水井 | 🖊半承压畜水层中抽水井 | 🖊两畜水层底部抽水的井。

   2. 🎯均匀流情况下的孔：🖊单孔 | 🖊涡流孔 | 🖊沿河流域的单孔 | 🖊具有渗流层的沿河流域单孔 | 🖊沿海流域的单孔。

   3. 🎯解析元建模 | 🎯瞬态流 | 🎯垂直水流​。

3. 🎯Python统计模型处理流体数据 | 🎯Julia模拟逼近非饱和区水域有限差分传输模型。

🍇Python计算畜水层补给成本

畜水层管理补给是有意、有管理性地对畜水层进行补给，以供日后使用或环境效益。

• 可以根据不同的目标进行操作，例如平衡季节性需求和供应或促进水循环利用。

• 还可用于增强长期水安全和抗旱能力。

• 含水层储存实际上是一个不需要建造的水库，主要的基础设施成本仅限于进入含水层进行补给和恢复。

根据基础设施项目经济评估指南，通常需要测试输入参数和假设的敏感性，以确保该过程对潜在变化和情景范围具有鲁棒性。通过概率建模和单变量敏感性分析来量化畜水层不确定性因素。

分析框架和评估模型：

分析框架如下图所示。使用针对一系列计划规模和运行条件的水平衡和财务评估函数，将输入变量传递到迭代模型中。 产出包括平准化成本分布、敏感性分析以及分类资本和运营成本。

用于计算不同规模和运营条件下的畜水层计划成本分布的子模型和可变范围场景：

$$\begin{array}{|l|l|l|l|l|} \hline \text { 子模型 } & \text { 情况 } & \text { 补给 }( MCM / y )^* & \text { 井注入或渗透率 } & \text { 回收效率 }(\%) \\ \hline \text { 井注入 } & 1 & 0.25,0.5,1,2,3,4,5 & 25-32.5 L / s & 95-100 \\ \hline & 2 & 0.25,0.5,1,2,3,4,5 & 20-26 L / s & 90-94.5 \\ \hline & 3 & 0.25,0.5,1,2,3,4,5 & 15-19.5 L / s & 85-89 \\ \hline & 4 & 0.25,0.5,1,2,3,4,5 & 10-13 L / s & 80-84 \\ \hline \text { 渗透盆地 } & 1 & 0.25,0.5,1,2,3,4,5 & 0.8-1 m / d & 95-100 \\ \hline & 2 & 0.25,0.5,1,2,3,4,5 & 0.6-0.78 m / d & 90-94.5 \\ \hline & 3 & 0.25,0.5,1,2,3,4,5 & 0.4-0.52 m / d & 85-89 \\ \hline & 4 & 0.25,0.5,1,2,3,4,5 & 0.2-0.26 m / d & 80-84 \\ \hline \end{array}$$

计算畜水层补给成本

 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")