基于既往優化知識的穩健性優化方法

呂國臣，程遠勝，易家祥，劉均

華中科技大學 船舶與海洋工程學院，湖北 武漢 430074

0 引 言

產品設計會存在很多不確定性因素[1]，例如船舶結構的材料參數、尺寸、載荷。這些不確定性因素的波動可能使產品性能參數發生變化，導致產品性能的波動量超過允許范圍而使產品不滿足設計要求。穩健性設計優化問題最早由Taguchi[2]提出，目的是在存在不確定性的情況下，使目標/約束函數的值保持在可接受的范圍內，同時獲得盡可能好的設計問題的解。學者們也對不同問題做 了 工 程 應 用 研 究[3-4]。

不確定性問題主要包括模型不確定性、參數不確定性以及預測不確定性等[5]。本文針對設計變量或參數在區間上存在的不確定性進行研究。傳統解決非線性區間不確定性穩健優化問題的方法通常為嵌套優化方法，這會導致大量的資源消耗。Siddiqui 等[6]提出了一種修改的benders 分解方法來降低計算量，即通過離散化不確定性區間來解耦嵌套優化。但由于該方法對不確定性區間進行了離散，并不能保證優化解的穩健性。Zhou等[7]提出了一種區間不確定性魯棒優化的序列二次規劃方法，該方法兼具目標魯棒性和可行性魯棒性。Hu 等[8]提出了一種基于系統性能對不確定性的靈敏度雅可比矩陣的有界約束問題穩健性指標計算方法，雅可比矩陣可以定量地測量系統允許的最大變化幅度。辛鵬飛等[9]運用基于切比雪夫多項式的區間擴張函數，將含區間參數的微分代數方程轉化為切比雪夫多項式插值點處確定參數的動力學方程，減少了評估穩健性的仿真工作量。王瓊等[10]利用區間可能度處理不確定性約束，將穩健性優化模型轉換為確定性多目標優化模型。這些方法雖然能夠評估穩健性，但都需要利用函數的梯度信息，對于需要采用數值仿真方法計算的實際工程問題，很難利用其梯度信息，所以基于梯度的求解方法不適用。Gunawan 等[11-12]根據逆向思維和最壞情況分析，提出了反映設計方案穩健性的指標，對其施加約束，并作為約束條件加入到優化數學模型中，構成穩健性設計優化問題；但該方法采用了逆向思維來反向映射自變量和函數區間，難以保證穩健優化解的精度。所以對于變化趨勢不確定的黑箱工程優化問題，基于嵌套進化算法能夠有效解決此問題，但會消耗大量的計算資源。Cheng 等[13]將局部搜索算法應用到嵌套優化中加速其收斂，但所需的計算資源仍然巨大。

根據前人文獻的算例經驗，若直接利用基于嵌套進化算法來求解穩健性優化問題，通常會調用百萬次甚至千萬次的原函數，在進行產品設計如船舶結構優化時會消耗大量的計算資源和時間。故本文提出一種基于既往優化知識的穩健性優化方法，以大幅度減少計算資源，提高優化效率。

1 問題描述及求解

1.1 穩健性優化數學模型

不失一般性地，區間穩健性設計優化問題的數學模型可以由式(1)表示：

其中：

式中：x為D維設計變量，xd為 確定性變量；xu為不確定性變量；x和x分別為確定性變量的下界和上界；x和x分別為不確定性變量的下界和上界； δxu為 不確定性變量的波動值； Δxu為不確定性變量的最大波動范圍； ηf為目標穩健性指標； ηg為約束穩健性指標； Δf0為目標函數最大允許波動量；f為目標函數；gm為小于等于形式的約束函數。

需要指出的是，對于確定性優化，其約束條件可以統一寫成gm≤0 的 形式，當 δxu=0時，式(1)中約束穩健性指標 ηg≤0就表示確定性優化的約束條件。目標穩健性指標為不確定性變量變化下的目標函數的最大波動量與最大允許波動量的比值，其約束穩健性指標為不確定性變量變化下的所有約束函數的最大值。

1.2 基于差分進化的區間不確定性穩健性優化問題求解

1.2.1 差分進化算法

差分進化算法（DE）是工程優化領域應用較多的一種智能進化算法，它由Storn 等[14]提出。該算法能夠較好地對函數進行全局搜索，被廣泛應用到工程優化問題中。差分進化算法存在多種進化格式，本文采用DE-rand/1/bin 格式。其主要步驟如下：

通常，對于一個種群數為N的種群，第G代的第i個體可以表示為=(p,···,p)，i=1，2，···，N，其下一代對應的變異個體由式(2)產生：

式中：v為（G+1）代變異個體；F為放大因子，是定 義 在 [0,2]范 圍 內 的 常 數 因 子；k1，k2和k3為 從{1,2,···,N}中隨機產生的3 個整數，且要滿足k1≠k2≠k3≠i。

然后，對種群個體進行交叉操作，遵循式（3）的機制將變異個體與當前代的個體交叉產生下一代的新個體。

式中：CR為 在 [0,1]范 圍內的交叉因子；rand(j)為[ 0,1]區 間內產生的隨機數；jRand為從{1,2,···,N}中隨機產生的一個整數。

最后，選擇x與x中目標函數較好（或其他準則）的一個進化到下一代。

1.2.2 求解方法

針對式（1）表征的區間不確定性穩健性優化數學模型，一種常規的求解方法就是直接應用差分進化算法進行嵌套優化求解（DE-RO）。外層優化過程是基于差分進化求解最優目標函數值，內層優化是差分進化求解2 個穩健性指標時需要進行2 個無約束的優化，即目標穩健性指標和約束穩健性指標。本文內外層優化均采用相同的進化格式。

2 基 于 既 往 優 化 知 識 的 穩 健 優 化方法

若應用差分進化算法求解方法，在嵌套優化過程中會進行內層優化，從而大量增加目標函數和約束函數的調用次數。本文提出基于既往優化知識的穩健優化方法（DE-BPOK-RO），該方法涉及穩健性指標的近似評估以及基于穩健性指標重新評估的臨界距離調整2 個核心策略。

2.1 穩健性指標的近似評估

假設一個穩健性優化問題的不確定性變量的數目為2，畫出不確定性變量變化區間的示意圖，如圖1 所示。圖中，p1和p2為2 個外層個體，藍色點和綠色點分別為p1和p2進行內層優化生成的個體。

由圖1 可見，若p1和p2的空間距離較近，兩者變化區間的重疊區域較大，則其內層優化生成的個體之間會出現距離極近甚至重合的情況。在內層優化過程中，若能夠利用已經精確計算過的個體響應值的信息來近似預測其他個體響應值，則可在很大程度上降低原函數的調用次數。

圖1 不確定性變量變化區間示意圖Fig. 1 Schematic diagram of the variation interval of uncertain variables

建立精確計算過樣本點庫Dp={p1,p2,···}及對應的響應值庫Dt={t(p1),t(p2),···}。近似評估穩健性指標策略的核心是近似評估內層樣本點的響應值，需要定義2 個參數，分別為近似距離閾值R及上限近似點數Ns，Ns的 引入是為了降低R的取值對近似的影響。圖2 為近似評估內層個體響應值的示意圖，圖中藍色點表示精確計算過的樣本點，A，B為2 個個體，其響應值計算公式為：

圖2 近似評估內層個體響應值示意圖Fig. 2 Schematic diagram of approximate evaluation of inner individual response values

式中：t?(p)為 被評估個體的響應值；q為近似距離閾值內精確計算的樣本點數；pl為近似距離閾值內樣本點庫中經距離排序后的樣本點。需要指出的是，用來近似評估內層個體響應值的目標函數庫和約束函數庫均通過外層和內層在既往優化過程中的信息進行更新，是一個動態增加過程。

2.2 臨界距離調整策略

在近似評估穩健性指標時，若進化前期樣本點庫的信息不足，使內層樣本點近似得到的響應值誤差較大，這可能導致實際滿足穩健性約束條件的部分外層個體被誤評估為不穩健，從而誤導了整體進化過程。而隨著代數增加，計算的樣本點增加后，重新近似評估其穩健性指標會在一定程度上修正誤判解，基于此提出了一種基于穩健性指標重新評估的臨界距離調整策略。

在進行第G代外層優化時，建立第G代修正庫D={p|f(p)＜f∧gj(p)≤0 ∧(ηf(p)＞1‖ηg(p)＞0)，,p∈WG}， 其中f為第G-1 代的最優目標值，WG為截止到G代每代外層個體的集合。在判斷當前代最優之前，隨機選取修正庫中N個點重新近似評估穩健性指標，定義參數修正率閾值RC0，重新評價后穩健性變化（即由不穩健變為穩?。┑膫€數為Nc，將其與當前種群中min(Nc,N·RC0)（取整）數目的隨機個體進行替換。若RC=＞RC0，則近似距離閾值R=R/2 。需要指出的是，用于臨界距離調整的修正庫在優化的過程中不斷變化，其建立和更新是基于外層既往優化過程中可能搜索到的真正最優解附近的信息。

2.3 改進的可行性準則

在應用基于既往優化知識的穩健性優化方法求解穩健性優化問題時，新種群個體的選擇策略采用改進的可行性準則，用新生成的個體代替上一代個體需要符合以下幾種情況：

1） 新生成的個體是可行解，而原個體是非可行解；

2） 新生成的個體與原個體均是非可行解，而新生成個體約束違反的總體程度更??；

3） 新生成的個體與原個體均是可行解，且新生成的個體滿足穩健性約束條件而原個體不滿足；

4） 新生成的個體與原個體均是可行解，且新生成的個體和原個體均滿足穩健性約束條件而新生成的個體的目標函數值更??；

5） 新生成的個體與原個體均是可行解，且新生成的個體和原個體均不滿足穩健性約束條件而新生成的個體的穩健性指標的和更小。

2.4 算法求解框架及流程

算法的基本框架是基于差分進化求解策略的雙層循環模式，其中內層優化用于求解外層穩健性指標約束。具體流程圖如圖3 所示。

圖3 基于既往優化知識的穩健性優化方法流程圖Fig. 3 Flow chart of robust optimization method based on previous optimization knowledge

步驟1：初始化外層優化的種群。

步驟2：調用原函數精確計算初始外層個體的目標函數f(x)和 約束函數gj(x)。

步驟3：評估初始外層個體的目標穩健性指標 ηf和約束穩健性指標 ηg，即用差分進化算法求解外層個體在不確定區間內目標函數的最大變化量和最大約束函數值（本質是分別以 ηf和 ηg為內層目標函數的無約束優化問題）。

步驟4：用所有精確計算過的點（設計方案）建立初始目標函數庫和約束函數庫。

步驟5：根據DE-rand/1/bin 的差分進化格式生成新一代外層個體。

步驟6：調用原函數精確計算新一代外層個體的目標函數f(x)和 約束函數gj(x)。

步驟7：根據2.1 節提出的策略近似評估新一代外層個體的目標穩健性指標 ηf和約束穩健性指標 ηg。

步驟8：更新目標函數庫和約束函數庫。

步驟9：結合2.2 節提出的臨界距離調整策略對當代外層個體進行選擇并判斷當前最優。

步驟10：精確評估當前最優個體的穩健性指標。

步驟11：更新目標函數庫和約束函數庫。

步驟12：若滿足停止條件，則結束外層優化；否則，返回步驟5。

需要指出的是：

1） 只需要對滿足gj≤0 的個體進行穩健性指標評價，這在一定程度上減少了計算資源的消耗。

2） 精確評估穩健性指標就是在內層優化中計算所有內層個體的響應值。

3） 對于外層優化過程，其停止準則為進化代數達到了最大進化代數；對于內層優化過程，其停止準則為進化代數達到了最大進化代數或最優值在 μ （一般取20～30）代的差值小于等于10e-3。

3 算例測試結果及分析

采用2 個數學算例和1 個工程算例驗證DEBPOK-RO 算法的可行性與有效性。

3.1 算例測試結果

3.1.1 算例1

算例1 是一個目標函數比較復雜的測試函數，其數學表達式為

式中：x1為 不確定性變量，其變化區間為[-0.05,0.05]；x2為確定性變量；目標函數允許的最大波動值 Δf取為0.02。

分別運用DE-BPOK-RO 和DE-RO 算法求解以上算例，并與Cheng 等[13]提出的混合進化方法（DE-SQP-RO）得到的結果比較。對提出的DEBPOK-RO 算法中3 個關鍵參數進行研究發現，采用表1 所示的參數能得到最佳結果。

表1 DE-BPOK-RO 和DE-RO 算法參數Table 1 Parameters of DE-BPOK-RO and DE-RO algorithms

其中，R0=Δx1/5=0.01。為了降低智能優化算法的不確定性，DE-BPOK-RO 和DE-RO 算法都獨立運行了10 次，統計結果如表2 所示，其中FE為所有原函數調用的次數。根據表2，在目標函數最優的結果中，DE-BPOK-RO 算法求解得到的目標函數值與DE-RO 算法得到的目標函數值的誤差約為0.072%，與DE-SQP-RO 算法得到的誤差約為0.05%，而調用的函數次數僅為另外2 種算法的2.393% 和13.385%；10 次運行目標函數的中位數值要優于DE-RO 算法得到的結果，而原函數的平均調用次數僅為差分進化求解策略的2.620%；最差一次目標函數值雖然劣于DERO 算法得到的結果，但2 種算法得到的最優解均與全局的穩健最優解差距很大，可能由初始種群的隨機性導致。

表2 測試函數1 統計結果Table 2 Statistical results of test function-1

單次優化過程的收斂情況如圖4 所示。由圖可見，2 種算法都收斂到最優值附近。圖4 還展示了算法收斂到最優的局部情況。次優化中DEBPOK-RO 算法在67 代時收斂，最優值為-5.952 9，此時原函數調用次數為191 484；而DE-RO 算法要在78 代才能達到相同程度，原函數調用次數為5 259 690，直觀展示了DE-BPOK-RO 算法的優勢。

圖4 算例1 的DE-BPOK-RO 和DE-RO 算法收斂曲線Fig. 4 Convergence curves of DE-BPOK-RO and DE-RO algorithms of numerical example-1

通過以上結果分析，基于既往知識的穩健性優化方法在最優目標誤差小于2.4%或更優的情況下，大大減少了真實函數調用的次數，節省了97%以上的計算資源。

3.1.2 算例2

算例2 是一個約束函數比較復雜的測試函數，其數學表達式為

式中：不確定性變量x1和x2的變化區間均為[-0.4,0.4]； 目標函數允許的最大波動值 Δf=2.5。

采用與算例1 相同的參數，其中，R0=min(Δx1,Δx2)/5=0.08。采用DE-RO 算法和DE-BPOK-RO算法對此問題進行10 次獨立的優化求解，得到的統計結果如表3 所示，單次優化過程的收斂情況如圖5 所示。

表3 測試函數2 的統計結果Table 3 Statistical results of test function-2

圖5 算例2 的DE-BPOK-RO 和DE-RO 算法收斂曲線Fig. 5 Convergence curves of DE-RO-DE and DE-RO algorithms of numerical example-2

根據表3，在目標函數最優的結果中, DEBPOK-RO 算法求解得到的目標函數值與DE-RO算法得到的誤差約為0.243%，而調用函數次數僅為DE-RO 算法的0.956%；10 次運行的目標函數的中位數值的誤差約為2.273%，平均調用原函數的次數僅為DE-RO 算法的0.870%；最差一次目標函數值雖然劣于用DE-RO 算法得到的結果，但2 種算法得到的最優解均與全局的穩健最優解差距很大，表明可能是初始種群的分布對進化的影響較大。

從圖5 中可以看出，2 種算法都收斂到最優值附近。根據收斂曲線的局部圖，DE-BPOK-RO算法在28 代時收斂，最優值為-1.768 5，原函數調用次數為128 609；而DE-RO 算法要在38 代才能收斂到相同程度，而原函數調用次數為4 759 410，同樣證明了DE-BPOK-RO 算法的優勢。

由以上分析可知，基于既往知識的穩健性優化方法在最優目標誤差小于2.3%的情況下，節省了99%以上的計算資源。

3.2 壓力容器設計

選用的工程算例來源于經典的壓力容器設計（Kannan 等[15]），如圖6 所示。該壓力容器也可作為水下結構物內置液艙，獲得對設計變量波動不敏感的設計，即穩健性設計對確保這種內置液艙的安全性是非常重要的。壓力容器設計經過修改后成為一個穩健性優化問題，本文將比較差分進化求解策略和基于既往優化知識算法對此問題的優化結果，并結合混合優化方法進行對比分析。

圖6 壓力容器的中心和端部示意圖Fig. 6 Schematic diagram of the center and end of pressure vessel

原問題的目標是最小化包括圓柱容器材料、成型和焊接成本的總成本，1 個設計變量分別是x1(Ts，殼體厚度)、x2(Th， 頭部厚度)、x3(R，內徑)和x4(L，容器圓柱形部分的長度，不包括頭部)。修改后的優化問題如式(8)所示：

式中：不確定性變量x4為容器圓柱形部分的長度，其變化區間為 [-0.05,0.05]；其他變量為確定性變量，目標函數（總體設計成本）允許的最大波動值Δf=50。

采用與數學算例相同的算法參數，其中R0=Δx4/5=0.01。用差分進化求解策略和基于既往優化知識算法對此問題進行10 次獨立的優化計算，得到的統計結果如表4 所示。

表4 工程算例統計結果Table 4 Statistical results of engineering example

根據表4，最好的優化結果中，DE-BPOK-RO算法求解得到的目標函數值優于DE-RO 算法，與DE-SQP-RO 算法得到的最優值誤差僅為0.008%，而調用的函數次數僅為另外2 種算法的5.655%和8.229%；10 次運行的目標函數的中位數值也要優于用DE-RO 算法得到的結果，原函數的平均調用次數也僅僅為DE-RO 算法的5.936%；最差一次目標函數值雖然劣于用DE-RO 算法得到的結果，但出現的概率較低，且誤差約為4.193%，在工程優化可接受范圍內。通過以上分析可知，基于既往知識的穩健性優化方法在優于或者極少程度犧牲最優解的條件下，大幅減少了真實函數調用次數，對于工程優化問題具有較大意義。

在實際的船舶優化設計中，不確定性因素的存在可能會導致船舶性能達不到預定的設計目標。直接應用嵌套優化算法求解基于區間不確定性的船舶結構穩健性優化問題顯然不可行，本文提出的基于既往知識的穩健性優化方法雖大幅度節省了計算資源，但幾十萬次的有限元計算對船舶結構優化來說仍存在一定壓力。

劉婧等[16]將動態代理模型應用到船舶結構可靠性優化中，在未來工作中將把代理模型及本文提出的方法結合進一步降低計算成本，研究提出船舶結構強度或穩定性穩健優化設計方法。

4 結 語

本文針對區間不確定性穩健性優化問題提出了一種基于既往優化知識的穩健性優化方法。該方法基于嵌套差分進化算法的格式，利用臨界距離內已經精確計算過的個體響應值近似評估內層個體響應值，以此近似評估個體穩健性指標來大幅度地減少原函數的調用次數，降低計算成本。同時利用進化過程中逐步擴充的精確個體響應值的信息，選擇性地重新評估已往個體的穩健性，并根據評估誤判率自適應地調整臨界距離來降低穩健性誤判率。

本文通過2 個數學算例和1 個工程算例驗證了算法的有效性。通過與DE-RO 及DE-SQP-RO算法的綜合統計結果對比得出，DE-BPOK-RO 算法在最優目標誤差為2.5% 以下的情況下，節省了94%以上的計算資源，對工程優化問題也取得了比較好的效果。因此，提出的基于既往優化知識的穩健性優化方法，為解決區間不確定性的產品穩健性優化提供了新思路，大幅節省了基于雙層嵌套產品穩健性優化所需要的計算資源。