MATLAB的lsqnonlin函数拟合非线性式子中的未知参数

往往在实际问题中都存在exp(x)、lnx、sinx等多种函数组合的非线性经验公式。对此我们就可以通过lsqnonlin函数进行求解，该函数的方法被称为非线性最小二乘，损失函数一样，只不过类似于优化算法，给定参数初始值，然后优化参数，非线性最小二乘模型如下，即目标函数。

1. lsqnonlin语法

lsqnonlin函数用于求解非线性最小二乘（非线性数据拟合）问题。

语法
x=lsqnonlin(fun,x0)
x=lsqnonlin(fun,x0,lb,ub)
x=lsqnonlin(fun,x0,lb,ub,options)
x=lsqnonlin(problem)
[x,resnorm]=lsqnonlin(___)
[x,resnorm,residual,exitflag,output]=lsqnonlin(___)
[x,resnorm,residual,exitflag,output,lambda,jacobian]=lsqnonlin(___)
lb、ub为解向量的下界和上界lb≤a≤ub，若没有指定界，则lb=[ ]，ub=[ ]；
options为指定的优化参数；
fun为待拟合函数，计算x处拟合函数值，其定义为 function F = myfun(a,xdata)
resnorm=sum ((fun(a,xdata)-ydata).^2)，即在a处残差的平方和；
residual=fun(a,xdata)-ydata，即在x处的残差；
exitflag为终止迭代的条件；
output为输出的优化信息；
lambda为解x处的Lagrange乘子；jacobian为解x处拟合函数fun的jacobian矩阵。

2. 拟合简单指数

对数据进行简单的指数衰减曲线拟合。

从添加了噪声的指数衰减模型生成数据。模型是：

y=exp(?1.3t)+ε

其中 t 的范围是从 0 到 3，ε 是均值为 0、标准差为 0.05 的正态分布噪声。

问题表述为：给定数据（d、y），求出与数据拟合最佳的指数衰减率。

创建一个匿名函数，该函数接受指数衰减率 r 的值作为输入，并返回采用该衰减率的模型与数据之差组成的向量。

程序

clc;
clear all;
close all;
rng default
% rng('default') 将 rand、randi 和 randn 使用的随机数生成器的设置重置为其默认值。
%这样，会生成相同的随机数，就好像您重新启动了 MATLAB。默认设置是种子为 0 的梅森旋转生成器。
d = linspace(0,3);
y = exp(-1.3*d) + 0.05*randn(size(d));
fun = @(r) exp(-d*r)-y;
%找到最佳衰减率的值。任意选择一个初始估计值 x0 = 4。
x0 = 4;
x = lsqnonlin(fun,x0)
%绘制数据和最佳拟合指数曲线。
plot(d,y,'ko',d,exp(-x*d),'b-')
legend('原始数据','拟合数据')
xlabel('t')
ylabel('exp(-tx)')

运行结果

Local minimum possible.

lsqnonlin stopped because the final change in the sum of squares relative to 
its initial value is less than the value of the function tolerance.

<stopping criteria details>

x =

    1.2645

3. 拟合具有边界约束的问题

当某些拟合参数有边界时，找到最佳拟合模型。找到合适的中心化参数 b 和缩放参数 a，以便拟合以下函数：a*exp(?t)exp(?exp(?(t?b)))。

程序

clc;
clear all;
close all;
%最好地拟合以下标准正态密度分布：
%创建数据点组成的向量 t，以及在这些点上的对应正态密度。
t = linspace(-4,4);
y = 1/sqrt(2*pi)*exp(-t.^2/2);
%创建一个缩放参数 a 为 x(1) 且中心化参数 b 为 x(2) 的函数，来计算中心化并缩放的函数与正态 y 之间的差。
fun = @(x)x(1)*exp(-t).*exp(-exp(-(t-x(2)))) - y;
%从 x0 = [1/2,0] 开始寻找最佳拟合，缩放参数 a 的值介于 1/2 和 3/2 之间，中心化参数 b 的值介于 -1 和 3 之间。
lb = [1/2,-1];
ub = [3/2,3];
x0 = [1/2,0];
x = lsqnonlin(fun,x0,lb,ub)
%绘制这两个函数来查看拟合的质量。
plot(t,y,'r-',t,fun(x)+y,'b-')
xlabel('t')
legend('原始数据函数曲线','拟合函数曲线')

运行结果

Local minimum possible.

lsqnonlin stopped because the final change in the sum of squares relative to 
its initial value is less than the value of the function tolerance.

<stopping criteria details>

x =

    0.8231   -0.2444

4. 检查求解过程

在求解过程发生时和发生后都对其进行检查（通过将 Display 选项设置为 'iter' 在发生时进行检查，通过检查 output 结构体在发生后进行检查）。

假设您有观测时间数据 xdata 和观测响应数据 ydata，并且要求得参数 x(1) 和 x(2) 以拟合以下形式的模型：ydata=x(1)exp(x(2)xdata)。

程序

clc;
clear all;
close all;
% 输入观测时间和响应。
xdata = ...
 [0.9 1.5 13.8 19.8 24.1 28.2 35.2 60.3 74.6 81.3];
ydata = ...
 [455.2 428.6 124.1 67.3 43.2 28.1 13.1 -0.4 -1.3 -1.5];
% 创建简单的指数衰减模型。该模型计算预测值和观测值之间的差组成的向量。

fun = @(x)x(1)*exp(x(2)*xdata)-ydata;
% 以 x0 = [100,-1] 为起点拟合模型。通过将 Display 选项设置为 'iter' 来检查求解过程。获取 output 结构体以获取有关求解过程的详细信息。
x0 = [100,-1];
options = optimoptions('lsqnonlin','Display','iter','PlotFcn','optimplotx');
[x,resnorm,residual,exitflag,output] = lsqnonlin(fun,x0,[],[],options)
%Algorithm在 'trust-region-reflective'（默认值）和 'levenberg-marquardt' 之间进行选择。
options.Algorithm = 'levenberg-marquardt';
[x1,resnorm1,residual1,exitflag1,output1] = lsqnonlin(fun,x0,[],[],options)

运行结果

                                      Norm of      First-order 
 Iteration  Func-count     f(x)          step          optimality
     0          3          359677                      2.88e+04
Objective function returned Inf; trying a new point...
     1          6          359677        11.6976       2.88e+04      
     2          9          321395            0.5       4.97e+04      
     3         12          321395              1       4.97e+04      
     4         15          292253           0.25       7.06e+04      
     5         18          292253            0.5       7.06e+04      
     6         21          270350          0.125       1.15e+05      
     7         24          270350           0.25       1.15e+05      
     8         27          252777         0.0625       1.63e+05      
     9         30          252777          0.125       1.63e+05      
    10         33          243877        0.03125       7.48e+04      
    11         36          243660         0.0625        8.7e+04      
    12         39          243276         0.0625          2e+04      
    13         42          243174         0.0625       1.14e+04      
    14         45          242999          0.125        5.1e+03      
    15         48          242661           0.25       2.04e+03      
    16         51          241987            0.5       1.91e+03      
    17         54          240643              1       1.04e+03      
    18         57          237971              2       3.36e+03      
    19         60          232686              4       6.04e+03      
    20         63          222354              8        1.2e+04      
    21         66          202592             16       2.25e+04      
    22         69          166443             32       4.05e+04      
    23         72          106320             64       6.68e+04      
    24         75         28704.7            128       8.31e+04      
    25         78         89.7947        140.674       2.22e+04      
    26         81         9.57381        2.02599            684      
    27         84         9.50489      0.0619926           2.27      
    28         87         9.50489    0.000462262         0.0114      
Local minimum possible.
lsqnonlin stopped because the final change in the sum of squares relative to 
its initial value is less than the value of the function tolerance.
<stopping criteria details>
x =
  498.8309   -0.1013
resnorm =
    9.5049
residual =
    0.1817   -0.0610   -0.7628   -0.1196    0.2659    0.5979    1.0261    1.5124    1.5615    1.6327
exitflag =

     3
output = 

  包含以下字段的 struct:

    firstorderopt: 0.0114
       iterations: 28
        funcCount: 87
     cgiterations: 0
        algorithm: 'trust-region-reflective'
         stepsize: 4.6226e-04
          message: '?Local minimum possible.??lsqnonlin stopped because the final change in the sum of squares relative to ?its initial value is less than the value of the function tolerance.??<stopping criteria details>??Optimization stopped because the relative sum of squares (r) is changing?by less than options.FunctionTolerance = 1.000000e-06.??'

                                        First-Order                    Norm of 
 Iteration  Func-count    Residual       optimality      Lambda           step
     0           3          359677        2.88e+04         0.01
Objective function returned Inf; trying a new point...
     1          13          340761        3.91e+04       100000       0.280777
     2          16          304661        5.97e+04        10000       0.373146
     3          21          297292        6.55e+04        1e+06      0.0589933
     4          24          288240        7.57e+04       100000      0.0645444
     5          28          275407        1.01e+05        1e+06      0.0741266
     6          31          249954        1.62e+05       100000       0.094571
     7          36          245896        1.35e+05        1e+07      0.0133606
     8          39          243846        7.26e+04        1e+06     0.00944311
     9          42          243568        5.66e+04       100000     0.00821622
    10          45          243424        1.61e+04        10000     0.00777936
    11          48          243322         8.8e+03         1000      0.0673933
    12          51          242408         5.1e+03          100       0.675209
    13          54          233628        1.05e+04           10        6.59804
    14          57          169089        8.51e+04            1        54.6992
    15          60         30814.7        1.54e+05          0.1        196.939
    16          63         147.496           8e+03         0.01        129.795
    17          66         9.51503             117        0.001        9.96069
    18          69         9.50489          0.0714       0.0001       0.080486
    19          72         9.50489        4.96e-05        1e-05    5.07028e-05
Local minimum possible.
lsqnonlin stopped because the relative size of the current step is less than
the value of the step size tolerance.
<stopping criteria details>
x1 =
  498.8309   -0.1013
resnorm1 =
    9.5049
residual1 =
    0.1817   -0.0610   -0.7628   -0.1196    0.2659    0.5979    1.0261    1.5124    1.5615    1.6327
exitflag1 =

     4
output1 = 
  包含以下字段的 struct:
       iterations: 19
        funcCount: 72
         stepsize: 5.0703e-05
     cgiterations: []
    firstorderopt: 4.9629e-05
        algorithm: 'levenberg-marquardt'
          message: '?Local minimum possible.?lsqnonlin stopped because the relative size of the current step is less than?the value of the step size tolerance.??<stopping criteria details>??Optimization stopped because the relative norm of the current step, 1.016433e-07,?is less than options.StepTolerance = 1.000000e-06.??'

5.属性设置

优化选项，指定为 optimoptions 的输出或 optimset 返回的结构体。
MaxIterations
允许的迭代最大次数，为正整数。默认值为 400。
OptimalityTolerance
一阶最优性的终止容差（正标量）。默认值为 1e-6。
PlotFcn
对算法执行过程中的各种进度测量值绘图，可以选择预定义的绘图，也可以自行编写绘图函数。传递名称、函数句柄或者由名称或函数句柄组成的元胞数组。对于自定义绘图函数，传递函数句柄。默认值是“无”([])：
'optimplotx' 绘制当前点。
'optimplotfunccount' 绘制函数计数。
'optimplotfval' 绘制函数值。
'optimplotresnorm' 绘制残差范数。
'optimplotstepsize' 绘制步长大小。
'optimplotfirstorderopt' 绘制一阶最优性度量。

参考内容

[1] https://ww2.mathworks.cn/help/optim/ug/lsqnonlin.html

作者：郭志龙
编辑：郭志龙
校对：郭志龙