彭細(xì)榮， 隋允康

(1.湖南城市學(xué)院土木工程學(xué)院， 湖南 益陽(yáng) 413000；2.北京工業(yè)大學(xué)工程數(shù)值模擬中心， 北京 100124)

?

應(yīng)該為結(jié)構(gòu)拓?fù)鋬?yōu)化的設(shè)計(jì)變量正名和處理方法順言

彭細(xì)榮1， 隋允康2

(1.湖南城市學(xué)院土木工程學(xué)院， 湖南 益陽(yáng) 413000；2.北京工業(yè)大學(xué)工程數(shù)值模擬中心， 北京 100124)

為了澄清拓?fù)浠拍罴按龠M(jìn)發(fā)展，從基本概念入手，闡述了為拓?fù)渥兞空?、為有關(guān)映射順言的必要性和迫切性. 連續(xù)體結(jié)構(gòu)拓?fù)鋬?yōu)化方面的研究，以往學(xué)者們不注意厘清拓?fù)渥兞康母拍?，沒(méi)有探討其獨(dú)立地位，而是依附于結(jié)構(gòu)優(yōu)化低層次，采用懲罰函數(shù)使變量趨于0/1，有研究者于2004年意識(shí)到應(yīng)當(dāng)引進(jìn)Heaviside函數(shù)取代懲罰函數(shù). 其實(shí)，筆者在1996年提出的獨(dú)立連續(xù)映射(independent continuous mapping，ICM)方法中已經(jīng)徹底解決了這些問(wèn)題：定義獨(dú)立層次的拓?fù)渥兞?；采用階躍函數(shù)及其逆函數(shù)的近似逼近. 本文闡述了ICM方法的概念，并以頻率約束拓?fù)鋬?yōu)化為例詳述了建模及求解. 算例表明ICM方法因概念清晰顯示出明顯優(yōu)勢(shì). 《論語(yǔ)·子路篇》曰：“名不正則言不順，言不順則事不成.”正名方可順言，是時(shí)候該為結(jié)構(gòu)拓?fù)鋬?yōu)化的設(shè)計(jì)變量正名和處理方法順言了.

連續(xù)體結(jié)構(gòu)拓?fù)鋬?yōu)化；拓?fù)湓O(shè)計(jì)變量的獨(dú)立化；近似逼近函數(shù)；懲罰函數(shù)；階躍函數(shù)；Heaviside函數(shù)；獨(dú)立連續(xù)映射(independent continuous mapping，ICM)方法；變密度法

連續(xù)體結(jié)構(gòu)拓?fù)鋬?yōu)化是指通過(guò)子域的“有”或“無(wú)”的拓?fù)洳季郑跐M足一定約束條件下，求得使某目標(biāo)最優(yōu)的結(jié)構(gòu)構(gòu)型[1]. 約束及目標(biāo)可取為結(jié)構(gòu)強(qiáng)度、剛度等性能指標(biāo)或材料用量等經(jīng)濟(jì)指標(biāo).

目前使用的連續(xù)體結(jié)構(gòu)拓?fù)鋬?yōu)化大部分是構(gòu)建在基結(jié)構(gòu)(ground structure)之上的[2]，其拓?fù)鋬?yōu)化模型在本質(zhì)上屬于大型離散規(guī)劃問(wèn)題. 對(duì)離散模型的求解目前沒(méi)有高效算法，智能算法，如模擬退火算法(simulated annealing，SA)[3]、遺傳算法(genetic algorithm，GA)[4]、粒子群優(yōu)化算法(partical swarm optimization algorithm，PSOA)[5]等，其普遍特點(diǎn)是收斂速度都較慢. 啟發(fā)式的進(jìn)化結(jié)構(gòu)優(yōu)化(evolutionary structural optimization，ESO )[6]法按某種準(zhǔn)則，對(duì)單元進(jìn)行刪除或增加，其收斂速度也較慢，且最優(yōu)構(gòu)型受準(zhǔn)則參數(shù)影響大等嚴(yán)重困擾[2,7].

為了克服離散變量?jī)?yōu)化問(wèn)題的求解困難，出現(xiàn)了用連續(xù)變量“代替”或“逼近”離散變量的做法. 連續(xù)變量“代替”離散變量的大部分做法是將拓?fù)鋵哟蝺?yōu)化問(wèn)題降低為低層次尺寸或形狀優(yōu)化問(wèn)題，用連續(xù)的材料性能變量、截面變量或形狀變量取代離散的拓?fù)渥兞? 主要可分為兩大類(lèi)方法：一類(lèi)方法是依據(jù)設(shè)計(jì)變量在設(shè)計(jì)區(qū)域內(nèi)的分布場(chǎng)來(lái)確定，典型方法有均勻化方法[8]、變密度法[9]、變厚度法[10]等；另一類(lèi)方法是在設(shè)計(jì)區(qū)域內(nèi)通過(guò)孔洞的形成、合并及孔洞與邊界的形狀改變來(lái)確定，典型方法有水平集法[11]、相場(chǎng)法[12]及拓?fù)鋵?dǎo)數(shù)法[13-14]等.

由于不同方法間的參考與借鑒，不同方法間呈現(xiàn)交叉、融合的趨勢(shì). 如變密度法的密度過(guò)濾所采用的投影技術(shù)類(lèi)似于水平集方法中的水平集函數(shù)[15]，可看作一種參數(shù)化的水平集方法[16]. 使用人工假想材料描述邊界的水平集方法，“空”域內(nèi)單元指定一個(gè)很小的人工密度，“實(shí)”域內(nèi)單元指定人工密度為1，邊界附近單元人工密度依據(jù)Heaviside 投影確定取適當(dāng)?shù)闹虚g密度值，其實(shí)質(zhì)是在用水平集函數(shù)控制人工密度場(chǎng)分布，而非真實(shí)的有限元網(wǎng)格邊界[17]，可看作一種變密度法[16]. 基于敏度分析及過(guò)濾技術(shù)改進(jìn)的雙向進(jìn)化結(jié)構(gòu)優(yōu)化(bi-directional evolutionary structural optimization，BESO)方法，已不完全是初始的基于生物進(jìn)化啟發(fā)而提出的完全離散的優(yōu)化算法[6,18]，可看作一種采取離散更新策略[16]的固體各向同性材料懲罰(solid isotropic material with penalization，SIMP)[9]法.

從連續(xù)體結(jié)構(gòu)拓?fù)鋬?yōu)化方法的發(fā)展來(lái)看，無(wú)論是基于投影技術(shù)的變密度法，還是基于人工材料描述邊界的水平集方法，材料方式描述拓?fù)浠蛐螤罘绞矫枋鐾負(fù)湟呀?jīng)開(kāi)始融合，基于設(shè)計(jì)變量分布場(chǎng)來(lái)確定結(jié)構(gòu)拓?fù)洳季值母拍钜呀?jīng)形成. 然而，不幸的是，由于歷史發(fā)展的原因，這個(gè)設(shè)計(jì)變量分布場(chǎng)卻生硬地依附在人工材料密度或水平集函數(shù)所描述的邊界形狀等一些低層次材料或形狀優(yōu)化變量上，使得在問(wèn)題描述及建模上出現(xiàn)了各種所謂的雜交方法. 既然各種方法開(kāi)始融合，是時(shí)候應(yīng)努力尋求各種方法歸一的可能，即各類(lèi)初始從不視角提出的方法，是否有可能在某種理論框架下進(jìn)行統(tǒng)一.

突破口應(yīng)回歸到連續(xù)體結(jié)構(gòu)拓?fù)鋬?yōu)化問(wèn)題的本源上來(lái). 一是選取什么量作為設(shè)計(jì)變量？二是如何利用基于連續(xù)變量的優(yōu)化算法來(lái)解決大規(guī)劃離散優(yōu)化的困難？

本質(zhì)上，上述概念和處理方法相關(guān)的問(wèn)題，早在1996年由隋允康提出的拓?fù)鋬?yōu)化獨(dú)立連續(xù)映射(independent contimous mapping，ICM)方法已經(jīng)解決：一是ICM方法定義了獨(dú)立層次的拓?fù)渥兞浚欢峭ㄟ^(guò)階躍函數(shù)及其逆函數(shù)的近似逼近，在連續(xù)拓?fù)渥兞康幕A(chǔ)上建立及求解拓?fù)鋬?yōu)化模型[19].

連續(xù)體結(jié)構(gòu)拓?fù)鋬?yōu)化發(fā)展之初，對(duì)于這2個(gè)問(wèn)題，大多方法將設(shè)計(jì)變量依附于低層次的變量上，如板的厚度、微結(jié)構(gòu)尺寸、人工材料密度、水平集函數(shù)描述的隱式邊界形狀等，通過(guò)這樣處理，以達(dá)到將一個(gè)大規(guī)劃離散優(yōu)化問(wèn)題轉(zhuǎn)化為連續(xù)變量?jī)?yōu)化問(wèn)題求解的目的. Sigmund等[16]提到基于密度的各類(lèi)拓?fù)鋬?yōu)化方法中，設(shè)計(jì)變量場(chǎng)為純數(shù)學(xué)意義上的概念. 以人工密度來(lái)表征子域的“有”或“無(wú)”，采用懲罰函數(shù)的概念來(lái)實(shí)現(xiàn)大部分設(shè)計(jì)變量取0/1，其實(shí)，這種處理方式是算不上純數(shù)學(xué)的拓?fù)鋬?yōu)化變量. Heasivide投影法[15,20]用于變密度法中過(guò)濾處理以得到無(wú)棋盤(pán)格及網(wǎng)格依賴(lài)現(xiàn)象、拓?fù)淝逦?即過(guò)渡單元很少)的最優(yōu)拓?fù)?，其用連續(xù)函數(shù)近似Heasivide函數(shù)(其為階躍函數(shù)的別名)類(lèi)似于ICM法中磨光函數(shù)逼近于階躍函數(shù)，但其拓?fù)鋬?yōu)化建模仍是變密度法，設(shè)計(jì)變量仍是人工密度. 到目前為止，綜合國(guó)內(nèi)外相關(guān)研究文獻(xiàn)來(lái)看，除ICM方法外，還沒(méi)有其他方法提出過(guò)獨(dú)立層次的拓?fù)渥兞考跋鄳?yīng)的建模處理方法.

近年來(lái)，一些學(xué)者敏銳地意識(shí)到還原拓?fù)渥兞开?dú)立地位的重要意義，ICM方法得到他們的關(guān)注，并且做出了可喜的研究工作，如龍凱等[21-22]基于ICM方法提出了基于節(jié)點(diǎn)獨(dú)立變量的方法和混合插值建模方法，后來(lái)定義物質(zhì)點(diǎn)拓?fù)渥兞繛閇0 , 1]上用于表征物質(zhì)點(diǎn)及其領(lǐng)域存在與否的實(shí)數(shù)[23]，以拓?fù)渥兞縼?lái)表征一個(gè)“物質(zhì)點(diǎn)”的“有”或“無(wú)”，從概念上來(lái)講還是ICM方法所提出的獨(dú)立連續(xù)拓?fù)渥兞? 鄧果等[24]基于ICM方法，研究了位移線性近似式和應(yīng)力約束轉(zhuǎn)換表達(dá)式以及有效結(jié)構(gòu)信息到結(jié)構(gòu)最大設(shè)計(jì)域信息的映射轉(zhuǎn)換方法.

連續(xù)體結(jié)構(gòu)拓?fù)鋬?yōu)化發(fā)展到今天，可謂方興未艾，卻在最基本概念上如此良莠紛繁而未能趨同. 然而，在建筑結(jié)構(gòu)形態(tài)創(chuàng)構(gòu)領(lǐng)域，研究者卻自然應(yīng)用了拓?fù)渥兞康母拍顏?lái)描述拓?fù)鋬?yōu)化問(wèn)題[25]，作為結(jié)構(gòu)拓?fù)鋬?yōu)化的一個(gè)延展應(yīng)用的領(lǐng)域，卻在基本概念的使用上更為自覺(jué)及合理，這不得不引起結(jié)構(gòu)拓?fù)鋬?yōu)化領(lǐng)域研究者的反思乃至自省：是時(shí)候該為結(jié)構(gòu)拓?fù)鋬?yōu)化設(shè)計(jì)變量正名了，同時(shí)也該為相應(yīng)處理方法順言了.

1 獨(dú)立層次拓?fù)湓O(shè)計(jì)變量定義及拓?fù)鋬?yōu)化模型

求解連續(xù)體結(jié)構(gòu)拓?fù)鋬?yōu)化問(wèn)題，旨在確定設(shè)計(jì)區(qū)域內(nèi)各子域的“有”或“無(wú)”，應(yīng)該怎樣定義拓?fù)湓O(shè)計(jì)變量呢？

ICM方法解決這一問(wèn)題時(shí)，首先考察了低層次的截面及形狀優(yōu)化設(shè)計(jì)變量的取法：在截面優(yōu)化中，需要設(shè)計(jì)的是桿件或板的截面，其設(shè)計(jì)變量通常取桿件的截面尺寸或面積、板的厚度等；在形狀優(yōu)化中，需要設(shè)計(jì)的是結(jié)構(gòu)的形狀，其設(shè)計(jì)變量通常取骨架類(lèi)結(jié)構(gòu)的連接節(jié)點(diǎn)坐標(biāo)或連續(xù)體結(jié)構(gòu)的邊界或孔洞的有限元網(wǎng)格節(jié)點(diǎn)坐標(biāo)或形狀描述函數(shù)參數(shù)，亦即取“形狀”作為設(shè)計(jì)變量. 到了拓?fù)鋬?yōu)化層次，設(shè)計(jì)變量不再是具體的幾何量或物理量，而是各子域的“有”或“無(wú)”，很自然，各子域“有”或“無(wú)”的拓?fù)錉顟B(tài)就可以分別用0和1兩個(gè)數(shù)來(lái)表征，即：

子域——一個(gè)幾何點(diǎn)鄰域內(nèi)足夠小的區(qū)域.

子域的拓?fù)渥兞渴请x散量0或1，分別代表此子域?yàn)椤盁o(wú)”與“有”.

不難看出，上述定義的拓?fù)湓O(shè)計(jì)變量是一個(gè)純粹的數(shù)學(xué)量. 對(duì)連續(xù)體結(jié)構(gòu)拓?fù)鋬?yōu)化問(wèn)題，作為子域“有”或“無(wú)”的表征，也可以稱(chēng)為子域拓?fù)浞植紙?chǎng)的平均值. 可稱(chēng)為離散拓?fù)湓O(shè)計(jì)變量或拓?fù)浞植紙?chǎng)的離散平均值. 這里強(qiáng)調(diào)離散，是為了區(qū)分ICM方法定義的連續(xù)拓?fù)渥兞?

可以看到，由于定義了表征子域“有”或“無(wú)”的拓?fù)錉顟B(tài)量，結(jié)構(gòu)拓?fù)涞拿枋鼍妥兊梅浅Ｇ逦白匀? 選取其他物理量作為拓?fù)鋬?yōu)化模型的設(shè)計(jì)變量，如變厚度法中的板厚度、均勻化方法中的微結(jié)構(gòu)尺寸、變密度法中的人工密度等均沒(méi)有直接反映拓?fù)鋬?yōu)化的本質(zhì)，應(yīng)是拓?fù)鋬?yōu)化發(fā)展初期，為解決大規(guī)模離散優(yōu)化模型的“組合爆炸”困難，將問(wèn)題轉(zhuǎn)化為基于連續(xù)變量的規(guī)劃模型而采取的不得已的一些物理轉(zhuǎn)化方式.

ICM方法則不依賴(lài)物理轉(zhuǎn)化，而是把離散的0/1變量延拓到[0,1]，0與1之間的中間實(shí)數(shù)值類(lèi)似于模糊集合論中的隸屬度，表示對(duì)于0或1的靠近程度. 這就是連續(xù)的拓?fù)渥兞?

另一方面，發(fā)展到現(xiàn)在，變密度法由于其廣泛的研究及應(yīng)用，“密度”一詞一直沿用，但其所表達(dá)的意義已由初始的人工密度的起始概念悄悄發(fā)生了變化，其意義已非常接近于ICM方法中的連續(xù)拓?fù)渥兞? 對(duì)比一下變密度法中用“密度”這個(gè)詞匯來(lái)表征子域“有”或“無(wú)”的拓?fù)錉顟B(tài)，其值取值在[0,1]，利用此人工密度變量，建立單元體積及單元彈性模量的函數(shù)關(guān)系. 很明顯，變密度法在此處使用密度這一詞匯是極易與真正的材料密度概念產(chǎn)生混亂的. 為什么要用“密度”來(lái)說(shuō)一個(gè)子域的“有”或“無(wú)”，讓人感覺(jué)非常奇怪，難道真正的材料密度必須用人造的、人工的、假的或偽的“密度”表達(dá)嗎？

更重要的是，伴隨著純數(shù)學(xué)的拓?fù)湓O(shè)計(jì)變量由1變成0，單元所具備的所有的幾何量、物理量到從有變成無(wú)，不僅僅是密度從有變成無(wú)，因此，單獨(dú)關(guān)注密度，并不是公正的. 可見(jiàn)，最合理的做法是采用純數(shù)學(xué)的拓?fù)湓O(shè)計(jì)變量，ICM方法稱(chēng)之為獨(dú)立連續(xù)的拓?fù)渥兞?

選定了設(shè)計(jì)變量為離散拓?fù)渥兞浚B續(xù)體結(jié)構(gòu)拓?fù)鋬?yōu)化的提法可表達(dá)為

(1)

式中：ti為子域i對(duì)應(yīng)的離散拓?fù)湓O(shè)計(jì)變量，取值0或1，t為由tj組成的設(shè)計(jì)變量向量；c(t)為目標(biāo)函數(shù)；Aj(t)及Bk(t)為等式及不等式約束條件；J及K分別為等式及不等式約束數(shù)目；N為設(shè)計(jì)變量數(shù)目.

這是個(gè)離散數(shù)學(xué)規(guī)劃問(wèn)題，在基于有限單元法的結(jié)構(gòu)拓?fù)鋬?yōu)化中，通常取一個(gè)單元對(duì)應(yīng)于一個(gè)子域，即對(duì)應(yīng)于一個(gè)拓?fù)湓O(shè)計(jì)變量，由于有限元單元數(shù)量通常成千上萬(wàn)，因此，此離散數(shù)學(xué)規(guī)劃的設(shè)計(jì)變量數(shù)通常很大，而求解大規(guī)模離散數(shù)學(xué)規(guī)劃問(wèn)題遭遇極為棘手的“組合爆炸”困難. 如引言中所述，為了克服離散變量?jī)?yōu)化問(wèn)題的求解困難，通常用連續(xù)變量“代替”或“逼近”離散變量的做法，均勻化方法、變厚度法、變密度法及水平集法等采用“代替”的做法，唯有ICM法從函數(shù)逼近的視角采用“逼近”的做法.

2 離散及連續(xù)拓?fù)渥兞哭D(zhuǎn)化——階躍函數(shù)及其逆函數(shù)近似逼近

2.1 離散拓?fù)渥兞颗c階躍函數(shù)

考察變厚度法、均勻化方法或變密度法等中的設(shè)計(jì)變量如板厚度、微結(jié)構(gòu)尺寸或人工材料密度等，這些在結(jié)構(gòu)拓?fù)鋬?yōu)化過(guò)程中變化的物理量，無(wú)論其如何小，對(duì)應(yīng)的拓?fù)渥兞繛?，而當(dāng)其等于0時(shí)，拓?fù)渥兞客蛔優(yōu)?；物理量是連續(xù)變化的，但拓?fù)渥兞繀s是在零點(diǎn)跳躍的. ICM法引入階躍函數(shù)描述上述關(guān)系，函數(shù)圖形如圖1所示.

階躍函數(shù)為

(2)

2.2 連續(xù)拓?fù)渥兞颗c磨光函數(shù)

(3)

磨光函數(shù)逼近階躍函數(shù)，使拓?fù)渥兞坑扇≈?/1的離散變量擴(kuò)展為到取值在[0, 1]的連續(xù)變量. 此函數(shù)逼近過(guò)程稱(chēng)為磨光近似映射. 連續(xù)拓?fù)渥兞咳?0,1)值時(shí)反映了對(duì)應(yīng)子域“有”與“無(wú)”的程度.

函數(shù)圖形如圖2所示.

2.3 基于連續(xù)拓?fù)渥兞康耐負(fù)浠Ｅc求解

取連續(xù)拓?fù)渥兞孔鳛樵O(shè)計(jì)變量，連續(xù)體結(jié)構(gòu)拓?fù)鋬?yōu)化的提法可表達(dá)為

(4)

式中ti為子域i對(duì)應(yīng)的連續(xù)拓?fù)湓O(shè)計(jì)變量，取值[0, 1]. 這是個(gè)連續(xù)變量數(shù)學(xué)規(guī)劃問(wèn)題，可應(yīng)用基于連續(xù)函數(shù)的高效求解算法求解，如對(duì)偶序列二次規(guī)劃法、移動(dòng)漸近線法(method of moving asymptotes,MMA)等.

2.4 離散連續(xù)拓?fù)渥兞块g的映射與求解

求解基于連續(xù)拓?fù)湓O(shè)計(jì)變量建立的結(jié)構(gòu)拓?fù)鋬?yōu)化模型(如式(4)所示)，得到最優(yōu)的設(shè)計(jì)變量分布場(chǎng)t*，稱(chēng)之為連續(xù)拓?fù)浞植紙?chǎng). 最后取連續(xù)值的拓?fù)渥兞啃枰貧w到取0/1離散值的拓?fù)渥兞浚陔A躍函數(shù)的逆函數(shù)(為引用方便，以下均稱(chēng)之為跨欄函數(shù))中，拓?fù)渥兞恐蝗‰x散值0/1，如圖3所示，因此，要使得拓?fù)渥兞咳≈虚g值，需要一個(gè)連續(xù)可導(dǎo)函數(shù)(稱(chēng)為過(guò)濾函數(shù)，如圖4所示)逼近跨欄函數(shù)

(5)

過(guò)濾函數(shù)實(shí)現(xiàn)將連續(xù)拓?fù)渥兞侩x散化，同時(shí)在建模時(shí)，起到了對(duì)相應(yīng)的子域物理量進(jìn)行識(shí)別的作用：

(6)

過(guò)濾函數(shù)為磨光函數(shù)的逆函數(shù)，跨欄函數(shù)為階躍函數(shù)的逆函數(shù)，過(guò)濾函數(shù)近似逼近跨欄函數(shù)，而磨光函數(shù)近似逼近階躍函數(shù)，4種函數(shù)的關(guān)系如圖5所示.

過(guò)濾函數(shù)或磨光函數(shù)可以為多樣函數(shù)形式，到目前為止，ICM法研究過(guò)的過(guò)濾函數(shù)大致有3類(lèi)：冪函數(shù)、修正的Sigmoid函數(shù)和指數(shù)函數(shù)[7].

綜上所述，ICM法定義了獨(dú)立拓?fù)渥兞?，?yīng)用函數(shù)逼近理論，對(duì)階躍函數(shù)及其逆函數(shù)進(jìn)行連續(xù)可導(dǎo)化的函數(shù)逼近近似，從而實(shí)現(xiàn)用基于連續(xù)變量的優(yōu)化模型求解原本離散的大規(guī)劃優(yōu)化問(wèn)題.

(7)

所示Heaviside函數(shù)表示這種關(guān)系，然后，用一個(gè)連續(xù)光滑的函數(shù)近似逼近此函數(shù)，如

(8)

所示.

β取不同值時(shí)的近似程度如圖6所示，可以看到，此與ICM法中磨光函數(shù)對(duì)階躍函數(shù)的近似是類(lèi)似的(見(jiàn)圖2)，然而其應(yīng)用卻是完全不同的. 兩者之間的共同點(diǎn)是均用了階躍函數(shù)來(lái)表征單元的“有”或“無(wú)”，均用了連續(xù)光滑函數(shù)對(duì)階躍函數(shù)逼近. 但I(xiàn)CM法是用階躍函數(shù)的近似逼近來(lái)將離散拓?fù)渥兞窟B續(xù)化，將不可導(dǎo)函數(shù)變換為連續(xù)光滑的可導(dǎo)函數(shù)，從而用連續(xù)變量規(guī)劃模型來(lái)求解原離散規(guī)劃模型，而Heaviside投影法是應(yīng)用在拓?fù)溥吔绲募怃J化處理上，迫使過(guò)渡單元的人工密度盡量取值為0或1.

從上述對(duì)Heaviside投影法及ICM法的對(duì)比分析可以看到，發(fā)展到現(xiàn)在，人們逐漸意識(shí)到階躍函數(shù)與單元“有”或“元”表達(dá)間的關(guān)系及應(yīng)用光滑連續(xù)函數(shù)進(jìn)行建模處理的好處，但還處于朦朧階段，在此發(fā)展階段，深入闡述ICM法的獨(dú)立層次拓?fù)渥兞考捌浠陔A躍函數(shù)及其逆函數(shù)近似的建模處理方法，對(duì)加深對(duì)拓?fù)鋬?yōu)化問(wèn)題本質(zhì)的認(rèn)識(shí)及促進(jìn)其發(fā)展是有非常重要的意義.

3 ICM法合理化建模及求解

選取連續(xù)拓?fù)渥兞孔鳛樵O(shè)計(jì)變量，可以建立如式(4)形式的連續(xù)體結(jié)構(gòu)拓?fù)鋬?yōu)化模型. 這里還有2個(gè)問(wèn)題：

問(wèn)題1 式(4)所示模型中目標(biāo)函數(shù)c(x)，等式及不等式約束條件Aj(x)及Bk(x)的選取問(wèn)題.

問(wèn)題2 如何將目標(biāo)函數(shù)及約束條件表達(dá)為連續(xù)拓?fù)湓O(shè)計(jì)變量的近似顯式函數(shù)的問(wèn)題.

對(duì)問(wèn)題1，將結(jié)構(gòu)優(yōu)化問(wèn)題涉及的物理量劃分為兩大類(lèi)指標(biāo)：經(jīng)濟(jì)指標(biāo)，如結(jié)構(gòu)總重量(或總體積)、總造價(jià)等；性能指標(biāo)，如結(jié)構(gòu)應(yīng)力、位移、頻率等. 則結(jié)構(gòu)優(yōu)化模型可分為兩大類(lèi)：1) 在滿足性能指標(biāo)下結(jié)構(gòu)經(jīng)濟(jì)指標(biāo)最小化問(wèn)題；2) 在結(jié)構(gòu)經(jīng)濟(jì)指標(biāo)限制下的結(jié)構(gòu)性能最優(yōu)化問(wèn)題.

連續(xù)體結(jié)構(gòu)拓?fù)鋬?yōu)化中研究得最多的模型是以結(jié)構(gòu)總體積(即對(duì)應(yīng)于結(jié)構(gòu)總重量)約束下結(jié)構(gòu)柔順度極小化問(wèn)題. Sigmund[16]指出：應(yīng)實(shí)際工程問(wèn)題的需要，可以加入其他性能約束條件，但體積約束即使不重要，也作為一個(gè)基本約束條件加入在優(yōu)化模型中以提高求解過(guò)程收斂的穩(wěn)定性.

考察低層次的截面及形狀優(yōu)化問(wèn)題，可以發(fā)現(xiàn)通常建立的優(yōu)化模型是第1類(lèi)以經(jīng)濟(jì)指標(biāo)最小化為目標(biāo)的模型，而發(fā)展到拓?fù)鋬?yōu)化層次，模型的主流卻發(fā)生了變化，轉(zhuǎn)變成了第2類(lèi)模型. 拓?fù)鋬?yōu)化以體積約束作為一個(gè)基本約束，建立第2類(lèi)模型， 存在如下不足：

1) 體積約束值的確定是沒(méi)有依據(jù)的、先驗(yàn)的，指定不同的體積約束值會(huì)得到不同的最優(yōu)拓?fù)錁?gòu)型.

2) 取性能為目標(biāo)，在很多優(yōu)化問(wèn)題中不得不建立多目標(biāo)優(yōu)化模型，而多目標(biāo)間的權(quán)重系數(shù)也是先驗(yàn)的.

3) 與低層次的截面及形狀優(yōu)化的模型不一致，不便于建立截面、形狀及拓?fù)浼蓛?yōu)化模型.

隋允康等[26-27]研究表明：第1種模型比第2種模型更加合理. 在ICM法中，均是建立第1類(lèi)優(yōu)化模型.

對(duì)問(wèn)題2，在過(guò)濾函數(shù)逼近跨欄函數(shù)的過(guò)程中，建立了連續(xù)拓?fù)湓O(shè)計(jì)變量與各物理量的顯式函數(shù)關(guān)系，如式(6)所示，在此基礎(chǔ)上，通過(guò)敏度分析求得各物理量對(duì)拓?fù)湓O(shè)計(jì)變量的偏導(dǎo)數(shù)，通過(guò)一階Taylor展式，可以建立目標(biāo)函數(shù)或約束與連續(xù)拓?fù)湓O(shè)計(jì)變量間的近似顯式函數(shù).

以重量極小受結(jié)構(gòu)頻率約束的拓?fù)鋬?yōu)化問(wèn)題為例，其優(yōu)化模型為

(9)

上述3個(gè)過(guò)濾函數(shù)取為冪函數(shù)形式，分別為

fw(ti)=tαw,fk(ti)=tαk,fm(ti)=tαm

(10)

為避免當(dāng)拓?fù)渥兞縯i→0時(shí)，質(zhì)量剛度比趨于無(wú)窮大，取αw=a,αk=αm=b.

各單元重量、單元?jiǎng)偠染仃嚰皢卧|(zhì)量矩陣與拓?fù)湓O(shè)計(jì)變量間的關(guān)系分別用對(duì)應(yīng)的過(guò)濾函數(shù)進(jìn)行識(shí)別：

(11)

式中上標(biāo)k指第k次迭代的值.

記

令

則所有類(lèi)型頻率約束可統(tǒng)一表述為

由此得到頻率約束下、重量最輕為目標(biāo)的連續(xù)體結(jié)構(gòu)拓?fù)鋬?yōu)化模型為

(12)

式(12)可分離變量的非線性規(guī)劃問(wèn)題，可采用對(duì)偶系列二次規(guī)劃法求解.

迭代求解的收斂準(zhǔn)則為

ΔW=|(W(k+1)-W(k))/W(k+1)|≤ε

(13)

式中：W(k)及W(k+1)為前輪與本輪迭代的結(jié)構(gòu)總重量；ε為收斂精度，本文取ε=0.001.

其他類(lèi)型性能約束的拓?fù)鋬?yōu)化問(wèn)題建模及求解的過(guò)程類(lèi)似[8].

4 數(shù)值算例

算例1 多點(diǎn)位移約束拓?fù)鋬?yōu)化

如圖7所示，基本結(jié)構(gòu)為16 mm×10 mm的平面，初始厚度為1 mm，材料彈性模量為200 GPa，泊松比為0.3. 劃分為64×40個(gè)矩形單元，載荷F=1 000 N. 左、右下腳點(diǎn)采用固定約束.A點(diǎn)處豎直向下位移約束為0.25 mm，B、C點(diǎn)處豎直向下位移約束為0.2 mm.

結(jié)構(gòu)初始總體積為160 mm3，A、B及C三點(diǎn)豎向的位移分別為-0.100 6、-0.940 2及-0.940 2 mm. 經(jīng)39次迭代后得到的拓?fù)鋱D形，如圖8所示，目標(biāo)及約束的歷史迭代曲線如圖9及10所示. 最優(yōu)點(diǎn)總體積35.852 1 mm3，A、B及C三點(diǎn)最優(yōu)點(diǎn)位移分別為-0.249 7、-0.199 9及-0.199 9 mm.

算例2 多個(gè)頻率約束拓?fù)鋬?yōu)化

如圖11所示，60 mm×20 mm×2 mm的矩形板，2個(gè)角點(diǎn)固定，在下邊界布置3個(gè)集中質(zhì)量塊，1/2跨處的集中質(zhì)量塊為8 g，1/4及3/4跨處的集中質(zhì)量塊為4 g，集中質(zhì)量塊只有Y方向的慣性，彈性模量E=1×106MPa，泊松比為0.3， 材料密度ρ=1 mg/mm3. 用矩形單元，劃分為120×40=4 800個(gè)單元.

為比較動(dòng)態(tài)加入頻率約束的方式解決模態(tài)交換的效果，計(jì)算2種情況：

1) 只考慮指定的頻率約束.

2) 除考慮指定的頻率約束外，動(dòng)態(tài)加入防止模態(tài)交換的頻率約束.

計(jì)算情況1)不收斂，得到的部分迭代過(guò)程目標(biāo)及約束的歷史曲線如圖12、13所示，計(jì)算情況2)經(jīng)過(guò)51次迭代后收斂，得到的最優(yōu)拓?fù)鋱D形如圖14，目標(biāo)及約束的迭代歷史曲線如圖15、16所示，最優(yōu)點(diǎn)處體積為512.12 mm3， 前三階頻率為15 007.78、15 770.05、30 058.07 Hz.

從迭代歷史曲線中可以看出：在沒(méi)有加入防止模態(tài)交換的頻率約束的情況下，在迭代進(jìn)行到一定程度后，第一階頻率值與第二階頻率值非常接近，從而發(fā)生模態(tài)交換，使得敏度計(jì)算得到錯(cuò)誤的信息，從而引起迭代過(guò)程中出現(xiàn)振蕩，最終不收斂；加入防止模態(tài)交換的頻率約束后，迭代過(guò)程非常平穩(wěn). 由此可見(jiàn)，動(dòng)態(tài)加入頻率約束的方式確實(shí)可以有效地防止迭代過(guò)程中模態(tài)交換現(xiàn)象的發(fā)生，從而達(dá)到穩(wěn)定收斂過(guò)程的目的.

算例3 多點(diǎn)頻率響應(yīng)位移幅值約束拓?fù)鋬?yōu)化

如圖17所示，基結(jié)構(gòu)尺寸為520 mm×260 mm×6 mm, 材料彈性模量為68.890 GPa，泊松比為0.3，密度為7.8×103kg/m3. 集中載荷F=7 000×sin (2πft)N分別作用于下邊界A、B及C三點(diǎn)，激勵(lì)頻率為500 Hz. 底部2個(gè)角點(diǎn)固定. 位移約束A點(diǎn)許用Y方向位移為0.5 mm，B及C點(diǎn)許用Y方向位移為0.45 mm. 計(jì)算無(wú)阻尼及結(jié)構(gòu)阻尼為0.3兩種情況.

在無(wú)阻尼情況下，41次迭代后計(jì)算得到的最優(yōu)拓?fù)鋱D形如圖18所示，有阻尼情況下，43次迭代后計(jì)算得到的最優(yōu)拓?fù)鋱D形如圖19所示. 初始基結(jié)構(gòu)的位移及最優(yōu)結(jié)構(gòu)的位移比較如表1所示，位移約束迭代歷史迭代曲線及目標(biāo)迭代歷史曲線如圖20～22所示. 從最優(yōu)拓?fù)鋱D形中可以看出：在相同約束條件下，無(wú)阻尼結(jié)構(gòu)及有阻尼結(jié)構(gòu)得到的最優(yōu)拓?fù)鋱D形不同.

結(jié)構(gòu)A點(diǎn)B點(diǎn)C點(diǎn)基結(jié)構(gòu)(無(wú)阻尼)0.488260.444540.44454最優(yōu)結(jié)構(gòu)(無(wú)阻尼)0.499860.449150.44915基結(jié)構(gòu)(有阻尼)0.418480.381090.38108最優(yōu)結(jié)構(gòu)(有阻尼)0.499010.449570.44957

5 結(jié)論

1) 以往連續(xù)體結(jié)構(gòu)拓?fù)鋬?yōu)化沒(méi)有能夠厘清拓?fù)渥兞炕靖拍?，而是用依附于低層次結(jié)構(gòu)優(yōu)化上，而ICM方法定義拓?fù)渥兞勘碚髯佑颉坝小被颉盁o(wú)”的拓?fù)錉顟B(tài)，自然而清晰，不再需要將拓?fù)鋬?yōu)化問(wèn)題降格為材料優(yōu)化、尺寸優(yōu)化或形狀優(yōu)化問(wèn)題. 可惜，連續(xù)體結(jié)構(gòu)拓?fù)鋬?yōu)化研究領(lǐng)域內(nèi)還沒(méi)有對(duì)此引起足夠重視，例如變密度方法中“密度”一詞極易與真實(shí)的材料密度概念混淆.

2) 不同于采用懲罰函數(shù)進(jìn)行“有”或“無(wú)”拓?fù)錉顟B(tài)的變密度方法處理，ICM方法通過(guò)對(duì)階躍函數(shù)及其逆函數(shù)的逼近，定義了磨光函數(shù)及過(guò)濾函數(shù)，識(shí)別出各物理量與連續(xù)拓?fù)渥兞块g的函數(shù)關(guān)系，解決了目標(biāo)函數(shù)與約束條件按拓?fù)湓O(shè)計(jì)變量近似顯式化的建模問(wèn)題. 而Heaviside投影法盡管有了用階躍函數(shù)表征子域的“有”或“無(wú)”的意識(shí)，也應(yīng)用到了用連續(xù)光滑函數(shù)逼近階躍函數(shù)，但其應(yīng)用卻不是用于建立連續(xù)變量規(guī)劃模型以解決大規(guī)模離散規(guī)劃模型的求解困難，而是用在人工密度的過(guò)濾處理方面，僅起到對(duì)存在過(guò)渡單元的拓?fù)溥吔邕M(jìn)行銳化的作用.

3) ICM方法以經(jīng)濟(jì)指標(biāo)為目標(biāo)函數(shù)、性能指標(biāo)為約束函數(shù)建立的拓?fù)鋬?yōu)化模型可以避免通常的以體積為約束的拓?fù)鋬?yōu)化模型的不足，與低層次的截面及形狀優(yōu)化在目標(biāo)取法上保持了一致，更有利于各層次優(yōu)化間的集成.

[1] ROZVANY G I N. A critical review of established methods of structural topology optimization [J]. Structural Multidisciplinary Optimization, 2009, 37(3): 217-237.

[2] BENDSOE M P, SIGMUND O. Topology optimization: theory, methods and applications [M]. Berlin: Springer, 2003.

[3] RAO S S. Simulated annealing and parallel processing: an implementation for constrained global design optimization [J]. Engineering Optimization, 2000, 32(5): 659-685.

[4] CHEN S Y, RAJAN S D. Robust genetic algorithm for structural optimization [J]. Structural Engineering and Mechanics, 2000, 10(4): 313-336.

[5] FOURIE P C, GROENWOLD A A. Particle swarms in topology optimization [C]∥Proceedingsof Fouthworld Congress of Structural and Multidisciplinary Optimization. Dalian: WCSMO_4, 2001.

[6] HUANG X, XIE Y M. A further review of ESO type methods for topology optimization [J]. Structural and Multidisciplinary Optimization, 2010, 41: 671-683.

[7] 隋允康, 葉紅玲. 連續(xù)體結(jié)構(gòu)拓?fù)鋬?yōu)化的ICM方法[M]. 北京: 科學(xué)出版社, 2013: 50-80.

[8] BENDSOE M P, KIKUCHI N. Generating optimal topologies in structural design using a homogenization method [J]. Computer Methods in Applied Mechanics and Engineering, 1988, 69: 197-224.

[9] MLEJNEK H P. Some aspects of the genesis of structures [J]. Structural Optimization, 1992, 5: 64-69.

[10] 程耿東, 張東旭. 受應(yīng)力約束的平面彈性體的拓?fù)鋬?yōu)化[J]. 大連理工大學(xué)學(xué)報(bào), 1995, 35(1): 1-9. CHENG G D, ZHANG D X. Topology optimization of plane elastic continuum with stress constraints [J]. Journal of Dalian University of Technology, 1995, 35(1): 1-9. (in Chinese)

[11] OSHER S, SETHIAN J. Fronts propagating with curvature dependent speed-algorithms based on Hamilton-Jacobi formulations [J]. J Comput Phys, 1988, 79(1): 12-49.

[12] BOURDIN B, CHAMBOLLE A. Design-dependent loads in topology optimization [J]. ESAIM: Control, Optimisation and Calculus of Variations, 2003, 9(8): 19-48.

[13] ESCHENAUER H A, KOBELEV V V, SCHUMACHER A. Bubble method for topology and shape optimization of structures [J]. Structural Optimization, 1994, 8(1): 42-51.

[14] NORATO J, BENDSOE O, HABER R, et al. A topological derivative method for topology optimization [J]. Structural and Multidisciplinary Optimization, 2007, 33(4/5): 375-386.

[15] GUEST J, PREVOST J, BELYTSCHKO T. Achieving minimum length scale in topology optimization using nodal design variables and projection functions [J]. Int J Numer Methods Eng, 2004, 61(2): 238-254.

[16] SIGMUND O, MAUTE K. Topology optimization approaches — a comparative review [J]. Structural Multidisciplinary Optimization, 2013, 48(6): 1031-1055.

[17] PINGEN G, WAIDMANN M, EVGRAFOV A, et al. A parametric level-set approach for topology optimization of flow domains [J]. Structural and Multidisciplinary Optimization, 2010, 41: 117-131.

[18] HUANG Y, XIE Y M. Evolutionary topology optimization of geometrically and materially nonlinear structures under prescribed design load [J]. Struct Eng Mech, 2010, 34(5): 581-595.

[19] 隋允康. 建?！ぷ儞Q·優(yōu)化——結(jié)構(gòu)綜合方法新進(jìn)展[M]. 大連：大連理工大學(xué)出版社, 1996: 90-100.

[20] KAWAMOTO A, MATSUMORI T, YAMASAKI S, et al. Heaviside projection based topology optimization by a PDE-filtered scalar function [J]. Structural and Multidisciplinary Optimization, 2011, 44: 19-24.

[21] 龍凱, 左正興, 閆清東. 靜動(dòng)態(tài)多目標(biāo)下的連續(xù)體結(jié)構(gòu)拓?fù)鋬?yōu)化[J]. 宇航學(xué)報(bào), 2008, 29(2): 456-460. LONG K, ZUO Z X, YAN Q D. Static and dynamic multi-objective topological optimization of continuum structure [J]. Journal of Astronautics, 2008, 29(2): 456-460. ( in Chinese)

[22] 龍凱, 左正興. 基于不同過(guò)濾函數(shù)的ICM拓?fù)鋬?yōu)化混合建模方法[J]. 中國(guó)機(jī)械工程, 2007, 18(19): 2303- 2306. LONG K, ZUO Z X. ICM method for topological Optimization of continuum structure based on different filter functions [J]. Chinese Journal of Mechanical Engineering, 2007, 18(19): 2303-2306. (in Chinese)

[23] 龍凱, 左正興. 物質(zhì)點(diǎn)拓?fù)渥兞糠ㄔ谌嵝詸C(jī)構(gòu)設(shè)計(jì)中的應(yīng)用[J]. 計(jì)算機(jī)輔助設(shè)計(jì)與圖形學(xué)學(xué)報(bào), 2010, 22(1): 158-164. LONG K, ZUO Z X. Topology optimization for compliant mechanism using material point topological variable[J]. Journal of Computer-Aided Design & Computer Graphics, 2010, 22(1): 158-164. (in Chinese)

[24] 鄧果, 榮見(jiàn)華, 邢曉娟. 一種考慮位移和應(yīng)力約束的結(jié)構(gòu)拓?fù)鋬?yōu)化方法[J]. 長(zhǎng)沙理工大學(xué)學(xué)報(bào), 2008, 5(3): 26-33. DENG G, RONG J H, XING X J. A structural topological optimization method considering displacement and stress constraints [J]. Journal of Changsha University of Science and Technology, 2008, 5(3): 26-33. (in Chinese)

[25] 吳志海. 基于遺傳算法的網(wǎng)殼結(jié)構(gòu)形態(tài)創(chuàng)構(gòu)方法研究[D]. 哈爾濱: 哈爾濱工業(yè)大學(xué), 2013. WU Z H. Research on form finding method of mesh shell structures based ongenetic algorithm [D]. Harbin: Harbin Institute of Technology, 2013. (in Chinese)

[26] YI G L, SUI Y K. Different effects of economic and structural performance indicators on model construction of structural topology optimization[J]. Acta Mechanica Sinia, 2015, 9: 1-12.

[27] 彭細(xì)榮, 隋允康. 對(duì)連續(xù)體結(jié)構(gòu)拓?fù)鋬?yōu)化合理模型的再探討[J]. 固體力學(xué)學(xué)報(bào), 2016, 37(1): 1-11. PENG X R, SUI Y K. A further discussion on rational topology optimization models for continuum structures [J]. Chinese Journal of Solid Mechanics, 2016, 37(1): 1-11. (in Chinese)

(責(zé)任編輯 楊開(kāi)英)

Name Correction for Design Variable of Structural Topology Optimization and Presentation of Its Corresponding Treatment Method

PENG Xirong1, SUI Yunkang2

(1.School of Civil Engineering, Hunan City University, Yiyang 413000, Hunan, China;2.Numerical Simulation Center for Engineering, Beijing University of Technology, Beijing 100124, China)

To clarify the concepts and promote developments of the structural topology optimization by analyzing the basic concepts, the necessity and urgency of name correction of the design variable of structural topology optimization and presenting its corresponding mapping method are discussed. For the topology optimization of continuum structures, the concept of topology variable was not clarified. The independent level of the topology variable was not discussed and the topology variable was attached to the variables of the structural optimization in lower levels. A penalty function method was adopted to assure the design variables of lower levels to approach to 0/1. In 2004, the Heaviside function was introduced to replace the penalty function by some researchers. Actually, the above problems and their treatments were solved completely in the independent continuous mapping (ICM) method, which was put forward in 1996. The one was the definition of the topology variable with independent level. The other was the approximations of the step function and its inverse function. In this paper, the basic concepts in the ICM method were first introduced. The process of modeling and solution was described by taking the topology optimization problem with frequency constraints as an example. Some examples show that the ICM method has obvious advantages in the aspects of modeling and solution because of its clear concepts. The Zi lu article in the Analects said that a speech will not be in right order if a name is not correct, and nothing will be successful if a speech is not in right order. It is time to correct the name of design variable of structural topology optimization and present its corresponding treatment method.

topology optimization of continuum structures; Independent topology variable; approximation functions; penalty function; step function; heaviside function; independent continuous mapping (ICM) method; variable density method

2016- 07- 24

國(guó)家自然科學(xué)基金資助項(xiàng)目(11672103)；湖南省自然科學(xué)基金資助項(xiàng)目(2016JJ6016)

彭細(xì)榮(1972—)， 男， 副教授， 主要從事結(jié)構(gòu)優(yōu)化方面的研究， E-mail：pxr568@163.com

隋允康(1943—)，男，教授，博士生導(dǎo)師，主要從事結(jié)構(gòu)優(yōu)化方面的研究，E-mail：ysui@bjut.edu.cn

O 343.1

A

0254-0037(2016)12-1787-11

10.11936/bjutxb2016070072