基于中間變量的乘子法

李 亮，孫 秦

LI Liang,SUN Qin

西北工業大學 航空學院，西安 710072

School of Aviation,Northwestern Polytechnic University,Xi’an 710072,China

1 引言

乘子法是求解約束優化問題的一類重要優化算法，該法最早由Powell[1]和Hestenes[2]于1969年針對等式約束優化問題同時獨立提出，后又于1973年由Rockfellar[3]推廣到求解不等式約束優化問題。該法的基本思想是在原約束優化問題的拉格朗日函數的基礎上再加上適當的罰函數，從而將原問題轉化為一系列的無約束優化子問題，并通過求解序列子問題的解來逐次逼近原問題的解[4-5]。

結構優化是一類典型的不等式約束優化問題，可以很好地用乘子法求解。在結構優化中，作為約束函數的應力、應變、位移、屈曲特征值等響應與作為設計變量的結構尺寸、坐標位置之間呈復雜的非線性關系，且變量往往存在于分母中[6-7]。為了降低約束函數的非線性程度，在對近似概念的研究中，Schmit[8]于1974年引入了倒變量作為中間變量，后來Fleury[9]又于1986年引入了混合變量作為中間變量。考慮到結構優化中設計變量值一般為正，本文將倒變量和混合變量這兩種中間變量引入乘子法，提出基于中間變量的乘子法。中間變量可以降低無約束子問題的非線性程度，可以提高乘子法求解結構優化問題的效率和精度。

2 理論基礎

2.1 乘子法

一般的結構優化問題可以表示為如下數學形式：

其中，x為設計變量的集合，f(x)為目標函數，cj(x)為約束函數，和分別為設計變量xi在設計空間中的下界和上界。此外，f(x)和cj(x)都是連續可微的，且問題（1）的可行域非空。

根據乘子法，可以將問題（1）轉化為如下形式的增廣拉格朗日函數：

其中，λ=(λ1，λ2，…，λm)T為非負的拉格朗日乘子向量，σ為正的罰參數，這兩個量的初值由人為給定。邊界約束可以在線性搜索中得到滿足，因此未加入到增廣拉格朗日函數中。

由式（2）可在點xκ處構造如下無約束子問題：

有多種方法可用于求解子問題（3），本文采用擬牛頓法。將求得的解記作xκ+1，若xκ+1滿足如下收斂條件：

則將xκ+1作為原問題（1）的近似極小點。其中，ε是人為設定的收斂閾值，0＜ε＜1。否則，則按如下形式更新拉格朗日乘子向量λ和罰參數σ的值：

并在點xκ+1處繼續構造無約束子問題，依此迭代，直到收斂。

2.2 擬牛頓法

擬牛頓法[10]是求解如下形式的無約束優化問題的一種重要方法：

該方法的其基本思路是，在當前迭代點xk處構造如下二次近似子問題來近似原問題（7）：

再求解子問題（8）并用其解dk作為搜索方向，即有：

其中，gk是 f(x)在點 xk處的一階偏導數，Gk是 f(x)在點xk處的Hesse陣，Hk是Hesse陣逆的近似。Hk可以用校正公式求得，本文采用如下BFGS校正公式[11]：

其中，yk=gk-gk-1，sk=xk-xk-1。

再沿dk做線性搜索，求得步長αk，得到xk+1=xk+αkdk。本文采用Armijo線性搜索技術[12]，并對其進行了改進，以保證搜索一開始就在設計空間內進行，算法步驟如下：

步驟1 給定0＜β＜1，0＜μ＜0.5，令 m1:=0，m2=20。

步驟2 求 p(i)=xk(i)+βm1d(i)，i=1，2，…，n。

若線性搜索求得的xk+1滿足如下收斂條件：

則將xk+1作為問題（7）的近似解；否則，在xk+1處繼續構造搜索方向，并求步長，以得到下一個點，依此迭代，直到收斂。

3 基于中間變量的乘子法

3.1 引入中間變量

本文引入了兩種中間變量，即倒變量和混合變量。對于原變量x，若變量t滿足：

則稱其為x的倒變量；若變量t滿足：

則稱其為x的混合變量。

3.2 構造無約束子問題

引入了中間變量t，則原問題（1）可以化為：

其中，若ti是倒變量，滿足式（12）中的關系，則有：

若ti是混合變量，滿足式（13）中的關系，則有：

將問題（14）轉化為增廣拉格朗日函數形式為：

由式（17）可在點tκ處構造如下無約束子問題：

3.3 求解無約束子問題

在求解子問題（18）時，用擬牛頓法求搜索方向dk，用改進的Armijo線性搜索技術求搜索步長αk，以及用BFGS校正法構造Hesse陣逆的近似Hk，都是對中間變量t進行的。在求解中需要用到 fM(tk)和cMj(tk)的函數值和梯度值，求解方式為，若t是倒變量，滿足式（12）中的關系，則

若t是混合變量，滿足式（13）中的關系，則

將求得的子問題（18）的解記作 tκ+1。

需要注意的是，若采用的是混合變量時，由于不同迭代點處?L/?xi取值的正負情況并不一定相同，使得各迭代點處所用的混合變量的形式也會不同。因此，在用BFGS校正法構造Hesse陣逆的近似Hk時，為了獲得前后兩迭代點的梯度差，需要把兩點處所用的混合變量的形式進行統一。本文的做法是，按當前點處的?L/?xi取值的正負情況來確定當前點和上一迭代點處所用的混合變量的形式。

3.4 算法的迭代和收斂

若tκ+1滿足如下收斂條件：

則將tk+1作為原問題（14）的近似極小點；否則，則按如下形式更新拉格朗日乘子向量λ和罰參數σ的值：

并在點tκ+1處繼續構造無約束子問題，依此迭代，直到收斂。

4 數值算例

為了驗證本文提出的基于中間變量的乘子法的計算效率和收斂性能，用Fortran語言編寫了乘子法（MM）、基于倒變量的乘子法（RVMM）和基于混合變量的乘子法（MVMM）的計算程序，并用這三種方法分別求解了兩個有一定非線性程度的算例[13-15]。為了保證算法的可比性，在求解同一個算例時，三種方法都采用相同的初始值和收斂閾值。

4.1 兩桿桁架的形狀優化

考慮一個如圖1所示的兩桿桁架的形狀優化。該優化問題有兩個類型不同的設計變量，其中x1定義為兩桿的截面面積，x2定義為兩桿下端點距離的一半。設計目標是桿的重量最小。設計約束施加在每個桿的應力上。

圖1 兩桿桁架

該結構優化問題可表示為如下數學形式：

以 x=(1.5，0.5)為起始點，用三種方法求解該問題，結果見表1。

表1 兩桿桁架形狀優化結果

4.2 懸臂梁的重量優化

考慮如圖2所示的懸臂梁的重量優化。這個梁分為5段，且每段的截面都是一個有給定均勻厚度的空心四方體。該優化問題有5個設計變量，分別定義為每個梁段的高度。設計目標是梁的重量最小。設計約束施加在梁頂端的位移上。

圖2 懸臂梁

該結構優化問題可表示為如下數學形式：

以 xj=5.0(j=1，2，…，5)為起始點，用三種方法求解該問題，結果見表2。

表2 懸臂梁重量優化結果

4.3 討論

在兩個算例中，三種方法都達到了收斂，且求得的近似最小點和近似最小值都很接近，這表明基于中間變量的乘子法具有良好的收斂性。

在兩個算例中，RVMM和MVMM的迭代次數都少于MM的，MVMM的迭代次數都略少于RVMM的，這表明在求解結構優化問題時，基于中間變量的乘子法比一般的乘子法具有更快的收斂速度，而基于混合變量的乘子法的收斂速度比基于倒變量的乘子法略有提高。

在第一個算例中，RVMM和MVMM的迭代次數比MM的減少得較少，在第二個算例中則減少得較多。經分析，發現這個差別是由兩個算例中目標函數和約束函數的結構形式的不同引起的。從兩個問題的數學形式可見，與第二個算例相比，第一個算例中增廣函數中變量在分子上占的比重更大、次方數更高；與第一個算例相比，第二個算例中增廣函數中變量在分母上占的比重更大、次方數更高。事實上，若一個函數中變量在分母上占的比重越大、次方數越高，則引入中間變量代替原變量后，函數的非線性程度降低得也就越多。相反，若一個函數中變量在分子上占的比重越大、次方數越高，則引入中間變量后，函數的非線性程度降低得就沒那么明顯。如果一個函數中所有變量都只出現在分子上，則引入中間變量后，函數的非線性程度不但沒有降低，反而會增大。因此，在結構優化問題中，若增廣函數中變量在分子上占的比重較大、次方數比較高，則RVMM和MVMM的收斂速度與MM相比提高得就較少；而若增廣函數中變量在分母上占的比重較大、次方數比較高，則提高得就較大。

5 結論

為了更有效地用乘子法求解結構優化問題，本文將乘子法與中間變量的概念相結合，提出了基于中間變量的乘子法。引入兩種中間變量，即倒變量和混合變量，并詳述了構造基于增廣拉格朗日函數的無約束子問題，并用擬牛頓法、Armijo線性搜索法和BFGS校正法求解該子問題，以及實現拉格朗日乘子和罰參數的更新的算法流程。

兩個數值算例的測試表明，基于中間變量的乘子法具有良好的收斂性，并比一般的乘子法具有更快的收斂速度。基于混合變量的乘子法的收斂速度比基于倒變量的乘子法有提高，但提高不顯著。基于中間變量的乘子法的收斂速度，會受到優化問題中目標函數和約束函數的結構形式的影響。若增廣函數中變量在分子上占的比重較大、次方數比較高，收斂速度就低些，而若增廣函數中變量在分母上占的比重較大、次方數比較高，收斂速度就高些。

[1]Powell M J D.A method for nonlinear constraints in minimization problems[M].New York：Academic Press，1969：283-298.

[2]Hestenes M R.Multiplier and gradient methods[J].Journal of Optimization Theory and Applications，1969，4：303-320.

[3]Rockafellar R T.Augmented Lagrange multiplier functions and duality in nonconver programming[J].SIAM Journal on Control，1974，12：268-285.

[4]Bertsekas D P.Constrained optimization and lagrange multiplier methods[M].New York：Academic Press，1982：104-205.

[5]Nocedal J，Wright S J.Numerical optimization[M].New York：Springer-Verlag，1999：511-520.

[6]Venkayya V B.Structural optimization：a review and some recomm-endations[J].International Journal for Numerical Methods in Engineering，1978，13：203-228.

[7]Vanderplaats G N.Structural design optimization status and directon[J].AIAA Journal，1999，36：11-20.

[8]Schmit L A，Farshi B.Some approximation concepts for structural synthesis[J].AIAA Journal，1974，12：692-699.

[9]Fleury C.Structural optimization：a new dual method using mixed variables[J].International Journal for Numerical Methods in Engineer，1986，23：409-428.

[10]Wei Z X，Li G Y，Qi L Q.New quasi-Newton methods for unconstrained optimization[J].Applied Mathematics and Computation，2006，175：1156-1188.

[11]Yuan Y X.A modified BFGS algorithm for unconstrained optimization[J].IAM Journal of Numerical Analysis，1991，11：325-332.

[12]Shi Z J，Wang S Q.Modified nonmonotone Armijo line search for descent method[J].Numerical Algorithms，2011，57：1-25.

[13]Svanberg K.The method of moving asymptotes—a new method for structural optimization[J].International Journal for Numerical Methods in Engineer，1987，24：359-373.

[14]Yang D X，Yang P X.Numerical instabilities and convergence control for convex approximation methods[J].Nonlinear Dynamics，2010，61：605-622.

[15]Zhang W H，Fleury C.A modification of convex approximation methods for structural optimization[J].Computer&Structures，1997，64：89-95.