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

🎯要点

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

🍇Python计算畜水层补给成本

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

  • 可以根据不同的目标进行操作，例如平衡季节性需求和供应或促进水循环利用。
  • 还可用于增强长期水安全和抗旱能力。
  • 含水层储存实际上是一个不需要建造的水库，主要的基础设施成本仅限于进入含水层进行补给和恢复。

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

分析框架和评估模型：

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

用于计算不同规模和运营条件下的畜水层计划成本分布的子模型和可变范围场景：
子模型  情况  补给  ( M C M / y ) ∗ 井注入或渗透率  回收效率  ( % ) 井注入  1 0.25 , 0.5 , 1 , 2 , 3 , 4 , 5 25 − 32.5 L / s 95 − 100 2 0.25 , 0.5 , 1 , 2 , 3 , 4 , 5 20 − 26 L / s 90 − 94.5 3 0.25 , 0.5 , 1 , 2 , 3 , 4 , 5 15 − 19.5 L / s 85 − 89 4 0.25 , 0.5 , 1 , 2 , 3 , 4 , 5 10 − 13 L / s 80 − 84 渗透盆地  1 0.25 , 0.5 , 1 , 2 , 3 , 4 , 5 0.8 − 1 m / d 95 − 100 2 0.25 , 0.5 , 1 , 2 , 3 , 4 , 5 0.6 − 0.78 m / d 90 − 94.5 3 0.25 , 0.5 , 1 , 2 , 3 , 4 , 5 0.4 − 0.52 m / d 85 − 89 4 0.25 , 0.5 , 1 , 2 , 3 , 4 , 5 0.2 − 0.26 m / d 80 − 84 \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}  子模型  井注入  渗透盆地 ​ 情况 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 pltp = 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, pf = 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 = 10000rechMLy = [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] = 1rech[0] = 0 recov[trig < thres] = 1aquifML = 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] = aquifMLelse: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] = aquifMLelif c<0:d[j] = 0elif c<aquifML:d[j] = (d[j-1] * re) + a[j] + b[j]# calc rech recov vols from de = 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 invo = e * recov # vols outvi_tot = np.sum(vi.round(1))vo_tot = np.sum(vo.round(1)*-1)vi_vo = vi + vovi_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 *100def rech_capex (rechMLy, injwell, pumpinstrcap, boreyieldLs, dy, rate):Nbores = math.ceil(rechMLy / (boreyieldLs * 86400 * dy * 0.000001))Wellconstr = Nbores * injwellPumpinstrcap = npf.npv(rate, (pumpinstrcap * Nbores) * t25) Rech_cap = Pumpinstrcap + WellconstrRech_capML = Rech_cap / rechMLyreturn Rech_cap, Rech_capML, Nbores, Wellconstr, Pumpinstrcapdef 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_capMonitor_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-03pipeop = Pipe_cap * 0.02 Pipe_op = npf.npv(rate, ([pipeop] * 50)) - pipeop # no pipe maint in 1st yrPipe_tot = Pipe_cap + Pipe_opreturn 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_costWater_alloc_op = npf.npv(rate, water_alloc) - alloc_cost return Water_alloc_opWater_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 * -1Recov_op = npf.npv(rate, recov)return Recov_opRecov_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_inRiver_pump_op = npf.npv(rate, riverpump)river_pump_cap = math.ceil(rechMLy / (dy * 1.296)) * transpumpRiver_pump_cap = npf.npv(rate, river_pump_cap * t25)return River_pump_op, River_pump_capRiver_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 * rechMLytreat =  (-465.3 * np.log(rechMLy) + 4302.2) * rechMLy * rf * rech_yrsPretreat_op = npf.npv(rate, treat)return Pretreat_cap, Pretreat_opPretreat_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_MLPretreat_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_m3Pretreat_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_m3dp = 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)): # volfor m in range (int(xn)): # scendp['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)): # volfor m in range (int(xn)): # scend['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")

参阅一：计算思维

参阅二：亚图跨际

本文来自互联网用户投稿，该文观点仅代表作者本人，不代表本站立场。本站仅提供信息存储空间服务，不拥有所有权，不承担相关法律责任。如若转载，请注明出处：http://www.mzph.cn/web/5640.shtml

如若内容造成侵权/违法违规/事实不符，请联系多彩编程网进行投诉反馈email:809451989@qq.com，一经查实，立即删除！

相关文章

Flask简介

Flask简介

Flask简介 安装概述使用PyCharm创建一个Flask程序 Flask程序的基本结构初始化路由和视图函数启动服务器请求-响应循环 安装 概述 Flask算是小型框架&#xff0c;小到可以称为“微框架”。Flask 非常小&#xff0c;因此你一旦能够熟练使用它&#xff0c;很可能就能读懂它所有的…
阅读更多...
《MySQL对库的基本操作》

《MySQL对库的基本操作》

文章目录 一、查看数据库列表查看数据库中的所有表想知道当前处于哪个数据库里 二、创建一个数据库三、删除一个数据库知道两个集1.字符集2.校验集修改数据库的字符集和编码集 不同的校验码对数据库的影响四、数据库的备份与恢复注意事项&#xff1a;备份数据库中的表 总结 一、…
阅读更多...
现代神经网络总结(AlexNet VGG GoogleNet ResNet的区别与改进)

现代神经网络总结(AlexNet VGG GoogleNet ResNet的区别与改进)

VGG NIN GoogleNet 1.VGG&#xff0c;NIN&#xff0c;GoogleNet的块结构图对比(注意:无AlexNet) 这些块带来的区别与细节 AlexNet未使用块,主要对各个层进行了解: 卷积:捕捉特征 relu:增强非线性 池化层:减少计算量 norm:规范数据分布 全连接层:分类VGG块的改善(对比AlexNe…
阅读更多...
开源博客项目Blog .NET Core源码学习（18：App.Hosting项目结构分析-6）

开源博客项目Blog .NET Core源码学习（18：App.Hosting项目结构分析-6）

本文学习并分析App.Hosting项目中后台管理页面的_AminLayout.cshtml模版页面和登录页面。 _AminLayout.cshtml模版页面 后台管理页面中的大部分页面都使用_AminLayout.cshtml作为模板页面&#xff0c;如下图所示&#xff0c;后台页面的视图内容放置在表单中&#xff0c;使用la…
阅读更多...
wps用js宏给文档增加用户名密码验证

wps用js宏给文档增加用户名密码验证

偶然看见别人一个给office文档增加验证的vba教程&#xff0c;觉得很有意思&#xff0c;就打开wps验证了一下 wps在windows版本支持vba和js&#xff0c;js默认&#xff0c;vba需要自己下载相关插件折腾&#xff0c;linux版本wps个人版默认没有支持2次开发的高级功能&#xff0c…
阅读更多...
unity想让方法带一个默认参数怎么写

unity想让方法带一个默认参数怎么写

在C#中&#xff0c;包括Unity使用的C#版本&#xff0c;你可以为方法参数提供默认值。这允许你在调用方法时省略某些参数&#xff0c;并使用这些参数的默认值。以下是如何为一个方法参数设置默认值的示例&#xff1a; using UnityEngine; public class MyClass : MonoBehaviou…
阅读更多...
【C++航海王：追寻罗杰的编程之路】C++11(四)

【C++航海王：追寻罗杰的编程之路】C++11(四)

目录 1 -> 相关文章 【C航海王&#xff1a;追寻罗杰的编程之路】C11(一) 【C航海王&#xff1a;追寻罗杰的编程之路】C11(二) 【C航海王&#xff1a;追寻罗杰的编程之路】C11(三) 2 -> lambda表达式 2.1 -> C98中的一个例子 2.2 -> lambda表达式 2.3 ->…
阅读更多...
Python 与 TensorFlow2 生成式 AI（三）

Python 与 TensorFlow2 生成式 AI（三）

原文&#xff1a;zh.annas-archive.org/md5/d06d282ea0d9c23c57f0ce31225acf76 译者&#xff1a;飞龙 协议&#xff1a;CC BY-NC-SA 4.0 第七章&#xff1a;使用 GAN 进行风格转移 神经网络在涉及分析和语言技能的各种任务中正在取得进步。创造力是人类一直占有优势的领域&…
阅读更多...
关于ChatGPT的论文Demo

关于ChatGPT的论文Demo

ChatGPT: 解锁人工智能的无限可能 引言 随着科技的飞速发展&#xff0c;人工智能&#xff08;AI&#xff09;逐渐成为改变世界的强大力量。而在众多AI技术中&#xff0c;ChatGPT以其独特的魅力和广泛的应用领域&#xff0c;吸引了全球的关注。本文将深入探讨ChatGPT的技术原理…
阅读更多...
MetaGPT初体验之HelloWorld-Git教程编写

MetaGPT初体验之HelloWorld-Git教程编写

[目录] 1.环境准备 2.效果预览 3.总结 4.智能体完整输出 5.源码及教程点我去AIGIS公众号查看本文 前言 5.1假期坚持研究智能体的玩法可以说非常敬业了。今天我们来小试一把目前GitHub最火爆智能体框架MetaGPT,让它给我们写一篇Git教程&#xff0c;看看是不是像传说中的那么神奇…
阅读更多...
如何使用KCF算法。

如何使用KCF算法。

KCF&#xff08;Kernelized Correlation Filters&#xff09;算法是一种高效的目标跟踪算法&#xff0c;它结合了核技巧和相关滤波器的思想。以下是使用KCF算法进行目标跟踪的一般步骤&#xff1a; 初始化&#xff1a; 在视频的第一帧中&#xff0c;手动选择或自动检测要跟踪的…
阅读更多...
VSCode编译C++连接lib文件

VSCode编译C++连接lib文件

期货CTP在Windows上需要连接静态链接库&#xff0c;在VS2022一切正常&#xff0c;在VSCode却始终失败。 原因是Windows系统的dll用的vs编译器&#xff0c;导出的Dll没有用extern c &#xff0c;gcc 编译各种坑。 最后通过在VSCode中配置VS2022的编译器&#xff0c;才终于成功。…
阅读更多...
判断循环链表以及其环入口

判断循环链表以及其环入口

文章目录 题目题目链接题目要求 解题思路方法一&#xff1a;哈希表方法二&#xff1a;双指针 进阶思考快指针一次走三步 进阶问题&#xff08;入口点&#xff09;题目链接题目要求问题思路 总结 题目 题目链接 环形链表 题目要求 解题思路 显而易见的是&#xff0c;单纯的遍…
阅读更多...
day7 c++

day7 c++

整理代码 1、unique_ptr 指针 #include <iostream> #include <memory> using namespace std; class Demo {public:Demo(){cout<<"无参构造"<<endl;}~Demo(){cout<<"Demo的析构函数"<<endl;} };int main() {//unique…
阅读更多...
python 笔记：cls VS self

python 笔记：cls VS self

cls&#xff1a; 用于类方法&#xff1a; cls 通常作为类方法&#xff08;用 classmethod 装饰&#xff09;中的第一个参数。它指代调用该方法的类本身&#xff0c;无论该类有没有被实例化访问类级别的属性和方法 通过 cls&#xff0c;可以访问类级别的属性和方法。可以通过 c…
阅读更多...
【IDEA】IDEA自带Maven/JDK，不需要下载

【IDEA】IDEA自带Maven/JDK，不需要下载

IDEA是由Java编写的&#xff0c;为了保证其运行&#xff0c;内部是自带JDK的。IDEA 2021 及 之后的版本是自带Maven的&#xff1a; 视频连接&#xff1a; https://www.bilibili.com/video/BV1Cs4y1b7JC?p4&spm_id_frompageDriver&vd_source5534adbd427e3b01c725714cd…
阅读更多...
理解Linux文件系统

理解Linux文件系统

文章目录 一、引言二、Linux文件系统概述1、文件系统的结构2、文件系统目录树的逻辑结构 二、文件系统的特性1、super block&#xff1a;文件系统的超级块2、inode&#xff1a;文件系统的索引节点3、inode table4、block&#xff1a;文件系统的数据块5、块组描述符表&#xff0…
阅读更多...
链表例题（分割链表）

链表例题（分割链表）

链接&#xff1a; . - 力扣&#xff08;LeetCode&#xff09; 题目描述&#xff1a; 即将小于特定值的节点放在前面&#xff0c;大于特定值的节点接在后面 思路&#xff1a; 我们可以创建两个链表分别存放大于的值和小于的值。5个变量&#xff08;记录链表当前位置的指针st…
阅读更多...
Spring Kafka——基于 Spring Kafka 实现动态管理 Kafka 连接和 topic 的监听

Spring Kafka——基于 Spring Kafka 实现动态管理 Kafka 连接和 topic 的监听

文章目录 使用 Spring Kafka 动态管理 Kafka 连接和主题监听1. 前言2. 简单的消费程序配置3. Spring Kafka 主要的相关类的说明4. KafkaListener 注解的加载执行流程解析5. 动态监听消费订阅的设计与实现 使用 Spring Kafka 动态管理 Kafka 连接和主题监听 文章内容较长&#x…
阅读更多...
云托管 代码 项目

云托管 代码 项目

代码(云)托管平台 国内&#xff1a; Coding&#xff08;https://coding.net&#xff09;&#xff1a;Coding是国内知名的代码托管平台&#xff0c;支持Git和SVN版本控制系统。它提供了代码托管、项目管理、协作开发、持续集成等功能&#xff0c;并且支持与其他开发工具的集成&a…
阅读更多...
最新文章

Copyright @ 2022~2023