Main Content

基于问题的非线性最小二乘

此示例说明如何使用 基于问题的优化工作流 执行非线性最小二乘曲线拟合。

模型

此问题的模型方程是

y(t)=A1exp(r1t)+A2exp(r2t),

其中 A1A2r1r2 是未知参数,y 是响应,t 是时间。该问题需要时间 tdata 和(含噪)响应测量 ydata 的数据。目标是找到最佳 Ar,也就是最小化下式的那些值:

ttdata(y(t)-ydata)2.

样本数据

通常,您有针对某个问题的数据。在这种情况下,请为该问题生成人为含噪数据。使用 A = [1,2]r = [-1,-3] 作为基础值,使用从 0 到 3 的 200 个随机值作为时间数据。绘制生成的数据点。

rng default % For reproducibility
A = [1,2];
r = [-1,-3];
tdata = 3*rand(200,1);
tdata = sort(tdata); % Increasing times for easier plotting
noisedata = 0.05*randn(size(tdata)); % Artificial noise
ydata = A(1)*exp(r(1)*tdata) + A(2)*exp(r(2)*tdata) + noisedata;
plot(tdata,ydata,'r*')
xlabel 't'
ylabel 'Response'

Figure contains an axes object. The axes object contains an object of type line.

数据是含噪数据。因此,该解可能不会很好地匹配原始参数 Ar

基于问题的方法

要找到最佳拟合参数 Ar,首先用这些名称定义优化变量。

A = optimvar('A',2);
r = optimvar('r',2);

为目标函数创建一个表达式,它是要最小化的平方和。

fun = A(1)*exp(r(1)*tdata) + A(2)*exp(r(2)*tdata);
obj = sum((fun - ydata).^2);

用目标函数 obj 创建一个优化问题。

lsqproblem = optimproblem("Objective",obj);

对于基于问题的方法,将初始点指定为结构体,并将变量名称作为结构体的字段。指定初始 A = [1/2,3/2] 和初始 r = [-1/2,-3/2]

x0.A = [1/2,3/2];
x0.r = [-1/2,-3/2];

检查问题表示。

show(lsqproblem)
  OptimizationProblem : 

	Solve for:
       A, r

	minimize :
       sum(arg6)

       where:

           arg5 = extraParams{3};
           arg6 = (((A(1) .* exp((r(1) .* extraParams{1})))
     + (A(2) .* exp((r(2) .* extraParams{2})))) - arg5).^2;

         extraParams{1}:

         0.0139
         0.0357
         0.0462
         0.0955
         0.1033
         0.1071
         0.1291
         0.1385
         0.1490
         0.1619
         0.1793
         0.2276
         0.2279
         0.2345
         0.2434
         0.2515
         0.2533
         0.2894
         0.2914
         0.2926

         :
         :

       extraParams{2}:

         0.0139
         0.0357
         0.0462
         0.0955
         0.1033
         0.1071
         0.1291
         0.1385
         0.1490
         0.1619
         0.1793
         0.2276
         0.2279
         0.2345
         0.2434
         0.2515
         0.2533
         0.2894
         0.2914
         0.2926

         :
         :

       extraParams{3}:

         2.9278
         2.7513
         2.7272
         2.4199
         2.3172
         2.3961
         2.2522
         2.1974
         2.1666
         2.0944
         1.9566
         1.7989
         1.7984
         1.7540
         1.8318
         1.6745
         1.6874
         1.5526
         1.5229
         1.5680

         :
         :

基于问题的解

求解。

[sol,fval] = solve(lsqproblem,x0)
Solving problem using lsqnonlin.

Local minimum found.

Optimization completed because the size of the gradient is less than
the value of the optimality tolerance.
sol = struct with fields:
    A: [2x1 double]
    r: [2x1 double]

fval = 0.4724

绘制得到的解和原始数据。

figure
responsedata = evaluate(fun,sol);
plot(tdata,ydata,'r*',tdata,responsedata,'b-')
legend('Original Data','Fitted Curve')
xlabel 't'
ylabel 'Response'
title("Fitted Response")

Figure contains an axes object. The axes object with title Fitted Response contains 2 objects of type line. These objects represent Original Data, Fitted Curve.

该图显示拟合数据与原始含噪数据匹配良好。

查看拟合参数与原始参数 A = [1,2]r = [-1,-3] 的匹配程度。

disp(sol.A)
    1.1615
    1.8629
disp(sol.r)
   -1.0882
   -3.2256

A 中的拟合参数偏差约为 15%,r 中的拟合参数偏差约为 8%。

不受支持的函数要求 fcn2optimexpr

如果您的目标函数不是由初等函数组成,您必须使用 fcn2optimexpr 将该函数转换为优化表达式。请参阅Convert Nonlinear Function to Optimization Expression。对于本示例:

fun = @(A,r) A(1)*exp(r(1)*tdata) + A(2)*exp(r(2)*tdata);
response = fcn2optimexpr(fun,A,r);
obj = sum((response - ydata).^2);

求解此问题的其余步骤是相同的。唯一不同是在绘图例程中调用 response 而不是 fun

responsedata = evaluate(response,sol);

有关支持的函数列表,请参阅Supported Operations for Optimization Variables and Expressions

另请参阅

相关主题