CN107218019B - 一种基于蒙特卡洛算法的聚合物驱生产优化方法和系统 - Google Patents

一种基于蒙特卡洛算法的聚合物驱生产优化方法和系统 Download PDF

Info

Publication number
CN107218019B
CN107218019B CN201710056212.5A CN201710056212A CN107218019B CN 107218019 B CN107218019 B CN 107218019B CN 201710056212 A CN201710056212 A CN 201710056212A CN 107218019 B CN107218019 B CN 107218019B
Authority
CN
China
Prior art keywords
variable
monte carlo
vector
control
optimization
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN201710056212.5A
Other languages
English (en)
Other versions
CN107218019A (zh
Inventor
赵辉
谢晓庆
张贤松
康晓东
石爻
曾杨
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
China National Offshore Oil Corp CNOOC
CNOOC Research Institute Co Ltd
Original Assignee
China National Offshore Oil Corp CNOOC
CNOOC Research Institute Co Ltd
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by China National Offshore Oil Corp CNOOC, CNOOC Research Institute Co Ltd filed Critical China National Offshore Oil Corp CNOOC
Priority to CN201710056212.5A priority Critical patent/CN107218019B/zh
Publication of CN107218019A publication Critical patent/CN107218019A/zh
Application granted granted Critical
Publication of CN107218019B publication Critical patent/CN107218019B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • EFIXED CONSTRUCTIONS
    • E21EARTH OR ROCK DRILLING; MINING
    • E21BEARTH OR ROCK DRILLING; OBTAINING OIL, GAS, WATER, SOLUBLE OR MELTABLE MATERIALS OR A SLURRY OF MINERALS FROM WELLS
    • E21B43/00Methods or apparatus for obtaining oil, gas, water, soluble or meltable materials or a slurry of minerals from wells
    • E21B43/16Enhanced recovery methods for obtaining hydrocarbons
    • GPHYSICS
    • G05CONTROLLING; REGULATING
    • G05BCONTROL OR REGULATING SYSTEMS IN GENERAL; FUNCTIONAL ELEMENTS OF SUCH SYSTEMS; MONITORING OR TESTING ARRANGEMENTS FOR SUCH SYSTEMS OR ELEMENTS
    • G05B13/00Adaptive control systems, i.e. systems automatically adjusting themselves to have a performance which is optimum according to some preassigned criterion
    • G05B13/02Adaptive control systems, i.e. systems automatically adjusting themselves to have a performance which is optimum according to some preassigned criterion electric
    • G05B13/04Adaptive control systems, i.e. systems automatically adjusting themselves to have a performance which is optimum according to some preassigned criterion electric involving the use of models or simulators
    • G05B13/042Adaptive control systems, i.e. systems automatically adjusting themselves to have a performance which is optimum according to some preassigned criterion electric involving the use of models or simulators in which a parameter or coefficient is automatically adjusted to optimise the performance
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00Computer-aided design [CAD]
    • G06F30/20Design optimisation, verification or simulation
    • G06F30/23Design optimisation, verification or simulation using finite element methods [FEM] or finite difference methods [FDM]
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F2111/00Details relating to CAD techniques
    • G06F2111/04Constraint-based CAD
    • YGENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02ATECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
    • Y02A10/00TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE at coastal zones; at river basins
    • Y02A10/40Controlling or monitoring, e.g. of flood or hurricane; Forecasting, e.g. risk assessment or mapping

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Mining & Mineral Resources (AREA)
  • Evolutionary Computation (AREA)
  • Geology (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Fluid Mechanics (AREA)
  • Health & Medical Sciences (AREA)
  • General Engineering & Computer Science (AREA)
  • Geometry (AREA)
  • Computer Hardware Design (AREA)
  • General Life Sciences & Earth Sciences (AREA)
  • Geochemistry & Mineralogy (AREA)
  • Environmental & Geological Engineering (AREA)
  • Artificial Intelligence (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Medical Informatics (AREA)
  • Software Systems (AREA)
  • Automation & Control Theory (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

本发明涉及一种基于蒙特卡洛算法的聚合物驱生产优化方法和系统,其步骤:根据油藏数值模拟器生成油藏优化的最优控制函数;由求模型生成协方差矩阵,利用蒙特卡洛逼近算法在最优控制变量附近生成初始扰动向量,并采用协方差矩阵对初始扰动向量进行处理,生成目标扰动向量;采用包括目标扰动向量的蒙特卡洛算法对最优控制函数进行求解,生成最优开发效益对应的最优注采参数。本发明计算过程简单,容易实现,而且优化效率高,收敛速度快,易于和任意油藏数值模拟器相结合进行聚合物驱生产优化的计算,能够同时优化注采井生产参数和注入聚合物浓度,能显著改善油藏开发效果,为科学合理的进行油田高效开发提供依据。本发明能广泛在油气田开发领域中应用。

Description

一种基于蒙特卡洛算法的聚合物驱生产优化方法和系统
技术领域
本发明涉及一种油气田开发领域,特别是关于一种基于蒙特卡洛算法的聚合物驱生产优化方法和系统。
背景技术
国内外油田经过多年水驱开发,已进入高含水阶段,水驱驱油效率明显下降,严重影响了开发效果。因此,为了对老油田进行再次开发,三次采油逐步推广,其中聚合物驱是最为重要的一种改善波及效率进而提高开发效果的三次采油技术。尽管聚合物驱成本较高,但因其能够实现明显降水增油或者控水增油的效果,也广被工程师们所接受。正是由于该技术具有成本较高的缺点,因此,需要通过油藏开发实时生产优化技术,节约成本,改善聚合物驱开发,最大化开发收益和效果。
油藏开发生产实时优化是近年来新兴的一种开发方案预测方法,其目前主要应用于水驱开发油藏,与聚合物驱应用结合研究目前较少。该方法与聚合物驱结合是基于对当前油藏地质和生产条件的认识,并从油藏长期开发效益出发,利用油藏数值模拟和最优控制理论自动优化计算油水井不同阶段的注采参数(如注入聚合物浓度、注入聚合物时机、井底流压、油水井流量等),确定最优的生产开发方案,使油藏开发尽可能处于最佳状态,从而节约聚合物驱生产成本,改善开发效果。油藏生产优化属于大系统最优化控制问题,涉及到的变量多,维数高,实现求解十分困难。目前主要使用的伴随梯度方法进行优化计算,该方法需要了解油藏模拟器代码层面实现,通过编写伴随矩阵嵌入油藏数值模拟计算来获取伴随梯度,求解过程异常复杂,仅限于一些概念油藏模型进行理论研究与应用,难以推广应用于成熟的商业模拟器。因此,无梯度类方法成为生产优化算法的热点,EnOpt算法、NEWUOA算法、QIM-AG算法等被应用水驱生产优化,还未涉及聚合物驱生产优化。现有技术中利用随机扰动逼近算法(SPSA)结合有限差分算法(FDSA)对聚合物用量进行了优化实现了优化后的调控方案在使用了同样的聚合物用量的情况下,降水增油效果更加明显,改善了聚合物驱的效果。但作为随机扰动类逼近算法所生成的生产调控方案在时间上变化十分剧烈不够平滑,难以在现场实际应用。
发明内容
针对上述问题,本发明的目的是提供一种基于蒙特卡洛算法的聚合物驱生产优化方法和系统,其容易实现,而且优化效率高,收敛速度快,易于和任意油藏数值模拟器相结合进行聚合物驱生产优化的计算,能够同时优化注采井生产参数和注入聚合物浓度。
为实现上述目的,本发明采取以下技术方案:一种基于蒙特卡洛算法的聚合物驱生产优化方法,其特征在于包括以下步骤:1)根据油藏数值模拟器生成油藏优化的最优控制函数;2)由球形模型生成协方差矩阵,利用蒙特卡洛逼近算法在最优控制变量附近生成初始扰动向量,并采用协方差矩阵对初始扰动向量进行处理,生成目标扰动向量;3)采用包括目标扰动向量的蒙特卡洛算法对最优控制函数进行求解,生成最优开发效益对应的最优注采参数。
进一步,所述步骤1)中,以开发期内经济净现值NPV来评价聚合物驱开发效益,得到最优控制函数:
第一约束条件为:ei(u)=0,i∈Ii
第二约束条件为:ej(u)≤0,j∈Ij
第三约束条件为:ulow≤u≤uup
其中,J(u)表示聚合物驱的开发效益,u为注采参数对应的控制向量;L为控制时间;NP为生产井数;NI为注水井数,NPI为注聚井数;分别为产油、产水、注水速度及注聚量;ro、rw、rwi、rp、b分别表示原油价格、产水成本、注水成本、聚合物价格和年利率,Δtn为n时刻模拟计算时间步,tn为n时刻累积计算时间;ei(u)表示等式约束方程;Ii表示等式约束个数;ej(u)表示不等式约束方程;Ij表示不等式约束个数;ulow表示控制变量下边界值;uup表示控制变量上边界值。
进一步,所述步骤2)中,采用协方差矩阵对所述初始扰动向量进行光滑处理和搜索方向筛选,生成目标扰动向量。
进一步,所述步骤3)中,对最优控制函数进行求解具体过程如下:3.1)生成优化参数以及迭代步l等于0对应的第一控制变量ul;3.2)采用对数变换方法将第一控制变量ul转换为第三约束条件对应的无约束优化条件下的第二控制变量sl;3.3)将第二控制变量sl带入到最优控制函数中,生成对应的第一函数值J(sl);3.4)利用包括目标扰动向量和优化参数的蒙特卡洛算法求取最优控制函数在第二控制变量sl处的蒙特卡洛梯度;3.5)根据投影梯度优化方法,采用优化参数和步骤3.4)计算出的蒙特卡洛梯度,将第二控制变量sl转换为满足第一约束条件和第二约束条件的第三控制变量sl+1;3.6)计算第三控制变量sl+1对应的第二函数值,并比较第二函数值和第一函数值,若第二函数值大于第一函数值,则将第三控制变量sl+1设定为目标控制变量;否则,将迭代步长减半,并返回至步骤3.5);3.7)判断第二函数值是否满足预设的收敛条件,若满足,则输出第二函数值为最优开发效益,输出目标控制变量为最优注采参数对应的控制向量;若不满足,则令迭代步加1,并返回至步骤3.2);或者判断迭代次数是否达到预设的最大迭代次数,若是,则输出第二函数值为最优开发效益,输出目标控制变量为最优注采参数对应的控制向量,若否,则令迭代步加1,并返回至步骤3.2)。
进一步,所述步骤3.2)中,采用第一公式将第一控制变量ul转换为第三约束条件对应的无约束优化条件下的第二控制变量sl,第一公式为:
其中,表示变换后的第i个控制变量;表示变量u在第lth迭代步下的第ith个控制分量,表示上界,表示下界。
进一步,所述步骤3.4)中,蒙特卡洛梯度求解过程为:3.4.1)在第l个迭代步,生成了Nr个用于计算近似梯度的目标函数值设ΔJj为第j个扰动向量sl ,j生成的目标函数值与当前最优控制变量sl所对应的目标函数值之差;3.4.2)在第l个迭代步,令第j个扰动向量的第i个元素为则蒙特卡洛梯度为:
进一步,所述步骤3.5)中,采用第二公式将第二控制变量sl转换为满足第一约束条件和第二约束条件的第三控制变量sl+1
进一步,所述第二公式包括无需归位移动时对应的第二公式:以及需要归位移动时对应的第二公式:
一种基于蒙特卡洛算法的聚合物驱生产优化系统,其特征在于:该系统包括模型建立模块、扰动向量生成模块和优化模块;所述模型建立模块用于根据油藏数值模拟器生成油藏优化的最优控制函数;所述扰动向量生成模块用于生成协方差矩阵和初始扰动向量,并采用协方差矩阵对初始扰动向量进行处理,生成目标扰动向量;所述优化模块用于采用包括目标扰动向量的蒙特卡洛算法对最优控制函数进行求解,生成最优开发效益对应的最优注采参数。
进一步,所述优化模块包括第一生成单元、第一变换单元、第一计算单元、第二计算单元、第二变换单元、比较单元和收敛判断单元;所述第一生成单元,用于生成优化参数以及迭代步l等于0对应的第一控制变量ul,优化参数包括迭代步长和扰动常数;所述第一变换单元,用于采用对数变换方法将第一控制变量ul转换为第一约束条件对应的无约束优化条件下的第二控制变量sl;所述第一计算单元,用于将第二控制变量sl带入到最优控制函数中,生成对应的第一函数值;所述第二计算单元,用于利用包括目标扰动向量和优化参数的蒙特卡洛算法求取最优控制函数在第二控制变量sl处的蒙特卡洛梯度;所述第二变换单元,用于根据投影梯度优化方法,采用优化参数和所述蒙特卡洛梯度,将第二控制变量sl转换为满足第二约束条件和第三约束条件的第三控制变量sl+1;所述比较单元,用于将第三控制变量sl+1带入到最优控制函数中,生成对应的第二函数值;并比较第二函数值和第一函数值,若第二函数值大于第一函数值,则将第三控制变量sl+1设定为目标控制变量,否则,将迭代步长减半,并驱动第二变换单元生成新的第三控制变量sl+1;所述收敛判断单元,用于判断第二函数值是否满足预设的收敛条件,若是,则输出第二函数值为最优开发效益,输出目标控制变量为最优注采参数对应的控制向量,若否,则令迭代步加1,并驱动第一变换单元;或者收敛单元用于判断迭代次数是否达到预设的最大迭代次数,若是,则输出第二函数值为最优开发效益,输出目标控制变量为最优注采参数对应的控制向量,若否,则令迭代步加1,并驱动第一变换单元开始下一次迭代过程。
本发明由于采取以上技术方案,其具有以下优点:1、本发明对现有的油藏开发方法进行了优化,不仅计算过程简单,容易实现,而且优化效率高,收敛速度快,易于和任意油藏数值模拟器相结合进行聚合物驱生产优化的计算,能够同时优化注采井生产参数和注入聚合物浓度,并且本发明的生产优化调控方案考虑了在时间上的相关性,更为平滑易于现场实施并能显著改善油藏开发效果,为科学合理的进行油田高效开发提供了重要的决策依据。2、由于不同的最优控制函数可以得到不同的优化结果,针对油田开发的实际情况,本发明以开发期内经济净现值(NPV)来评价聚合物驱开发效益,同时采用了新的最优控制函数,可以同时对注采量、注聚量和段塞大小进行优化,具有较高的优化效率且优化结果准确直观。3、本发明引入协方差矩阵对现有的蒙特卡洛算法进行改进,具有使梯度光滑、对梯度方向进行筛选的作用,进而使得近似梯度与真实梯度夹角更小、更接近于真实梯度,从而达到加快收敛速度、提高算法的收敛速度的目的。4、本发明将相应的约束优化问题转换为无约束优化问题,简化了优化过程,提高了优化速度和效率。同时,采用投影梯度优化方法,即采用投影移动和归位移动两个过程,使得到的第三控制变量满足预设的优化条件,提高了优化结果的准确性。本发明可以广泛在油气田开发领域中应用。
附图说明
图1是本发明基于蒙特卡洛算法的聚合物驱生产优化方法的流程示意图;
图2a是本发明实施例中采用现有技术聚合物驱生产优化方法的归一化梯度图;
图2b是本发明实施例中采用聚合物驱生产优化方法优化后的归一化梯度图;
图3是本发明实施例中采用聚合物驱生产优化方法优化后的渗透率场分布图;
图4是本发明实施例中采用本发明方法和采用现有技术聚合物驱生产优化方法的净现值优化结果对比图;
图5是本发明实施例中采用聚合物驱生产优化方法和采用现有技术聚合物驱生产优化方法的区块含水率优化结果对比图;
图6a是本发明实施例中优化前饱和度场图;
图6b是本发明实施例中采用现有技术聚合物驱生产优化方法的饱和度场图;
图6c是本发明实施例中采用聚合物驱生产优化方法优化后的饱和度场图;
图7是本发明实施例中采用聚合物驱生产优化方法和采用现有技术聚合物驱生产优化方法的聚合物浓度优化结果对比图;
图8是本发明实施例中采用聚合物驱生产优化方法和采用现有技术聚合物驱生产优化方法的聚合物段塞优化结果对比图;
图9是本发明实施例中采用聚合物驱生产优化方法和采用现有技术聚合物驱生产优化方法的注水量优化结果对比图;
图10a是本发明实施例中采用现有技术聚合物驱生产优化方法的产液量优化结果图;
图10b是本发明实施例中采用聚合物驱生产优化方法优化后的产液量优化结果图;
图11是本发明基于蒙特卡洛算法的聚合物驱生产优化系统的结构示意图。
具体实施方式
本发明的基本思想是首先把对油藏生产体系的控制描述成一个最优化问题,借助油藏数值模拟技术建立了油藏生产最优控制模型;在考虑油藏实际注采参数特点的基础上,通过引入控制参数协方差矩阵,并采用蒙特卡洛梯度逼近算法对所建控制模型进行求解,获得最优油藏注采参数及注入聚合物浓度。所得最优调控方案考虑了油水井注采参数间的相关性,便于实际操作并能显著改善油藏开发效果,为科学合理的进行油田高效开发提供了重要的决策依据。下面结合附图和实施例对本发明进行详细的描述。
如图1所示,本发明提供一种基于蒙特卡洛算法的聚合物驱生产优化方法,其包括以下步骤:
1)根据油藏数值模拟器生成油藏优化的最优控制函数;最优控制函数的输入为注采参数对应的控制向量,输出为聚合物驱的开发效益;
2)由球形模型生成协方差矩阵,利用蒙特卡洛逼近算法在最优控制变量附近生成初始扰动向量,并采用协方差矩阵对初始扰动向量进行处理,生成目标扰动向量;
3)采用包括目标扰动向量的蒙特卡洛算法对最优控制函数进行求解,生成最优开发效益对应的最优注采参数。
可选地,作为本发明的一个实施例,步骤1)中基于油藏数值模拟器来描述油藏开发生产系统,建立了油藏生产优化的最优控制模型,即最优控制函数,不同的最优控制模型会得到不同的优化结果。针对当前油田开发的实际情况,该实施例以开发期内经济净现值(NPV)来评价聚合物驱开发效益,以NPV为性能指标函数结合约束条件得出如下最优控制函数:
约束条件为:
第一约束条件为:ei(u)=0,i∈Ii (2)
第二约束条件为:ej(u)≤0,j∈Ij (3)
第三约束条件为:ulow≤u≤uup (4)
其中,J(u)表示聚合物驱的开发效益,u为注采参数对应的控制向量;L为控制时间;NP为生产井数;NI为注水井数,NPI为注聚井数;分别为产油、产水、注水速度及注聚量;ro、rw、rwi、rp、b分别表示原油价格、产水成本、注水成本、聚合物价格和年利率,Δtn为n时刻模拟计算时间步,tn为n时刻累积计算时间;ei(u)表示等式约束方程;Ii表示等式约束个数;ej(u)表示不等式约束方程;Ij表示不等式约束个数;ulow表示控制变量下边界值;uup表示控制变量上边界值。
油藏生产优化问题就是控制变量在满足约束条件下,求取性能指标的最大值及其对应的最优控制。对于最优控制模型的求解,由于J计算公式中不显含u,无法获得解析梯度,而伴随法计算又过于复杂,而蒙特卡洛梯度逼近算法(MCGA)是一种简单的有限差分梯度近似算法,每个迭代步最少仅需进行两次数值计算即可获得优化方向,因此,本发明实施例结合实际油藏注采参数特点,采用计算更为高效的蒙特卡洛梯度逼近算法对上述最优控制模型求解。
可选地,作为本发明的一个实施例,步骤2)中,采用协方差矩阵对所述初始扰动向量进行光滑处理和搜索方向筛选,生成目标扰动向量。
步骤2)中的蒙特卡洛梯度逼近算法的基本思想是:
首先,利用蒙特卡洛法在当前控制变量周围生成若干个随机变量的实现;
然后,分别求取各随即变量实现所对应的目标函数值;
最后,利用各实现及其目标函数值来估计目标函数的梯度。
在利用蒙特卡洛逼近算法时,首先会在第l迭代步的控制变量附近生成扰动向量,现有技术的蒙特卡洛逼近算法通常直接采用这个扰动向量进行求解,而本实施例则采用协方差矩阵对所述扰动向量进行光滑处理和搜索方向筛选后,采用生成的目标扰动向量进行求解,从而使近似梯度与真实梯度夹角更小、更接近于真实梯度,达到加快收敛速度、提高算法的收敛速度的目的。
具体的,协方差矩阵由球形模型进行定义,表达式为:
其中,σ是扰动向量的标准方差,表示时间关联程度;m和n分别表示时间和步数;T表示时间步数预设阈值。协方差矩阵是一个方块对角矩阵,其特殊的形式是由于考虑了单井的时间关联性,其表达式为:
接下来对协方差矩阵进行乔里斯基分解,分解后满足等式:
式中,为下三角矩阵,表示的转置;进而定义目标扰动向量Δ为:
其中,Z为服从标准正态分布的Nu维初始扰动向量。
可选地,作为本发明的一个实施例,步骤3)中对最优控制函数进行求解具体过程如下:
3.1)生成优化参数以及迭代步l等于0对应的第一控制变量ul;其中,优化参数包括迭代步长和扰动常数。
3.2)采用对数变换方法将第一控制变量ul转换为第三约束条件对应的无约束优化条件下的第二控制变量sl
其中,具体转换过程为:
主要是对于边界约束公式(4),采用对数变换方法,将边界约束转换为无约束优化问题,具体公式如下。
以上公式即为第一公式,其中,表示变换后的第i个控制变量;;表示变量u在第lth迭代步下的第ith个控制分量,表示上界,表示下界。利用对数变换方法,值域从正无穷变为负无穷,相应的约束优化问题被转换成了无约束优化问题,接下来的优化以为对象,且变换回公式如下:
3.3)将第二控制变量sl带入到最优控制函数中,生成对应的第一函数值J(sl);
3.4)利用包括目标扰动向量和优化参数的蒙特卡洛算法求取最优控制函数在第二控制变量sl处的蒙特卡洛梯度,具体的求解过程为:
3.4.1)在蒙特卡洛梯度逼近算法的计算过程中,在第l个迭代步,生成了Nr个目标函数值,即用于计算近似梯度。设ΔJj为第j个扰动向量sl,j生成的目标函数值与当前最优控制变量sl所对应的目标函数值之差,则ΔJj的表达式为:
ΔJj=J(sl,j)-J(sl)=J(sl,j+ε·Δl,j)-J(sl) j=1,2...,Nr (11)
3.4.2)在第l个迭代步,令第j个扰动向量的第i(i=1,2,...,Nu)个元素为则蒙特卡洛梯度为:
3.5)根据投影梯度优化方法,采用优化参数和步骤3.4)计算出的蒙特卡洛梯度,将第二控制变量sl转换为满足第一约束条件和第二约束条件的第三控制变量sl+1。该步骤具体的计算过程如下:
针对油藏开发实际情况而言,现场的注水量、采出量变化不宜过大,因此对控制变量进一步采用了等式和不等式约束的约束条件进行约束,即公式(2)~(3)。对于等式约束和不等式约束,采用了投影梯度优化方法进行处理。投影梯度优化方法是:当迭代点在可行域的内部时,以该点的负梯度方向为下降可行方向;而当迭代点位于可行域的边界上且其梯度方向指向可行域外部时,则取它的负梯度方向在边界上的投影为下降可行方向,若这个投影为零向量,则停止迭代,得到问题的极小点。
投影梯度优化方法的优化过程包括投影移动和归位移动两个过程。设在第l个迭代步,对于sl有Ns个激活的约束条件,其所组成的激活约束向量为ea,设其元素ea,i(i=1,2,…,Ns)为优化控制模型中的第s个约束条件,即
ea,i=es(ul)≥0,s∈Ii∪Ij (13)
其中,ul为sl对应的真实控制变量。那么,对于sl进行更新的投影移动计算公式为:
式中,α为迭代步长;P为Nu维投影矩阵:
P=I-▽ea((▽ea)T▽ea)-1(▽ea)T (15)
式中,I为Nu维单位阵;▽ea为由激活约束对sl的梯度所组成矩阵,其第i列向量表达式为:
对于非线性约束问题,进行投影移动后变量sl+1通常会越过边界,因此需要进行归位移动计算。但是在sl处进行归位计算的方法在实际应用中进行一次归位计算可能难以满足全部约束,为此这里采用迭代归位移动的求解方法。
设在sl+1处有Ns'个激活的约束,其组成的激活约束向量为e′a,将向量e′a(sl+1+β▽e′a)在sl+1处进行线性展开可得,
e′a(sl+1+β▽e′a)=e′a(sl+1)+(▽e′a)T▽e′aβ (17)
令e′a(sl+1)+(▽e′a)T▽e′aβ=0,求取β为
β=((▽e′a)T▽e′a)-1e′a (18)
从而sl+1进行归位可得,
其中,▽e′a为e′a对sl+1的梯度所组成的矩阵,β为Ns'维向量。根据公式(19),对sl+1进行迭代归位移动计算的表达式为:
式中,n是归位移动中迭代的步数。以上公式(14)即为无需归位移动时对应的第二公式,公式(20)为需要归位移动时对应的第二公式。
3.6)计算第三控制变量sl+1对应的第二函数值,并比较第二函数值和第一函数值,若第二函数值大于第一函数值,则将第三控制变量sl+1设定为目标控制变量,然后执行步骤3.7);否则,将迭代步长减半,并返回至步骤3.5)。
3.7)判断第二函数值是否满足预设的收敛条件,若满足,则输出第二函数值为最优开发效益,输出目标控制变量为最优注采参数对应的控制向量;若不满足,则令迭代步加1,并返回至步骤3.2);或者判断迭代次数是否达到预设的最大迭代次数,若是,则输出第二函数值为最优开发效益,输出目标控制变量为最优注采参数对应的控制向量,若否,则令迭代步加1,并返回至步骤3.2)。其中,收敛条件为:
公式(21)和(22)中,ξ为预设极小常数。
实施例:
本实施例的模型区块是一个典型的含高渗带地层,网格数为21×21×1,网格长度均为50m,厚度为20m。平均孔隙度和渗透率分别为949mD和0.25,初始含油饱和度为0.67,初始油藏压力为10.426MPa。如图2a所示,为采用现有技术聚合物驱生产优化方法,即采用标准蒙特卡洛逼近方法优化后的归一化梯度图,而图2b为本发明采用改进后的蒙特卡洛逼近方法优化后的归一化梯度图,结合图2a和图2b可以看到,采用本发明的优化方法得到的梯度更接近真实梯度,因此优化效率更高,收敛速度更快。如图3所示,展示了该模型区块的渗透率场,可见该模型为一注四采5点法井网,其中分别有两口井分布在高渗带和低渗区。原生产方案在高产井含水率约91%时开始注入聚合物,注聚时间为900天,总模拟时间为4200天。对原生产方案进行生产优化时,其中设置的边界约束条件为:高产井为0到200m3/天,低产井为0到70m3/天,注水井为0到600m3/天,聚合物浓度为0到4kg/m3
图4为NPV的优化结果,横坐标表示迭代次数。如图4所示,在约20次迭代后,本发明改进后的蒙特卡洛梯度逼近算法优化结果收敛,而现有技术的标准蒙特卡洛梯度逼近算法的NPV仍然在缓慢增加,需要约30次迭代,且其优化后的值始终低于本发明算法。经过优化,NPV值在原始生产制度的基础上增加了约21%。由此可见本发明的蒙特卡洛梯度逼近算法具有较高的计算效率。
图5为优化前后的区块含水率对比结果,由图5可见,本发明改进后的蒙特卡洛梯度逼近算法的区块含水率下降至77%,比优化前下降了约5%;而现有技术的标准蒙特卡洛梯度逼近算法的区块含水率只下降到了80%,由此可见本发明的算法具有更好的优化表现。
图6a~6c为优化前后的饱和度场对比结果图,可以看到经过本发明蒙特卡洛梯度逼近算法的优化,在注水井周围的网格的注入水流动到了更广的区域,同时从注水井到低产井方向上也有更多的网格含油饱和度降低,高产井周边的网格含油饱和度也进一步下降。与现有技术的标准蒙特卡洛梯度逼近算法相比,现有技术在注水井附近网格含水饱和度更高,说明本发明算法的波及效果更好。因此说明了本发明的蒙特卡洛梯度逼近算法具有较好的优化效果。
图7为注聚浓度的优化结果,图8为段塞大小的优化结果,从图7和图8中可以看到,本发明改进后蒙特卡洛梯度逼近算法的优化结果比现有算法的优化结果具有更强的波动性,也得到了更好的优化结果。
图9为注水量的优化结果,图10a、图10b为产液量的优化结果,如图9、图10a、图10b所示,与注聚浓度和段塞大小的优化结果类似,本发明蒙特卡洛梯度逼近算法的优化结果具有更强的规律性,而现有算法,特别是从图10b可以看出,其值均在某个值附近随机增减,导致最终的优化效果不佳,由此也说明了本发明的蒙特卡洛梯度逼近算法具有更好的优化效果。
如图11所示,本发明还提供一种基于蒙特卡洛算法的聚合物驱生产优化系统,其包括模型建立模块、扰动向量生成模块和优化模块。其中:
模型建立模块用于根据油藏数值模拟器生成油藏优化的最优控制函数,最优控制函数的输入为注采参数对应的控制向量,输出为聚合物驱的开发效益;
扰动向量生成模块用于生成协方差矩阵和初始扰动向量,并采用协方差矩阵对初始扰动向量进行处理,生成目标扰动向量;
优化模块用于采用包括目标扰动向量的蒙特卡洛算法对最优控制函数进行求解,生成最优开发效益对应的最优注采参数。
可选的,在本发明的一个优选实施例中,采用开发期内经济净现值NPV表示聚合物驱的开发效益,最优控制函数具体为:
第一约束条件为:ei(u)=0,i∈Ii
第二约束条件为:ej(u)≤0,j∈Ij
第三约束条件为:ulow≤u≤uup
其中,J(u)表示聚合物驱的开发效益,u为注采参数对应的控制向量;L为控制时间;NP为生产井数;NI为注水井数,NPI为注聚井数;分别为产油、产水、注水速度及注聚量;ro、rw、rwi、rp、b分别表示原油价格、产水成本、注水成本、聚合物价格和年利率,Δtn为n时刻模拟计算时间步,tn为n时刻累积计算时间。
在另一优选实施例中,扰动向量生成模块中采用协方差矩阵对初始扰动向量进行光滑处理和搜索方向筛选,生成目标扰动向量。
可选的,在本发明的一个优选实施例中,优化模块包括第一生成单元、第一变换单元、第一计算单元、第二计算单元、第二变换单元、比较单元和收敛判断单元。其中:
第一生成单元,用于生成优化参数以及迭代步l等于0对应的第一控制变量ul,优化参数包括迭代步长和扰动常数;
第一变换单元,用于采用对数变换方法将第一控制变量ul转换为第一约束条件对应的无约束优化条件下的第二控制变量sl
第一计算单元,用于将第二控制变量sl带入到最优控制函数中,生成对应的第一函数值;
第二计算单元,用于利用包括目标扰动向量和优化参数的蒙特卡洛算法求取最优控制函数在第二控制变量sl处的蒙特卡洛梯度;
第二变换单元,用于根据投影梯度优化方法,采用优化参数和所述蒙特卡洛梯度,将第二控制变量sl转换为满足第二约束条件和第三约束条件的第三控制变量sl+1
比较单元,用于将第三控制变量sl+1带入到最优控制函数中,生成对应的第二函数值;并比较第二函数值和第一函数值,若第二函数值大于第一函数值,则将第三控制变量sl +1设定为目标控制变量,否则,将迭代步长减半,并驱动第二变换单元生成新的第三控制变量sl+1
收敛判断单元,用于判断第二函数值是否满足预设的收敛条件,若是,则输出第二函数值为最优开发效益,输出目标控制变量为最优注采参数对应的控制向量,若否,则令迭代步加1,并驱动第一变换单元;或者收敛单元用于判断迭代次数是否达到预设的最大迭代次数,若是,则输出第二函数值为最优开发效益,输出目标控制变量为最优注采参数对应的控制向量,若否,则令迭代步加1,并驱动第一变换单元开始下一次迭代过程。
在本实施例中,第一变换单元具体用于利用第一预设公式将第一控制变量ul转换为所述第一约束条件对应的无约束优化条件下的第二控制变量sl
第二变换单元具体用于利用第二预设公式将第二控制变量sl转换为满足第二约束条件和第三约束条件的第三控制变量sl+1
综上所述,本发明建立了一种考虑注采量相关性的基于蒙特卡洛梯度逼近算法的油藏生产优化方法。现阶段,各大油田均进入开发的中后期,含水率不断上升,单纯的使用水驱已经不能进一步提高采收率;同时作为三次采油的一种重要方法,聚合物驱在制定方案时多利用正交设计等方法,方案往往不是最优。本发明针对以上的一些问题,提出了一套新的可应用于聚合物驱的生产优化方法。本发明首先基于聚合物驱方法建立了新的优化模型,同时引入了协方差矩阵对蒙特卡洛梯度逼近算法进行了改进。蒙特卡洛梯度逼近算法是一种无梯度类算法,具有算法步骤简单,容易实现,同时具有易于与商业油藏模拟器结合的特点。最后在实施例中,使用了一含高渗条带的一注四采模型对本方法进行了验证,其结果显示改进后的蒙特卡洛梯度逼近算法比其标准算法效率更高,因而证明了本发明的方法具有更高的计算效率与实用价值。由于水驱可以视为聚合物驱的一种特殊情况,即注入聚合物浓度为0的情况,因此本发明依然可用于解决水驱优化问题,具有较广泛的适用性,对本发明进行适当改进,可以进一步推广到其他的三次采油方法中。
在本发明的描述中,需要理解的是,术语“第一”、“第二”仅用于描述目的,而不能理解为指示或暗示相对重要性或者隐含指明所指示的技术特征的数量。由此,限定有“第一”、“第二”的特征可以明示或者隐含地包括至少一个该特征。在本发明的描述中,“多个”的含义是至少两个,例如两个,三个等,除非另有明确具体的限定。
在本说明书的描述中,参考术语“一个实施例”、“一些实施例”、“示例”、“具体示例”、或“一些示例”等的描述意指结合该实施例或示例描述的具体特征、结构、材料或者特点包含于本发明的至少一个实施例或示例中。在本说明书中,对上述术语的示意性表述不必须针对的是相同的实施例或示例。而且,描述的具体特征、结构、材料或者特点可以在任一个或多个实施例或示例中以合适的方式结合。
此外,在不相互矛盾的情况下,本领域的技术人员可以将本说明书中描述的不同实施例或示例以及不同实施例或示例的特征进行结合和组合。凡根据本发明原理进行的改进和等同变换,均不应排除在本发明的保护范围之外。

Claims (7)

1.一种基于蒙特卡洛算法的聚合物驱生产优化方法,其特征在于包括以下步骤:
1)根据油藏数值模拟器生成油藏优化的最优控制函数;
所述步骤1)中,以开发期内经济净现值NPV来评价聚合物驱开发效益,得到最优控制函数:
第一约束条件为:ei(u)=0,i∈Ii
第二约束条件为:ej(u)≤0,j∈Ij
第三约束条件为:ulow≤u≤uup
其中,J(u)表示聚合物驱的开发效益,u为注采参数对应的控制向量;L为控制时间;NP为生产井数;NI为注水井数,NPI为注聚井数;分别为产油、产水、注水速度及注聚量;ro、rw、rwi、rp、b分别表示原油价格、产水成本、注水成本、聚合物价格和年利率,Δtn为n时刻模拟计算时间步,tn为n时刻累积计算时间;ei(u)表示等式约束方程;Ii表示等式约束个数;ej(u)表示不等式约束方程;Ij表示不等式约束个数;ulow表示控制变量下边界值;uup表示控制变量上边界值;
2)由球形模型生成协方差矩阵,利用蒙特卡洛逼近算法在最优控制变量附近生成初始扰动向量,并采用协方差矩阵对初始扰动向量进行处理,生成目标扰动向量;
3)采用包括目标扰动向量的蒙特卡洛算法对最优控制函数进行求解,生成最优开发效益对应的最优注采参数;
所述步骤3)中,对最优控制函数进行求解具体过程如下:
3.1)生成优化参数以及迭代步l等于0对应的第一控制变量ul
3.2)采用对数变换方法将第一控制变量ul转换为第三约束条件对应的无约束优化条件下的第二控制变量sl
3.3)将第二控制变量sl带入到最优控制函数中,生成对应的第一函数值J(sl);
3.4)利用包括目标扰动向量和优化参数的蒙特卡洛算法求取最优控制函数在第二控制变量sl处的蒙特卡洛梯度;
3.5)根据投影梯度优化方法,采用优化参数和步骤3.4)计算出的蒙特卡洛梯度,将第二控制变量sl转换为满足第一约束条件和第二约束条件的第三控制变量sl+1
3.6)计算第三控制变量sl+1对应的第二函数值,并比较第二函数值和第一函数值,若第二函数值大于第一函数值,则将第三控制变量sl+1设定为目标控制变量;否则,将迭代步长减半,并返回至步骤3.5);
3.7)判断第二函数值是否满足预设的收敛条件,若满足,则输出第二函数值为最优开发效益,输出目标控制变量为最优注采参数对应的控制向量;若不满足,则令迭代步加1,并返回至步骤3.2);或者判断迭代次数是否达到预设的最大迭代次数,若是,则输出第二函数值为最优开发效益,输出目标控制变量为最优注采参数对应的控制向量,若否,则令迭代步加1,并返回至步骤3.2)。
2.如权利要求1所述的一种基于蒙特卡洛算法的聚合物驱生产优化方法,其特征在于:所述步骤2)中,采用协方差矩阵对所述初始扰动向量进行光滑处理和搜索方向筛选,生成目标扰动向量。
3.如权利要求1所述的一种基于蒙特卡洛算法的聚合物驱生产优化方法,其特征在于:所述步骤3.2)中,采用第一公式将第一控制变量ul转换为第三约束条件对应的无约束优化条件下的第二控制变量sl,第一公式为:
其中,表示变换后的第i个控制变量;表示变量u在第lth迭代步下的第ith个控制分量,表示上界,表示下界。
4.如权利要求1所述的一种基于蒙特卡洛算法的聚合物驱生产优化方法,其特征在于:所述步骤3.4)中,蒙特卡洛梯度求解过程为:
3.4.1)在第l个迭代步,生成了Nr个用于计算近似梯度的目标函数值设△Jj为第j个扰动向量sl,j生成的目标函数值与当前最优控制变量sl所对应的目标函数值之差;
3.4.2)在第l个迭代步,令第j个扰动向量的第i个元素为i=1,2,...,Nu,则蒙特卡洛梯度为:
5.如权利要求1所述的一种基于蒙特卡洛算法的聚合物驱生产优化方法,其特征在于:所述步骤3.5)中,采用第二公式将第二控制变量sl转换为满足第一约束条件和第二约束条件的第三控制变量sl+1
6.如权利要求5所述的一种基于蒙特卡洛算法的聚合物驱生产优化方法,其特征在于:所述第二公式包括无需归位移动时对应的第二公式:
以及需要归位移动时对应的第二公式:
7.一种基于蒙特卡洛算法的聚合物驱生产优化系统,其特征在于:该系统包括模型建立模块、扰动向量生成模块和优化模块;
所述模型建立模块用于根据油藏数值模拟器生成油藏优化的最优控制函数;
所述扰动向量生成模块用于生成协方差矩阵和初始扰动向量,并采用协方差矩阵对初始扰动向量进行处理,生成目标扰动向量;
所述优化模块用于采用包括目标扰动向量的蒙特卡洛算法对最优控制函数进行求解,生成最优开发效益对应的最优注采参数;所述优化模块包括第一生成单元、第一变换单元、第一计算单元、第二计算单元、第二变换单元、比较单元和收敛判断单元;
所述第一生成单元,用于生成优化参数以及迭代步l等于0对应的第一控制变量ul,优化参数包括迭代步长和扰动常数;
所述第一变换单元,用于采用对数变换方法将第一控制变量ul转换为第一约束条件对应的无约束优化条件下的第二控制变量sl
所述第一计算单元,用于将第二控制变量sl带入到最优控制函数中,生成对应的第一函数值;
所述第二计算单元,用于利用包括目标扰动向量和优化参数的蒙特卡洛算法求取最优控制函数在第二控制变量sl处的蒙特卡洛梯度;
所述第二变换单元,用于根据投影梯度优化方法,采用优化参数和所述蒙特卡洛梯度,将第二控制变量sl转换为满足第二约束条件和第三约束条件的第三控制变量sl+1
所述比较单元,用于将第三控制变量sl+1带入到最优控制函数中,生成对应的第二函数值;并比较第二函数值和第一函数值,若第二函数值大于第一函数值,则将第三控制变量sl +1设定为目标控制变量,否则,将迭代步长减半,并驱动第二变换单元生成新的第三控制变量sl+1
所述收敛判断单元,用于判断第二函数值是否满足预设的收敛条件,若是,则输出第二函数值为最优开发效益,输出目标控制变量为最优注采参数对应的控制向量,若否,则令迭代步加1,并驱动第一变换单元;或者收敛单元用于判断迭代次数是否达到预设的最大迭代次数,若是,则输出第二函数值为最优开发效益,输出目标控制变量为最优注采参数对应的控制向量,若否,则令迭代步加1,并驱动第一变换单元开始下一次迭代过程。
CN201710056212.5A 2017-01-25 2017-01-25 一种基于蒙特卡洛算法的聚合物驱生产优化方法和系统 Active CN107218019B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201710056212.5A CN107218019B (zh) 2017-01-25 2017-01-25 一种基于蒙特卡洛算法的聚合物驱生产优化方法和系统

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201710056212.5A CN107218019B (zh) 2017-01-25 2017-01-25 一种基于蒙特卡洛算法的聚合物驱生产优化方法和系统

Publications (2)

Publication Number Publication Date
CN107218019A CN107218019A (zh) 2017-09-29
CN107218019B true CN107218019B (zh) 2019-05-21

Family

ID=59928156

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201710056212.5A Active CN107218019B (zh) 2017-01-25 2017-01-25 一种基于蒙特卡洛算法的聚合物驱生产优化方法和系统

Country Status (1)

Country Link
CN (1) CN107218019B (zh)

Families Citing this family (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN107939372B (zh) * 2017-10-23 2020-03-06 南京特雷西能源科技有限公司 针对小断块油藏的最优井位部署方法和装置
CN109145251B (zh) * 2018-08-22 2023-03-24 合肥工业大学 一种改进型同步扰动随机逼近算法的大气参数求解方法
CN109899043B (zh) * 2019-03-20 2021-06-22 中国海洋石油集团有限公司 一种周期注聚提高原油采收率幅度的定量预测方法
CN110941890B (zh) * 2019-09-27 2022-11-01 中国海洋石油集团有限公司 基于最优控制理论的海上油藏动态实时生产优化方法
CN112984919A (zh) * 2021-05-12 2021-06-18 创新奇智(北京)科技有限公司 制冷系统能效优化方法、装置、电子设备及存储介质

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102289559A (zh) * 2011-05-30 2011-12-21 复旦大学 蒙特卡洛模拟预测自由基共聚体系中共聚物序列分布的方法
CN104216341A (zh) * 2013-05-31 2014-12-17 中国石油化工股份有限公司 一种基于改进随机扰动近似算法的油藏生产实时优化方法
EP2892974A1 (en) * 2012-09-03 2015-07-15 Poweltec Use of thermo-thickening polymers in the gas- and oilfield industry
CN105019894A (zh) * 2015-07-29 2015-11-04 长江大学 一种多层油藏井间连通性模型建立方法及系统

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102289559A (zh) * 2011-05-30 2011-12-21 复旦大学 蒙特卡洛模拟预测自由基共聚体系中共聚物序列分布的方法
EP2892974A1 (en) * 2012-09-03 2015-07-15 Poweltec Use of thermo-thickening polymers in the gas- and oilfield industry
CN104216341A (zh) * 2013-05-31 2014-12-17 中国石油化工股份有限公司 一种基于改进随机扰动近似算法的油藏生产实时优化方法
CN105019894A (zh) * 2015-07-29 2015-11-04 长江大学 一种多层油藏井间连通性模型建立方法及系统

Also Published As

Publication number Publication date
CN107218019A (zh) 2017-09-29

Similar Documents

Publication Publication Date Title
CN107218019B (zh) 一种基于蒙特卡洛算法的聚合物驱生产优化方法和系统
CN108868712B (zh) 一种基于连通性方法的油藏开发生产优化方法和系统
CN105740563B (zh) 一种成熟油田二次开发之优势通道识别方法
CN104216341A (zh) 一种基于改进随机扰动近似算法的油藏生产实时优化方法
CA2514516A1 (en) Performance prediction method for hydrocarbon recovery processes
CN106351624B (zh) 特高含水期断块油藏分区调控提高采收率方法
CN109948272A (zh) 基于井间连通性的调堵动态预测方法和系统
CN109426672B (zh) 基于不确定地质模型的油藏注采参数优化方法
CN106779211A (zh) 一种用于油田注采井网驱替开发的产能预测方法
Sarma Efficient closed-loop optimal control of petroleum reservoirs under uncertainty
CN105868508A (zh) 一种基于气测录井信息的产能定量预测方法
CN107038516A (zh) 一种中渗复杂断块油藏水驱开发效果定量评价方法
CN104060985A (zh) 一种层状油藏调剖堵水堵剂进入深度测试方法及系统
Zhang et al. Water flooding optimization with adjoint model under control constraints
CN109958413A (zh) 一种特高含水期油藏动态流动单元划分方法
Chen et al. Water flooding performance prediction in layered reservoir using big data and artificial intelligence algorithms
CN111350498B (zh) 一种中高渗油藏特高含水开发期弱驱分布特征描述的方法
CN109339777A (zh) 基于改进qfd的低渗透老油田开发经济潜力评估方法
CN109829586A (zh) 一种低渗透油藏聚合物微球深部调驱技术增油量评价方法
De Montleau et al. Production optimization under constraints using adjoint gradients
CN107491568A (zh) 基于实数编码遗传算法的复杂结构井优化设计方法
US11501043B2 (en) Graph network fluid flow modeling
CN111625925B (zh) 一种基于色谱分离的三元复合驱注采优化方法
CN108491625A (zh) 一种三元复合驱油体系提高采收率的预测方法
CN111091476A (zh) 针对强非均质地层的油藏数值模拟方法

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
CB02 Change of applicant information
CB02 Change of applicant information

Address after: 100010 Beijing, Chaoyangmen, North Street, No. 25, No.

Applicant after: China Offshore Oil Group Co., Ltd.

Applicant after: CNOOC research institute limited liability company

Address before: 100010 Beijing, Chaoyangmen, North Street, No. 25, No.

Applicant before: China National Offshore Oil Corporation

Applicant before: CNOOC Research Institute

GR01 Patent grant
GR01 Patent grant