如何利用Matlab的lsqnonlin函數(shù)處理復(fù)雜函數(shù)的擬合問(wèn)題
1. 簡(jiǎn)介
在工程問(wèn)題的計(jì)算中，我們經(jīng)常需要處理一些離散數(shù)據(jù)的擬合問(wèn)題，而最小二乘法是處理曲線擬合問(wèn)題的常用方法。目前，許多軟件都提供有基于最小二乘法進(jìn)行曲線擬合的功能，例如在Origin和Excel中均可直接利用離散數(shù)據(jù)進(jìn)行曲線擬合。然而，這些軟件只能處理一些簡(jiǎn)單函數(shù)的擬合問(wèn)題，當(dāng)需要擬合的函數(shù)較為復(fù)雜時(shí)，或者無(wú)法用簡(jiǎn)單的表達(dá)式來(lái)表述時(shí)，則往往無(wú)法直接進(jìn)行擬合。為此，本文將對(duì)最小二乘法的基本原理做簡(jiǎn)單介紹，隨后介紹如何利用Matlab的lsqnonlin函數(shù)處理復(fù)雜函數(shù)的擬合問(wèn)題。
1.1 曲線擬合的最小二乘法原理
利用最小二乘法進(jìn)行曲線擬合的本質(zhì)為尋找某個(gè)近似函數(shù) φ ( x )，使得該函數(shù)與離散點(diǎn)之間盡可能地逼近。若將偏差定義為近似函數(shù)的近似值 φ ( xi )與離散點(diǎn) yi *之間的差值：
求解上述線性方程組即可得到擬合多項(xiàng)式的系數(shù)。
2. 利用Matlab處理曲線擬合問(wèn)題
基于上述計(jì)算原理，Matlab提供了polyfit函數(shù)用于處理多項(xiàng)式曲線的擬合問(wèn)題，對(duì)于一些較為復(fù)雜但仍可通過(guò)簡(jiǎn)單表達(dá)式進(jìn)行表述的函數(shù)，也可以利用Matlab的擬合工具箱（Curve fitting）進(jìn)行擬合。但在某些情況，當(dāng)擬合函數(shù)非常復(fù)雜，以致于無(wú)法用簡(jiǎn)單表達(dá)式進(jìn)行表述時(shí)（例如分段函數(shù)以及涉及到條件語(yǔ)句），則無(wú)法使用擬合工具箱進(jìn)行擬合。對(duì)于此類(lèi)問(wèn)題，可以使用Matlab優(yōu)化工具箱中的lsqnonlin函數(shù)進(jìn)行解決。
2.1 lsqnonlin函數(shù)
lsqnonlin函數(shù)用于求解以下述形式表示的非線性最小二乘法擬合問(wèn)題：
在使用該函數(shù)進(jìn)行最小二乘法擬合時(shí)，lsqnonlin函數(shù)并不需要用戶(hù)提供min || f ( x )||（平方和），而是需要用戶(hù)提供自定義函數(shù)fun，用于計(jì)算矢量形式表示的 f ( x )：
lsqnonlin函數(shù)常用語(yǔ)法為：
x = lsqnonlin(fun,x0)
x = lsqnonlin(fun,x0,lb,ub)
其中fun為用戶(hù)自定義函數(shù)，x0為計(jì)算采用的初始值，lsqnonlin函數(shù)首先利用x0通過(guò)自定義函數(shù)fun計(jì)算 fi (x)的取值并計(jì)算平方和，隨后通過(guò)優(yōu)化算法調(diào)整x的取值直至得到平方和的最小值。此外，lb和ub還可以用于定義x的取值范圍，使得x滿(mǎn)足lb≤x≤ub。
例如，對(duì)于節(jié)1.1中所述的多項(xiàng)式，根據(jù)最小二乘法的定義，則自定義函數(shù) f ( x )應(yīng)表示為：
注意此時(shí) f ( x )中的x為以向量形式表示的多項(xiàng)式 P ( x )的系數(shù)：
在計(jì)算時(shí)，用戶(hù)需要指定多項(xiàng)式系數(shù)的初始值，則lsqnonlin函數(shù)將利用最小二乘法計(jì)算多項(xiàng)式系數(shù)。
下面，本文將以筆者所在領(lǐng)域常用的NASGRO方程為例，介紹如何利用lsqnonlin函數(shù)處理此類(lèi)復(fù)雜函數(shù)的曲線擬合問(wèn)題。
2.2 NASGRO方程簡(jiǎn)介
在進(jìn)行基于斷裂力學(xué)的損傷容限分析時(shí)，應(yīng)力強(qiáng)度因子和裂紋擴(kuò)展速率模型是最為重要的輸入。一般來(lái)說(shuō)，應(yīng)力強(qiáng)度因子可以通過(guò)經(jīng)驗(yàn)公式或數(shù)值方法進(jìn)行計(jì)算，而裂紋擴(kuò)展速率模型則需要通過(guò)裂紋擴(kuò)展速率試驗(yàn)獲得的試驗(yàn)數(shù)據(jù)擬合得到。例如，大量的試驗(yàn)結(jié)果表明，在裂紋擴(kuò)展的中速率區(qū)域，應(yīng)力強(qiáng)度因子幅值ΔK和裂紋擴(kuò)展速率d a /dN滿(mǎn)足良好的對(duì)數(shù)線性關(guān)系，可以通過(guò)Paris公式進(jìn)行描述：
其中C和m為材料常數(shù)。
盡管Paris公式已經(jīng)得到廣泛的應(yīng)用，但是Paris公式僅僅描述了裂紋在中速率區(qū)域的擴(kuò)展行為，沒(méi)有描述近門(mén)檻區(qū)域和接近斷裂的高速率區(qū)域的擴(kuò)展行為，也沒(méi)有考慮應(yīng)力比R和裂紋閉合效應(yīng)對(duì)裂紋擴(kuò)展速率的影響，因此給出的計(jì)算結(jié)果將過(guò)于保守。另一個(gè)常用的裂紋擴(kuò)展速率模型為Newman提出的NASGRO模型，該模型基于Forman模型改進(jìn)了裂紋擴(kuò)展速率模型，同時(shí)比Paris和Walker模型更加全面，不僅考慮了應(yīng)力強(qiáng)度因子門(mén)檻值和斷裂韌度，還體現(xiàn)了應(yīng)力比以及裂紋閉合效應(yīng)對(duì)裂紋擴(kuò)展速率d a /dN的影響，如圖2.1所示，其表達(dá)式如下：
其中R為應(yīng)力比，ΔK為應(yīng)力強(qiáng)度因子幅值，ΔKth為應(yīng)力強(qiáng)度因子幅值門(mén)檻值，Kmax為最大應(yīng)力強(qiáng)度因子，可表示為：
圖2.1 NASGRO方程
NASGRO方程中的應(yīng)力強(qiáng)度因子門(mén)檻值ΔKth可采用下面的經(jīng)驗(yàn)公式進(jìn)行估算：
其中A0為裂紋張開(kāi)函數(shù)中的多項(xiàng)式系數(shù)，ΔK1是 R =1時(shí)的應(yīng)力強(qiáng)度因子門(mén)檻值，Cth是對(duì)于正應(yīng)力(上標(biāo)為p=positive)和負(fù)應(yīng)力比(上標(biāo)為m=minus, negative)取不同值的材料常數(shù)，*a*~0~是內(nèi)在小裂紋尺寸(典型值為0.0381mm)。在基于NASGRO方程開(kāi)發(fā)的疲勞裂紋擴(kuò)展分析軟件NASGRO中，正應(yīng)力比下*C*~th~^p^和Δ*K*~1~是保存在數(shù)據(jù)庫(kù)里的值，負(fù)應(yīng)力比下*C*~th~^m^的默認(rèn)值為0.1。
2.3 NASGRO方程擬合
圖2.2為疲勞裂紋擴(kuò)展分析軟件NASGRO材料庫(kù)中某鋁合金材料的裂紋擴(kuò)展速率數(shù)據(jù)，已知試驗(yàn)時(shí)采用的試樣為中心平板試樣(M(T))，σmax和σF的比值為0.3，塑性約束因子α為2.0，材料斷裂韌度Kc為65.7，應(yīng)力比R為1時(shí)的門(mén)檻值ΔK1為1.23，C th ^p^為1.06，C th ^m^為0.1，下面需要通過(guò)擬合試驗(yàn)數(shù)據(jù)獲得NASGRO方程的參數(shù) C , m , p , q 。
圖2.2 疲勞裂紋擴(kuò)展數(shù)據(jù)
擬合NASGRO方程的難點(diǎn)主要有以下幾點(diǎn)：
(1)裂紋擴(kuò)展速率d a /dN不僅與應(yīng)力強(qiáng)度因子幅值ΔK有關(guān)，還與使用的應(yīng)力比R有關(guān)，因此實(shí)際上為多變量的擬合問(wèn)題；
(2)裂紋張開(kāi)函數(shù)f為分段函數(shù)，并且使用了計(jì)算最大值的max函數(shù)，該函數(shù)在擬合時(shí)無(wú)法用簡(jiǎn)單函數(shù)進(jìn)行表述。
針對(duì)以上問(wèn)題，NASGRO軟件給出的擬合方法為首先給參數(shù)p和q確定一個(gè)初始值，并利用最小二乘法確定參數(shù)C和 m ，隨后根據(jù)工程經(jīng)驗(yàn)來(lái)獲得可接受的結(jié)果，如果對(duì)擬合效果不滿(mǎn)意，可以調(diào)整任意參數(shù)，直至獲得滿(mǎn)意的結(jié)果。
顯然，這樣的擬合策略具有很大的隨意性，如果參數(shù)p和q選取不當(dāng)，很可能對(duì)擬合效果有很大的影響。下面，本文將介紹如何利用lsqnonlin函數(shù)在不提前定義參數(shù)p和q的情況下對(duì)NASGRO方程進(jìn)行擬合。
根據(jù)lsqnonlin函數(shù)的介紹，首先需要構(gòu)造自定義函數(shù) f ( x )使其滿(mǎn)足最小二乘法計(jì)算的基本原理，由于Paris公式具有對(duì)數(shù)線性的關(guān)系，因此嘗試將NASGRO方程兩邊取對(duì)數(shù)，可得：
上式可以用如下所示的通式表示：
系數(shù)bj為與 C , n , p和q有關(guān)(b 0 =log( C ), b 1 = n , b 2 = p , b 3 =- q )的系數(shù)，而gj為與Δ K 、R和NASGRO中所有剩余參數(shù)有關(guān)的函數(shù)。
根據(jù)最小二乘法的定義，應(yīng)選取參數(shù)bj使得：
參考lsqnonlin函數(shù)對(duì)目標(biāo)函數(shù)的定義，則自定義函數(shù) f ( x )應(yīng)表示為：
y= R .^2 * (R >0) + R * (R <= 0)
而裂紋張開(kāi)函數(shù)f中涉及到求取最大值的計(jì)算以及分段函數(shù)的處理，也可以通過(guò)上述語(yǔ)法實(shí)現(xiàn)，具體的計(jì)算過(guò)程可參見(jiàn)程序代碼（參見(jiàn)附錄）。
此外，由于自定義函數(shù) f ( x )為關(guān)于系數(shù) bj ( j =1,2,3,4)的函數(shù)，為了將試驗(yàn)數(shù)據(jù)（不同應(yīng)力比R下的應(yīng)力強(qiáng)度因子幅值ΔK和裂紋擴(kuò)展速率d a /d N ）傳遞到函數(shù) f ( x )中進(jìn)行計(jì)算，可以將試驗(yàn)數(shù)據(jù)定義為全局變量，以便被 f ( x )調(diào)用。
通過(guò)編寫(xiě)程序，可以計(jì)算得到NASGRO方程的系數(shù)如表1所示。
擬合曲線與試驗(yàn)數(shù)據(jù)如圖2.3所示。
圖2.3 試驗(yàn)數(shù)據(jù)及擬合曲線對(duì)比
附錄1 NASGRO方程曲線擬合程序
NASGRO_LSQ.m
NASGRO_LSQ用于定義采用最小二乘法擬合NASGRO方程時(shí)的自定義函數(shù) f ( x )，輸入?yún)?shù)Coeff為NASGRO方程系數(shù) bj ，輸出參數(shù)為擬合函數(shù)與試驗(yàn)數(shù)據(jù)誤差的平方和。
function F=NASGRO_LSQ(Coeff)
%程序用于計(jì)算最小二乘法擬合NASGRO方程的目標(biāo)函數(shù)
%程序返回一個(gè)N×1的數(shù)值，其中N為數(shù)據(jù)對(duì)的個(gè)數(shù)
%Coeff為擬合時(shí)待求的系數(shù)（共4個(gè)系數(shù)）
%4個(gè)系數(shù)分別為log(C)、n、p和-q
%DataX(:,1)為應(yīng)力強(qiáng)度因子幅值
%DataX(:,2)為應(yīng)力比
%DataY為裂紋擴(kuò)展速率
%*********全局變量傳遞**************
global S_max_flow alpha DKth1 a0 Cth_p Cth_m a Kc
global DataX DataY
%S_max_flow為施加的最大應(yīng)力與流動(dòng)應(yīng)力的比值
%alpha為塑性約束因子
%DKth1為應(yīng)力比為1時(shí)對(duì)應(yīng)的門(mén)檻值
%a0為與門(mén)檻值有關(guān)的常數(shù)
%Cth_p為正應(yīng)力比下與門(mén)檻值有關(guān)的常數(shù)
%Cth_m為負(fù)應(yīng)力比下與門(mén)檻值有關(guān)的常數(shù)
%a為計(jì)算門(mén)檻值時(shí)使用的裂紋長(zhǎng)度，建議取為遠(yuǎn)大于a0的值
%Kc為材料斷裂韌度
%DataX為應(yīng)力強(qiáng)度因子幅值
%DataY為裂紋擴(kuò)展速率
%***********************************
%********Newman裂紋張開(kāi)函數(shù)計(jì)算**********
R=DataX(:,2); %應(yīng)力比
DK=DataX(:,1); %應(yīng)力強(qiáng)度因子幅值
%計(jì)算系數(shù)A0(與應(yīng)力比和應(yīng)力強(qiáng)度因子幅值無(wú)關(guān))
A0=(0.825-0.34*alpha+0.05*alpha^2)*...
(cos(pi/2*S_max_flow))^(1/alpha);
A1=(0.415-0.071*alpha)*S_max_flow;
A3=2*A0+A1-1;
A2=1-A0-A1-A3;
%計(jì)算向量形式的裂紋張開(kāi)函數(shù)
f1=max(A0+A1*R+A2*R.^2+A3*R.^3,R);
f2=A0-2*A1;
f3=A0+A1*R;
f=f1.*(R >=0)+f2.*(R< -2)+...
f3.*(R >=-2&R< 0); %裂紋張開(kāi)函數(shù)
%****************************************
%********應(yīng)力強(qiáng)度因子門(mén)檻值計(jì)算**********
DKth_p1=DKth1*sqrt(a/(a+a0))*((1-R)./(1-f)).^(1+R*Cth_p)./...
(1-A0).^((1-R)*Cth_p); %正應(yīng)力比下的門(mén)檻值
DKth_p2=DKth1*sqrt(a/(a+a0))*((1-R)./(1-f)).^(1+R*Cth_m)./...
(1-A0).^(Cth_p-R*Cth_m); %負(fù)應(yīng)力比下的門(mén)檻值
DKth=DKth_p1.*(R >=0)+...
DKth_p2.*(R< 0); %應(yīng)力強(qiáng)度因子門(mén)檻值
%****************************************
%******根據(jù)NASGRO方程計(jì)算函數(shù)F***********
F1=log10((1-f)./(1-R).*DK); %DataX(1,:)為應(yīng)力強(qiáng)度因子幅值
F2=log10(1-DKth./DK);
F3=log10(1-1./(1-R).*(DK./Kc));
%****************************************
%******根據(jù)NASGRO方程計(jì)算裂紋擴(kuò)展速率***********
y=Coeff(1)+Coeff(2)*F1+Coeff(3)*F2+Coeff(4)*F3;
%***********************************************
%*****構(gòu)造基于最小二乘法的目標(biāo)函數(shù)F**************
%最小二乘法應(yīng)保證目標(biāo)函數(shù)F中所有原始之和達(dá)到最小
F=(y-log10(DataY)).^2;
%***********************************************
end