前前处理，把后台需要的参数都读出来。写进name.txt

生成新的sch文件，需要在data里追加新sch名字

没改变的井就不用重新卸载关键字里了。

重启动模型：制作出来是空的，得自己加别的东西

模型要准确，否则不好和历史模型对比，历史模型就是从优化开始的那个时间步全都不处理，按照原来的进行。

sch关键字随便加

跑到几月算几月。不改变就按照之前的井跑。软件加载到sch的时间步，未计算就是灰色的。
最好直接利用别人建立好的重启动模型去做，自己不要轻易做模型，一定保证模型不能出错，不能自己乱加，这样结果也好对比。

也可以不用重启动模型，直接使用命令重启动原来的模型。
不管是eclipse还是tnavigater他们的结果几乎是相同的。

1设置在solution关键字，RESTARTDATE可以设置日期，而不是时间步

2. 要和已有的历史制度对比
   从十月优化，十月以后就都是按照原有进行，利用好别人给的重启动模型

3.重启后的sch时间步都要删掉，否则会重复跑出多出自己设置的时间步数。

RESTARTDATE
‘…/HM_FINAL/HM_FINAL’ 1* 01 SEP 2022
/

模型里追加新sch重启动文件，直接跑即可

输出格式转换,按照时间步输出

RPTRST
‘BASIC=3’ ‘FREQ=1’ /

F:\\tNavigator-con-Win-64\\tNavigator-con-Win-64\\tNavigator-con --use-gpu  -e -i -s -r -u –m --server-url=http://103.90.189.57:5805/tnav_cgi1/tnav_license_server.cgi D:\\moxing\\HM_FINAL_RESTART_CASE5_1\\HM_FINAL_RESTART_CASE5_1.datafor %%i in (D:\Injection_production_optimization\item\HM_FINAL_RESTART_CASE5_1\RESULTS\*.S*) do ( 
F:\\tNavigator-con-Win-64\\tNavigator-con-Win-64\\tNavigator-con --use-gpu --server-url=http://103.90.189.57:5805/tnav_cgi1/tnav_license_server.cgi  --convert-ecl-bin-to-text %%i ) for %%i in (D:\Injection_production_optimization\item\HM_FINAL_RESTART_CASE5_1\RESULTS\*.X0*) do ( 
F:\\tNavigator-con-Win-64\\tNavigator-con-Win-64\\tNavigator-con --use-gpu --server-url=http://103.90.189.57:5805/tnav_cgi1/tnav_license_server.cgi  --convert-ecl-bin-to-text %%i ) 

降维方法

每个时间步的变量进行统计，都加起来就是最后的维数。
1.按照井组优化。
2.转注前变量不优化，专注后才算维数。
在随机生成种子的时候就200次130维度。130=10*13口井。根据上下限和拉丁分布算法生成种子。也可以根据变异生成种子，局部效果好，二次规划。

注水井的注水量

只有注水井有

生产井的累积产量

生产井可以产水和产油

对比两个模型的累计产油量

含油饱和度图

是重要指标

所有模式都带入输出变量维度，并装入字典，相当于java map

装入字典map中if horimode != 'None':#initialmode= 'process_needed'sampdic['hw_new_well_name'] = hw_new_well_namesampdic['hw_new_group_name'] = hw_new_group_namesampdic['hw_new_well_r'] = hw_new_well_r#sampdic['hw_new_ninj'] = hw_new_ninjsampdic['hw_new_npro'] = hw_new_nprosampdic['hw_new_nwell'] = sampdic['hw_new_ninj'] + sampdic['hw_new_npro']#sampdic['hw_consmethod'] = hw_consmethod  # 约束处理方法sampdic['hw_layermode'] = hw_layermode#

Xrange, varlen = preRan_Len(Xsampmode, sampdic)

目标函数为npv

python二维数组

井位优化spec

加井，spec里面加。

1.水平井井位优化参数输入，水平井范围

和直井井位优化的前处理区别，水平井优化过程
一口井有多个井位，在射孔那里看
射孔那里水平井多了很多层

整体优化里选择水平井井位优化模式

    # ====================================================================if func == 'Joint':#vertmode = 'None'  # 直井模式选择horimode = 'cooropt-Lfixed-alpopt-4'  # 水平井模式选择transinjmode = 'None'  # 转注模式选择sealmode = 'None'  # 封堵模式选择prodmode = 'None'  # 注采模式选择extramode = 'None'  # 提液模式选择#水平井范围设置if horimode != 'None':# hw_new_well_name = ['New-1', 'New-2', 'New-3']hw_new_group_name = ['G1', 'G1', 'G1']hw_new_well_r = [0.178, 0.178, 0.178]hw_new_ninj = 1  # 新加水平井中注水井数hw_new_npro = 2  # 新加水平井中产油井数# ----与COMPDAT关键字写入时有关的接口（方便新加水平井井位与新加水平井的层位进行综合优化）----hw_layermode = 'all-open'# ----与WCONINJE和WCONPROD关键字写入时有关的接口（方便新加水平井井位与新加水平井的注采进行综合优化）----hw_valuemode = 'allwell-opt'# [1] allwell-same 所有井都是一样的注采量或者压力（此处输入为float数字，后面转化为list形式）# [2] allwell-diff 所有井可以是不一样的注采量或者压力（可以①分井指定，list形式；也可以是②优化获取分井数据，如果与注采优化耦合需要选择[2]的②模式）if hw_valuemode == 'allwell-same':hw_injvalue = 0hw_provalue = 900elif hw_valuemode == 'allwell-opt':hw_value_low = [10, 10, 10]hw_value_up = [1000, 1000, 1000]elif hw_valuemode == 'allwell-diff':hw_injvalue = [100]hw_provalue = [185, 200]# ----与井的选取范围有关的关键字----hw_regionmode = 'allwell-diff'# [1] allwell_same 所有井都是在同一个选定的区域内进行位置优化（矩形区域、需在活网格内）# [2] allwell_diff 每口井一个选的的区域进行位置优化，每口井的位置待选区域可以不同（矩形区域，需在活网格内）if hw_regionmode == 'allwell-same':hw_new_wxlow = 10  # 新加水平井井眼的X坐标范围的下限hw_new_wxup = 80  # 新加水平井井眼的X坐标范围的上限hw_new_wylow = 10hw_new_wyup = 80elif hw_regionmode == 'allwell-diff':hw_new_wxlow = [10, 20, 30]hw_new_wxup = [40, 50, 60]hw_new_wylow = [35, 36, 37]hw_new_wyup = [41, 42, 43]

2.封堵—层位优化参数

open shut关键字

#封堵优化用油藏模型参数 ！！！！！！！！！！！！！
history_timestep = 5
well_name = ['PRO-1', 'PRO-2', 'PRO-3', 'PRO-4', 'PRO-5', 'PRO-6', 'PRO-7', 'PRO-8', 'PRO-9', 'INJ-01', 'INJ-02','INJ-03', 'INJ-04']井类型 1油井0水井
initial_welltype = [1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0]
timestep = 10
max_layer = 12
seal_timestep = [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]
init_status = np.ones((1, max_layer * len(well_name)))//封堵井组分组
seal_list = ['PRO-1', 'PRO-2', 'PRO-3', 'PRO-4', 'PRO-5', 'PRO-6', 'PRO-7', 'PRO-8', 'PRO-9', 'INJ-01', 'INJ-02', 'INJ-03', 'INJ-04']
seal_list_group = [['PRO-1', 'PRO-2'], ['PRO-3', 'PRO-4', 'PRO-5', 'PRO-6'], ['PRO-7', 'PRO-8', 'PRO-9'], ['INJ-01','INJ-02','INJ-03','INJ-04']]
seal_list_group_num = len(seal_list_group)
maxnum_seal = 2
well_group_name = ['G1', 'G1', 'G1', 'G1', 'G1', 'G1', 'G1', 'G1', 'G1', 'G2', 'G2', 'G2', 'G2']
well_r = [0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5]
井位
well_location = [[1, 1], [13, 1], [25, 1], [1, 13], [13, 13], [25, 13], [1, 25], [13, 25], [25, 25], [7, 7], [19, 7],[7, 19], [19, 19]]
well_hori_location = np.tile(np.array(well_location), (1, max_layer))
transinj_list = ['PRO-1', 'PRO-2', 'PRO-3', 'PRO-4', 'PRO-5', 'PRO-6', 'PRO-7', 'PRO-8', 'PRO-9']
transinj_list_group = [['PRO-1', 'PRO-2', 'PRO-3', 'PRO-4'], ['PRO-5', 'PRO-6', 'PRO-7', 'PRO-8', 'PRO-9']]
pro_opt_list = ['PRO-1', 'PRO-2', 'PRO-3', 'PRO-4', 'INJ-01']
history_pro = np.ones((len(well_name)))*350  ##跟well_name对应，np
low = np.zeros((len(well_name)))
up = np.ones((len(well_name))) * 700
up[9] = 1500
up[10] = 1500
up[11] = 1500
up[12] = 1500
prod_low = 0  ##跟well_name对应，np
prod_up = 700
inj_low = 0
inj_up = 1500
extra_low = 50
extra_up = 100
//分为五层
segment = [[1, 2], [3, 4], [5], [6, 7, 8, 9, 10], [11, 12]]
specific_seal = [[[0, 0, 1, 1, 0], [1, 1, 0, 1, 0], [0, 1, 0, 1, 0]],[[0, 0, 1, 1, 0], [1, 1, 0, 1, 0], [0, 1, 0, 1, 0]],[[0, 0, 1, 1, 0], [1, 1, 0, 1, 0], [0, 1, 0, 1, 0]],[[0, 0, 1, 1, 0], [1, 1, 0, 1, 0], [0, 1, 0, 1, 0]],[[0, 0, 1, 1, 0], [1, 1, 0, 1, 0], [0, 1, 0, 1, 0]],[[0, 0, 1, 1, 0], [0, 1, 0, 1, 0]],[[0, 0, 1, 1, 0], [1, 1, 0, 1, 0]]]
segment_num = 5
timestep_len = 180
maxnum_tran = 2
objmode = 'NPV'

井位*层数



特殊封堵处理

if max_layer > 10:for well_num in range(len(specific_seal)):for i in range(len(specific_seal[well_num])):temp_scheme = np.ones((max_layer))count = 0for l in range(segment_num):for num in range(len(segment[l])):temp_scheme[count] = specific_seal[well_num][i][l]count = count + 1specific_seal[well_num][i] = temp_scheme

3.转注优化过程参数输入

water oil关键字，采油井转为注水井

参数if transinjmode == 'None':transinj_list = []if transinjmode == 'specific_max':maxnum_tran = 3转注+封堵一起操作elif transinjmode == 'specific_joint':specific_seal = []  #list,指定的封堵方案，与transinj_list的顺序对应maxnum_tran = 2seal_timestep = np.ones((len(well_name), 1))大于10段，分段来减少维数if max_layer > 10:segment = []  #list，分段segment_num = []  #段数

4.注采优化过程参数输入

BHP控制和液体控制

本示例为两个井设置控制。1043井处于注入速率控制，速率等于253.1。没有底部压力限制，所以模拟器将设置注入速率为253.1。如果用户指定的注入速率如此之高，以至于必要的底部孔压力超过默认极限（6801patm），该井将自动切换到必和必拓控制。1054井也有注入速率控制，速率等于253.1。该井的井底孔压力限制为400；因此，如果计算压力超过400，该井将自动切换到必和必拓控制。

                if prodmode != 'None':low = []  # list，注采下界，与well_name顺序对应up = []  # list，注采上界，与well_name顺序对应if prodmode == 'specific':history_pro = [] #历史注采量/压力，与well_name顺序对应pro_opt_list = [] #进行注采优化的井

选择注采模式，层位和注采、转注等可以一起优化。

13口井10个时间步

生成变量

            Xrange, varlen = preRan_Len(Xsampmode, sampdic)sampdic['varlen'] = varlen

5.提液模式

注采上下限特殊处理，直接上限设置好就行

6. 直井井位优化参数输入

wellspec关键字，新加井就是井位优化

            if vertmode != 'None':# vw_new_well_name = ['NEW-1', 'NEW-2', 'NEW-3']vw_new_group_name = ['GROUP1', 'GROUP1', 'GROUP1']vw_new_well_r = [0.178, 0.178, 0.178]vw_new_ninj = 1  # 新加直井中注水井数vw_new_npro = 2  # 新加直井中产油井数# ----与COMPDAT关键字写入时有关的接口（方便新加直井井位与新加直井的层位进行综合优化）----vw_layermode = 'all-open'  # 选择不同模式，对应ecl_joint中的xp处的生成传递# [1] all-open 默认新加井的所有层打开，从第1层到最后一层（max_layer）# ----与WCONINJE和WCONPROD关键字写入时有关的接口（方便新加直井井位与新加直井的注采进行综合优化）----vw_valuemode = 'allwell-same'# [1] allwell-same 所有井都是一样的注采量或者压力（此处输入为float数字，后面转化为list形式）# [2] allwell-diff 所有井可以是不一样的注采量或者压力（可以①分井指定，list形式；也可以是②优化获取分井数据，如果与注采优化耦合需要选择[2]的②模式）#if vw_valuemode == 'allwell-same':vw_injvalue = 100vw_provalue = 500elif vw_valuemode == 'allwell-diff':vw_injvalue = [100, 120, 300]vw_provalue = [125, 185, 200]# vw_inj_controlmode = 'RATE'    # 新加的注水直井的控制模式选择# [1] RATE 流量控制# [2] BHP 井底流压控制 vw_pro_controlmode = 'LRAT'    # 新加的采油直井的控制模式选择# [1] LRAT 液量控制# [2] BHP 井底流压控制 指定了一个矩形范围内优化if vertmode == 'specific_region':# specific_region 指定优化区域，单独优化加井的井位。这种模式下直井井眼的选择范围是预先选定的一个规则矩形（需要在活网格内）vw_new_wxlow = 1    # 新加直井的X坐标范围的下限vw_new_wxup = 50   # 新加直井的X坐标范围的上限vw_new_wylow = 20vw_new_wyup = 100

新加水平井#hw_consmethod = 'notcons-del-relhsget' # 约束处理方法# [1] notcons-del-relhsget 不满足约束的样本删除，并利用拉丁超立方生成horimode = 'cooropt-Lfixed-alpopt-4'if horimode == 'cooropt-Lfixed-alpopt-4':#hw_new_wxlow = 25    # 新加水平井井眼的X坐标范围的下限hw_new_wxup = 55   # 新加水平井井眼的X坐标范围的上限hw_new_wylow = 35hw_new_wyup = 65#hw_new_wdim = 3  # 单口井的优化维数hw_new_wL = 500 / 100  # 水平段长度hw_new_walplow = 0.51  # 1hw_new_walpup = 4.49  # 4#elif horimode == 'cooropt-Lopt-alpopt-4':#hw_new_wxlow = 25    # 新加水平井井眼的X坐标范围的下限hw_new_wxup = 55   # 新加水平井井眼的X坐标范围的上限hw_new_wylow = 35hw_new_wyup = 65hw_new_wdim = 4  # 单口井的优化维数hw_new_wLlow = 251 / 100  # 3hw_new_wLup = 549 / 100  # 5hw_new_walplow = 0.51  # 1hw_new_walpup = 4.49  # 4

python 报错

策略优化模式选择：优化主要分类6类

除了公共参数输入，还有各种模式的自己的参数。井的参数都随着时间步变化。

        # ====================================================================if func == 'Joint':# vertmode = 'None'  # 直井井位优化horimode = 'cooropt-Lfixed-alpopt-4'  # 水平井井位优化transinjmode = 'None'  # 转注模式选择sealmode = 'None'  # 层位封堵模式选择prodmode = 'None'  # 注采模式选择extramode = 'None'  # 提液模式选择#if vertmode != 'None':# 

每种模式的变量数，注采是一起的

油藏参数处理：面向井对象，井单位（井位（xy上下限）、井类型（01就是上下限）、注采量（上下限）、射孔（井数总层数01上下限））+分优化模式（针对不同的井组）+是否随时间步变（井数时间步）

优化必有上下界，一个上下界确定一个变量

如果不涉及优化就给定一个固定值特殊处理，或者不操作。
一般不根据时间步优化，维数太大，调控一次即可。specific_time和这是specific的差别

目标函数，公共参数和各种模式的参数，算法
map
所有要优化的变量都是输入参数
timestep_history要历史注采量干嘛，水平井的模式，
1.先选目标函数
分为模拟和函数，模拟里面有水平井，联合优化，tnavigater优化

2.选输入公共参数和各种模式的参数赋值
直井综合优化里面有油藏模型基本信息和很多模式，每个模式的输入参数不同

timestep_history = 0
well_name = [‘PRO-1’, ‘PRO-2’, ‘PRO-3’, ‘PRO-4’, ‘PRO-5’, ‘PRO-6’, ‘PRO-7’, ‘PRO-8’, ‘PRO-9’, ‘INJ-1’, ‘INJ-2’,
‘INJ-3’, ‘INJ-4’]
1是油井
initial_welltype = [1, 1, 1, 1, 1, 1, 1, 1, 1,0, 0, 0, 0]
max_layer = 1
init_status = np.ones((1, max_layer * len(well_name)))
well_group_name = [‘G1’, ‘G1’, ‘G1’, ‘G1’, ‘G1’, ‘G1’, ‘G1’, ‘G1’, ‘G1’, ‘G1’, ‘G1’, ‘G1’, ‘G1’]
well_r = [0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.5, 0.5, 0.5, 0.5]
well_location = [[1, 1], [1, 51], [1, 101], [51, 1], [51, 51], [51, 101], [101, 1], [101, 51], [101, 101],
[77, 81], [74, 18],
[23, 94], [71, 5]]
up = np.ones((len(well_name), 1)) * 1800
low = np.ones((len(well_name), 1)) * 10
up[9] = 4050
up[10] = 4050
up[11] = 4050
up[12] = 4050

    # ----优化的信息----时间步变量timestep = 10timestep_len = 180objmode = 'NPV'# ----路径设置有关----casefolder_name = 'HORI'rundata_name = 'HORI'

字典套字典Xrange

solvepara[‘realfunc’] = func # 适应度函数的具体选择 ---- 字符串形式
solvepara[‘Xrange’] = Xrange # 样本的上下界范围 ---- np二维数组，size为（2，ndim）
solvepara[‘sampdic’] = sampdic

综合优化和tnavigater优化的直井输入一样。

3.选择算法

优化变量处理 1.确定上下界2.变量归一化0-1之间的随机值 2.（上界-下界）*归一随机值

明确要优化的变量是谁
拉丁分布前处理

    #上下限确定，维数长度确定Xrange, varlen = preRan_Len(Xsampmode, sampdic)sampdic['varlen'] = varlen根据时间步变化和不根据时间步变化
elif prodmode == 'specific':pro_opt_list = sampdic['pro_opt_list']len_prod = len(pro_opt_list) * timestepran_prod = np.vstack([np.zeros((1, len_prod)), np.ones((1, len_prod))])
elif prodmode == 'specific_time':specific_time = sampdic['specific_time']len_prod = len(specific_time)ran_prod = np.vstack([np.zeros((1, len_prod)), np.ones((1, len_prod))])

归一化

乘以上界减去下界

预计算：计算各个模式加起来的总共维度、上下限

统计各个模式加起来的总共维度

    if prodmode == 'None':len_prod = 0ran_prod = np.zeros((2, len_prod))elif prodmode == 'wgiven-optTgiven':len_prod = nwell * timestepran_prod = np.vstack([np.zeros((1, len_prod)), np.ones((1, len_prod))])elif prodmode == 'specific':pro_opt_list = sampdic['pro_opt_list']len_prod = len(pro_opt_list) * timestepran_prod = np.vstack([np.zeros((1, len_prod)), np.ones((1, len_prod))])# ----提液相关----#!!!!!!!!!!!!!!!这一部分还没写#varlen = {'len_vert' : len_vert,'len_hori' : len_hori,'len_transinj' : len_transinj,'len_seal' : len_seal,'len_prod' : len_prod# 'len_extra' : len_extra}#Xrange = np.hstack([ran_vert, ran_hori, ran_transinj, ran_seal, ran_prod])# Xrange = np.hstack([ran_vert, ran_hori, ran_transinj, ran_seal, ran_prod, ran_extra])return Xrange, varlen

 #生成随机值0-1Xrange, varlen = preRan_Len(Xsampmode, sampdic)sampdic['varlen'] = varlen# ----拉丁超立方初始种群模式initmode产生方式---- #!!!!!if initialmode == 'process_needless':initmode = 'Basic'elif initialmode == 'process_needed':# 如果lhs生成的初始解需要进行进一步处理（删除、重新生成等），这个时候需要根据不同的特定模式执行不同的处理操作if horimode != 'None': initmode = 'Joint-Hori'# elif #!!! 【留的接口】# 如果有其他模式（比如直井不规则的四个点判断井位区域、实际油藏活死网格不规则边界判断等），在这个接口后面加入elif，定义新的initmode，并在latin函数中增加与新定义的initmode对应的新模块

生成随机值
def wellprod_proce(X, well_type, optional):mode = optional['Xsampmode']sampdic = optional['sampdic']xp = copy.deepcopy(X)if mode['prodmode'] == 'wgiven-optTgiven':xp = xp.reshape(-1, len(sampdic['well_name']))for time in range(sampdic['timestep']):for well in range(len(sampdic['well_name'])):xp[time, well] = sampdic['low'][well] + xp[time, well] * (sampdic['up'][well] - sampdic['low'][well])elif mode['prodmode'] == 'specific':history_pro = sampdic['history_pro']pro_opt_list = sampdic['pro_opt_list']xp = xp.reshape(-1, len(sampdic['pro_opt_list']))temp = np.ones((sampdic['timestep'], len(sampdic['well_name'])))for time in range(sampdic['timestep']):for well in range(len(sampdic['well_name'])):for opt_well in range(len(pro_opt_list)):if pro_opt_list[opt_well] == sampdic['well_name'][well]:temp[time, well] = xp[time, opt_well]xp = tempfor time in range(sampdic['timestep']):for well in range(len(sampdic['well_name'])):flag = 0for opt_well in range(len(pro_opt_list)):if pro_opt_list[opt_well] == sampdic['well_name'][well]:flag = 1breakif flag == 1:xp[time, well] = sampdic['low'][well] + xp[time, well] * (sampdic['up'][well] - sampdic['low'][well])else:xp[time, well] = history_pro[well]xp = xp.reshape(1, -1)return xp

设置初始位置和优化位置，不同

在这里插入图片描述](https://img-blog.csdnimg.cn/7977d85e1b744a2db1afcb47b9518394.png)

选择算法，并行，跑单次

并行代码，i为并行的次序
def tnav_joint(X, optional, i, index=''):X = latin(Xrange, npop, initmode, solvepara)#Y = Parallel(n_jobs=1)(delayed(ecl_joint_parallel)(X, solvepara, i, str(i % (2 * cpu_count()))) for i in range(npop))Y = np.array(Y).reshape(-1, 1)
调用ecl_joint_parallel跑单次模拟

单次模拟：前处理 处理 后处理

循环写入井控关键字
复制文件夹，制作xp，写入文件夹。跑data、读结果、算npv。

  para = optional['sampdic']timestep = para['timestep']if os.path.exists(simupath):  # 判断括号里的文件夹路径或文件路径是否存在，存在返回Ture，不存在返回Falsermtree(simupath)  # 递归删除一个目录以及目录内的所有内容copytree(initpath,simupath)  # shutil.copytree( olddir, newdir, True/Flase) #把olddir拷贝一份newdir，如果第3个参数是True，则复制目录时将保持文件夹下的符号连接，如果第3个参数是False，则将在复制的目录下生成物理副本来替代符号连接#xp = {}#if vertmode != 'None':xs_vert = copy.deepcopy(XS['Xvert'][i, :])wvert = wellvert_proce(xs_vert, optional)xp['WX_vert'] = wvert['wx_vert']xp['WY_vert'] = wvert['wy_vert']xp['Wstatus_vert'] = wvert['wstatus_vert']xp['Wpro_vert'] = wvert['wpro_vert']xp['Winj_vert'] = wvert['winj_vert']
分时间步写writeSCH(simupath, optional)
井位优化WELSPECS(simupath, xp, optional)
层位优化（封堵）COMPDAT(simupath, xp, optional)
转注和注采优化、提液模式WCONINJE(simupath, xp, optional)WCONPROD(simupath, xp, optional)

根据时间步写入每个文件sch文件

判断开始时间步是否写入，判断井类型写油水井。
每个时间步用多个循环写。

6直井井位优化模式，优化xy变量，

每个井都有个x和y位置，可随着时间步变化。

前处理，数据格式的解释，二维数组的处理（三维数组）
井类型数组是13*10，13口井和10个时间步，每个时间步的值对应好了

循环写入井定义文件
这是一个时间步的，井类型根据时间步变化，其他不变
只写include文件，时间步在sch文件里面写了

2.封堵射孔（层位优化模式原理）优化well_status变量，每个井的每一层都有个状态，随着时间步变化，也可以不变化。

前处理三重循环：时间步、井名、层数
优化well_status变量，每个井的每一层都有个状态，随着时间步变化，也可以不变化。
首先，在重启动sch中，手动把data里面写进去新得sch文件，手动改summary

层位优化模式分为两种，按照时间步优化和不按照时间步优化。
在这56个方案里面随机选择层位方案。0-56方案就是上下限。有了上下限就可以生成具体的值了（单次维数*种群数）

种群*井数

井状态

well_status井状态随着时间步变化，其他的不变
well_status = xp[‘Wstatus’] Wstatus就是选中的方案
没有层位优化则默认全打开

 def COMPDAT(simupath, xp, optional):''':param well_hor_location: numpy,dim:井数*（最大层数*2）第一行为：[第一个井在第一层的横坐标，第一个井在第一层的纵坐标，第一个井在第二层的横坐标，第一个井在第二层的纵坐标，....]:param simulatuon_folder: 文件路径:param timestep: int,总时间步，一个数:param well_status: numpy,所有时间步内所有井的状态；dim：总时间步*(总井数*总层数），要求为numpy数组！:param well_r: numpy,所有井的井眼直径，油井在前水井在后；dim：1*总井数:param well_name: list,所有井的井名，与well_status和well_r相对应:return: 重新编写Compdat_ctrl+timestep .dat 文件(文件名里的timestep从1开始)注：所有有关井的变量井的顺序都与well_name相对应'''well_status = xp['Wstatus']sampdic = optional['sampdic']total_well_num = len(sampdic['well_name'])for time in range(sampdic['timestep']):wellcon_file = open(simupath + '/Compdat' + str(time + 1) + '.dat', 'w')contents = ['COMPDAT\n']for well in range(total_well_num):for layer in range(sampdic['max_layer']):if well_status[time, well * sampdic['max_layer'] + layer] == 1:status = '\'OPEN\''else:status = '\'SHUT\''contents.append(sampdic['well_name'][well] + '   ' + str(sampdic['well_hori_location'][well, layer * 2]) + '  '+ str(sampdic['well_hori_location'][well, layer * 2 + 1])+ '  ' + str(layer + 1) + '  ' + str(layer + 1) + '  ' +status + '  2*  ' + str(sampdic['well_r'][well]) + '  3*  \'Z\' 1*  / \n')contents.append('/\n')wellcon_file.writelines(contents)wellcon_file.close()

井定义和注采啥都不写，按照原来的跑
只写个WCONINJE

    if prodmode == 'None' and horimode == 'None' and vertmode == 'None':for time in range(sampdic['timestep']):iline = ['WCONINJE \n']iline.append('/ \n')# ----先获取这个时间步要写入的WCONPROD的地址----wconinje_path = os.path.join(simupath, 'Wconinje' + str(time + 1) + '.dat')#with open(wconinje_path, 'a+') as f:  # !!!!!!这里a+还是w+模式的选择f.writelines(iline)

而射孔是按照大层分段层都写一遍

先根据层上下限随机生成数据，预处理得到01射孔矩阵

    if sealmode != 'None':xs_seal = copy.deepcopy(XS['Xseal'][i, :])xp['Wstatus'] = wellseal_proce(xs_seal, xs_transinj, optional)  # 结合wellstatus_generate和范围

层状态wstaue表示，1是开open，10156，10是时间步，156=1312，一一对应了，这里可以写死。
init_status = np.ones((1, max_layer * len(well_name)))
写入射孔文件中

特殊模式：不优化的写死的，27口井需要特殊处理，分别对应，每口井的关闭层和打开层。以井为研究对象。一口井就是一个对象，它的属性有处理时间步、定义、射孔、注采、井位等处理。

1水平井井位优化，xy，长度、角度，四个变量*井数的优化

水平井参数

转注（优化模式）优化注采变量，每个井都有井类型，随着时间步变化

4采油（注采优化模式）优化注采变量，每个井都有个主材量，随着时间步变化

前处理井两层循环：时间步、井名

注采量Wprod维数是一维数组130=13口井*10个时间步，一一对应了。关键是上下界。

# ----注采相关----
if prodmode == 'None':len_prod = 0ran_prod = np.zeros((2, len_prod))
elif prodmode == 'wgiven-optTgiven':len_prod = nwell * timestepran_prod = np.vstack([np.zeros((1, len_prod)), np.ones((1, len_prod))])

把各个模式的变量合在一起，先确定上下界。varlen = {'len_vert' : len_vert,'len_hori' : len_hori,'len_transinj' : len_transinj,'len_seal' : len_seal,'len_prod' : len_prod# 'len_extra' : len_extra}#Xrange = np.hstack([ran_vert, ran_hori, ran_transinj, ran_seal, ran_prod])# Xrange = np.hstack([ran_vert, ran_hori, ran_transinj, ran_seal, ran_prod, ran_extra])return Xrange, varlen

根据上下界生产随机值wellprod_proce
def wellprod_proce(X, well_type, optional):mode = optional['Xsampmode']sampdic = optional['sampdic']xp = copy.deepcopy(X)if mode['prodmode'] == 'wgiven-optTgiven':xp = xp.reshape(-1, len(sampdic['well_name']))for time in range(sampdic['timestep']):for well in range(len(sampdic['well_name'])):xp[time, well] = sampdic['low'][well] + xp[time, well] * (sampdic['up'][well] - sampdic['low'][well])

只对特定井名进行优化elif mode['prodmode'] == 'specific':history_pro = sampdic['history_pro']pro_opt_list = sampdic['pro_opt_list']xp = xp.reshape(-1, len(sampdic['pro_opt_list']))temp = np.ones((sampdic['timestep'], len(sampdic['well_name'])))for time in range(sampdic['timestep']):for well in range(len(sampdic['well_name'])):for opt_well in range(len(pro_opt_list)):if pro_opt_list[opt_well] == sampdic['well_name'][well]:temp[time, well] = xp[time, opt_well]xp = tempfor time in range(sampdic['timestep']):for well in range(len(sampdic['well_name'])):flag = 0for opt_well in range(len(pro_opt_list)):if pro_opt_list[opt_well] == sampdic['well_name'][well]:flag = 1breakif flag == 1:xp[time, well] = sampdic['low'][well] + xp[time, well] * (sampdic['up'][well] - sampdic['low'][well])else:xp[time, well] = history_pro[well]

返回随机值=10时间步*井数

如果井类型是1油井
def WCONPROD(simupath, xp, optional):well_type = xp['Wtype']pro = copy.deepcopy(xp['Wprod'])sampdic = optional['sampdic']total_well_num = len(sampdic['well_name'])pro = pro.reshape(-1, len(sampdic['well_name']))for time in range(sampdic['timestep']):f = open(simupath + '/Wconprod' + str(time + 1) + '.dat', 'w')f.write('WCONPROD \n')for well in range(total_well_num):如果井类型是1油井if well_type[well][time] == 1:f.write('\'' + sampdic['well_name'][well] + '\'' + "    OPEN   LRAT    3*    " + str(pro[time, well]) + '/\n')f.write('/\n')f.close()

如果井类型是0水井
def WCONINJE(simupath, xp, optional):sampdic = optional['sampdic']well_type = xp['Wtype']pro = copy.deepcopy(xp['Wprod'])total_well_num = len(sampdic['well_name'])pro = pro.reshape(-1, len(sampdic['well_name']))for time in range(sampdic['timestep']):f = open(simupath + '/Wconinj' + str(time + 1) + '.dat', 'w')f.write('WCONINJE \n')for well in range(total_well_num):如果井类型是0水井if well_type[well][time] == 0:f.write('\'' + sampdic['well_name'][well] + '\'' + "     WATER    OPEN    RATE    " + str(pro[time, well]) + '/\n')f.write('/\n')f.close()

数值模拟，换模拟器，要有模拟器位置。模拟器限制很多

tnavigater不一样就把他定死summary顺序。
最新破解版关键字必须有
MULTOUT它
RPTRST
‘BASIC=3’ ‘FREQ=1’ /

SUMMARY   =============================================================
EXCEL
FWPR
FOPR
FWCT
FWIR
FOPT
FWPT
FWIT
FPR
/# ----进行数值模拟----DATA_path = os.path.join(simupath, rundata_name + '.DATA')command = 'eclrun eclipse' + ' ' + DATA_pathos.system(command)print('数值评估', i + 1, '次')# ----读取结果----output, outnum = read_sumkeywd(DATA_path)  # 读取DATA文件SUMMARY部分关键字SCH_path = fullfile(simupath, 'schedule.dat')# steps = get_schsteps(SCH_path)    # 查找schedule文件中的时间步数'''steps = para['timestep'] + para['history_timestep']'''## timestep_n = 30# timesteplen = sampdic['timestep_len']# tn = int(timesteplen/timestep_n)# steps = para['timestep']*tn + para['history_timestep']#steps = para['timestep'] + para['timestep_history']A_path = find_A(simupath, rundata_name, steps)  # 确定要读取的结果的文件名称（A文件）# 在确定的A文件中读取出相关结果(有两种思路，一种是和DATA文件读出类似，先全部读出来再查找和抽取，一种是利用指针，利用A文件特点，需要的结果在最后几行，直接指针到最后几行)with open(A_path, 'rb') as f:  # 要用seek()方法需要用二进制的读，所以用’rb‘而不是’r‘outnum = outnum + 2  # 在A文件输出的关键字数目加2（前面还有天和年两个关键字）# 计算需要将指针向前移动的字符距离movecha = (outnum // 4) * (17 * 4 + 2) + (outnum % 4) * 17 + 2  # //是取整符号，%是取余符号，计算的数目是根据A文件中每个SUMMARY关键字输出结果所占的字符数f.seek(-movecha, 2)  # seek(offset, whence)，offset是文件指针偏移量，whence为0是文件开头，1是文件当前位置，2是文件末尾outext = f.read().decode().split()  # 二进制形式read()出来需要用decode()进行解码，并用split()进行分开del outext[0:2]  # 删去前两个元素（不是SUMMARY中规定的关键字结果）outnum = outnum - 2  # 删去前两个元素的数量for j in range(outnum):order = output[j] + ' = ' + str(float(outext[j]))  # 把结果用order命令逐个写入，并用exec运行exec(order, globals())  # exec()函数的作用域问题，在自己定义的子函数里面使用exec()去define一个变量，后续再使用时会报错，需要加入globals()#objmode = para['objmode']if objmode == 'FOPT':SimRes[i] = FOPTelif objmode == 'NPV':#oil_price = 100  # 油价(默认单位为美元/桶)cost_prowat = 42  # 处理水成本(默认单位为美元/桶)cost_injwat = 2  # 注水成本(默认单位为美元/桶)#units = 'METRIC'  # Eclipse模型的单位进制if units == 'FIELD':npv = FOPT * oil_price - FWPT * cost_prowat - FWIT * cost_injwat

三种调用方式

    conexe = "F:\\tNavigator-con-Win-64\\tNavigator-con-Win-64\\tNavigator-con.exe"config = "F:\\tNavigator-con-Win-64\\tNavigator-con-Win-64\\license.conf"conexe_sever = "C:\\Users\\Administrator\\AppData\\Local\\RFD\\tNavigator\\22.3\\tNavigator-con"config_sever = 'http://103.90.189.57:5805/tnav_cgi1/tnav_license_server.cgi'

读取结果、历史时间步是为了读取到当前最后一个时间步

历史时间步是old里的时间步

bug

mpi并行环境安装

算法

自适应代理

针对维度和问题鲁棒性的
多代理的组合混合，但是会选择

并行变异

两种思路可以结合：
1.利用DE和代理模型并行产生多种x结果，再从x中选择几个进行并行跑。
2.并行多代理（多代理不组合，各跑各的），同时采用多种代理模型，全局和局部、截断学习，等等。