检查 lsqcurvefit 梯度
当您提供梯度计算函数时，非线性最小二乘求解器（如 lsqcurvefit）可以更快、更可靠地运行。然而，与其它求解器相比，lsqcurvefit 求解器在检查梯度时语法略有不同。
拟合问题和 checkGradient 语法
要检查 lsqcurvefit 梯度，不要将初始点 x0 作为数组传递，而是传递元胞数组 {x0,xdata}。例如，对于响应函数 ，从模型中创建添加了噪声的数据 ydata。响应函数为 fitfun。
function [F,J] = fitfun(x,xdata) F = x(1) + x(2)*exp(-x(3)*xdata); if nargout > 1 J = [ones(size(xdata)) exp(-x(3)*xdata) -xdata.*x(2).*exp(-x(3)*xdata)]; end end
生成 xdata 作为从 0 到 10 的随机点，并将 ydata 作为响应加上添加的噪声。
a = 2;
b = 5;
c = 1/15;
N = 100;
rng default
xdata = 10*rand(N,1);
fun = @fitfun;
ydata = fun([a,b,c],xdata) + randn(N,1)/10;检查响应函数在点 [a,b,c] 的梯度是否正确。
[valid,err] = checkGradients(@fitfun,{[a b c] xdata})valid = logical
1
err = struct with fields:
Objective: [100×3 double]
您可以安全地使用提供的目标梯度。将所有参数的下界设置为 0，无上界。
options = optimoptions("lsqcurvefit",SpecifyObjectiveGradient=true);
lb = zeros(1,3);
ub = [];
[sol,res,~,eflag,output] = lsqcurvefit(fun,[1 2 1],xdata,ydata,lb,ub,options)Local minimum possible. lsqcurvefit 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>
sol = 1×3
2.5872 4.4376 0.0802
res = 1.0096
eflag = 3
output = struct with fields:
firstorderopt: 4.4156e-06
iterations: 25
funcCount: 26
cgiterations: 0
algorithm: 'trust-region-reflective'
stepsize: 1.8029e-04
message: 'Local minimum possible.↵↵lsqcurvefit 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.'
bestfeasible: []
constrviolation: []
lsqcurvefit 中的非线性约束函数
非线性约束函数所需的语法与 lsqcurvefit 目标函数的语法不同。具有梯度的非线性约束函数具有以下形式：
function [c,ceq,gc,gceq] = ccon(x) c = ... ceq = ... if nargout > 2 gc = ... gceq = ... end end
梯度表达式的大小必须为 NxNc，其中 N 是问题变量的数量，Nc 是约束函数的数量。例如，下面的 ccon 函数只有一个非线性不等式约束。因此，对于这个包含三个变量的方程，该函数返回一个 3×1 的梯度。
function [c,ceq,gc,gceq] = ccon(x) ceq = []; c = x(1)^2 + x(2)^2 + 1/x(3)^2 - 50; if nargout > 2 gceq = []; gc = zeros(3,1); % Gradient is a column vector gc(1) = 2*x(1); gc(2) = 2*x(2); gc(3) = -2/x(3)^3; end end
检查函数 ccon 在点 [a,b,c] 返回的梯度是否正确。
[valid,err] = checkGradients(@ccon,[a b c],IsConstraint=true)
valid = 1×2 logical array
1 1
err = struct with fields:
Inequality: [3×1 double]
Equality: []
要使用约束梯度，请将选项设置为使用梯度函数，然后使用非线性约束再次求解问题。
options.SpecifyConstraintGradient = true; [sol2,res2,~,eflag2,output2] = lsqcurvefit(@fitfun,[1 2 1],xdata,ydata,lb,ub,[],[],[],[],@ccon,options)
Local minimum found that satisfies the constraints. Optimization completed because the objective function is non-decreasing in feasible directions, to within the value of the optimality tolerance, and constraints are satisfied to within the value of the constraint tolerance. <stopping criteria details>
sol2 = 1×3
4.4436 2.8548 0.2127
res2 = 2.2623
eflag2 = 1
output2 = struct with fields:
iterations: 15
funcCount: 22
constrviolation: 0
stepsize: 1.7914e-06
algorithm: 'interior-point'
firstorderopt: 3.3350e-06
cgiterations: 0
message: 'Local minimum found that satisfies the constraints.↵↵Optimization completed because the objective function is non-decreasing in ↵feasible directions, to within the value of the optimality tolerance,↵and constraints are satisfied to within the value of the constraint tolerance.↵↵<stopping criteria details>↵↵Optimization completed: The relative first-order optimality measure, 1.621619e-07,↵is less than options.OptimalityTolerance = 1.000000e-06, and the relative maximum constraint↵violation, 0.000000e+00, is less than options.ConstraintTolerance = 1.000000e-06.'
bestfeasible: [1×1 struct]
残差 res2 是残差 res 的两倍以上，而后者没有非线性约束。该结果表明，非线性约束使解远离无约束的最小值。
绘制有和没有非线性约束的解。
scatter(xdata,ydata) hold on scatter(xdata,fitfun(sol,xdata),"x") hold off xlabel("x") ylabel("y") legend("Data","Fitted") title("Data and Fitted Response, No Constraint")
figure scatter(xdata,ydata) hold on scatter(xdata,fitfun(sol2,xdata),"x") hold off xlabel("x") ylabel("y") legend("Data","Fitted with Constraint") title("Data and Fitted Response with Constraint")
