如何利用Matlab的lsqnonlin函数处理复杂函数的拟合问题
1. 简介
在工程问题的计算中，我们经常需要处理一些离散数据的拟合问题，而最小二乘法是处理曲线拟合问题的常用方法。目前，许多软件都提供有基于最小二乘法进行曲线拟合的功能，例如在Origin和Excel中均可直接利用离散数据进行曲线拟合。然而，这些软件只能处理一些简单函数的拟合问题，当需要拟合的函数较为复杂时，或者无法用简单的表达式来表述时，则往往无法直接进行拟合。为此，本文将对最小二乘法的基本原理做简单介绍，随后介绍如何利用Matlab的lsqnonlin函数处理复杂函数的拟合问题。
1.1 曲线拟合的最小二乘法原理
利用最小二乘法进行曲线拟合的本质为寻找某个近似函数 φ ( x )，使得该函数与离散点之间尽可能地逼近。若将偏差定义为近似函数的近似值 φ ( xi )与离散点 yi *之间的差值：
求解上述线性方程组即可得到拟合多项式的系数。
2. 利用Matlab处理曲线拟合问题
基于上述计算原理，Matlab提供了polyfit函数用于处理多项式曲线的拟合问题，对于一些较为复杂但仍可通过简单表达式进行表述的函数，也可以利用Matlab的拟合工具箱（Curve fitting）进行拟合。但在某些情况，当拟合函数非常复杂，以致于无法用简单表达式进行表述时（例如分段函数以及涉及到条件语句），则无法使用拟合工具箱进行拟合。对于此类问题，可以使用Matlab优化工具箱中的lsqnonlin函数进行解决。
2.1 lsqnonlin函数
lsqnonlin函数用于求解以下述形式表示的非线性最小二乘法拟合问题：
在使用该函数进行最小二乘法拟合时，lsqnonlin函数并不需要用户提供min || f ( x )||（平方和），而是需要用户提供自定义函数fun，用于计算矢量形式表示的 f ( x )：
lsqnonlin函数常用语法为：
x = lsqnonlin(fun,x0)
x = lsqnonlin(fun,x0,lb,ub)
其中fun为用户自定义函数，x0为计算采用的初始值，lsqnonlin函数首先利用x0通过自定义函数fun计算 fi (x)的取值并计算平方和，随后通过优化算法调整x的取值直至得到平方和的最小值。此外，lb和ub还可以用于定义x的取值范围，使得x满足lb≤x≤ub。
例如，对于节1.1中所述的多项式，根据最小二乘法的定义，则自定义函数 f ( x )应表示为：
注意此时 f ( x )中的x为以向量形式表示的多项式 P ( x )的系数：
在计算时，用户需要指定多项式系数的初始值，则lsqnonlin函数将利用最小二乘法计算多项式系数。
下面，本文将以笔者所在领域常用的NASGRO方程为例，介绍如何利用lsqnonlin函数处理此类复杂函数的曲线拟合问题。
2.2 NASGRO方程简介
在进行基于断裂力学的损伤容限分析时，应力强度因子和裂纹扩展速率模型是最为重要的输入。一般来说，应力强度因子可以通过经验公式或数值方法进行计算，而裂纹扩展速率模型则需要通过裂纹扩展速率试验获得的试验数据拟合得到。例如，大量的试验结果表明，在裂纹扩展的中速率区域，应力强度因子幅值ΔK和裂纹扩展速率d a /dN满足良好的对数线性关系，可以通过Paris公式进行描述：
其中C和m为材料常数。
尽管Paris公式已经得到广泛的应用，但是Paris公式仅仅描述了裂纹在中速率区域的扩展行为，没有描述近门槛区域和接近断裂的高速率区域的扩展行为，也没有考虑应力比R和裂纹闭合效应对裂纹扩展速率的影响，因此给出的计算结果将过于保守。另一个常用的裂纹扩展速率模型为Newman提出的NASGRO模型，该模型基于Forman模型改进了裂纹扩展速率模型，同时比Paris和Walker模型更加全面，不仅考虑了应力强度因子门槛值和断裂韧度，还体现了应力比以及裂纹闭合效应对裂纹扩展速率d a /dN的影响，如图2.1所示，其表达式如下：
其中R为应力比，ΔK为应力强度因子幅值，ΔKth为应力强度因子幅值门槛值，Kmax为最大应力强度因子，可表示为：
图2.1 NASGRO方程
NASGRO方程中的应力强度因子门槛值ΔKth可采用下面的经验公式进行估算：
其中A0为裂纹张开函数中的多项式系数，ΔK1是 R =1时的应力强度因子门槛值，Cth是对于正应力(上标为p=positive)和负应力比(上标为m=minus, negative)取不同值的材料常数，*a*~0~是内在小裂纹尺寸(典型值为0.0381mm)。在基于NASGRO方程开发的疲劳裂纹扩展分析软件NASGRO中，正应力比下*C*~th~^p^和Δ*K*~1~是保存在数据库里的值，负应力比下*C*~th~^m^的默认值为0.1。
2.3 NASGRO方程拟合
图2.2为疲劳裂纹扩展分析软件NASGRO材料库中某铝合金材料的裂纹扩展速率数据，已知试验时采用的试样为中心平板试样(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，下面需要通过拟合试验数据获得NASGRO方程的参数 C , m , p , q 。
图2.2 疲劳裂纹扩展数据
拟合NASGRO方程的难点主要有以下几点：
(1)裂纹扩展速率d a /dN不仅与应力强度因子幅值ΔK有关，还与使用的应力比R有关，因此实际上为多变量的拟合问题；
(2)裂纹张开函数f为分段函数，并且使用了计算最大值的max函数，该函数在拟合时无法用简单函数进行表述。
针对以上问题，NASGRO软件给出的拟合方法为首先给参数p和q确定一个初始值，并利用最小二乘法确定参数C和 m ，随后根据工程经验来获得可接受的结果，如果对拟合效果不满意，可以调整任意参数，直至获得满意的结果。
显然，这样的拟合策略具有很大的随意性，如果参数p和q选取不当，很可能对拟合效果有很大的影响。下面，本文将介绍如何利用lsqnonlin函数在不提前定义参数p和q的情况下对NASGRO方程进行拟合。
根据lsqnonlin函数的介绍，首先需要构造自定义函数 f ( x )使其满足最小二乘法计算的基本原理，由于Paris公式具有对数线性的关系，因此尝试将NASGRO方程两边取对数，可得：
上式可以用如下所示的通式表示：
系数bj为与 C , n , p和q有关(b 0 =log( C ), b 1 = n , b 2 = p , b 3 =- q )的系数，而gj为与Δ K 、R和NASGRO中所有剩余参数有关的函数。
根据最小二乘法的定义，应选取参数bj使得：
参考lsqnonlin函数对目标函数的定义，则自定义函数 f ( x )应表示为：
y= R .^2 * (R >0) + R * (R <= 0)
而裂纹张开函数f中涉及到求取最大值的计算以及分段函数的处理，也可以通过上述语法实现，具体的计算过程可参见程序代码（参见附录）。
此外，由于自定义函数 f ( x )为关于系数 bj ( j =1,2,3,4)的函数，为了将试验数据（不同应力比R下的应力强度因子幅值ΔK和裂纹扩展速率d a /d N ）传递到函数 f ( x )中进行计算，可以将试验数据定义为全局变量，以便被 f ( x )调用。
通过编写程序，可以计算得到NASGRO方程的系数如表1所示。
拟合曲线与试验数据如图2.3所示。
图2.3 试验数据及拟合曲线对比
附录1 NASGRO方程曲线拟合程序
NASGRO_LSQ.m
NASGRO_LSQ用于定义采用最小二乘法拟合NASGRO方程时的自定义函数 f ( x )，输入参数Coeff为NASGRO方程系数 bj ，输出参数为拟合函数与试验数据误差的平方和。
function F=NASGRO_LSQ(Coeff)
%程序用于计算最小二乘法拟合NASGRO方程的目标函数
%程序返回一个N×1的数值，其中N为数据对的个数
%Coeff为拟合时待求的系数（共4个系数）
%4个系数分别为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为与门槛值有关的常数
%Cth_p为正应力比下与门槛值有关的常数
%Cth_m为负应力比下与门槛值有关的常数
%a为计算门槛值时使用的裂纹长度，建议取为远大于a0的值
%Kc为材料断裂韧度
%DataX为应力强度因子幅值
%DataY为裂纹扩展速率
%***********************************
%********Newman裂纹张开函数计算**********
R=DataX(:,2); %应力比
DK=DataX(:,1); %应力强度因子幅值
%计算系数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;
%计算向量形式的裂纹张开函数
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); %裂纹张开函数
%****************************************
%********应力强度因子门槛值计算**********
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); %应力强度因子门槛值
%****************************************
%******根据NASGRO方程计算函数F***********
F1=log10((1-f)./(1-R).*DK); %DataX(1,:)为应力强度因子幅值
F2=log10(1-DKth./DK);
F3=log10(1-1./(1-R).*(DK./Kc));
%****************************************
%******根据NASGRO方程计算裂纹扩展速率***********
y=Coeff(1)+Coeff(2)*F1+Coeff(3)*F2+Coeff(4)*F3;
%***********************************************
%*****构造基于最小二乘法的目标函数F**************
%最小二乘法应保证目标函数F中所有原始之和达到最小
F=(y-log10(DataY)).^2;
%***********************************************
end