概率功能度量求解的共軛梯度步長調節法

易 平， 謝東赤

（大連理工大學 建設工程學部，大連116024）

1 引 言

工程結構中的材料屬性、幾何尺寸和外部荷載通常都具有隨機性，基于概率參數的優化設計不可避免，也更符合工程實際。概率結構優化設計PSDO（Probabilistic Structural Design Optimization）［1－5］通常表述為滿足概率約束條件下目標函數的最小化。直接的PSDO本質上是一個兩層次迭代問題，外部為確定性優化，內部為可靠性分析。內層可靠度分析可采用功能度量法PMA（Performance Measure Approach）［4］和 可 靠 指 標 法RIA（Reliability Index Approach）［5］，其中功能度量法更加穩定高效。

功能度量法通過在標準正態空間中搜索目標可靠指標球面上的最小功能目標點MPTP（Minimum Performance Target Point）獲得其概率功能度量。改進均值法 AMV（Advanced Mean Value）［6］常用于確定 MPTP，但對于一些非線性程度較高的功能函數，AMV迭代可能會陷入周期振蕩或混沌解。因此，相繼提出多種改善其收斂性的算法。Youn等［7－9］利用連續三次迭代點處的梯度信息和功能函數凹凸性，提出了混合均值法HMV（Hybrid Mean Value method），然而該算法對某些問題的求解仍因不收斂而失效［10］；Yang等［11］從混沌動力學的角度出發，利用穩定轉換法對AMV迭代格式實施收斂控制；Meng等［12］改進混沌控制方法CC（Chaos Control），提出了混合混沌控制策略MCC（Modified Chaos Control），兩種混沌控制方法雖然提高了效率，但在實際應用中，尤其是對于高維問題，很難選取合適的迭代因子及對合矩陣。

共軛梯度法是求解非線性無約束優化問題的常用方法，穩定高效的特點使其在優化領域得到了廣泛的發展。基于 RMIL（Rivaie Mustafa Ismail and Leong）共軛梯度方法［13］的思想并結合自適應步長調節策略［10，14，15］，本文提出了共軛梯度步長調節法 CGS（Conjugate gradient step length adjustment method）。多個數值算例表明，無論是凸函數還是凹函數，與其他求解方法相比，CGS方法更加高效且穩健。

2 概率約束評估的功能度量法

2.1 基于概率功能度量的概率結構優化設計

概率結構優化設計數學模型通常為［5］

式中d為設計變量向量，x為隨機變量向量，gj（d，x）為第j號功能函數，gj（d，x）≤0代表失效，失效概率定義為

式中 Fgj（＊，d）為功能函數gj（d，x）的累積分布函數，f（x）為x的聯合概率密度函數，Pj，t為給定的許可失效概率，與目標可靠指標βj，t的關系為Pj，t≈Φ（－βj，t）。

式（1）的概率約束可采用概率功能度量的方式表達，其形式為

式中 Gpj（d）為第j號功能函數gj（d，x）的概率功能度量。采用概率功能度量表達概率約束，即為功能度量法（PMA）。

2.2 概率功能度量求解

概率功能度量求解需要在標準正態空間中進行，若初始隨機向量不是標準正態分布向量，則需借助Rosenblatt／Nataf等變換［5］，即u＝T（x），相應功能函數的變換為g（d，x）＝g（d，T－1（u））＝G（d，u）。

概率功能度量的求解本質上為標準正態空間中的一個優化問題，其數學優化模型可表示為

式中u＊為在目標可靠指標球面上功能函數值最小的點，稱為最小功能目標點（MPTP）。從而概率功能度量Gp為MPTP處的功能函數值，也就是目標可靠指標球面上功能函數的最小值。

改進均值法（AMV）［6］常用于確定 MPTP，其迭代格式為

式中uk代表第k次迭代點，n（uk）為點uk處的最速下降方向。AMV的迭代過程如圖1所示。

從式（6）可以看出，AMV方法中通過uk處的最速下降方向來實現MPTP的更新。對于非線性程度較低的功能函數以及凸函數，AMV方法具有良好的收斂性能；然而，對于非線性程度較高的功能函數，AMV迭代格式可能出現周期振蕩或混沌等不收斂現象［10］。針對AMV方法求解存在的弊端，HMV，CC 和 MCC 等方法相繼提出［8，11，12］，雖然改善了收斂性能，但其操作復雜，難以滿足實際應用的需求。

3 共軛梯度步長調節法

20世紀50年代初，為求解線性方程組提出了最早的共軛梯度法［16］。隨后，Fletcher等［17］將共軛梯度法應用于求解非線性無約束優化問題，稱之為Fletcher-Reeves（FR）共軛梯度法；Keshtegar［14，15］將共軛搜索方向用于可靠指標的求解，相比最速下降搜索方向，采用共軛梯度方向更具有魯棒性。基于共軛搜索方向的思想，并結合步長調節策略，本文提出了共軛梯度步長調節法（CGS）以求解概率功能度量。

3.1 基本的迭代格式

CGS方法的迭代過程如圖1所示。可以看出，uk為第k次迭代點，Sk為功能函數在點uk處的共軛梯度方向。在點uk處沿著Sk方向移動距離λk（λk＞0），可得到，

所以CGS迭代格式為

CGS的共軛梯度方向Sk為

式中 共軛梯度系數θk用于控制前一次共軛搜索方向Sk－1對當前搜索方向Sk的影響，如圖2所示，通常不同共軛梯度算法將定義不同的θk。若θk的值較大，會使新的共軛搜索方向太接近前一次搜索方向，影響收斂性［18］。因此，必須設定θk的合理取值以保證算法的收斂性和高效性。本文基于（Rivaie Mustafa Ismail and Leong，RMIL［13］）確定共軛梯度系數θk，其中θkRMIL定義為

圖1 AMV和CGS方法的迭代Fig.1 Iterative procedure of AMV and CGS

圖2 CGS方法迭代參數Fig.2 Iterative parameters of CGS

從圖2可以看出，負梯度方向（－Ug（uk））和共軛搜索方向Sk之間的夾角為φk，φk的值主要由共軛梯度系數θk確定。為使迭代逐漸收斂，每次迭代時應滿足0≤φk＜π／2［14］。

經算例分析（見3.4節中算例1）發現，當共軛梯度系數θk采用時，迭代過程中出現θk大幅波動的情況。當θk＜0時，新的共軛搜索方向Sk與前一次搜索方向Sk－1幾乎反向，會延緩迭代進程，降低收斂效率；當θk太大時，則可能導致φk≥π／2，影響算法的收斂性。因此，為了避免θk大幅波動影響收斂效率，可以假定0≤θk≤‖Ug（uk）‖／‖Sk－1‖。結合，本文共軛梯度系數θk的取值為

3.4 節算例1表明，采用如上共軛梯度系數時，θk的取值相對穩定，保證了收斂穩定性，提高了效率。

3.2 確定步長的方法

除共軛搜索方向，步長調節是本文所提方法的另一重要措施。步長λ的選擇通常與功能函數的非線性程度有關。然而，在實際問題求解中，步長λ的取值與函數非線性程度的量化關系無法確立，因此通常先假定一個較大的步長λ，使迭代點快速向最優點附近靠近，隨后逐漸調小λ以精確逼近最優點。Keshtegar［15］在可靠指標求解中假定初始步長也即最大步長為

式中 常數M滿足10≤M≤100。Keshtegar提出的步長準則看似解決了初始步長選擇的難題，但實際求解過程中發現，根據Keshtegar步長準則給出的初始步長經常因為太小造成收斂緩慢的局面。因此，通過限定初始步長下限值并結合Keshtegar步長準則，最終得到初始步長取值為

同時，取λ1＝λ0。而在后續迭代過程中，若某一步出現了 ‖uk＋1－uk‖＞‖uk－uk－1‖，將步長λk調整為λk＋1＝ckλk，ck稱為步長調節因子［14］，ck根據每次迭代過程中梯度和共軛梯度的情況進行取值，取值如下。

綜上所述，當k＞0時，λk的取值如下：

3.3 迭代過程

CGS方法迭代過程如下。

（1）將原始隨機變量x轉化為標準正態變量u，設定常數M，迭代停止準則ε＝0.001。

（2）令k＝0，設定初始迭代點u0，算例中沒有明確說明則取u0為坐標原點，根據式（14）計算初始步長λ0，根據式（9，10）計算u1，取λ1＝λ0，k＝k＋1，式（14）的常數M本文均取35。

（3）采用有限差分法計算功能函數在點uk處的梯度Ug（uk）。

（4）根據式（12）確定θk，根據式（10）確定Sk，結合式（9）計算uk＋1，根據式（15，16）調整步長。

（5）若‖uk＋1－uk‖＜ε，停止迭代，uk＋1即為u＊，g（uk＋1）為所求Gp；否則，令k＝k＋1，返回步驟（3）。

3.4 算 例

通過多個算例驗證CGS方法求解概率功能度量時的有效性和穩健性，并與其他常見方法如AMV，HMV，CC和MCC進行比較，考證CGS方法的計算效率。參與對比的這些方法均在Matlab R2014a軟件中實現，均采用同一收斂準則‖uk＋1－uk‖＜0.001，靈敏度分析均采用有限差分法。

算例1 考慮如下凹功能函數［12］，

圖3 功能函數示意圖（算例1）Fig.3 Performance function diagram （example 1）

標準正態空間中，功能函數和最小功能目標點u＊如圖3所示，可以看出，在u＊附近的功能函數為凹函數。表1～表11中，Iters為迭代次數，NFE為功能函數計算次數。表1對比了分別采用和本文所提出的兩種共軛梯度系數的計算結果，初始迭代點分別取u0＝（0.0，0.0）和u0＝（3.0，3.0）。可以看出，采用兩種共軛搜索方向雖然均能得到收斂解，但采用時，效率明顯高于。表2給出了初始迭代點為u0＝（0.0，0.0）時，兩種不同共軛梯度系數的迭代歷程。可以看出，整個迭代過程中，共軛梯度系數的取值發生大幅波動，因而共軛搜索方向的波動也較大，影響計算效率；而的取值相對平穩，避免了連續迭代過程中共軛搜索方向大幅波動的情況；表1結果也驗證了采用的高效性和穩健性。

表3給出了初始迭代點為u0＝（0.0，0.0）時，CGS，AMV，HMV和CC等常見方法的計算結果。文獻［12］給出MCC方法計算結果Gp＝－21.672，但經計算發現該答案錯誤。可以看出，當采用CGS方法時，計算效率很高，僅經過12次迭代和37次函數調用就得到了收斂解，概率功能度量值為Gp＝－76.035。因此，針對非線性較高的凹函數，CGS方法具有良好的收斂性。

算例2 考慮凹功能函數［12］：

表4給出了各方法的計算結果，除AMV方法外，其余方法都能收斂，且與文獻［12］給出的求解結果Gp＝－2.2293一致。可以看出，本文提出的CGS方法不僅能夠得到收斂解，其計算效率也明顯高于其他方法，僅需11次迭代和34次功能函數調用，再次驗證了CGS方法求解凹函數的高效性。

CGS方法的迭代歷程列入表5。可以看出，通過步長準則自動選取的初始步長λ＝51.356，經兩次調節減小到3.340。由于CGS采用了自適應步長調節策略，因此無需估計功能函數非線性程度，也無需假定λ的合適初值，通過步長準則會自動選取初值λ，并在迭代過程中不斷調整，直至最終收斂。

算例3 考慮凸功能函數［7］：G（X）＝－exp（X1－7）－X2＋10

Xi～N（6.0，0.82），i＝1，2，βt＝3.0

文獻［7］給出了求解結果，最小功能目標點u＊＝（2.8963，0.7813），概率功能度量Gp＝－0.358。從表6可以看出，各種方法的求解結果與文獻所得結果均非常吻合。當采用AMV方法時，經過8次迭代和25次功能函數調用達到收斂，說明對于凸功能函數，AMV方法具有較好的收斂性能。同時，相對于其他方法，CGS方法的計算效率也較高，經過11次迭代和34次功能函數調用得到收斂解，即CGS方法對于凸函數的計算效率與AMV基本相當。

表1 采用不同共軛梯度系數的計算結果Tab.1 Results of different conjugate gradient coefficients

表2 不同共軛梯度系數的迭代歷程Tab.2 Iterative history of different conjugate gradient coefficients

表3 算例1中各方法的計算結果Tab.3 Results of various methods in example 1

表4 算例2中各方法的計算結果Tab.4 Results of various methods in example 2

算例4 考慮如下功能函數［9］：

該功能函數隨機變量較多，有7個隨機變量。表7給出了不同方法的計算結果，除了AMV方法以外，其他方法均能收斂，概率功能度量Gp＝0.075，與文獻［9］的結果一致。CGS方法經過13次迭代和105次函數調用，計算效率相對較高。

算例5 考慮如下高度非線性功能函數［11］：

表8給出了各常見方法的計算結果。采用初始迭代點u0＝（0.0，0.0）時，只有CGS方法收斂，且所得結果與文獻［11］結果吻合。若采用初始點u0＝（－2.0，－2.0），CGS方法僅經過15次迭代就收斂到最小功能目標點，與CC方法相比，CGS方法更加高效。由此可看出，CGS方法對初始點的依賴較小，且能顯著改善高度非線性功能函數概率功能度量求解的收斂性能。

算例6 如圖4所示，三跨十二層框架結構［10］由同種材料制成，楊氏模量為20GPa，P為水平荷載。梁柱單元截面面積分為5組，圖中數字編號表示不同單元截面面積Ai（i＝1，2，…，5）的組別，截面慣性矩為Ii＝αiA2i（i＝1，2，…，5），系數αi的值列入表9。Ai（i＝1，2，…，5）與P是相互獨立的隨機變量，其統計信息列入表9。當結構最大位移超過0.096m時認為結構失效，目標可靠指標為βt＝3.0。

文獻［10］給出概率功能度量為Gp＝－0.0134，最小功能目標點u＊＝（－0.6923，－0.2559，－1.0827，－1.1523，0.0002，2.4403）。從表10可以看出，僅AMV方法和CGS方法能夠收斂，且計算結果與文獻結果相近，其中CGS方法調用函數次數較少，效率較高。

表5 算例2中CGS方法的迭代歷程Tab.5 Iterative history of CGS in example 2

表6 算例3中各方法的計算結果Tab.6 Results of various methods in example 3

表7 算例4中各方法的計算結果Tab.7 Results of various methods in example 4

表8 算例5中各方法的計算結果Tab.8 Results of various methods in example 5

算例7 圖5所示十桿桁架，材料楊氏模量為1.0E7psi，施加兩個垂直向下荷載P1＝P2＝1.0E5lb，各桿單元橫截面面積是相互獨立的正態分布隨機變量，均值為10in2，變異系數為0.05。結構功能函數定義為桁架的最大豎向位移不超過4in，目標可靠指標βt＝3.0。

圖4 十二層框架結構Fig.4 Twelve-story frame structure

圖5 十桿桁架Fig.5 Ten-bar truss

表9 算例6中隨機變量信息Tab.9 Information of random variables in example 6

表10 算例6中各方法的計算結果Tab.10 Results of various methods in example 6

表11 算例7中各方法的計算結果Tab.11 Results of various methods in example 7

從表11可以看出，僅AMV方法和CGS方法能夠收斂，CGS方法經18次迭代收斂，其效率遠遠高于AMV方法，再次驗證了CGS方法的高效性。

4 結 論

針對AMV方法求解高度非線性功能函數的概率功能度量難以收斂的情況，本文提出了共軛梯度步長調節法（CGS），引入新的共軛梯度搜索方向和自適應步長調節策略，以求解概率功能度量。新的共軛搜索方向保證了收斂穩定性；自適應步長調節策略無需了解功能函數凹凸性及非線性程度等先驗信息，無需尋找步長λ的合適取值。通過限定步長準則自動選取初始步長λ，并不斷調節λ，直至最終收斂。多個數值算例表明，相比于AMV和HMV等方法，本文提出的CGS方法不僅顯著地改善了高度非線性功能函數以及凹功能函數概率功能度量求解的收斂性能，同時對于凸功能函數也保持了較高的計算效率。