彭細(xì)榮 隋允康 葉紅玲 鐵 軍

* (湖南城市學(xué)院土木工程學(xué)院,湖南益陽 413000)

? (北京工業(yè)大學(xué)材料與制造學(xué)部,北京 100124)

** (天津財經(jīng)大學(xué)理工學(xué)院,天津 300222)

引言

連續(xù)體結(jié)構(gòu)應(yīng)力拓?fù)鋬?yōu)化由于面對局部性能約束,較之全局性能約束問題,具有更多求解困難.對此,諸多學(xué)者例如Le 等[1]總結(jié)了主要困難,一是應(yīng)力奇異性,二是應(yīng)力局部特性導(dǎo)致敏度計算量大,約束數(shù)目多.

本文在贊同這一看法的同時,認(rèn)為首先應(yīng)當(dāng)區(qū)分困難的主次,且不可漏掉相關(guān)的困難.一是主要困難:約束數(shù)目多,應(yīng)力的局部性特點,要求每個單元的約束都不可違背,于是導(dǎo)致約束總數(shù)劇增,以靜力問題為例,約束數(shù)目等于單元總數(shù)與載荷工況數(shù)目的乘積,從而敏度計算量難以接受,造成建立優(yōu)化模型的極大困難;二是現(xiàn)有的有些解法存在個別單元應(yīng)力超限的現(xiàn)象,這緣于應(yīng)力較強非線性特征,在形狀與拓?fù)渫蛔儏^(qū)域出現(xiàn)應(yīng)力集中,也必須解決;三是應(yīng)力奇異性的困難,根本原因在于,在迭代中傾向于被刪除的子區(qū)域或單元集,其中的拓?fù)渥兞亢蛻?yīng)力分別趨于0 的速度決定了是否出現(xiàn)應(yīng)力奇異性現(xiàn)象.旨在解決應(yīng)力奇異性問題,Cheng 和Guo[2]提出了ε-松馳法,Bruggi[3]提出了qp-松馳法等.

為了克服約束數(shù)目多的困難,以往的研究提出了分部、化整和集成3 種解法.

分部方法[4-8]即滿應(yīng)力方法,或應(yīng)力的零階近似處理.隋允康等[4-5]提出了基于ICM 方法的應(yīng)力零階近似處理的方法.程耿東等[6-7]提出了修改的滿應(yīng)力方法求解二維連續(xù)結(jié)構(gòu)的應(yīng)力拓?fù)鋬?yōu)化問題.Guan 等[8]用雙向進化結(jié)構(gòu)優(yōu)化方法(BESO)處理應(yīng)力約束問題,依據(jù)單元應(yīng)力值進行刪除或加入.應(yīng)力零階近似的方法避免了求應(yīng)力敏度,約束則轉(zhuǎn)化為設(shè)計變量動態(tài)下限.但Zhou 和Sigmund[9]指出,滿應(yīng)力方法僅對設(shè)計區(qū)域無幾何限制的拓?fù)鋬?yōu)化問題有效,對有幾何限制的問題,無法得到合理拓?fù)?

化整方法[10-12]將應(yīng)力約束集合成為全局化的約束.隋允康等[10-11]將局部應(yīng)力約束轉(zhuǎn)化為結(jié)構(gòu)應(yīng)變能約束,或按結(jié)構(gòu)畸變比能進行轉(zhuǎn)化處理[12].但化整解法存在個別單元應(yīng)力超限的現(xiàn)象,一些迭代的算法被提出[10-12],但是解決方案屬于準(zhǔn)則方法,不能令人滿意.

集成方法[1,13-25]采用集成函數(shù)如模函數(shù)(或稱為P 范數(shù))或 Kreisselmeier-Steinhauser (K-S)函數(shù)等將諸多的局部應(yīng)力約束集成為一個或幾個全局化的函數(shù).K-S 函數(shù)最初應(yīng)用于控制系統(tǒng)的設(shè)計[13],后來被應(yīng)用于工程結(jié)構(gòu)優(yōu)化中,將多目標(biāo)或大量約束集成為單一目標(biāo)或約束.Yang 和Chen[14]最先提出了用K-S 函數(shù)將眾多應(yīng)力集成為一個應(yīng)力函數(shù),進行結(jié)構(gòu)拓?fù)鋬?yōu)化設(shè)計.隋允康等[15]為解決多工況下應(yīng)力約束的連續(xù)體結(jié)構(gòu)拓?fù)鋬?yōu)化問題,引入了K-S 函數(shù)對應(yīng)力約束進行集成化處理,并進而提出拋物型凝聚函數(shù)對應(yīng)力約束進行集成化處理[16].Le 等[1]及Luo 等[17]基于SIMP 法用集成處理方式求解應(yīng)力拓?fù)鋬?yōu)化問題.Picelli 等[18]采用集成函數(shù)及水平集方法解決應(yīng)力拓?fù)鋬?yōu)化問題.Xia 等[19]提出了一種基于集成函數(shù)的BESO解法.王選等[20]在BESO 方法中采用K-S 函數(shù)集成應(yīng)力,用拉格朗日乘子法求解,高云凱等[21]使用K-S 函數(shù)集成結(jié)構(gòu)畸變比能.Wares和Hideyuki[22]對最大應(yīng)力極小化的結(jié)構(gòu)拓?fù)鋬?yōu)化問題,應(yīng)用K-S 函數(shù)對應(yīng)力集成,并給出了二階近似,應(yīng)用牛頓梯度法進行求解.但集成的方法亦難以保證應(yīng)力處處滿足,增大集成函數(shù)的參數(shù)可以更好逼近應(yīng)力最大值,但又易引起優(yōu)化求解中的數(shù)值困難.于是一些解決此問題的算法被提出.Andrew 等[23]以翼盒結(jié)構(gòu)質(zhì)量極小化應(yīng)力約束優(yōu)化問題為例,比較研究了K-S 函數(shù)及模函數(shù)等不同集成函數(shù)及其參數(shù)對優(yōu)化結(jié)果的影響,算例結(jié)果表明采用誘導(dǎo)指數(shù)集成函數(shù)并分為多個子域凝集的效果最好.Zhang等[24]對機翼結(jié)構(gòu)使用代理模型進行全局優(yōu)化,對比了K-S 函數(shù)對應(yīng)力約束的不同凝集處理策略,比較了固定K-S 函數(shù)參數(shù)及自適應(yīng)參數(shù)、全部結(jié)構(gòu)應(yīng)力約束集成及劃分多個子域應(yīng)力約束凝集,結(jié)果表明采用較大的固定K-S 函數(shù)參數(shù)及多個子域集成的策略效果最好.Graeme[25]對應(yīng)力拓?fù)鋬?yōu)化問題提出了基于K-S 函數(shù)的自適應(yīng)方法,可自動調(diào)整K-S 函數(shù)的參數(shù)及應(yīng)力約束集成的子域劃分.但是以上所提出的諸多解決方案均不能令人非常滿意.

克服約束數(shù)目多的困難其他處理策略還有Guo等[26]提出的主動集策略,Zhang 等[27-28]提出了結(jié)構(gòu)邊界曲率信息及應(yīng)力梯度信息加權(quán)的應(yīng)力集成策略,Xia 等[29]提出懲罰加權(quán)的應(yīng)力集成策略,Fernando等[30]基于增加拉格朗日法的應(yīng)力懲罰處理策略等.

為了徹底解決應(yīng)力超限現(xiàn)象和同時提高求解效率.對應(yīng)力、局部穩(wěn)定及疲勞壽命等局部性能結(jié)構(gòu)優(yōu)化問題,隋允康等[31-33]提出了互逆規(guī)劃理論,闡述了單目標(biāo)多約束模型(s 方模型)及多目標(biāo)單約束模型(m 方模型)之間存在互逆關(guān)系,并應(yīng)用互逆規(guī)劃理論解決結(jié)構(gòu)拓?fù)鋬?yōu)化建模中的合理性問題,進而基于互逆規(guī)劃的理論提出了交融優(yōu)化的3 種解法——分部-集成、化整-集成、集成-集成解法.

本文的工作針對應(yīng)力拓?fù)鋬?yōu)化的交融優(yōu)化解法,由于分部-集成解法的前半部分屬于準(zhǔn)則方法,不必再討論,集成-集成解法的研究將另文給出,本文的研究聚焦在化整交融(化整-集成)方法的研究上,并且同化整方法進行比較.

首先,敘述集成方法即s 方模型的化整解法,然后,對m 方模型集成解法,提出了乘子法、一階近似及二階近似等不同解法,最后,敘述了基于化整處理的的s 方模型和K-S 函數(shù)集成處理的m 方模型交替交融的迭代解法,亦即化整交融解法.用數(shù)值算例對比了不同m 方模型集成解法的效率,選擇了求解效率最佳的二階近似法對m 方模型進行求解,在此基礎(chǔ)上,對比了化整解法與化整交融解法的求解效果.

1 s 方模型化整解法[34-35]

應(yīng)力約束下最小化結(jié)構(gòu)重量(或體積)化整解法的優(yōu)化模型,屬于s 方模型,即

式中,t=[t1,t2,···,tN] 為拓?fù)湓O(shè)計向量,W(t) 為結(jié)構(gòu)重量,wi(ti) 為單元重量,ψil(t) 為i 單元在l 號載荷工況下的局部性能如單元應(yīng)力,N 為載荷工況總數(shù),L為載荷工況總數(shù).t 為防止結(jié)構(gòu)分析時產(chǎn)生奇異而設(shè)置的拓?fù)渥兞肯孪拗?

當(dāng)應(yīng)力約束化整為全局化應(yīng)力約束時,有

其中,eil為i 單元在l 載荷工況下的應(yīng)變能,Vs為結(jié)構(gòu)體積,E,μ 和 σ 分別為材料彈性模量、泊松比和許用應(yīng)力.

式(2)的左端為結(jié)構(gòu)應(yīng)變能,其中單元應(yīng)變能為

其中 σil,εil與 Pil為i 單元在l 載荷工況下的應(yīng)力、應(yīng)變與節(jié)點力向量,Ki與 Ωi分別為i 單元的剛度矩陣及子域.

將拓?fù)渥兞?ti代入單元剛度矩陣中,得到

再代入式(3)得

在第 ν 次迭代中,由靜定化假設(shè),該輪中的內(nèi)力不變,于是得應(yīng)變能函數(shù)為

將式(6)代入式(7)得

將式(8)代入式(2)得

其中Ve為單元體積,N(ν)為第 ν 次迭代的主動拓?fù)渥兞靠倲?shù),取 N(ν-1)近似替代,為結(jié)構(gòu)許用應(yīng)變能.

由此,s 方優(yōu)化模型式(1)中的局部應(yīng)力約束轉(zhuǎn)化為L 個全局化的應(yīng)變能約束,得到應(yīng)力全局化的s 方優(yōu)化模型

2 m 方集成模型的不同逼近解法

化整交融解法的第二步是多目標(biāo)的m 方模型求解算法.m 方模型是以結(jié)構(gòu)重量(或體積)為約束,單元最大局部性能如應(yīng)力極小化為目標(biāo)的優(yōu)化模型,即

式中符號意義與式(1)相同.

對式(12)中的目標(biāo)函數(shù),采用K-S 函數(shù)進行集成化處理,即為m 方模型的集成解法[33],得到集成轉(zhuǎn)化后的優(yōu)化模型

以下給出m 方集成化處理模型式(13)的3 種解法,為后文敘述的方便,每種解法在括號內(nèi)給出了其簡稱,分別如下:

(1) Lagrange 乘子解法(乘子法);

(2) Taylor 一階近似逼近的解法,采用移動漸近線法[36]求解(MMA 法);

(3) Taylor 二階近似逼近的解法,采用序列二次規(guī)劃法求解(SQP 法).

2.1 Lagrange 乘子解法

對于應(yīng)力約束,式(13)中

為了求出規(guī)劃(13)的解,先列出它的Lagrange 函數(shù)

由式(15)得鞍點條件

從式(16)得

從而有

定義

按式(20)代入式(17)得

由式(22)得

由式(23)代入式(20)得

其中 βi見式(21).注意,式(24)的 ti取因而式(24)左端 ti取,于是式(24)就構(gòu)成了一個迭代式.迭代求解中當(dāng) t(v+1)與 t(v)充分接近后收斂,每次小循環(huán)迭代中都用 0 ≤ti≤1 修改迭代解.

2.2 一階近似的MMA 解法

約束函數(shù)為線性函數(shù),其函數(shù)及一階導(dǎo)數(shù)為

很多學(xué)者都在研究高管薪酬激勵與技術(shù)創(chuàng)新之間的關(guān)系，其理論各不相同，但核心基本相似：為獲超額利潤與競爭優(yōu)勢，企業(yè)傾向于高管做技術(shù)創(chuàng)新，而高管進行技術(shù)創(chuàng)新的程度與企業(yè)的薪酬激勵相關(guān)，因此，提高高管薪酬是企業(yè)技術(shù)創(chuàng)新的基本手段。技術(shù)創(chuàng)新是企業(yè)發(fā)展的根本，技術(shù)創(chuàng)新需要資金等條件的保障。高管是企業(yè)發(fā)展的關(guān)鍵因素，讓高管利益與企業(yè)利益進行掛鉤，是企業(yè)進行技術(shù)創(chuàng)新突破的關(guān)鍵。通過薪酬激勵，促使高管引領(lǐng)企業(yè)技術(shù)創(chuàng)新，是保障企業(yè)科學(xué)發(fā)展的重要手段[4]。

移動漸近線法(MMA)是基于一階導(dǎo)數(shù)求解的算法[36],有了目標(biāo)函數(shù)及約束函數(shù)的原函數(shù)及其一階導(dǎo)數(shù)如式(26)～ 式(29)所示,利用MMA 求解器可直接對模型式(25)進行求解.

2.3 二階近似的序列二次規(guī)劃解法

優(yōu)化模型仍為式(13),轉(zhuǎn)換到x 空間里的式(25).對式(25)中目標(biāo)函數(shù)進行二階Taylor 近似,記原函數(shù)和一階偏導(dǎo)數(shù)分別為式(26) 及式(27) 所示.后面將采用可分離規(guī)劃解法迭代逼近式(25)的準(zhǔn)確解,為此將式(27)所含求和項中的xi暫時取為于是原函數(shù)的Hessian 矩陣中只有對角線上的二階導(dǎo)數(shù)為非零值,其為

其中xi(k)的初始值取為下標(biāo)(k)為小循環(huán)標(biāo)志,代表一次分析之后,通過二階近似迭代逼近求解;上標(biāo)(v)為大循環(huán)標(biāo)志,代表結(jié)構(gòu)分析的第v 次,亦即一次大迭代尋優(yōu)過程.

對原函數(shù)式(26)的二階近似式略去對優(yōu)化沒有影響的常數(shù)項,為

為了求出規(guī)劃式(33) 的解,列出其Lagrange函數(shù)

由式(34)得鞍點條件

由式(35)可求得主動變量為

解得

將式(40)代入式(38)得主動變量為

把式(41)由主動集延伸到非主動集,則有下式

式(42)的計算伴隨著主動集的更新,直到主動集不變?yōu)橹?

3 化整交融解法的計算流程

化整交融模型的解法是通過交替求解s 方模型及m 方模型,以達到更好收斂速度的目標(biāo).s 方為化整處理法建立的優(yōu)化模型,m 方為集成處理法建立的優(yōu)化模型.一次結(jié)構(gòu)分析的大循環(huán)內(nèi),先進行s 方化整模型式(11)的求解,再進行m 方模型公式(13)的求解,流程如圖1 所示.

圖1 化整與集成解法流程圖Fig.1 The flow chart of globalization-aggregation method

4 數(shù)值算例

4.1 算例1:m 方模型K-S 函數(shù)集成應(yīng)力拓?fù)鋬?yōu)化

如圖2 所示,基結(jié)構(gòu)平面體的尺寸為50 mm×100 mm×10 mm,采用 50×100 網(wǎng)格.材料彈性模量為E=210 GPa,泊松比 ν=0.3 .左邊界固定,右邊界中點作用集中載荷 F=5000 N 作用于上邊界中點,為避免應(yīng)力集中,分散作用在鄰近的3 個節(jié)點上.體積比約束為0.2%,收斂精度取0.01.

圖2 兩工況算例分析及優(yōu)化的模型Fig.2 Model for analysis and optimization of the example with two load cases

單獨的m 方模型存在優(yōu)化模型提法的不合理性[34-35,37].此算例的設(shè)計僅為了比較求解m 方集成模型式(25)的3 種不同求解方法——乘子法、MMA法及SQP 法的效果.

3 種求解方法得到的最優(yōu)拓?fù)淙鐖D3 所示,對應(yīng)的Mises 應(yīng)力分布如圖4 所示.迭代次數(shù)、體積比約束值、目標(biāo)函數(shù)應(yīng)力值及收斂精度對比如表1中所示.目標(biāo)及約束函數(shù)的迭代歷史曲線如圖5 所示.從結(jié)果對比可以看出,在0.01 的收斂精度下,乘子法及SQP 法均能收斂到清晰的拓?fù)鋱D形,兩種方法下的最大Mises 應(yīng)力接近(見表1),拓?fù)鋱D形也相似(見圖3(a)及圖3(c)),但SQP 的收斂速度更快.MMA 方法在0.01 的收斂精度下得不到清晰的拓?fù)鋱D形(見圖3(b)),相比于乘子法及SQP 法,其收斂時的Mises 應(yīng)力也非常大(見表1 及圖5(b)).當(dāng)將MMA法的收斂精度提高到0.001 時,經(jīng)過217 次迭代收斂,其拓?fù)鋱D形及Mises 應(yīng)力分布如圖6 所示,此時的最優(yōu)拓?fù)渑c乘子法及SQP 法得到的最優(yōu)拓?fù)涫窍嗨频?對比圖6(a) 與圖3(a) 及圖3(c)),Mises應(yīng)力為191.192 MPa,與乘子法及SQP 法的對應(yīng)值也是接近的(見表1).可見,MMA 法在提高一個量級收斂精度的情況下,才能得到與乘子法及SQP 法相同的結(jié)果.

圖3 算例1 最優(yōu)拓?fù)銯ig.3 Optimized topology of the Example 1

圖4 算例1 最優(yōu)拓?fù)涞腗ises 應(yīng)力分布Fig.4 Mises stress distribution of of the optimized topology for the Example 1

表1 m 方模型不同求解方法對比Table 1 Comparison of different methods to solve m-model

圖5 算例1 目標(biāo)及約束迭代歷史Fig.5 The iteration history of the objectives and constraints for the Example 1

圖6 算例1 收斂精度0.001 時MMA 法得到的最優(yōu)拓?fù)浼癕ises 應(yīng)力分布Fig.6 Optimized topology and Mises stress distribution for the Example 1 by MMA method with the convergence accuracy 0.001

從此算例對比中可知:3 種方法中SQP 法的收斂速度最快,乘子法與之接近.而MMA 法的收斂速度遠低于乘子法及SQP 法,其原因是MMA 方法是一階近似的求解方法.由此,在本文在化整-集成方法中,對m 方集成模型的求解選用SQP 解法.

4.2 算例2:s 方模型化整交融應(yīng)力拓?fù)鋬?yōu)化

此算例的條件與算例1 相同.但是求解的優(yōu)化模型為s 方模型,即應(yīng)力約束下極小化結(jié)構(gòu)體積的優(yōu)化模型.結(jié)構(gòu)允許應(yīng)力為200 MPa.分別用單獨的化整解法及化整-集成解法兩種方法進行求解.

單獨的化整解法經(jīng)過19 次收斂,最大Mises 應(yīng)力為201.743 MPa,略微超出了允許應(yīng)力200 MPa.化整-集成解法經(jīng)過17 次收斂,最大Mises 應(yīng)力為199.704 MPa.兩種方法得到的最優(yōu)拓?fù)淙鐖D7 所示,對應(yīng)的Mises 應(yīng)力如圖8 所示,目標(biāo)及約束迭代歷史曲線如圖9 所示.兩種方法的最優(yōu)拓?fù)浼癕ises應(yīng)力分布均是接近的,化整-集成解法比單獨的化整解法收斂速度略快.

圖7 算例2 最優(yōu)拓?fù)銯ig.7 Optimized topology of the Example 2

圖8 算例2 最優(yōu)拓?fù)涞腗ises 應(yīng)力分布Fig.8 Mises stress distribution of the optimized topology for the Example 2

圖9 算例2 目標(biāo)及約束迭代歷史.Fig.9 The iteration history of the objectives and constraints for the Example 2

4.3 算例3:多工況s 方模型化整交融應(yīng)力拓?fù)鋬?yōu)化

如圖10 所示,基結(jié)構(gòu)平面體的尺寸為200 mm×100 mm×10 mm,采用 200×100 網(wǎng)格.材料彈性模量為 E=210 GPa,泊松比 ν=0.3 .左右邊界固定,工況1:集中載荷 F1=2600 N 作用于上邊界中點;工況2:集中載荷 F2=2600 N 作用于下邊界中點.劃分為200×100個矩形單元.應(yīng)力約束為200 MPa,收斂精度取0.001.

圖10 兩工況算例分析及優(yōu)化的模型Fig.10 Model for analysis and optimization of the example with two load cases

此兩工況應(yīng)力約束拓?fù)鋬?yōu)化問題在文獻[8-9]中采用單獨的化整解法進行過求解.在此用它來比較多工況情況下化整解法與化整-集成解法的效果.

化整解法經(jīng)過41 次迭代收斂,最優(yōu)體積比為35.895%,最大Mises 應(yīng)力為200.771 MPa.化整-集成解法經(jīng)過40 次迭代收斂,最優(yōu)體積比為24.055%,最大Mises 應(yīng)力為200.171 MPa.兩種解法得到的最優(yōu)拓?fù)淙鐖D11 所示,對應(yīng)的兩個工況下的Mises 應(yīng)力分布如圖12～ 圖13 所示.經(jīng)對比可以看到:兩種解法得到的最優(yōu)拓?fù)錁?gòu)型相差較遠.雖然兩種解法收斂速度差不多,但在同樣滿足應(yīng)力約束的條件下,化整交融解法得到的最優(yōu)拓?fù)淇傮w積更小,表現(xiàn)出更好的尋優(yōu)能力.

圖11 算例3 最優(yōu)拓?fù)銯ig.11 Optimized topology of the Example 3

圖12 算例3 化整方法最優(yōu)拓?fù)涞腗ises 應(yīng)力分布Fig.12 Mises stress distribution of of the optimized topology for the Example 3 by globalization approach

圖13 算例3 化整-集成方法方法最優(yōu)拓?fù)涞腗ises 應(yīng)力分布Fig.13 Mises stress distribution of the optimized topology for the Example 3 by globalization-integration approach

5 結(jié)論

對連續(xù)體結(jié)構(gòu)應(yīng)力約束拓?fù)鋬?yōu)化問題,本文提出了求解m 方模型的乘子法及SQP 法,并與MMA法進行了對比,結(jié)果表明乘子法及SQP 法求解效率遠高于MMA 法.在此基礎(chǔ)上,構(gòu)建了應(yīng)力約束拓?fù)鋬?yōu)化問題的化整交融解法,其中求解m 方集成模型選用效果最好的SQP 解法.與文獻[8-9]中原所提出的化整解法進行了對比.單工況及多工況算例均表明:化整交融解法求解效率與化整解法基本相當(dāng),但尋優(yōu)能力更好,得到的最優(yōu)拓?fù)浣Y(jié)構(gòu)在滿足應(yīng)力約束條件下更輕.

順便指出:本文的方法對于其他類型的局部性能例如疲勞或局部穩(wěn)定約束下的結(jié)構(gòu)拓?fù)鋬?yōu)化都是適用的,相關(guān)研究將另文發(fā)表.