基于主成分分析和代理模型輔助的多目標(biāo)生產(chǎn)優(yōu)化方法

姚為英， 馮高城*， 任宜偉， 尹彥君， 王中正, 李振宇

(1.中海油能源發(fā)展股份有限公司工程技術(shù)分公司， 天津 300452; 2.中國石油大學(xué)(華東)石油工程學(xué)院， 青島 266580)

基于油藏數(shù)值模型的生產(chǎn)優(yōu)化是油藏開發(fā)管理的關(guān)鍵，即尋找最優(yōu)井控方案，改進油藏開發(fā)管理策略。通過優(yōu)化每口井的控制方案，獲得最大的油藏開發(fā)經(jīng)濟效益，一般用經(jīng)濟凈現(xiàn)值 (net present value, NPV)來衡量[1-3]。傳統(tǒng)的油藏生產(chǎn)優(yōu)化算法大都使用單個精確的油藏數(shù)值模型，結(jié)合優(yōu)化算法調(diào)整油藏區(qū)塊內(nèi)的注采井的工作參數(shù)，使其維持在最優(yōu)井控的工作狀態(tài)，從而最大化生產(chǎn)效益[4]。然而，此類方法忽略了油藏地質(zhì)的不確定性，由于一些未知的儲層參數(shù)，如滲透率和孔隙度，可能會顯著影響儲層性質(zhì)，因此油藏地質(zhì)具有很高的不確定性[5]。如果不考慮油藏地質(zhì)的不確定性，采用單一模型優(yōu)化得到的開發(fā)方案在現(xiàn)場實施時可能會面臨不同程度的開發(fā)風(fēng)險。因此，開發(fā)方案在實施前必須考慮油藏地質(zhì)的不確定性并進行多目標(biāo)風(fēng)險分析。Brouwer等[6]首次提出通過采用歷史擬合得到的多個概率模型來表征地質(zhì)參數(shù)的不確定性，然后綜合考慮多個油藏模型進行生產(chǎn)優(yōu)化，這類生產(chǎn)優(yōu)化問題被稱為魯棒生產(chǎn)優(yōu)化問題。此外，多目標(biāo)優(yōu)化方法[3]也常用于考慮地質(zhì)不確定性時進行油藏魯棒生產(chǎn)優(yōu)化，該方法能夠同時優(yōu)化歷史擬合得到的多個概率模型的平均經(jīng)濟凈現(xiàn)值和經(jīng)濟凈現(xiàn)值的離散程度，從而得到一組權(quán)衡多個目標(biāo)的最優(yōu)生產(chǎn)方案。然而，直接將多目標(biāo)優(yōu)化算法應(yīng)用于魯棒生產(chǎn)優(yōu)化問題時，需要面臨決策變量規(guī)模大和目標(biāo)函數(shù)評估計算耗時長等問題[5]。對此，國內(nèi)外學(xué)者近年來提出了多種基于代理模型的多目標(biāo)優(yōu)化算法：Chang等[7]將多項式擬合引入到考慮地質(zhì)不確定性的多目標(biāo)井位魯棒優(yōu)化中，與傳統(tǒng)方法相比有效降低了計算成本；馬寧[8]提出了一種基于遺傳算法的多目標(biāo)優(yōu)化建設(shè)項目調(diào)度模型，大大優(yōu)化了建設(shè)項目的調(diào)度過程；王鵬宇等[9]使用分布式訓(xùn)練的多目標(biāo)算法，顯著提高了多目標(biāo)算法優(yōu)化的速度。然而，現(xiàn)有代理模型應(yīng)用于高維參數(shù)空間時性能會顯著下降，導(dǎo)致處理油藏多目標(biāo)生產(chǎn)優(yōu)化問題時效果欠佳[10]。為解決上述問題，現(xiàn)提出結(jié)合主成分分析降維算法和代理模型輔助的多目標(biāo)生產(chǎn)優(yōu)化新方法，將該方法應(yīng)用到標(biāo)準(zhǔn)測試模型，并將計算結(jié)果與現(xiàn)有進化算法進行對比分析。

1 油藏魯棒生產(chǎn)優(yōu)化問題

隨著油氣田開采程度的不斷加深，開采方式的經(jīng)濟效益越來越被油田生產(chǎn)管理者所重視，因此通常需要選擇最優(yōu)的工作措施和開發(fā)方案。例如，在水驅(qū)油田開發(fā)過程中，需要通過提高原油產(chǎn)量的收益去補償注水措施帶來的額外成本。油藏生產(chǎn)優(yōu)化的目的是最大化油藏的開發(fā)潛力，從而獲得最大的油田開發(fā)利潤。因此，采用經(jīng)濟凈現(xiàn)值作為油藏生產(chǎn)優(yōu)化問題的性能指標(biāo)，表示在整個油田生產(chǎn)周期里，項目預(yù)期實現(xiàn)的經(jīng)濟效益與增產(chǎn)措施所產(chǎn)生的投資成本之間的差額。此時，水驅(qū)注采優(yōu)化模型的目標(biāo)函數(shù)可被定義為

(1)

式(1)中：J(x,s)為經(jīng)濟凈現(xiàn)值，$；x為控制變量；s為油藏狀態(tài)向量；Nt為總控制時間步數(shù)；b為年折現(xiàn)率，%；pt為第n個控制步的時間尺度參數(shù)，年；Δtn為第n個調(diào)控步的時長，d；ro表示油價，$/桶(1桶=158.98 L)；Qo,n為第n個控制步的平均產(chǎn)油速度，1桶/d；rw為產(chǎn)出水的處理費用，$/桶；Qw,n為第n個控制步的平均產(chǎn)水速度，1桶/d；rwi表示注入水費用，$/桶；Qwi,n為第n個控制步的平均注水速度，1桶/d。

魯棒生產(chǎn)優(yōu)化問題本質(zhì)上是一個多目標(biāo)優(yōu)化問題，即整個油藏系統(tǒng)的平均性能指標(biāo)最大化，同時相關(guān)的開發(fā)風(fēng)險最小化。其中，相關(guān)的開發(fā)風(fēng)險通常采用目標(biāo)函數(shù)在所有可能模型上的方差來衡量。最后，它可以生成具有一組最優(yōu)解決方案的Pareto前沿[11]，供決策者在不同的風(fēng)險水平下分析實施各解決方案的優(yōu)劣情況。

因此，在考慮實際油藏生產(chǎn)中地質(zhì)不確定性因素的基礎(chǔ)上，建立以NPV最大化的水驅(qū)油藏魯棒生產(chǎn)模型為

(2)

式(2)中：目標(biāo)函數(shù)Je(x)為平均經(jīng)濟凈現(xiàn)值；Ne為基于概率生成的能夠反映油藏特征的模型個數(shù)。

需要注意的是，僅僅優(yōu)化性能指標(biāo)的數(shù)學(xué)期望并不能降低不確定性，反而會導(dǎo)致一個很大的方差區(qū)間。針對上述問題，魯棒生產(chǎn)優(yōu)化具有更好的穩(wěn)定性，且優(yōu)化得到的開發(fā)方案更符合油田現(xiàn)場實際開發(fā)需求。換句話說，需要降低開發(fā)方案對地質(zhì)不確定性參數(shù)的敏感性，從而提高優(yōu)化結(jié)果的可靠性。

從平均經(jīng)濟凈現(xiàn)值和最小經(jīng)濟凈現(xiàn)值角度出發(fā)，建立降低不確定性油藏開發(fā)風(fēng)險的多目標(biāo)優(yōu)化模型。假設(shè)井控只考慮最簡單的邊界約束，則該魯棒生產(chǎn)優(yōu)化問題的數(shù)學(xué)模型為

(3)

2 算法機制

代理模型方法通過構(gòu)建計算廉價的近似模型代替原始數(shù)值模型，可以大幅度降低原始模型的評估代價，因此能夠有效提高優(yōu)化效率。基于代理模型的方法通常建立在具有強大全局搜索能力的進化算法基礎(chǔ)上。采用目前流行的基于分解的多目標(biāo)進化算法(multi-objective evolutionary algorithm based on decomposition, MOEA/D) 作為優(yōu)化求解器[12]。

2.1 基于分解的多目標(biāo)進化算法

基于分解的多目標(biāo)進化算法(MOEA/D)是一種流行的啟發(fā)式算法，用于解決全局優(yōu)化問題，并廣泛應(yīng)用于工程優(yōu)化問題。經(jīng)典的MOEA/D算法如下：一般來說，多目標(biāo)優(yōu)化問題 (multi-objective optimization problems, MOPs) 由多個互相沖突的目標(biāo)函數(shù)組成，一個具有m個目標(biāo)函數(shù)的多目標(biāo)優(yōu)化問題的數(shù)學(xué)模型表達式為

(4)

式(4)中：x=(x1,x2,…,xd)∈X?Rd為優(yōu)化變量，X為已知的決策空間，d為設(shè)計空間的維數(shù)；m為目標(biāo)函數(shù)的個數(shù)；F(x):Rd→Rm為目標(biāo)向量，所有F(x)構(gòu)成了未知的目標(biāo)空間；gi為p個不等式約束；hj為q個等式約束。

在此給出多目標(biāo)優(yōu)化中Pareto最優(yōu)的相關(guān)概念。給定兩個可行解x1和x2，當(dāng)且僅當(dāng)滿足

(5)

時，稱x1支配x2，顯然，此時解x1優(yōu)于x2，并記作x1

給定解x*∈X，當(dāng)不存在x滿足x

基于分解的多目標(biāo)優(yōu)化算法首先需要在目標(biāo)空間對多目標(biāo)優(yōu)化問題進行分解，在本文中，使用了切比雪夫的分解方法，示意圖如圖1所示。

圖1 切比雪夫分解法示意圖Fig.1 Schematic diagram of Tchebycheff decomposition

假設(shè)一組均勻分布的權(quán)重向量{λ1,λ2,…,λN}，逼近真實Pareto前沿的優(yōu)化問題可以通過切比雪夫分解法分解為N個標(biāo)量化的子問題，其中每個子問題的表達式為

(6)

2.2 代理模型技術(shù)

代理模型是為了降低計算復(fù)雜度而建立的真實模型的簡化近似[13]。因此，代理模型可以使用較低的評估成本預(yù)測新解的適應(yīng)度值。如果真實的模擬值表示為f(x)，代理模型的預(yù)測值表示為f′(x)，則f′(x)=f′(x) +e(x)，其中e(x)為近似誤差。f(x)的內(nèi)部計算過程是不需要被完全了解的，只有輸入與輸出的映射關(guān)系很重要，這類統(tǒng)稱為黑盒建模。其中，廣泛應(yīng)用的代理技術(shù)有多項式回歸、克里金、徑向基函數(shù)以及支持向量回歸等等，主要通過數(shù)據(jù)插值或回歸的方法擬合決策變量與適應(yīng)度函數(shù)之間的映射關(guān)系。本文中采用克里金插值模型[14]來逼近目標(biāo)函數(shù)。目前為止，克里金是應(yīng)用最廣泛的代理技術(shù)，主要是因為它能夠提供預(yù)測點的不確定性信息，這在代理的模型管理策略中非常有價值。

為了構(gòu)建計算廉價的代理模型y=g(x)，對于任意的個體x，克里金模型假設(shè)其目標(biāo)函數(shù)值可被估計為

g(x)=μ+ε(x)

(7)

式(7)中：ε(x)～N(0,σ2)；μ和σ為兩個與x無關(guān)的常數(shù)。換句話說，g(x)的先驗分布為高斯分布，且均值為μ，方差為σ2。

此外，對于任意x，x′∈Xd，ε(x)與ε(x′)之間的關(guān)聯(lián)式c(x,x′)，取決于|x-x′|。更準(zhǔn)確地表達為

c(x,x′)=exp[-d(x,x′)]

(8)

(9)

式(9)中：θi>0 且 1≤pi≤2；當(dāng)d(x,x′)→0時，c(x,x′)→1，當(dāng)d(x,x′)→+∞時，c(x,x′)→0，這種特征在為連續(xù)函數(shù)f(x)建模時是符合要求的；pi與g(x)在xi處的平滑度有關(guān)；θi表示xi對g(x)的重要性。更多關(guān)于克里金建模的細節(jié)見文獻[14]。

給定K個點x1,x2,…,xK∈Xd，與其對應(yīng)的函數(shù)值y1,y2,…,yK，超參數(shù)μ,θ1,θ2,…,θd和p1,p2,…,pd可以通過g(x)=yi在x=xi(i=1,2,…,K) 處的最大似然估計得到，即

(10)

式(10)中：y=(y1,y2,…,yK)T；I為K維列向量；C為K×K的矩陣，表達式為

(11)

為了最大化式(10)，μ和σ2的值分別取為

(12)

(13)

將式(12)和式(13)代入式(10)中，消除了式(10)中的未知參數(shù)μ和σ2。因此，似然函數(shù)僅依賴于θi和pi。然后，可以使用優(yōu)化方法最大化式(10)來估計θi和pi，再通過式(12)和式(13)可以很容易地估算出μ和σ2。

(14)

其均方誤差為

(15)

(16)

2.3 數(shù)據(jù)預(yù)處理的降維技術(shù)

盡管克里金是一個很流行的能夠提供不確定性信息的代理模型，但由于訓(xùn)練模型的計算復(fù)雜度，它也存在潛在的嚴(yán)重缺陷。研究發(fā)現(xiàn)，大規(guī)模優(yōu)化問題的克里金模型訓(xùn)練時間是無法忽略的，因此，訓(xùn)練時間可能會比真實評估目標(biāo)函數(shù)的時間長，從而破壞了降低計算成本的整體目標(biāo)。

由于油藏生產(chǎn)優(yōu)化方案設(shè)計的要求不斷提高，從決策變量空間來講，設(shè)計的變量愈加復(fù)雜，由此帶來變量數(shù)目的劇增?！熬S數(shù)災(zāi)難”問題的限制導(dǎo)致直接使用克里金插值構(gòu)建代理模型近似誤差大，且近似模型訓(xùn)練參數(shù)時間長，失去了代理輔助的意義，因此本文中采用主成分分析方法[17](principal component analysis, PCA) 對優(yōu)化問題進行降維處理。

PCA的目標(biāo)就是在盡力保證信息量不損失的前提下，對原始的數(shù)據(jù)特征進行降維，從而在較小的維度下展示數(shù)據(jù)的特征。PCA的原理是利用正交變換將潛在線性相關(guān)的樣本轉(zhuǎn)變?yōu)橐唤M線性不相關(guān)的低維變量，轉(zhuǎn)變完之后的這組變量稱作主成分。假設(shè)原始樣本點xi在新的低維空間上的投影為WTxi，根據(jù)最大可分性，所有的樣本點的投影需盡量分開，即最大化樣本點投影后的方差，因此優(yōu)化的目標(biāo)函數(shù)為

(17)

將原始的樣本集X=(x1,x2,…,xd)轉(zhuǎn)變?yōu)樾碌膌維數(shù)據(jù)集，只需要將協(xié)方差矩陣XXT進行特征值分解，并排序取前l(fā)個值對應(yīng)的特征向量構(gòu)成W=(w1,w2,…,wl)，這就是PCA的解。

通過PCA對高維的設(shè)計空間進行降維處理，把較高維數(shù)的油藏設(shè)計變量在低維空間表征，從而提高代理模型的近似精度以及優(yōu)化效率。

2.4 K-MOEA/D-PCA算法優(yōu)化計算步驟

K-MOEA/D-PCA算法包含兩個循環(huán)，如圖2所示，主循環(huán)和次循環(huán)分別是使用油藏數(shù)值模擬器和代理模型的優(yōu)化過程。該算法有3個步驟：①初始化；②克里金代理模型輔助MOEA/D進化搜索；③更新代理模型。

本文使用拉丁超立方試驗設(shè)計方法[18]從可行域中采樣生成初始訓(xùn)練樣本點。然后，利用真實的油藏數(shù)值模擬器對目標(biāo)函數(shù)值進行計算。這些經(jīng)過真實評估的解決方案存儲在樣本庫D中。以最大的函數(shù)評估數(shù)作為進化優(yōu)化過程的終止條件。

隨后，為每個目標(biāo)函數(shù)訓(xùn)練克里金代理模型，以減少計算成本。由于克里金模型能夠在誤差范圍內(nèi)進行預(yù)測，因此得到了廣泛應(yīng)用。另一方面，處理高維優(yōu)化問題時，構(gòu)造克里金模型的計算時間過長又導(dǎo)致該方法使用受限。因此，在構(gòu)建代理模型之前，采用主成分分析克里金將樣本數(shù)據(jù)從高維映射到低維。因此，低維數(shù)據(jù)可以提高克里金模型的精度。代理模型的管理策略對代理輔助進化計算的性能是至關(guān)重要的，即選擇個體使用數(shù)值模擬器重新評估的策略，然后使用這些重新評估的解決方案來更新代理模型，并將其添加到樣本庫D中。因此，需要綜合考慮這些解決方案的收斂性和多樣性。

圖2 K-MOEA/D-PCA算法流程Fig.2 Flow diagram of K-MOEA/D-PCA

3 計算實例

利用三維油藏模型測試K-MOEA/D-PCA算法在處理考慮地質(zhì)不確定性條件下魯棒生產(chǎn)優(yōu)化問題時的有效性和收斂速度。采用基準(zhǔn)測試模型(Egg模型)的集合版本(包含100個不確定模型)，該版本是公開的，詳細地質(zhì)統(tǒng)計學(xué)參數(shù)可見文獻[19]。由于模型之間具有相似性，基于100個數(shù)值模型的評估耗時長且浪費計算資源，因此在魯棒生產(chǎn)優(yōu)化過程前，采用多維尺度變換和K-means聚類算法來減少計算量，但不影響結(jié)果的準(zhǔn)確性。利用聚類算法，將隨機模型按滲透率場的相似度分為6類，并從每類中選擇一個具有代表性的模型。如圖3所示，每個模型共有12口直井，其中8口注水井，4口生產(chǎn)井。每個模型包含60 × 60 × 7個網(wǎng)格，其中18 553個為活網(wǎng)格。

生產(chǎn)井井底壓力恒為395 bar(1 bar=0.1 MPa)，因此控制輸入變量為每口注水井的流量，其上下限分別設(shè)置為0、100 桶/d。經(jīng)濟參數(shù)中的油價、注水成本、生產(chǎn)成本分別設(shè)置為45、5、5 $/桶，假設(shè)折現(xiàn)系數(shù)為0。每個油藏模型的生產(chǎn)周期為3 600 d，時間步長為180 d，因此，優(yōu)化決策變量的個數(shù)為8 × 20 = 160。注意，每個解決方案的目標(biāo)函數(shù)評估是基于六個油藏數(shù)值模型的。

圖4顯示了在優(yōu)化過程中平均和最小經(jīng)濟凈現(xiàn)值的變化情況?？梢钥闯?，未進行優(yōu)化時，初始解構(gòu)成的種群平均和最小經(jīng)濟凈現(xiàn)值均較低，通過提出的K-MOEA/D-PCA算法進行優(yōu)化后，能夠找到位于高經(jīng)濟凈現(xiàn)值區(qū)域的候選解，最終得到了帶來最高回報的最優(yōu)解，顯著提高了優(yōu)化方案的經(jīng)濟凈現(xiàn)值。

圖5顯示了經(jīng)過本文提出的K-MOEA/D-PCA算法優(yōu)化后，6個不同模型的最優(yōu)注水開發(fā)方案。在每個圖中，橫坐標(biāo)為調(diào)控的時間步，縱坐標(biāo)為不同注水井的注水量。所有模型不同時間步的注水量被限制在0～100 桶/d。

使用圖5所示的優(yōu)化方案對相應(yīng)的模型進行模擬開發(fā)，生產(chǎn)過程中累積產(chǎn)油量、累積產(chǎn)水量、累積注水量和含水率的變化情況如圖6所示?？梢钥闯?，考慮到模型的不確定性，生產(chǎn)3 600 d后，該模型的累積產(chǎn)油量可以達到(4 ～ 4.5)×105m3，含水率約為0.9。

圖3 Egg模型對數(shù)滲透率分布Fig.3 Log-permeability distribution of Egg models

為了保證解的穩(wěn)定性，避免偶然性，為不同模型執(zhí)行10次獨立計算，最終的優(yōu)化結(jié)果如圖7所示，箱型圖中的紅線代表10次運行結(jié)果的中位數(shù)?？梢钥闯?，由于不同模型的滲透率場不同，開發(fā)過程實現(xiàn)的經(jīng)濟凈現(xiàn)值也不同，經(jīng)過本文提出的算法優(yōu)化后的生產(chǎn)方案可以實現(xiàn)滿意的經(jīng)濟效益。

圖4 優(yōu)化過程經(jīng)濟凈現(xiàn)值的變化Fig.4 The change of NPV during optimization process

圖5 不同概率模型中各注水井的最優(yōu)注水方案Fig.5 The optimal water-injection scheme for each injection well within different models

圖6 不同模型開發(fā)方案下各指標(biāo)的變化Fig.6 The change of each indicator for the model with different development schemes

圖7 不同模型的開發(fā)方案的經(jīng)濟凈現(xiàn)值Fig.7 Net present value of development schemes for different models

4 結(jié)論

針對油藏魯棒生產(chǎn)優(yōu)化問題的大規(guī)模與計算耗時的特性，提出了新的求解思路。由于復(fù)雜油藏注采井?dāng)?shù)多、非線性強，直接將代理模型引入進化算法會導(dǎo)致優(yōu)化結(jié)果變差，因此提出了K-MOEA/D-PCA方法，通過主成分分析對數(shù)據(jù)進行降維預(yù)處理，然后采用先分后治的求解策略，利用權(quán)重向量將目標(biāo)空間分解為多個子問題從而協(xié)同計算，并通過構(gòu)建代理模型動態(tài)規(guī)劃搜索區(qū)域，實現(xiàn)了高效全局優(yōu)化。應(yīng)用K-MOEA/D-PCA方法對Egg模型進行了實例計算，優(yōu)化得到的最優(yōu)開發(fā)方案顯著提高了油田開發(fā)的經(jīng)濟凈現(xiàn)值。