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