劉春剛

(中鐵第一勘察設(shè)計(jì)院集團(tuán)有限公司，西安 710043)

1 最大熵法概率密度函數(shù)的確定和擬合

1.1 概率密度函數(shù)的確定

在概率論中，熵的定義［1］如式(1)所示

式中，S為熵;f(x)為概率密度函數(shù)。最大熵法的理論基礎(chǔ)是最大信息熵原理:在給定所有滿足約束條件的概率密度函數(shù)中，信息熵最大的概率密度函數(shù)為最佳的概率密度函數(shù)。設(shè)x為一個(gè)連續(xù)型隨機(jī)變量，其概率密度函數(shù)f(x)滿足以下條件

式中，Mi為x的i階原點(diǎn)矩，可通過(guò)統(tǒng)計(jì)樣本計(jì)算確定。這樣最大熵法就轉(zhuǎn)化為在式(2)、式(3)約束下求式(1)最大值問(wèn)題。構(gòu)造拉格朗日函數(shù)，如式(4)、式(5)所示

式中，λ0，λ1，λ2，…，λn為待定系數(shù)，由式(2)、式(3)構(gòu)成的非線性方程組求解。

式(7)給出了n階矩法概率密度函數(shù)表達(dá)式，韋征等［3］對(duì)n值的選取問(wèn)題進(jìn)行的研究，以一個(gè)單層網(wǎng)殼結(jié)構(gòu)可靠度分析為驗(yàn)證算例，比較了最大熵理論三階矩算法、四階矩算法及五階矩算法的計(jì)算結(jié)果。通過(guò)對(duì)比發(fā)現(xiàn)三階矩算法得到的概率密度不能正確地反應(yīng)結(jié)構(gòu)位移響應(yīng)變化規(guī)律，與四階矩算法相比，五階矩算法尾部誤差較大，已經(jīng)不能有效地提高概率密度函數(shù)、超越概率和可靠指標(biāo)的精度，并且求解待定系數(shù)的方程組過(guò)于復(fù)雜，采用數(shù)值解法求解方程組時(shí)對(duì)初始值要求嚴(yán)格，不容易得到能使方程收斂的初始值。最終韋征等建議n值取4，即采用四階矩算法。

1.2 擬合檢驗(yàn)

根據(jù)最大熵原理求得函數(shù)只有通過(guò)擬合檢驗(yàn)才能確定其性能的好壞，即能否最大程度反映隨機(jī)變量本身概率特性。由于總體分布未知，在檢驗(yàn)過(guò)程中自然就不會(huì)牽扯到總體分布的參數(shù)，對(duì)于這種總體參數(shù)未知的檢驗(yàn)，被稱為“非參數(shù)檢驗(yàn)”，常用的非參數(shù)檢驗(yàn)方法有K-S檢驗(yàn)法、變量值隨機(jī)性檢驗(yàn)等，以K-S檢驗(yàn)為例對(duì)最大熵法求得概率密度函數(shù)進(jìn)行檢驗(yàn)，其具體過(guò)程如下:

(1)原假設(shè)根據(jù)最大熵法得到的最大概率密度函數(shù)積分求得概率分布函數(shù)F0(X);

(2)利用樣本數(shù)據(jù)計(jì)算個(gè)樣本數(shù)據(jù)點(diǎn)的累計(jì)概率得到檢驗(yàn)累計(jì)概率分布函數(shù)S0(X);

我國(guó)高鐵走出國(guó)家的形勢(shì)雖然良好，有著很大的潛力，但一些挑戰(zhàn)是不可避免的，激烈的競(jìng)爭(zhēng)、國(guó)際形勢(shì)、未知的不利因素都是阻礙高鐵發(fā)展的重要因素。

(3)計(jì)算F0(X)和S0(X)在相應(yīng)的變量值點(diǎn)x上的差D(x)，得到差最大值Dmax;

(4)對(duì)于給定的檢驗(yàn)水平α和樣本容量n，查Kolmogorov分布的分位數(shù)表，得到臨界值Dn;

(5)判斷:若Dmax≤Dn時(shí)，則接受原假設(shè)，否則否定原假設(shè)。

2 概率密度函數(shù)的求解

最大熵法擬合概率密度函數(shù)的核心和難點(diǎn)在于以下非線性方程組的求解

式中，Mi統(tǒng)計(jì)樣本 x的 i階原點(diǎn)矩;λi為待定系數(shù)。

方程組的求解需要解決三方面的問(wèn)題:①各方程積分域的選擇;②方程組簡(jiǎn)化;③方程組求解方法。

2.1 積分域選取

方程組式(8)中各方程都含有一個(gè)不定積分項(xiàng)，直接求解非常困難，需要先把各方程的不定積分計(jì)算出來(lái)，去掉方程中的積分符號(hào)，簡(jiǎn)化方程組。由于被積函數(shù)比較復(fù)雜，需要人為的選擇一個(gè)積分區(qū)域，用被積函數(shù)在該區(qū)域上的定積分代替原不定積分，然后采用數(shù)值方法進(jìn)行積分運(yùn)算。趙國(guó)藩［2］建議積分區(qū)域選為［μ-10σ，μ+10σ］并用梯形積分公式進(jìn)行積分運(yùn)算。通過(guò)算例對(duì)比發(fā)現(xiàn)積分區(qū)域的選取是否合適對(duì)擬合的好壞有重要影響。以均值為25，標(biāo)準(zhǔn)差為11.2極值Ⅰ型分布為例，當(dāng)積分區(qū)域分別選擇［0，400］和［0，80］時(shí)擬合的概率密度函數(shù)如圖(1)、圖(2)所示。

由圖1、圖2對(duì)比發(fā)現(xiàn)，積分區(qū)域?。?，400］時(shí)得不到正確的擬合函數(shù)，擬合的效果很差，最大誤差為0.76。而積分區(qū)域?。?，80］時(shí)，最大誤差量級(jí)為10-2，二者相差76倍。通過(guò)算例驗(yàn)證發(fā)現(xiàn)，積分區(qū)域應(yīng)控制在［μ-10σ，μ+10σ］以內(nèi)為宜，同時(shí)積分區(qū)間端點(diǎn)處對(duì)應(yīng)的概率密度函數(shù)值不宜超過(guò)10-4。

圖1 積分域［0，400］

圖2 積分域［0，80］

2.2 方程組簡(jiǎn)化

觀察第一個(gè)方程，作以下變換

即

即

式(12)代入式(13)得

此時(shí)，最大熵法系數(shù)方程組就由式(8)化簡(jiǎn)為式(14)，即將原五元的方程組簡(jiǎn)化為四元的方程組，提高了方程組收斂速度和穩(wěn)定性。由方程組(14)得到λ1、λ2、λ3、λ44 個(gè)系數(shù)后，按照式(10)求解 λ0，從而得到所有的待定系數(shù)。

2.3 方程組求解方法

非線性方程組(8)比較復(fù)雜，求解起來(lái)比較麻煩，已有許多文獻(xiàn)對(duì)方程組的求解方法進(jìn)行了探索研究。趙國(guó)藩等［2］利用循環(huán)中點(diǎn)求積Newton算法對(duì)非線性方程組牛頓算法進(jìn)行改進(jìn)，改進(jìn)后的新算法具有收斂快、穩(wěn)定、對(duì)初始值無(wú)限制，成功率高的優(yōu)點(diǎn)。此外，還有一些文獻(xiàn)［4-5］應(yīng)用遺傳算法、粒子群算法等優(yōu)化算法求解，都取得了一定的成果。但是這些計(jì)算方法對(duì)數(shù)學(xué)知識(shí)要求較高，并且需要編寫較為復(fù)雜的程序才能實(shí)現(xiàn)。

Matlab具有強(qiáng)大的數(shù)值計(jì)算能力，并且整合了許多優(yōu)化算法、非線性方程組數(shù)字解法程序作為可以直接調(diào)用的內(nèi)部函數(shù)，求解方程組比較方便。在求解實(shí)例時(shí)發(fā)現(xiàn)直接調(diào)用Matlab優(yōu)化工具箱中的遺傳算法程序去求解方程組(14)，不需要設(shè)置初始值，但是得到的解精度并不高，調(diào)用ga函數(shù)求解方程組一般上所能達(dá)到的精度為0.1～0.01，很難得到精度更高的解。在Matlab中還可以用Fsolve函數(shù)求解非線性方程組，該方法可以得到精度為10-7以及更好的解，但是使用Fsolve函數(shù)對(duì)初始解的要求非常嚴(yán)格，初始解選擇不恰當(dāng)將導(dǎo)致程序不收斂。比較上述兩種算法的優(yōu)缺點(diǎn)后，將兩種方法有機(jī)的結(jié)合起來(lái)，先使用遺傳算法得到一個(gè)精度為0.1—0.01的解，然后將這個(gè)解作為方程組的初始解調(diào)用Fsolve函數(shù)得到精度更好的解。這種方法前后兩步都是直接調(diào)用Matlab內(nèi)置函數(shù)，實(shí)現(xiàn)起來(lái)非常簡(jiǎn)單，并且求解成功率很高。本文分別用遺傳算法和遺傳算法與Fsolve相結(jié)合的方法對(duì)四種分布進(jìn)行了擬合，四種分布的參數(shù)見表1。

表1 分布參數(shù)

分別使用上述兩種計(jì)算方法對(duì)上述四種分布進(jìn)行計(jì)算，計(jì)算結(jié)果如圖3～圖10所示，其中圖3～圖6為遺傳算法計(jì)算結(jié)果;圖7～圖10為遺傳算法和Fslove函數(shù)相結(jié)合的計(jì)算結(jié)果。

對(duì)以上模擬結(jié)果的精度進(jìn)行對(duì)比分析，結(jié)果如表2所示。

圖3 分布1遺傳算法擬合結(jié)果

圖4 分布2遺傳算法擬合結(jié)果

圖5 分布3遺傳算法擬合結(jié)果

圖6 分布4遺傳算法擬合結(jié)果

圖7 分布1改進(jìn)算法擬合

圖8 分布2改進(jìn)算法擬合

圖9 分布3改進(jìn)算法擬合

圖10 分布4改進(jìn)算法擬合

表2 各分部計(jì)算精度對(duì)比

由表2可以看出，單純用遺傳算法計(jì)算的結(jié)果精度不太高，概率密度函數(shù)擬合的誤差在0.3～0.1之間，將這樣的擬合概率密度函數(shù)代替真實(shí)密度函數(shù)進(jìn)行可靠度計(jì)算會(huì)帶來(lái)較大的誤差。而采用遺傳算法和Fsolve函數(shù)相結(jié)合的算法得到的擬合函數(shù)最大擬合誤差均在0.03以下，精度較高。其中正態(tài)分布(分布1)擬合精度由原來(lái)的0.25提高至10-5，提高了2 400倍;分布3的最大擬合誤差為0.02，代數(shù)與遺傳算法相比精度提高了9倍。因此用遺傳算法得到擬合函數(shù)的待定系數(shù)后，再以此結(jié)果為初始解用Fsolve進(jìn)行一次求解，可以有效地提高計(jì)算精度。解決方程組(8)求解困難的問(wèn)題，并且該方法直接調(diào)用Matlab內(nèi)置函數(shù)，不需要編寫實(shí)現(xiàn)遺傳算法和非線性方程組數(shù)值解法的程序，實(shí)際操作起來(lái)簡(jiǎn)便易行。

3 樁板墻的極限狀態(tài)方程與優(yōu)化模型

3.1 極限狀態(tài)方程

以路塹式樁板墻為研究對(duì)象，建立其可靠度分析方程，極限狀態(tài)方程的一般形式為

式中，R為抗力;S為荷載。

樁板墻破壞主要是樁發(fā)生破壞，其破壞形式主要有受彎破壞、受剪破壞和錨固段內(nèi)側(cè)土體受壓破壞三種形式，因此根據(jù)這三種形式求其極限方程。

按照式(15)極限狀態(tài)方程的形式把受彎和受剪的極限方程分別表示為

式中，Ωpw、Ωpj分別為受彎和受剪計(jì)算模式的隨機(jī)變量;Mmax、Qmax分別為樁板墻受到的最大彎矩和最大剪力;Mu、Qu分別為樁板墻極限彎矩和極限剪力。

以“m”法為例

式中，A3、A4、B3、B4、C3、C4、D3、D4為關(guān)于樁埋深的函數(shù);x0為滑面處的位移;M0、Q0分別為樁錨固段頂點(diǎn)的彎矩和剪力。對(duì)于式(20)和式(21)分別求一階偏導(dǎo)為零，即可得到Mmax和Qmax表達(dá)式［6］。

3.2 隨機(jī)變量的確定

綜上所述，式(16)、式(17)為失效模式下的偏微分方程。對(duì)于抗力項(xiàng)隨機(jī)變量，在式(16)中主要考慮縱向鋼筋受拉強(qiáng)度f(wàn)yu、縱向鋼筋截面面積As、混凝土軸心抗壓強(qiáng)度f(wàn)cu、樁的截面寬度b和有效高度h0以及分項(xiàng)系數(shù)Ωpw;在式(17)中主要考慮箍筋受拉強(qiáng)度f(wàn)yvu、混凝土抗拉強(qiáng)度f(wàn)tu、樁的截面寬度b和有效高度h0以及分項(xiàng)系數(shù)Ωpj。對(duì)于荷載項(xiàng)由于計(jì)算復(fù)雜性，暫時(shí)直接按照分項(xiàng)系數(shù)法直接計(jì)算。

4 算例驗(yàn)證

以路塹式樁板墻為例，樁長(zhǎng)23 m，樁間距6 m，樁采用2 m×3 m的T形樁，設(shè)計(jì)懸臂端長(zhǎng)8 m，錨固段長(zhǎng)15 m，懸臂端主要位于膨脹土和粉質(zhì)黏土，錨固段主要位于圓粒土、強(qiáng)風(fēng)化片巖和弱風(fēng)化片巖中，局部位于粉質(zhì)黏土中，其中設(shè)計(jì)土體容重取γ=18 kN/m3，綜合φ=28°，內(nèi)力計(jì)算按照“m”法，其中膨脹土取值3 000 kPa/m2，圓粒土和強(qiáng)風(fēng)化片巖取值12 000 kPa/m2，弱風(fēng)化片巖120 000 kPa/m，土壓力分布按照三角形計(jì)算。樁底支撐采用自由端。

對(duì)于荷載的計(jì)算，利用鐵一院研發(fā)的路基助手計(jì)算，其荷載分項(xiàng)系數(shù)按照K=1.35考慮，對(duì)于材料和樁的結(jié)構(gòu)相關(guān)參數(shù)如表3所示。

表3 計(jì)算資料

對(duì)于彎矩利用最大熵原理求解得到，γ0=-0.919，γ1=0.020，γ2=-0.500，γ3=-3.03，γ4=0，失效概率為3.34e-8，可靠指標(biāo) β=5.399;同樣對(duì)于剪力 γ0=-0.743，γ1=0.031，γ2=-0.459，γ3=-2.975，γ4=0.04，失效概率為7.95e-7，可靠指標(biāo)β=4.864。

5 結(jié)論

(1)通過(guò)選擇方程的積分區(qū)域和對(duì)方程組的簡(jiǎn)化，然后利用matlab中Fsolve函數(shù)和ga工具箱相結(jié)合的方式，實(shí)現(xiàn)了對(duì)基于最大熵函數(shù)概率密度函數(shù)方程組簡(jiǎn)單和精確求解。

(2)通過(guò)對(duì)路塹式樁板墻計(jì)算，然后利用熵函數(shù)分別得到受彎、受剪和錨固段樁側(cè)土體可靠指標(biāo)均值，而且一般只需要較少樣本即可得到相對(duì)精確的概率密度函數(shù)和可靠指標(biāo)。

(3)樁板墻的破壞屬于延性破壞，按照文獻(xiàn)［7］的規(guī)定，可靠度標(biāo)準(zhǔn)不低于3.2，所以計(jì)算出來(lái)的彎矩可靠指標(biāo)5.399和剪力可靠指標(biāo)4.864均滿足要求。

［1］ E.T.JaYNES.Information Theory and Stastical Mechanics［J］.The Physical Review，1957，106(104):620-630.

［2］ 趙國(guó)番.工程結(jié)構(gòu)可靠性理論與應(yīng)用［M］.大連理工大學(xué)出版社，1996:137-142.

［3］ 韋征，葉繼紅，沈士釗.最大熵法可靠度理論在工程中的應(yīng)用［J］.振動(dòng)與沖擊，2007，26(6):146-151.

［4］ 謝世坤，段芳，李強(qiáng)征，等.非線性方程組求解的三種Newton法比較［J］.井岡山學(xué)院學(xué)報(bào)(自然科學(xué))，2006，27(8):8-11.

［5］ 陳長(zhǎng)憶，葉永春.基于粒子群算法的非線性方程組求解［J］.計(jì)算機(jī)應(yīng)用與軟件，2006，23(5):137-139.

［6］ 新型支擋結(jié)構(gòu)設(shè)計(jì)與工程實(shí)例［M］.北京:人民交通出版社，2010:275-296.

［7］ 中華人民共和國(guó)建設(shè)部.GB50068—2001 建筑結(jié)構(gòu)可靠度統(tǒng)一標(biāo)準(zhǔn)［S］.北京:中國(guó)建筑工業(yè)出版社，2001.