基于敏度分析的框架抗震優(yōu)化

劉齊茂，燕柳斌

（1.廣西工學(xué)院 土木建筑系，廣西 柳州 545006；2.廣西大學(xué) 土木建筑學(xué)院，南寧 530004）

結(jié)構(gòu)的抗震優(yōu)化設(shè)計屬于結(jié)構(gòu)動力響應(yīng)優(yōu)化設(shè)計的范疇，具有計算量巨大、優(yōu)化設(shè)計問題建模非常復(fù)雜等特點(diǎn)，因此該領(lǐng)域的研究進(jìn)展非常緩慢。根據(jù)各國規(guī)范的規(guī)定、震害經(jīng)驗(yàn)和實(shí)驗(yàn)研究結(jié)果及工程實(shí)例分析，在地震作用下采用多高層建筑結(jié)構(gòu)層間相對位移作為衡量結(jié)構(gòu)變形能力從而判斷結(jié)構(gòu)是否滿足建筑功能要求的指標(biāo)是合理的。因此Zou等［1］建立以結(jié)構(gòu)質(zhì)量最小化為目標(biāo)，同時滿足層間相對位移和設(shè)計變量約束的優(yōu)化數(shù)學(xué)模型，分別以振型分解反應(yīng)譜法和時程分析法為基礎(chǔ)推導(dǎo)出抗震優(yōu)化設(shè)計準(zhǔn)則，根據(jù)優(yōu)化設(shè)計準(zhǔn)則提出優(yōu)化設(shè)計方法，是一種有效和實(shí)用的抗震優(yōu)化設(shè)計方法，但該方法也存在不足，如只要求地震響應(yīng)過程中層間相對位移在某些時點(diǎn)（加速度峰值及鄰近的一些時點(diǎn)）滿足約束條件，但實(shí)際上，當(dāng)在優(yōu)化過程中出現(xiàn)比較柔的結(jié)構(gòu)時，層間相對位移響應(yīng)的最大值并不一定出現(xiàn)在加速度峰值及鄰近的一些時點(diǎn)，故該方法得到的優(yōu)化設(shè)計在地震響應(yīng)過程中不一定能滿足約束條件，可能得到不可行的解。Lagarosa等［2］以時程分析法為基礎(chǔ)，并考慮材料的非線性，采用進(jìn)化算法求解了以結(jié)構(gòu)質(zhì)量最小化為目標(biāo)，同時滿足層間相對位移和設(shè)計變量約束的優(yōu)化數(shù)學(xué)模型，但該方法由于采用隨機(jī)搜索方法，存在計算工作量過大的缺點(diǎn)。黃冀卓等［3］以鋼框架質(zhì)量最輕和結(jié)構(gòu)總動應(yīng)變能最小為目標(biāo)，基于相關(guān)的設(shè)計規(guī)范，給出了抗震鋼框架多目標(biāo)優(yōu)化問題的一種合理提法，采用Pareto遺傳算法求解抗震鋼框架多目標(biāo)優(yōu)化問題，得到多目標(biāo)優(yōu)化問題的妥協(xié)解。結(jié)構(gòu)抗震優(yōu)化設(shè)計屬于結(jié)構(gòu)動力優(yōu)化設(shè)計領(lǐng)域，結(jié)構(gòu)動力優(yōu)化設(shè)計方法從數(shù)學(xué)的角度看可分為零階優(yōu)化方法（非梯度類算法）［4］、一階優(yōu)化方法［5］和二階優(yōu)化方法。零階優(yōu)化方法不需要計算結(jié)構(gòu)動力響應(yīng)對設(shè)計變量的一階導(dǎo)數(shù)，如粒子群方法［6］、遺傳算法［7-8］和模擬退火算法［9］。一階優(yōu)化方法需要計算結(jié)構(gòu)動力響應(yīng)對設(shè)計變量的一階導(dǎo)數(shù)如共軛梯度法［10］，一階優(yōu)化方法關(guān)鍵是計算結(jié)構(gòu)動力響應(yīng)對設(shè)計變量的一階導(dǎo)數(shù)（也稱敏度分析），國內(nèi)外的學(xué)者已經(jīng)發(fā)展了多種敏度分析方法［11-12］。二階優(yōu)化方法不僅需要計算結(jié)構(gòu)動力響應(yīng)對設(shè)計變量的一階導(dǎo)數(shù)，而且需要計算結(jié)構(gòu)動力響應(yīng)對設(shè)計變量的二階導(dǎo)數(shù)（也稱Hessian矩陣分析），如牛頓類算法［13］，結(jié)構(gòu)動力響應(yīng)對設(shè)計變量的二階導(dǎo)數(shù)的計算難度非常大，因此在國內(nèi)外很少見到報道。通常，二階優(yōu)化方法的計算效率高于一階優(yōu)化方法，一階優(yōu)化方法的計算效率高于零階設(shè)計方法。

本文的主要工作是Newmark-β法的基礎(chǔ)上推導(dǎo)了層間相對位移對設(shè)計變量一階導(dǎo)數(shù)和二階導(dǎo)數(shù)的計算公式，建立結(jié)構(gòu)抗震優(yōu)化設(shè)計的數(shù)學(xué)模型，構(gòu)造求解結(jié)構(gòu)抗震優(yōu)化設(shè)計問題二階優(yōu)化算法。最后演示了一個三層兩跨的平面框架進(jìn)行抗震優(yōu)化設(shè)計，優(yōu)化結(jié)果表明本文的抗震優(yōu)化設(shè)計方法能獲得結(jié)構(gòu)的最優(yōu)設(shè)計。

1 優(yōu)化數(shù)學(xué)模型

框架具有i=1，2，3，…，N個單元組，單元的橫截面為矩形，如圖1所示，橫截面的寬度（be）和高度（he）為設(shè)計變量。

為了討論的方便，設(shè)計變量向量定義為：

框架抗震優(yōu)化設(shè)計的目標(biāo)函數(shù)為結(jié)構(gòu)的質(zhì)量，由設(shè)計變量表示為：

其中：ρ為材料的密度；li為第i個單元組中單元的長度；ni為第i個單元組所含單元的數(shù)目。

圖1 平面梁單元Fig.1 Plane beam element

對高層建筑結(jié)構(gòu)而言，結(jié)構(gòu)的側(cè)移通常成為設(shè)計的控制指標(biāo)。高層建筑結(jié)構(gòu)抗震設(shè)計必須限制結(jié)構(gòu)的層間相對位移，因?yàn)檫^大的層間相對位移會使主體結(jié)構(gòu)出現(xiàn)裂縫甚至破損，使結(jié)構(gòu)產(chǎn)生附加內(nèi)力，嚴(yán)重使會加速結(jié)構(gòu)的倒塌，使填充墻及建筑裝修出現(xiàn)裂縫或破損，使電梯軌道變形過大，使人不舒服，影響使用。《建筑抗震設(shè)計規(guī)范》［14］規(guī)定，以彈性方法計算的層間相對位移應(yīng)滿足以下條件：

式中，ΔuJ（t）為在t時刻第J層的彈性層間相對位移；uJ（t）和uJ-1（t）分別為在t時刻相鄰的第J和J-1樓層的水平位移；HJ為第J層的高度；θJ為抗震設(shè)計規(guī)范規(guī)定的第J層層間轉(zhuǎn)角；TT為地震的持續(xù)時間。

在抗震優(yōu)化設(shè)計中，除了考慮層間相對位移的約束外，構(gòu)件的強(qiáng)度和延性約束也是非常重要的，但由于問題本身很復(fù)雜，在本文中只考慮層間相對位移的約束（假設(shè)結(jié)構(gòu)的側(cè)移是結(jié)構(gòu)設(shè)計的控制指標(biāo)），但文中的優(yōu)化設(shè)計方法將可推廣到考慮構(gòu)件的強(qiáng)度和延性約束、同時考慮受到靜力和動力約束的情況，抗震優(yōu)化設(shè)計問題通過式（1）～式（4）描述為：

本文將采用高效的牛頓類優(yōu)化算法求解優(yōu)化數(shù)學(xué)模型式（5），優(yōu)化數(shù)學(xué)模型可通過積分型罰函數(shù)轉(zhuǎn)變?yōu)橐幌盗袩o約束優(yōu)化設(shè)計問題，然后計算罰函數(shù)對設(shè)計變量的梯度和海森矩陣，從而可構(gòu)造高效的牛頓類優(yōu)化算法。由式（5）可見，要計算罰函數(shù)對設(shè)計變量的梯度和海森矩陣，必須計算層間相對位移和質(zhì)量對設(shè)計變量的一階和二階導(dǎo)數(shù)。質(zhì)量對設(shè)計變量的一階和二階導(dǎo)數(shù)的計算是很容易解決的問題，因此層間相對位移對設(shè)計變量的一階和二階導(dǎo)數(shù)的計算是解決問題的關(guān)鍵。

2 一、二階導(dǎo)數(shù)計算

2.1 層間位移及其對設(shè)計變量的一、二階導(dǎo)數(shù)

2.1.1 層間相對位移

初始條件：

其中：M、C、K分別為結(jié)構(gòu)的質(zhì)量矩陣、阻尼矩陣和剛度矩陣；x（t）、（t）、（t）分別為相對于地面的結(jié)點(diǎn)位移、速度和加速度；Eu為一向量，在節(jié)點(diǎn)線位移u方向、v方向和節(jié)點(diǎn)角位移θ方向上分別為1、0、0。

結(jié)構(gòu)的阻尼矩陣采用瑞利阻尼如下：

其中，

式中，ω1、ω2分別為結(jié)構(gòu)的第一、二階固有頻率；ζ1、ζ2為阻尼比，本文取 ζ1=ζ2=0.05。

將采用Newmark-β法求解（6）、式（7），框架在地震作用下層間相對位移響應(yīng)的計算步驟如下：

（1）初始計算

其中，積分參數(shù)δ≥0.5，β=0.25（0.5+δ）2；⑤ 計算有效的剛度矩陣K*：

（2）計算t+Δt時刻的響應(yīng)

① 計算t+Δt時刻有效的荷載向量F*：

② 計算t+Δt時刻的節(jié)點(diǎn)位移x（t+Δt）：

③ 計算t+Δt時刻的節(jié)點(diǎn)加速度（t+Δt）：

⑤ 用式（4）計算t+Δt時刻的層間相對位移ΔuJ（t+Δt）；

（3）重復(fù)（2）步驟①～⑤即可獲得各個時刻的層間相對位移。

2.1.2層間位移對設(shè)計變量的一階和二階導(dǎo)數(shù)

式（13）兩邊分別對設(shè)計變量di求導(dǎo)一次，得：

其中，

式（17）和式（18）分別通過式（11）和式（12）對設(shè)計變量di求導(dǎo)一次獲得。

通過式（16）計算t+Δt時刻位移對設(shè)計變量di的一階導(dǎo)數(shù)后，式（14）兩邊分別對設(shè)計變量di求導(dǎo)一次，得：

通過式（19）計算t+Δt時刻加速度對設(shè)計變量di的一階導(dǎo)數(shù)后，式（15）兩邊分別對設(shè)計變量di求導(dǎo)一次，得：

式（4）兩邊分別對設(shè)計變量di求導(dǎo)一次，得到層間相對位移對設(shè)計變量的一階導(dǎo)數(shù)：

式（16）兩邊分別對設(shè)計變量dj進(jìn)一步求導(dǎo)，得：

其中：

式（23）和式（24）分別由式（17）和式（18）對設(shè)計變量di求導(dǎo)一次獲得。

由式（22）求得?2x（t+ Δt）/（?di?dj）后，式（19）兩邊分別對設(shè)計變量dj求導(dǎo)，得：

式（21）兩邊分別對設(shè)計變量dj求導(dǎo)，得：

（1）初始計算

其中，積分參數(shù)δ≥0.5，β=0.25（0.5+δ）2；

⑥ 由式（11）計算K*。

（2）計算t+Δt時刻的響應(yīng)

① 由式（13）計算t+Δt時刻的節(jié)點(diǎn)位移x（t+Δt）；

② 由式（14）計算t+Δt時刻 的節(jié)點(diǎn)加速度（t+Δt）；

③ 由式（15）計算（t+Δt）時刻的節(jié)點(diǎn)速度（t+Δt）；

④由式（4）計算t+Δt時刻的層間相對位移ΔuJ（t+Δt）；

⑤ 由式（16）計算t+Δt時刻的節(jié)點(diǎn)位移對設(shè)計變量的一階導(dǎo)數(shù)?x（t+Δt）/?di；

⑦ 由式（20）計算t+Δt時刻的節(jié)點(diǎn)速度對設(shè)計變量的一階導(dǎo)數(shù)（t+Δt）/?di；

⑧ 由式（21）計算t+Δt時刻的層間相對位移對設(shè)計變量的一階導(dǎo)數(shù)?ΔuJ（t+Δt）/?di；

⑨ 由式（22）計算t+Δt時刻的節(jié)點(diǎn)位移對設(shè)計變量的二階導(dǎo)數(shù)?2x（t+ Δt）/（?di?dj）；

⑩ 由式（25）計算t+Δt時刻的節(jié)點(diǎn)加速度對設(shè)計變量的二階導(dǎo)數(shù)?2（t+Δt）/（?di?dj）；

(11) 由式（26）計算t+Δt時刻的節(jié)點(diǎn)速度對設(shè)計變量的二階導(dǎo)數(shù)?2（t+Δt）/（?di?dj）；

(12) 由式（27）計算t+Δt時刻的層間相對位移對設(shè)計變量的二階導(dǎo)數(shù)?2ΔuJ（t+Δt）/（?di?dj）。

（3）重復(fù)執(zhí)行①～(12)步即可求得地震過程中各個時刻的響應(yīng)。

2.2 結(jié)構(gòu)質(zhì)量對設(shè)計變量的一階和二階導(dǎo)數(shù)

式（2）兩邊分別對設(shè)計變量求導(dǎo)一次得結(jié)構(gòu)質(zhì)量對設(shè)計變量的一階導(dǎo)數(shù)如下：

式（28）兩邊分別對設(shè)計變量進(jìn)一步求導(dǎo)得結(jié)構(gòu)質(zhì)量對設(shè)計變量的二階導(dǎo)數(shù)如下：

3 優(yōu)化數(shù)學(xué)模型的求解

將數(shù)學(xué)模型（5）的約束條件規(guī)則化，得到等效的數(shù)學(xué)模型如下：

3.1 優(yōu)化設(shè)計問題的轉(zhuǎn)變

不等式優(yōu)化數(shù)學(xué)模型（30）可通過罰函數(shù)轉(zhuǎn)變?yōu)橐幌盗袩o約束的積分形式的優(yōu)化數(shù)學(xué)模型，為了求導(dǎo)的方便本文采用內(nèi)點(diǎn)罰函數(shù)，如下：

其中：罰系數(shù)rk是一遞減數(shù)列，當(dāng)rk→0，罰函數(shù)P（d，rk）的最小值逼近約束優(yōu)化問題的最小值。因此，求解不等式約束優(yōu)化問題（30）轉(zhuǎn)變?yōu)榍蠼庖幌盗袩o約束優(yōu)化問題如下：

罰系數(shù)的初始值r1由下式計算：

其中：d0為初始設(shè)計點(diǎn)，p0=1～50。在本文中p0=10。

罰系數(shù)rk根據(jù)下式遞減：

其中：c=10～50，本文取c=10；k為搜索過程中使用罰函數(shù)的個數(shù)。

3.2 罰函數(shù)的梯度和海森矩陣

罰函數(shù)的梯度和海森矩陣即罰函數(shù)對設(shè)計變量的一階和二階導(dǎo)數(shù)，式（31）兩邊分別對設(shè)計變量di求導(dǎo)一次得罰函數(shù)對設(shè)計變量一階導(dǎo)數(shù)如下：

式（35）兩邊對設(shè)計變量dj進(jìn)一步求導(dǎo)得罰函數(shù)對設(shè)計變量的二階導(dǎo)數(shù)如下：

通過式（35）和式（36）計算罰函數(shù)對設(shè)計變量的一階和二階導(dǎo)數(shù)后即可獲得罰函數(shù)的梯度和海森矩陣。

3.3 優(yōu)化設(shè)計問題的求解

利用梯度和海森矩陣信息，構(gòu)造求解不等式約束優(yōu)化數(shù)學(xué)模型（30）的優(yōu)化設(shè)計算法，其步驟如下：

第1步 選擇可行的初始設(shè)計點(diǎn)d0，由式（33）計算罰系數(shù)的初始值，定義收斂容差ε1=10，令k=1；

第2步 從設(shè)計點(diǎn)dk-1出發(fā)，采用 Marquardt方法［15］求解優(yōu)化數(shù)學(xué)模型（32）獲得最優(yōu)設(shè)計dk；用Marquardt方法求解優(yōu)化數(shù)學(xué)模型（32）的步驟：

② 令i=0，λ（0）=105；

⑤i≥MI?是，執(zhí)行第(11)步；否，執(zhí)行第⑥步；

⑥計算

⑩ 令λ（i）=2λ（i），返回第⑥步；

第3 步P（dk-1，rk）-P（dk，rk）≤ε1?是，dk為最優(yōu)設(shè)計，輸出結(jié)果，停止計算；否，執(zhí)行第④步；

第4步 由式（34）計算罰系數(shù)rk+1，k=k+1，回到第②步。

4 數(shù)值算例

一個三層兩跨的平面框架如圖2所示，假設(shè)每個樓層的質(zhì)量為30 t，框架材料的彈性模量為20 GPa，阻尼比為5%，梁和柱的截面尺寸空間見圖2的上方，層間轉(zhuǎn)角限值為1/450，則第一層層間相對位移的限值為0.012 m，第二、三層的層間相對位移限值為0.008 m。

支座1，5，9受到El Centro水平地震加速度時程的激勵，El Centro水平加速度峰值為0.349 g，記錄時長為 8 s，步長 Δt=0.02 s。

從初始可行設(shè)計出發(fā)，通過4.3節(jié)的優(yōu)化設(shè)計算法獲得一系列優(yōu)化設(shè)計如表1所示。

由表1和表2可見，隨著罰系數(shù)逐步趨向于0，目標(biāo)函數(shù)即結(jié)構(gòu)的質(zhì)量越來越小，優(yōu)化設(shè)計逐步趨向于最優(yōu)設(shè)計（r4=0.327和r5=0.032 7獲得的優(yōu)化設(shè)計），收斂時前后兩次最優(yōu)設(shè)計的質(zhì)量不變，說明收斂精度非常高。初始可行設(shè)計的質(zhì)量為156 735 kg，最優(yōu)設(shè)計的質(zhì)量為57 043 kg，結(jié)構(gòu)的質(zhì)量減輕了99 692 kg，最優(yōu)設(shè)計的質(zhì)量僅為初始可行設(shè)計的36.39%，優(yōu)化效果明顯。在優(yōu)化過程中僅僅用了5個罰系數(shù)，說明優(yōu)化算法的效率非常高。由表2還可看出，即使僅使用1個罰系數(shù)（r1=327），得到的優(yōu)化設(shè)計也能得到極大的改善（質(zhì)量由156 735 kg下降為82 876 kg，優(yōu)化設(shè)計的質(zhì)量為初始可行設(shè)計的52.88%）。

圖2 三層兩跨的平面框架Fig.2 A three-storey，two-bay planar frame

表1 初始可行設(shè)計和優(yōu)化設(shè)計（單位：mm）Tab.1 Initial feasible design and optimal design（Units：mm）

表2 初始可行設(shè)計和優(yōu)化設(shè)計的質(zhì)量（單位：kg）Tab.2 Weight of the initial feasible design and optimal design（Units：kg）

初始可行設(shè)計和一系列優(yōu)化設(shè)計的質(zhì)量如表2所示。

5 結(jié)論

（1）提出一種計算層間相對位移對設(shè)計變量的一階導(dǎo)數(shù)和二階導(dǎo)數(shù)的方法，該方法只需要進(jìn)行一次動力響應(yīng)分析就可完成，具有非常高的計算效率。

（2）提出框架的抗震優(yōu)化設(shè)計數(shù)學(xué)模型適用于結(jié)構(gòu)的側(cè)移是結(jié)構(gòu)設(shè)計的控制指標(biāo)的情況（如高層建筑）。當(dāng)還要考慮其他設(shè)計指標(biāo)時，需要在優(yōu)化數(shù)學(xué)模型中增加響應(yīng)的約束條件。

（3）采用積分型罰函數(shù)將顯含時間參數(shù)的不等式約束優(yōu)化問題轉(zhuǎn)變?yōu)橐幌盗胁缓瑫r間參數(shù)的無約束優(yōu)化問題，并且提出積分型罰函數(shù)梯度和海森矩陣的計算方法。但由于采用內(nèi)點(diǎn)罰函數(shù)，因此初始設(shè)計必須為可行的設(shè)計。

（4）提出的優(yōu)化算法充分利用梯度和海森矩陣的信息確定搜索的方向和步長，不需要進(jìn)行一維搜索，下降速度快，具有良好的收斂性，能獲得最優(yōu)設(shè)計。算例表明該優(yōu)化算法具有非常高的效率，并且在優(yōu)化過程中能獲得一系列可行的優(yōu)化設(shè)計供工程設(shè)計人員選用。

（5）本文的抗震優(yōu)化設(shè)計方法雖然只輸入了一條地震波，但實(shí)際上可以同時輸入多條地震波，只需要在優(yōu)化數(shù)學(xué)模型中增加多個時程響應(yīng)約束即可，計算方法基本不變。

［1］ Zou X K，Chan C M.An optimal resizing technique for seismic drift design of concrete buildings subjected to response spectrum and time history loadings［J］.Computers and Structures，2005，83：1689-1704.

［21］Lagarosa N D，F(xiàn)ragiadakisa M，Papadrakakisa M.Structural optimization：a tool for evaluating seismic design procedures［J］.Engineering Structures，2006，28：1623-1633.

［3］黃冀卓，王 湛.基于遺傳算法的抗震鋼框架多目標(biāo)優(yōu)化設(shè)計［J］.力學(xué)學(xué)報，2007，39（3）：389-397.

HUANG Ji-zhuo，WANG Zhan.Multiobjective optimization design of seismic steel frames using genetic algorithm［J］.Chinese Journal of Theoretical and Applied Mechanics，2007，39（3）：389-397.

［4］Kocer F Y，Arora J S.Optimal design of h-frame transmission poles for earthquake loading［J］.Journal of Structure Engineering，1999，125：1299-1308.

［5］ Arora J S，Cardoso J E B.A design sensitivity analysis principle and its implementation into ADINA［J］.Computer and Structure，1989，32：691-705.

［6］ Pereza R E，Behdinan K.Particle swarm approach for structural design optimization［J］. Computers and Structures，2007，85：1579-1588.

［7］ Togan V，Daloglu A T.Optimization of 3d trusses with adaptive approach in genetic algorithms［J］.Engineering Structures，2006，28：1019-1027.

［8］ Venanzi I，Materazzi A L.Multi-objective optimization of wind-excited structures［J］.Engineering Structures，2007，29（6）：983-990.

［9］ Lamberti L.An efficient simulated annealing algorithm for design optimization of truss structures［J］.Computers and Structures，2008，86：1936-1953.

［10］Fletcher R，Reeves C M.Function minimization by conjugate gradients［J］.Computer Journal，1964，6：163-168.

［11］Van Keulen F，Haftka R T，Kim N H.Review of options for structural design sensitivity analysis.part 1：linear systems［J］. ComputerMethodsin Applied Mechanicsand Engineering，2005，194：3213-3243.

［12］ Bogomolni M，Kirsch U，Sheinman I.Efficient design sensitivities of structures subjected to dynamic loading［J］.International Journal of Solids and Structures，2006，43：5485-5500.

［13］ Kirsch U，Bogomolni M，Van Keulen F.Efficient finitedifference design-sensitivities［J］.AIAA Journal，2005，43：399-405.

［14］ Broydon C G.Quasi-newton methods and their application to function minimization［J］.Mathematics of Computation，1967，21：368-381.

［15］ GB50011-2001，建筑抗震設(shè)計規(guī)范［S］.北京：中國建筑工業(yè)出版社，2001.

［16］ Reklaitis G V，Ravindran A，Ragsdell K M.Engineering optimization methods and applications［M］.New York：John Wiley＆ Sons，1983.