基于最优控制理论的聚合物驱方案优化设计方法 一、 技术领域
本技术涉及工业优化控制技术， 尤其是一种基于最优控制理论的聚合物驱方案优 化设计方法。 二、 背景技术
石油作为一种不可再生能源， 其地下蕴藏的储量是有限的。随着经济的发展， 对 石油的需求与日俱增。聚合物驱 (Polymer Flooding) 是指在注入的水中加入水溶性高分 子量聚合物， 通过增加水相粘度和降低水相渗透率， 改善流度比， 进而提高原油采收率的方 法。 聚合物驱只是在原来水驱的基础上添加了聚合物， 因此它又被称为改性水驱 (Modified Water Flooding)。聚合物驱始于 20 世纪 50 年代末和 60 年代初。在美国、 加拿大、 德国及 中国， 聚合物在提高采收率方法中占有重要的地位。
尽管大多数聚合物驱项目都取得了成功， 但是也有一些没有取得预期的效果。由 于聚合物的价格昂贵， 聚合物驱技术的应用存在风险高、 投资大的特点， 只有在经济预测表 明其投资可以获得丰厚的回报时， 才有可能进行商业应用。 因此， 科学地设计聚合物驱开发 方案以及根据矿场实际情况及时调整聚合物驱注采过程就显得尤为迫切和重要。
Colorado 大学的 Ramirez 等人最早利用最优控制理论研究了表面活性剂驱提高 采收率的最优注采策略。最优控制的性能指标为表面活性剂的注入费用减去产出原油的 价值， 控制变量为表面活性剂的注入浓度， 状态方程为一维油水两相的 Buckley-Leverett 方程和表面活性剂浓度的对流扩散和吸附方程。在求解时先用一维分布参数系统的极大 值原理推导了最优控制的必要条件， 在推导伴随方程时使用了拟线性假设， 最后采用梯度 法求解。Ramirez 的工作为推动最优控制理论在提高采收率技术中的应用开了一个先例， 随后这种方法被推广到二氧化碳驱、 氮气驱、 表面活性剂—聚合物驱等提高采收率方法中。 Ramirez 在研究二维模型的二氧化碳驱时采用了离散极大值原理， 把二维分布参数系统的 最优控制问题转化为离散最优控制问题， 推导了离散的伴随方程和性能指标的梯度， 最后 采用离散的梯度法进行求解。
在国内， 北京石油勘探开发科学研究院的齐与峰、 叶继根等人应用最优控制理论 研究了注气提高采收率的最优决策。齐与峰建立了蒸汽注入过程的最优控制模型， 支配方 程为描述流体和热能运动的偏微分方程， 控制变量为随时间变化的注采量和蒸汽的干度。 在推导该分布参数系统最优控制的必要条件时， 直接对最优控制模型进行变分， 并令增广 性能指标的一阶变分为零， 获得在伴随方程及梯度的表达式。叶继根在认真吸收多组分相 平衡求解技术和多组分油气渗流模型求解技术的基础上， 将最优方法应用于凝析气藏循环 注气的开发决策， 在求解时采用了离散极大值原理， 并设计出了优化决策软件， 在实际注气 采气速度和 实验中取得了较好效果。 其研究结果表明， 通过对气藏开发过程中的注气速度、 注入剂组成等进行优化， 不仅能有效地提高凝析油采收率， 而且还明显地改善了气田开发 的总体经济效益。
西南石油学院的刘昌贵等人根据系统控制论的观点和方法， 提出了油藏最优控制问题， 建立和求解了可用于注气提高采收率技术的油藏生产制度最优控制数学模型， 并完 成了方法的有效性检验和有关现场应用的实现方法。 他们在建立油藏生产制度最优控制数 学模型时， 以油藏开采系统为研究对象， 考虑了地层流体在多孔介质中的微分渗流规律。 在 求解模型时， 采用了变分和非线性规划的系统控制论方法， 并充分应用了油藏数值模拟软 件的优秀成果， 最终将该方法应用于现场二氧化碳吞吐注采参数的优化设计过程。
聚合物驱提高采收率问题中采用最优控制技术具有重要的理论与实际意义。 由于 聚合物驱模型为强非线性的三维分布参数模型， 模型求解十分复杂， 基于庞特里亚金极大 值原理推导伴随方程， 获得优化问题的梯度， 与差分方法和灵敏度法相比较， 具有对优化变 量个数不敏感的优点。子空间截断牛顿法具有局部二阶的收敛速度， 因此适于处理聚合物 驱这类复杂模型约束的最优控制问题。 目前国内关于聚合物驱提高采收率的最优控制方法 的专利尚无。 三、 发明内容
为解决聚合物驱的最优注入问题， 本发明提出了基于最优控制理论的聚合物驱方 案优化设计方法。 本发明解决其技术问题所采用的技术方案是 ：
基于最优控制理论的聚合物驱方案优化设计方法研究， 所述的聚合物驱最优控制 方法包括以下步骤 ：
1、 聚合物驱模型基本方程
设油藏所在的三维区域为 Ω ∈ R3， (x， y， z) ∈ Ω， Γ 为区域的边界，表示边界 的法线方向， 则聚合物驱的数学模型如下所示 ：
基本方程 ：
初始和边界条件 ：其中方程 (1)(2)(3) 分别为油相方程、 水相方程和聚合物的扩散吸附方程， 聚合 物驱模型中的辅助方程及辅助变量定义如下 ：
饱和度方程 ： So+Sw ＝ 1.
势能的定义 ： Φo ＝ p-ρogh， Φw ＝ p-pcow-ρwgh.
4流动系数的定义 ： 上述模型中各符号的物理意义如下 ：102352743 A CN 102352750 K 为绝对渗透率 ； Φw 为油相和水相的势位 ； Φo， p 为油相压力 ； Sw 分别为油相和水相的饱和度 ； So， pcow 为毛管力， Sw 是水相饱和度 ； g 为重力加速度 ； h 为深度， 向下为正 ； cp 为油藏水相中聚合物的浓度 ； ρw ， ρR 分别油、 水和岩石的密度 ； ρo， Bw 分别油和水的体积系数 ； Bo，说明书3/13 页φ 为岩石的孔隙度和聚合物可及孔隙度， φp 为聚合物可及孔隙度， fa 为可及孔隙系数 ； φp ＝ φfa， Pk 为相对渗透率下降系数 ； Crp 为单位质量岩石吸附聚合物的质量 ； ro， rw 分别为油相和水相的流动系数 ； rc 为聚合物方程对流项的流动系数 ； rd 为聚合物方程扩散项的流动系数 ； qo 为油相在标准状态下的流速， 流入为正 ； qw 为水相在标准状态下的流速， 流入为下 ； cpw 为井筒中聚合物溶液的浓度 ； D 为聚合物的扩散系数 ； x， y， z 为直角坐标系的三个方向， t 为时间 ；krw 分别为油相和水相的相对渗透率 ； kro， μp 为油相和含聚合物的水相粘度 ； μo，
2、 聚合物驱模型求解
聚合物驱模型中含有油水两相的渗流方程和聚合物的浓度方程， 是一个弱耦合的 非线性偏微分方程组。现在通常采用的方法是将模型的求解分为两部分 ： 一是油水两相方 程的求解 ； 二是聚合物对流扩散吸附方程的求解。 在求解油水两相渗流方程时， 可以把油水 两相方程化整理成压力方程和水相饱和度方程的形式， 然后采用隐压显饱 (IMPES) 的方法 求解， 即隐式求解压力方程， 显式求解饱和度方程。在每个时间步求解油水方程之后， 再利 用隐式或显式的方法求解聚合物的浓度方法。
1)、 全隐式差分格式
利用全隐式方差分方法将聚合物驱模型式 (1)(2)(3) 离散化， 设 τn 表示离散的 时间步长， n ＝ 0， 1， ...， N-1， x， y， z 三个方向的网格数分别为 nx， ny， nz， 将所有网格按一 定顺序排列， 记网格序号为 i(i ＝ 1， 2， ...， nb， nb ＝ nx×ny×nz)， 记第 i 个网格的坐标序号 T 为 (ix， iy， iz)。令 ui ＝ (pi， Swi， cpi) ， pi， Swi， cpi 分别表示网格 i 上的油藏压力、 水饱和度 和水相中的聚合物浓度， 则令 表示离散后所有网 格的求解变量， 则聚合物驱模型的全隐式差分格式可以表示成如下的一般形式
gn ＝ Fn+1(un+1)+Wn+1(un+1， vn)-[An+1(un+1)-An(un)] ＝ 0. (4)
其中 v 表示控制变量， 为聚合物的注入浓度等。上式是一个 nu 维的非线性代数方 程组， 其中 Fn(un+1) 表示流动项， Wn(un+1， vn) 表示井的流量项， An+1(un+1) 和 An(un) 分别表示n+1 和 n 时刻的累积项， 这两项只有求解变量取值的时刻不同。 设离散后的聚合物驱模型在 第 i 个网格上的节点方程为
其中式 (5)， (6) 和 (7) 分别对应连续模型中的式 (1)(2)(3) 表示离散的油相方 程、 水相方程和聚合物的对流扩散吸附方程。以下给出上述表达式中各项的具体定义。设 x， y， z 三个方向上与网格 i 相邻的网格数为 nsi， 记所有与 i 相邻网格的序号集合为 ηi， 令 k 表示与 i 相邻的网格序号， 则流动项的差分格式为
其中 SKi， SDi， Ki， k ＝ Ki， kAri， k/di， k， k ＝ D·Ari， k/di， k 为空间系数， k 为 i 与 k 节点之 分别表示 i 与 k 节点之间 处的系数， 按照上游权原则取值上述表达 设 Vi 表示第 i 个网格的体积， 则节点方程的累积项定义为间 处的绝对渗透率， 按调和平均数取值， Ari， di， k， k 分别为 i 与 k 节点的网格截面面积与距 离。
式中的所有待求变量 p， Sw， cp 都取 n+1 时刻的值， 即按照全隐式处理。
对于井的流量项， 分注入和产出两种情况， 同时考虑三维模型时一口井在纵向穿 过多个网格的情况。设油藏内注入井和生产井的数量分别为 Nw 和 Np， 记注入井和生产井序 号的集合分别为 κw 和 κp， 第 j 口注入井的共有 nwj 井段， 每个井段位于一个网格， 记这些
网格序号的集合为， 同样设生产井 j 所有井段所在的网格序号的集合为。所有的井均采用定总液量的工作方式， 按照流动系数分配各井段的流量， 并且采用全隐式方法计算流 动系数。
若注入井 j ∈ κw 的一个井段穿过第 i 个网格， 其总注水量为 ， 则第 i 组节点方程中井的流量项定义为聚合物注入浓度为
其中为流量分配系数， 令 WIi 为第 i 个网格的井指数， 则注入井的流量分配系数为
同样， 若生产井 j ∈ κp 的一个井段穿过第 i 个网格， 其总产液量则井的流量项定义为
其中生产井的流量分配系数定义为2)、 方程组的求解
利用全隐式差分格式离散后的聚合物驱模型式 (4) 是一个非线性代数方程组， 可 以写成如下一般形式
gn(un+1， un， vn) ＝ 0. (18)
对上述非线性代数方程组， 采用牛顿—拉斐森 (Newton-Raphson) 法求解， 其迭代公式为
k+1 k k un+1， ＝ un+1， +δun+1， .(20) 表示第 k 步迭代时系统的其中， un+1，k 表示第 k 步迭代系统状态的值，k k Jacobian 矩阵， 以下用 Jacn+1， 表示， δun+1， 为第 k 次迭代状态的变化量。方程 (19) 是一 k 个线性代数方程组， 左端的系数矩阵和方程的右端项已知， 求解该方程组可得出 δun+1， ， 从 n+1， k+1 而根据 (20) 计算出新的状态 u 。其具体过程为 ： n+1，
(1) 令 k ＝ 0， 选取初值 u 0 ＝ un ； k k
(2) 计算 Jacobian 矩阵 Jacn+1， 以及方程的右端项 -gn(un+1， ， un， vn) ； k+1
(3) 求解线性代数方程组 (19)， 根据式 (20) 计算出 un+1， ； n n+1， k n n n+1， k
(4) 若 ‖g (u ， u， v ‖ ≤ ε1 或 ‖δu ‖ ≤ ε2 则迭代停止， 否则返回第 (2) 步继续迭代。其中 ε1， ε2 为很小的正数， 分别为方程和被求变量所应满足的精度。3)、 Jacobian 矩阵的计算
在利用 Newton-Raphson 法求解全隐式方法离散的聚合物驱模型时， 需要计算系 统的 Jacobian 矩阵， 该矩阵是一个 nu×nu 维的稀疏矩阵， 可以表示为
以第 i 个节点方程为例， 计算节点方程 对 un+1 的导数。 在计算导数时先计算节点方程对中间变量的导数， 然后根据求导的链式法则计算中间变量对 un+1 的导数。
3、 离散最优控制模型
在离散的最优控制模型中， 设整个油藏被划分为 nb 个网格， 记网格序号为 i， i＝1， 2， ...nb。这样原连续性能指标中的注入项 qin 和产油项 fwqout 均可以用井的实际注入量 Qin 和产油量 Qo 表示。设油藏中注入井和生产井的数量分别为 Nw 和 Np， 记注入井和生产井 序号的集合分别为 κw 和 κp。考虑每口井穿越多个网格的情况， 设第 j 口注入井的共有 nwj 井段， 每个井段位于一个网格， 记这些网格序号的集合为 的网格序号的集合为 。 同样设生产井 j 所有井段所在 。由于每个井段位于一个网格， 则井段中的油相和水相中的流量即为网格节点的油相和水相流量项， 记第 i 个网格的油相和水相的流量分别为 Woi 和 Wwi， 流入 为正。则对于注入井的注入量 Qin 和生产井的产油量 Qo 可以表示为
其中网格节点流量和中的上标表示流量与 n+1 时刻的状态有关。设控制向量表示各水井的聚合物注入浓度， 其中 vj 表示第 j 口注入井聚合物的注入浓度 cpin j， j ＝ 1， 2， ...， Nw， 则聚合物驱的离散最优控制模型可以表示为 ：
s.t.gn(un+1， un， vn， n) ＝ 0，(23)0 ≤ v ≤ vmax， un|n ＝ 0 ＝ u0. 其中 τn， n ＝ 0， 1， ...， N-1 表示离散的时间步长 (d)， ξp 和 ξo 分别为聚合物的r 为折现率， tn 为时间 ( 年 )， 成本系数 ( 万元 /kg) 和原油的价格系数 ( 万元 /m3)， 为折现系数， r 为年折现率。性能指标 J 表示负的累积净现值。
4、 离散聚合物驱最优控制的必要条件
(1) 伴随方程
根据隐式离散极大值原理， 可得离散聚合物驱最优控制模型的伴随问题为 ：
伴随方程的具体表达式需要依照聚合物驱模型的离散化方法而定， 对于不同的离 散化方法， 其伴随方程的具体形式也不同。
(2) 梯度
根据隐式离散极大值原理， 可以获得性能指标的梯度为，
性能指标对第 j 个控制分量的导数为设注入井 j 的某个井段所在网格的序号为8， 由于控制变量只存在于井的流102352743 A CN 102352750说明书7/13 页量项中， 只考虑 的伴随向量为
即可。 对于第 i 个节点有该节点对应则 (24) 中的第二项对第 j 个控制分量的导数为则性能指标对 的梯度为5、 伴随方程的求解
在聚合物驱的最优控制问题中， 除了聚合物驱模型外还要求解由最优控制问题必 要条件获得的伴随方程， 这也是一组偏微分方程， 由于聚合物驱模型的复杂性， 这两组方程 的求解工作量很大。
根据离散最优控制问题的必要条件， 离散伴随方程的具体形式与支配方程的离散 化方法有关， 需根据不同的离散化方法推导伴随方程的具体形式。根据离散极大值原理得 到的伴随方程为关于 λn 的线性代数方程组， 求解时需要计算伴随向量前的两个系数矩阵
依据 Newton-Raphson 法的求解过程， 在由 un 计算 un+1 的迭代过程中， 第 k 步迭代时系统的 Jacobian 矩阵为
0 当 k ＝ 0 时， 若取 un+1， ＝ un， 则在计算 Jacobian 矩阵的过程中需要计算累积项 An+1 0 对 un+1， 的导数。根据累积项的定义 (10) ～ (13) 可得
设迭代过程在 k ＝ kf 步算法收敛， 取则此时的 Jacobian 矩阵为(31) 式即为伴随方程的一个系数矩阵的转置， 而 (32) 式与另一个系数矩阵的转 置相差一个时间步。实际上， 在利用 Newton-Raphson 法求解全隐式的聚合物驱模型时已经 计算出伴随方程的两个系数矩阵， 因此只要在求解过程中保存这两个系数矩阵， 就可以在 解伴随方程时直接调用。对于伴随方程的最后一项 由于性能指标 Jn-1 中只包 根据 (30) 式， 在 即可容易地求出含井的流量项 Wn， 因而只要计算出计算第 kf 步 Jacobian 矩阵时已经计算出流量项对状态变量的偏导数， 因而可以存储矩阵 供求解伴随方程使用， 这里仍相差一个时间步。 以上分析表明， 在求解聚合物驱模型的过程已经获得了计算伴随方程所需的导数 运算， 因而在求解过程中保存这些数据可以减少计算量， 并简化伴随方程的算法实现。
6、 最优控制求解
通过控制向量参数化， 聚合物驱最优控制问题的求解最终被转化为一个非线性规 划问题。将所有参数记为向量 x， 则转化后的非线性规划可表示为如下形式 ：
min F(x)，
s.t.gc(x) ≤ mpmax， (33)
xa ≤ x ≤ xb.
对于上述非线性规划问题采用惩罚函数法求解， 并采用子空间截断牛顿法来求解 其中的边界约束优化问题。
下面给出目标函数梯度的计算方法。令 xi 为在 t ∈ [ai，j， bi，j] 时间段内的确定 控制变量 vj 的某个参数， 记 ai， bi， nbi， 则目标函数对 j， j 对应的仿真时刻序号分别为 nai， j， j， xi 的导数为
通过求解上述非线性规划问题， 最终可以获得聚合驱最优控制问题的解。
本发明的有益效果主要表现在 ： 本发明方法可以在无人为干预下， 自动获得最优 开发方案， 且较传统的数模试凑方法有无人为干预， 更好的优化结果， 优化速度快等优点。
四、 附图说明
图 1 是本发明的方案求解框图 ； 图 2 是孤东 2 区块井位图 ； 图 3 是孤东 2 区综合含水率对比 ； 图 4 是孤东 2 区增油效果图。五、 具体实施方式
下面结合附图对本发明作进一步描述。
参照图 1， 聚合物驱提高原油采收率的最优控制方法， 所述方法包括如下步骤 ：
1)、 利用全隐式方法将聚合物驱模型离散化， 将求解过程离散为 N 个时间步， 则支 配方程的计算过程为 ：
(1) 令 n ＝ 0， 模型初始化， 初值为 u0 ；
(2) 已知 un， 利用 Newton-Raphson 法求解方程式 (4) 获得 un+1， 保存
un+1，τn 和 vn ；(3) 令 n ＝ n+1， 若 n ＜ N 返回第 (2) 步继续执行， 否则结束。 2)、 离散的伴随方程需要逆时间求解， 共需要求解 N-1 步。其计算过程为 ： (1) 令 n ＝ N， 计算出伴随向量的终值 λN， 保存 λN ； (2) 令 n ＝ n-1， 已 知 λn+1， 读取伴随方程的系数矩阵 求解线性代数方程组得到 λn， 保存 λn ； 和(3) 若 n ＞ 1， 返回第 (2) 步继续执行， 否则结束。
3)、 基于子空间的梯度优化方法， 包括子空间有限存储拟牛顿法 (SLMQN)、 有效集 截断牛顿法 (ASTN) 和子空间截断牛顿法 (SSTN)。此类方法将决策变量划分为自由变量和 有效约束变量， 对于自由变量采用拟牛顿法或截断牛顿法计算搜索方向， 对于有效约束变 量采用修正的梯度法计算搜索方向， 从而能够有效地处理边界约束并且具有较快的收敛速
度。
定义有效约束集合 A(x) 及其补集 B(x) 为 A(x) ＝ {i ： ai ≤ xi ≤ ai+εb or bi-εb ≤ xi ≤ bi}，其中 A(x) 中的变量为有效约束变量， 按照修正的梯度法计算搜索方向。B(x) 中 的变量为自由变量， 对应的搜索方向为子空间拟牛顿方向或子空间截断牛顿方向。在计算 A(x) 变量的搜索方向时， 把 A(x) 划分为 A1(x)， A2(x)， A3(x)3 个子集， 其定义如下 ：
A1(x) ＝ {i ： xi ＝ ai or xi ＝ bi， gi(x)(bi+ai-2xi) ≥ 0}，
A2(x) ＝ {i ： ai ≤ xi ≤ ai+εb or bi-εb ≤ xi ≤ bi， gi(x)(bi+ai-2xi) ＜ 0}，
A3(x) ＝ {i ： ai ＜ xi ≤ ai+εb or bi-εb ≤ xi ＜ bi， gi(x)(bi+ai-2xi) ≥ 0}.
其中 g(x) ＝▽ f(x)， εb 为很小的整数， 如 10-6。A1(x) 为负梯度方向指向边界外 的变量序号集合， A2(x) 为负梯度方向指向可行域内部的序号集合， A3(x) 为负梯度指向可 行域边界的序号集合。对于这 3 个子集分别计算方向搜索方向为
上述方法在 i ∈ A2(xk) 时采用了负梯度方向。由于其指向可行域内部， 并且在一 k 定搜索步长下新生成的 xi 仍在可行域内部。因此可以把 i ∈ A2(x ) 的变量归类于自由变 量， 通过拟牛顿法或截断牛顿法计算方向， 从而加快收敛速度。 改进后的有效边界约束集合 A(x) 及其自由变量集合 B(x) 的定义为
A(x) ＝ A1(x) ∪ A2(x)，
A1(x) ＝ {i ： li ≤ xi ≤ li+εb， g(xi) ≥ 0}， (35)
A2(x) ＝ {i ： ui-εb ≤ xi ≤ ui， g(xi) ≤ 0}，
下 标 在 A(x) 中 的 决 策 变 量 称 为 有 效 约 束 变 量， 其最速下降方向指向可行 域 的 外 侧， 下 标 在 B(x) 中 的 变 量 为 自 由 变 量。 对 于 有 效 约 束 变 量， 定义搜索方向
为
对于自由变量的搜索方向利用截断牛顿法计算。需要注意的是， 这里仅针对 B(x) 子空间中的变量进行计算， 对应的牛顿方程为 k k
HB(x)(x )pB(x) ＝ -gB(x)(x ). (37) k k
其中 HB(x)(x )， gB(x)(x ) 分别表示 B(x) 子空间的 Hessian 阵和梯度， 并且在计算子 空间截断牛顿方向时也可以使用有限差分来计算 Hessian 阵与向量的乘积， 该近似公式为
其中 diB(x) 为确定子空间截断牛顿方向过程中第 i 步的共轭方向。 最终获得的搜索方向 pk 为有了搜索方向， 在进行一维搜索时需要保证新的 xk+1 仍在可行域内， 需要采用投影 搜索法， 即定义 φk 为
φk(ω) ＝ f(PΩ(xk+ωdk)). (40) n
其中 Ω ＝ {x ∈ R ， a ≤ x ≤ b} 为可行域， PΩ(x) 是 x 在 Ω 上的投影。
算法 1、 控制受限最优控制的子空间截断牛顿法
(1) 通过控制向量参数化将最优控制问题转化为带有边界约束的优化问题， 记参 n 数向量为 x ∈ R ；
(2) 令 k ＝ 0， 选择初始参数向量 xk ＝ x0， 求解状态方程与伴随方程， 计算目标函 数 f(xk) 与梯度▽ f(xk) ；
(3) 根据式 (35) 划分边界约束集合， 按照式 (41) 计算搜索方向 pk， 求解子空间截 断牛顿方向的内部迭代中按照式 (38) 计算 Hessian 阵与向量的乘积， 通过求解状态方程和 伴随方程来计算计算梯度▽ f(xk+εdi) ；
(4) 采用投影搜索法进行一维搜索， 确定步长 ωk， 使得 Armijo 条件成立。
(5) 令 xk+1 ＝ PΩ(xk+ωkpk)， 若 满 足 收 敛 条 件 则 迭 代 结 束， 否 则 令 k ＝ k+1， 返 回 第 (3) 步。 迭 代 结 束 的 收 敛 条 件 为 |f(xk+1)-f(xk)| ＜ εf， 或 ‖xk+1-xk‖ ＜ εx， 或 ‖ ▽ f(xk+1)‖ ＜ εg， 其中 εf， εx， εg 为很小的正数。
对最优控制问题中含有一般形式控制约束的情形， 参数化转换后变为一个非线性 规划问题， 可通过惩罚函数法转换为一系列只有边界约束的非线性规划问题求解， 并利用 子空间截断牛顿法求解边界约束问题。
算法 2、 状态约束最优控制问题的罚函数法
(1) 令 k ＝ 0， 选择初始惩罚因子 μk ＝ μ0， 初始控制 (2) 以 为初值， 求解问题 (42)， 获得最优控制为 uk ；对于状态约束 gc， 先构造惩罚函数把含有状态约束的问题转化为只含有控制变量 边界约束的最优控制问题， 定义惩罚性能指标如下
其中 [y]- ＝ max(0， -y)， μ 为惩罚因子， i 表示第 i 个状态约束。然后， 采用算法 1 求解只含有控制变量边界约束的最优控制问题。
(3) 若 cCmax ≤ ε， ε ＞ 0， 则停止迭代 ；
(4) 选取控制的初值为-令 μk+1 ＝ ρμk， ρ ∈ (0， 1)， k ＝ k+1， 返回第 (2) 步。其中 cCmax ＝ max([gci] )， ε 为很小的整数， 为约束条件所要满足的精度。
实例 ： 孤东 2 区筛选出注聚区面积 4.2km2( 见图 2)， 有效厚度 12.4m， 入选储量 4 4 3 1039×10 t， 孔隙体积 1710.0×10 m 。模型设计水井 73 口， 油井 159 口。平面网格划分采 用 50m×50m 的均匀矩形网格， X 方向网格数 80， Y 方向网格数 110。在矿场实施注聚过程中， 聚合物溶液均是由配注站配置的， 各井的注入浓度与配 置聚合物溶液浓度相同， 难以实现对各井注聚浓度的分别调节。 因此在方案优化过程中， 各 井采用三段塞的注入方式， 取三段塞的浓度及各段塞的注入时间为优化变量。原注入方案 为单段塞 1.8g/L， 聚合物干粉用量为 15000t， 优化后注入段塞浓度分别为 2.5g/L， 1.93g/ L， 1.48g/L， 各井的注入方案如下表 1 所示。较水驱和原始方案的增油效果统计见表 2。后 续增油效果预测见表 3。图 3 为水驱， 初始方案， 优化方案， 实际注入的综合含水率对比图。 图 4 为增油效果图。
表 1 单井注入方案
表 2 增油效果统计14102352743 A CN 102352750 水驱 累积净现值 ( 万元 ) 增加净现值 ( 万元 ) 累积产油量 ( 万吨 ) 累积增油量 ( 万吨 ) 提高采收率 (％ ) 注聚量 ( 吨 )
说明书原方案 398195.9 156810.18 157.84 67.15 6.46 15000 优化方案 408957.83 167572.14 161.89 71.19 6.85 1500013/13 页241385.7 0.00 90.70 0.0 0.0 0表 3 增油效果预测表以上阐述的是本发明给出的实例表现出的优良效果， 显然本发明不只是限于上述 实例， 在不偏离本发明基本精神及不超出本发明实质内容所涉及范围的前提下对其可作种 种变形加以实施。本发明对求解其他行业过程优化问题有一定借鉴意义， 可推广应用于能 源、 交通、 冶金、 石化、 轻工、 医药、 建材、 纺织等行业的优化问题。