數(shù)學(xué)規(guī)劃法在有限元模型修正中的應(yīng)用

張 欣，于 瀾，張 淼

(1.長春工程學(xué)院電氣與信息工程學(xué)院，長春 130012；2.長春工程學(xué)院理學(xué)院，長春 130012)

0 引言

有限元模型修正是20世紀(jì)90年代以來模態(tài)分析與試驗(yàn)所面臨的一大挑戰(zhàn)。工程和科學(xué)系統(tǒng)的現(xiàn)代分析中，人們已經(jīng)投入了大量的努力用于修正被曲解的計(jì)算機(jī)模型。這種模型修正的主要目的是預(yù)測系統(tǒng)對干擾的響應(yīng)和預(yù)測從修正后系統(tǒng)中所獲得的可能的設(shè)計(jì)優(yōu)勢。從當(dāng)前的實(shí)踐來看，這一目的可以通過修正系統(tǒng)的外形結(jié)構(gòu)來獲得。在實(shí)踐中已經(jīng)能夠從離散模型中產(chǎn)生數(shù)值的預(yù)測，但當(dāng)預(yù)測結(jié)果與實(shí)驗(yàn)結(jié)果相比較時(shí)，人們發(fā)現(xiàn)沒有自信把這個(gè)模型用于實(shí)際中，因?yàn)轭A(yù)測結(jié)果與實(shí)驗(yàn)結(jié)果的相關(guān)度并不理想。

為了解釋在預(yù)測和觀察之間的相關(guān)度缺乏問題，有必要考慮導(dǎo)致解析模型不精確的幾種可能的原因。除了實(shí)驗(yàn)測量過程會(huì)產(chǎn)生誤差，還有3個(gè)一般會(huì)遇到的模型誤差的存在形式，它們也可能會(huì)引起模型預(yù)測的不準(zhǔn)確性：1)模型結(jié)構(gòu)誤差：在建模過程中發(fā)生在控制物理方程有不確定性的時(shí)候，或者在工程系統(tǒng)具有強(qiáng)非線性行為的時(shí)候；2)模型參數(shù)誤差：它通常包含于應(yīng)用了不適當(dāng)?shù)倪吔鐥l件和為了簡化模型所采用的不精確的假定的過程；3)模型階數(shù)誤差：它來自于復(fù)雜系統(tǒng)的離散化，它會(huì)導(dǎo)致產(chǎn)生了一個(gè)階數(shù)不足的模型，尤其重要的是模型階數(shù)可以被認(rèn)為是模型結(jié)構(gòu)的一部分[1]。

利用實(shí)測的結(jié)構(gòu)模態(tài)信息來修正結(jié)構(gòu)的有限元解析模型，使修正后的有限元模型的輸出信息與實(shí)測信息最大程度相符合，這就是有限元模型修正。修正后的有限元模型能夠更好地反映結(jié)構(gòu)的動(dòng)力特性[2]。目前，較為常見的有限元模型修正方法有矩陣元素型、靈敏度矩陣方程型和殘差型3種。矩陣元素型的模型修正方法的基本思路是：首先假定原始的有限元模型與真實(shí)結(jié)構(gòu)模型之間只存在微小差異，然后，除了保持剛度矩陣和質(zhì)量矩陣的對稱性、正定性和稀疏性以及振型正交性等性質(zhì)外，還要在滿足特征方程的條件下，利用最小二乘法求出符合實(shí)測結(jié)果的剛度和質(zhì)量矩陣[3]。靈敏度矩陣方程型的模型修正方法是為了使測量模態(tài)模型與分析模態(tài)模型之間的相關(guān)性達(dá)到最大化。這些方法可以廣泛地選擇修正參數(shù)，一般來說，這些方法都是基于模態(tài)數(shù)據(jù)的截?cái)嗵├占?jí)數(shù)展開作未知參數(shù)的函數(shù)建立靈敏度矩陣方程δz=Sδθ，其中δθ是參數(shù)擾動(dòng)，δz是測量輸出的擾動(dòng)，S是靈敏度矩陣。不同組的修正參數(shù)可能會(huì)使識(shí)別問題成為病態(tài)的或數(shù)值不穩(wěn)定的。殘差型模型修正方法是構(gòu)造含有解析的修正參數(shù)的函數(shù)與它的實(shí)測數(shù)據(jù)之間的殘差向量，利用各種尋優(yōu)算法，取得殘差為0時(shí)的設(shè)計(jì)參數(shù)的最優(yōu)解，其中修正參數(shù)關(guān)于設(shè)計(jì)參數(shù)的函數(shù)關(guān)系的顯化及優(yōu)化問題的求解是其主要難題。靈敏度矩陣方程型及殘差型模型修正的結(jié)果具有明確的物理意義，便于實(shí)際結(jié)構(gòu)的分析與計(jì)算，并與其他優(yōu)化設(shè)計(jì)過程兼容，實(shí)用性強(qiáng)。本文將對基于殘差法建立數(shù)學(xué)模型的模型修正方法進(jìn)行研究。

1 實(shí)模態(tài)參數(shù)與復(fù)模態(tài)參數(shù)[4-5]

描述自由度為N的線性阻尼離散系統(tǒng)的自由振動(dòng)方程為

(1)

相應(yīng)地其強(qiáng)迫振動(dòng)方程為

(2)

式中M，C和K分別為系統(tǒng)的質(zhì)量、阻尼和剛度矩陣。令C=0，則無阻尼固有頻率與實(shí)模態(tài)對(λi，vi)(i=1，2，…，N，λ=w2=-ω2)滿足方程

(K+λiM)vi=0。

(3)

可以計(jì)算得到實(shí)頻率λ1，λ2，…，λN及實(shí)模態(tài)v1，v2，…，vN。常用的規(guī)范正交化條件為

VTMV=E，VTKV=diag(-λ1，…，-λN)，

(4)

其中V=[v1，…，vN]。實(shí)頻率及實(shí)模態(tài)，有時(shí)也包括阻尼比，統(tǒng)稱為實(shí)模態(tài)參數(shù)。

考慮阻尼時(shí)的系統(tǒng)復(fù)頻率及復(fù)模態(tài)對(si，ui)(i=1，2，…，2N)滿足方程

(5)

它們的規(guī)范正交條件的形式有很多[6]，常用的形式為

ΦTAΦ=E，ΦTBΦ=diag(-s1，…，-s2N)，

(6)

其中Φ=[φ1，φ2，…，φ2N]，φi=[uiλiui]T(i=1，2，…，2N)，且

(7)

下文討論中用到的均為經(jīng)過式(4)和式(6)規(guī)范正交化后的實(shí)模態(tài)和復(fù)模態(tài)。

2 定義殘差法的目標(biāo)函數(shù)

minΔR=f(p)-f(p(o))，

(8)

尋找p的最優(yōu)解p*。

式(8)是一個(gè)無約束非線性優(yōu)化模型，從數(shù)學(xué)角度來說，它有眾多解法，人們除了使用最速下降法外，還使用并發(fā)展了Newton法及其改進(jìn)算法，例如阻尼Newton法，帶保護(hù)措施的阻尼Newton法，吉爾·默里穩(wěn)定Newton法等。盡管收斂速度較慢，但最速下降法最優(yōu)質(zhì)的特點(diǎn)是結(jié)構(gòu)簡單且計(jì)算代價(jià)較??；相較而言，Newton法及其改進(jìn)Newton法收斂速度較快，盡管在每一步迭代時(shí)，都需計(jì)算目標(biāo)函數(shù)的二階導(dǎo)數(shù)來構(gòu)成海森矩陣，確定搜索方向，然后再求解一個(gè)線性方程組，計(jì)算工作量大，從而抵消了Newton法收斂速度快的優(yōu)點(diǎn)。

由Daviden發(fā)明的變尺度法即擬Newton法，是無約束優(yōu)化計(jì)算方法中最杰出的、最富有創(chuàng)造性的工作。不同于Newton法，擬Newton法在構(gòu)造搜索方向時(shí)避免了計(jì)算海森矩陣，減少了計(jì)算量，且具有超線性收斂速度。Broyden、Powell和Fletcher等人沿著這一方向做了大量的工作，提出了PSB公式、BFGS公式以及DFP公式等一組著名的校正公式，形成了Broyden族擬Newton法。1970年Huang提出了一類比Broyden族擬Newton法更為廣泛的校正公式，即Huang族擬Newton法，它取消了Broyden族對修正矩陣的對稱性限制。Powell直接方法和共軛梯度法在無約束最優(yōu)化計(jì)算方法中也占有十分重要的地位，由于Broyden族擬Newton法對于中小規(guī)模的無約束優(yōu)化問題，是一類非常有效的方法，但在求解過程中，它要存儲(chǔ)一個(gè)n×n階矩陣，這對大規(guī)模無約束問題來說，會(huì)因?yàn)橛?jì)算機(jī)內(nèi)存的限制而無法使用。再則，它需要不停地求解線性方程組，因此，算法的效率受到了影響。高效且需較少存儲(chǔ)量的求解大規(guī)模無約束最優(yōu)化算法，即為共軛梯度法。例如，對正定二次函數(shù)的共軛梯度法有Polak-Ribiere-Rolyak共軛梯度法、Fletcher-Reeves共軛梯度法、Dixon共軛梯度法和Crowder-Wolfe共軛梯度法等。對于正定二次函數(shù)來說，如果以上幾種共軛梯度法均采用精確線性搜索，則它們是等價(jià)的。但在實(shí)際中，最常用的是FR方法和PR方法。由于FR方法中的βk計(jì)算最簡單，因此應(yīng)用得最多。而對于非二次函數(shù)的共軛梯度法，將共軛方向法應(yīng)用于非二次函數(shù)的極小化，每迭代n次就周期地取負(fù)梯度方向作為搜索方向，稱這種算法為重新開始的共軛梯度法。Beale-Powell提出在當(dāng)前點(diǎn)重新開始時(shí)，以算法產(chǎn)生的方向作為重新開始的搜索方向，即要求新一輪的搜索方向?qū)Χ魏瘮?shù)仍是共軛方向，即Beale-Powell三次重新開始共軛梯度法。

在一些無約束最優(yōu)化問題中，由于目標(biāo)函數(shù)的特殊性，它的海森矩陣G(x)是稀疏矩陣，若用通常的擬牛頓法求解這類問題，所構(gòu)造的用于逼近海森矩陣的矩陣Bk不具有與G(x)相同的稀疏結(jié)構(gòu)，因此，形成了一般稀疏擬牛頓法和基于BFGS校正公式的稀疏擬牛頓法。到目前為止，擬牛頓法中的BFGS算法是最有效的求解無約束優(yōu)化問題的方法。將其改進(jìn)，減少其存儲(chǔ)量，使之能適用于求解大規(guī)模無約束最優(yōu)化問題，即為無記憶擬牛頓法，它在計(jì)算過程中不需要存儲(chǔ)矩陣Bk，具有內(nèi)存小、收斂速度快的特點(diǎn)，特別適用于求解大規(guī)模無約束優(yōu)化問題。

但在實(shí)際問題中，由于具體建模所采用的修正參數(shù)不同，模型也不同，所選擇采用的算法也可能更加適用和靈活。

2.1 頻率殘差

模型修正的目的就是以t個(gè)修正對象的實(shí)測數(shù)據(jù)作為支撐來優(yōu)化得到設(shè)計(jì)參數(shù)的最優(yōu)值作為修正值。接下來，依據(jù)不同的修正參數(shù)，建立下列幾種殘差型的優(yōu)化模型。

將模型修正的目標(biāo)函數(shù)[8-9]定義為

(9)

式中：p1、p2為設(shè)計(jì)參數(shù)的上下限值；t為實(shí)測數(shù)據(jù)的個(gè)數(shù)，

D(p)=(D1(p)D2(p) …Dt(p))T

為t個(gè)頻率的殘差向量，其中

Di(p)=λi(p)-λi(p(o))。

將模型修正的目標(biāo)函數(shù)[2，10-12]定義為

(10)

式中：p1、p2為設(shè)計(jì)參數(shù)的上下限值；t為實(shí)測數(shù)據(jù)的個(gè)數(shù)。

Df(p)=(D1(p)D2(p) …Dt(p))

為t個(gè)頻率的殘差向量，其中

式中：λi(p)為第i階無阻尼固有頻率關(guān)于設(shè)計(jì)參數(shù)的顯式函數(shù)關(guān)系；λi(p(o))為第i階無阻尼固有頻率的實(shí)測數(shù)據(jù)值。

將模型修正的目標(biāo)函數(shù)[13]定義為

(11)

式中：W為權(quán)矩陣；p1、p2為設(shè)計(jì)參數(shù)的上下限值。

2.2 模態(tài)殘差

將模型修正的目標(biāo)函數(shù)[14]定義為

(12)

式中：W為權(quán)矩陣；p1、p2為設(shè)計(jì)參數(shù)的上下限值。

為t個(gè)頻率的殘差向量，且

將模型修正的目標(biāo)函數(shù)[15-16]定義為模態(tài)MAC值的殘差形式：

(13)

其中

式中：vi(p)為有限元模型第i階模態(tài)振型關(guān)于設(shè)計(jì)參數(shù)的顯式函數(shù)關(guān)系；vi(p(o))為第i階模態(tài)振型的實(shí)測數(shù)據(jù)值。

有時(shí)模態(tài)殘差也定義為如下形式：

(14)

但這種形式應(yīng)用得還不夠廣泛。

2.3 頻率和模態(tài)的綜合殘差

將模型修正的目標(biāo)函數(shù)[17-21]定義為頻率殘差與模態(tài)殘差的加權(quán)綜合式：

minf(p)=Df(p)WfDf(p)+r(p)TWmr(p)，

(15)

式中Wf、Wm為頻率殘差和模態(tài)殘差的權(quán)矩陣。

還可以將模型修正的目標(biāo)函數(shù)[22-23]定義為頻率殘差與模態(tài)殘差如式(16)的加權(quán)綜合式

(16)

式中W為權(quán)矩陣。

2.4 方程殘差

將模型修正的目標(biāo)函數(shù)定義為實(shí)測頻率數(shù)據(jù)和實(shí)測模態(tài)數(shù)據(jù)同時(shí)滿足特征方程，來獲得特征方程中的剛度或質(zhì)量矩陣的修正量。以剛度陣的修正為例，假設(shè)修正后的模型的剛度矩陣可以用初始模型中的各單元的剛度陣線性表示，即

式中：αi為第i個(gè)單元的剛度矩陣的修正系數(shù)；Koi為原始模型中第i個(gè)單元的剛度矩陣；Ne為結(jié)構(gòu)模型的單元數(shù)。目標(biāo)函數(shù)[24]定義為

(17)

3 對殘差法目標(biāo)函數(shù)的尋優(yōu)策略

在上文建立的10個(gè)優(yōu)化模型中，除了式(13)～(17)是無約束非線性規(guī)劃外，式(8)～(12)均為約束非線性規(guī)劃模型。針對約束優(yōu)化模型，求解約束最優(yōu)化問題的方法大致可分為3類：

第1類是將約束問題轉(zhuǎn)化為一系列無約束問題求解，用這一系列無約束問題的極小點(diǎn)去逼近原約束問題的最優(yōu)解——序列無約束極小化方法(Sequnential Unconstrained Minimization Techniqe)，簡稱SUMT方法，有代表性的是懲罰函數(shù)法(penalty function method)和廣義拉格朗日乘子法(Generalized Lagrange Multiplier Method)。這種方法主要是各種不同形式、不同特性的罰函數(shù)法，依賴于如何將目標(biāo)函數(shù)與約束函數(shù)進(jìn)行組合，導(dǎo)出不同形式的罰函數(shù)。只要恰當(dāng)?shù)貙?shí)施，這類方法，就可擁有較好的收斂速度與數(shù)值穩(wěn)定性，一直是求解約束優(yōu)化問題的主要方法。

第2類是在每個(gè)迭代點(diǎn)處構(gòu)造一個(gè)二次函數(shù)去逼近目標(biāo)函數(shù)，用線性函數(shù)逼近約束函數(shù)。將對原約束問題的求解轉(zhuǎn)化為對一系列二次規(guī)劃子問題的求解，以子問題的解作為本次迭代的搜索方向，沿搜索方向?qū)?yōu)得到新迭代點(diǎn)，迭代序列最終逼近原約束問題的最優(yōu)解——序列二次規(guī)劃法(Sequential Quadratic Programming，常簡記為SQP)方法。各種方法差異主要在于構(gòu)造二次函數(shù)的海森矩陣的方法不同。特別地，若以變尺度法(擬Newton法)構(gòu)造海森陣，形成二次規(guī)劃子問題——約束變尺度法，代表性方法為WHP (abbr.water horsepower)方法。后經(jīng)不斷改進(jìn)，SQP類方法已成為用于求解非線性約束最優(yōu)化問題最有效的算法，它不僅可以求解帶有等式約束的優(yōu)化問題，而且可以很容易地處理不等式約束，這類算法不僅可以具有全局收斂性，而且具有局部超線性收斂或更強(qiáng)的二次收斂性。

第3類是處理線性約束問題的方法在非線性約束時(shí)的推廣，如廣義消去法、可行方向法、廣義簡約梯度法和投影梯度法等。這些方法所產(chǎn)生的迭代點(diǎn)均是約束優(yōu)化問題的可行點(diǎn)。但由于約束函數(shù)的非線性屬性，這些方法的實(shí)施要比處理線性約束問題時(shí)所使用的相應(yīng)算法復(fù)雜得多，因此，算法的效果不算理想。代表性的是可行方向法(feasible direction method)，它是一種直接處理約束的方法。它的特點(diǎn)是在迭代過程中目標(biāo)函數(shù)要求沿著可行下降方向?qū)?yōu)，主要的方法是佐登第克(Zoutendijk)可行方向法、托克斯—凡努特(Topkis-Veinotl)可行方向法、投影(projection)梯度法、簡約梯度法等。

但在模型修正的應(yīng)用過程中，遠(yuǎn)沒有達(dá)到能廣泛使用約束優(yōu)化方法的程度，下面僅介紹在模型修正過程中建立殘差型目標(biāo)函數(shù)后的常用的優(yōu)化算法。

3.1 穩(wěn)定性條件方程求解

用穩(wěn)定性條件

(18)

求解式(17)，然后得到穩(wěn)定點(diǎn)，認(rèn)為是最優(yōu)點(diǎn)。顯然式(18)的求解變成一個(gè)超定或欠定方程組的求解，引入廣義逆技術(shù)求解[24]。

3.2 優(yōu)化算法求解

在不考慮約束條件，直接用共軛梯度法中的PR算法，求解形如式(9)的無約束優(yōu)化模型[8]，或自適應(yīng)模擬退火優(yōu)化算法求解[11]。對形如式(10)的非線性約束優(yōu)化問題，采用有限元分析軟件ANSYS先后用兩種優(yōu)化方法：零階和一階方法求解[2]，或用序列二次規(guī)劃法(SQP)求解[10]。對形如式(11)的優(yōu)化問題采用混合優(yōu)化算法求解[13]。對形如式(12)的約束優(yōu)化問題，采用罰函數(shù)法化為無約束優(yōu)化問題，對該無約束問題用最速下降法求解[14]。對形如式(13)的優(yōu)化問題，可以采用遺傳算法求解[16]。對于形如式(15)和式(16)的優(yōu)化問題，通常采用的優(yōu)化算法是基于靈敏度的迭代程序，在信賴域內(nèi)優(yōu)化目標(biāo)函數(shù)，利用目標(biāo)函數(shù)f(p)的泰勒級(jí)數(shù)截?cái)嘈问娇梢远x一個(gè)二次型

z(p)=f(p)+2f(p)]Δp。

Δp在迭代算法中代表一個(gè)步長向量，f(p)與[2f(p)]是f(p)的梯度與海森矩陣，經(jīng)過某次迭代過程，使f(p)≈0，從而得到近似的最優(yōu)值p*。在每次迭代中，優(yōu)化算法都在當(dāng)前點(diǎn)附近，構(gòu)建一個(gè)模型函數(shù)z(p)，并確定一個(gè)圍繞當(dāng)前點(diǎn)的信賴域，優(yōu)化值p*可以在信賴域中獲得，在每步迭代中搜索步長都限定在信賴域內(nèi)，避免了過大步長的出現(xiàn)。目前常用的算法有信賴域牛頓算法[19，21]或高斯牛頓算法[17]等。

4 結(jié)語

我們實(shí)際模型修正問題中遇到的目標(biāo)函數(shù)，常常利用非線性最小二乘法來建立殘差型優(yōu)化模型。對待非線性最小二乘問題，若采用牛頓類算法，則在迭代點(diǎn)取目標(biāo)函數(shù)的二次近似式的最優(yōu)解作為搜索方向或修正向量。由于是對殘差在平方和意義下取極小，極有可能殘差在最優(yōu)解處的值為零(零殘量問題)，或較小的值(小殘量問題)海森矩陣G(x)中的非線性項(xiàng)為0，或者相對線性項(xiàng)來說較小，因此，在最優(yōu)值的某個(gè)鄰域內(nèi)G(x)的線性項(xiàng)作為海森矩陣G(x)的一個(gè)很好的近似，即為著名的高斯—牛頓法。如果δ(k)為f(x)在點(diǎn)x(k)處的一個(gè)下降方向，取d(k)=δ(k)為搜索方向進(jìn)行線性搜索，并令x(k+1)=x(k)+αkd(k)即為阻尼高斯—牛頓法，其中αk為由線性搜索確定的步長。

布朗與丹尼斯(Brown，Dennis，1971)，貝茨(Betts，1976)，丹尼斯，蓋伊和韋爾舒(Dennis，Gay and Welsch，1981)以及吉爾—默里(Gill-Murray，1978)等提出用擬牛頓修正產(chǎn)生對海森矩陣G(x)中的非線性項(xiàng)的近似。上述幾種產(chǎn)生矩陣Ck去形成海森矩陣G(x)的近似矩陣Bk的方法，在實(shí)際計(jì)算時(shí)，有時(shí)效果不如預(yù)期的好。原因在于Ck的近似程序不夠理想，且隨著迭代過程的進(jìn)行，誤差會(huì)越來越大?？朔@一困難的一個(gè)可取辦法是采用調(diào)比技術(shù)，即在對矩陣Ck作修正之前，先選擇一個(gè)調(diào)比因子來建立修正公式?；旌纤惴?又稱雜交算法)是一類以高斯—牛頓法為基礎(chǔ)的算法?；旌纤惴ò迅咚埂ｎD法同一般的無約束最優(yōu)化方法進(jìn)行組合，并能在求解問題的過程中自動(dòng)地調(diào)節(jié)到適合問題特性的方法，因而也稱為開關(guān)算法。分解式擬牛頓法是另一類即利用非線性最小二乘問題的特點(diǎn)，又利用擬牛頓修正產(chǎn)生目標(biāo)函數(shù)海森矩陣近似的方法。以上非線性最小二乘法在模型修正中的應(yīng)用還不普遍，但它卻是最適合殘差型模型修正的算法之一。

另外，可行內(nèi)點(diǎn)弧算法(feasible arc interior point algorithm)[25]對非線性不等式和等式約束的優(yōu)化問題來說，是一種新的技術(shù)。它需要在不等約束的內(nèi)點(diǎn)中得到一個(gè)初始點(diǎn)，然后產(chǎn)生一系列的內(nèi)點(diǎn)。當(dāng)問題僅僅是含有不等約束，我們的模型用連接變量和等式約束來引入多學(xué)科帶來的影響，就像在典型的優(yōu)化模型中，狀態(tài)方程能得到隱式處理，從另一方面說，在同步分析和設(shè)計(jì)優(yōu)化問題中，狀態(tài)方程可以看成是那些規(guī)劃的一部分?？尚袃?nèi)點(diǎn)弧算法減少了狀態(tài)方程所給出的等式約束中的狀態(tài)變量，最新穎的一點(diǎn)是這種減少是通過不需要解狀態(tài)方程來完成的，它闡述了解以線性化平衡方程給出的等式的系統(tǒng)解法，而且它能有效地應(yīng)用在解決工程問題的規(guī)范的解算器上。