为什么 fsolve 似乎没有迭代解决方案？

Question

我正在尝试使用 fsolve 和 dogleg 方法求解两个非线性方程组。我的目标函数及其雅可比是这样的

function [F jacF]= objective(x)

    F(:,1) = ((((x(:,2)./10).*k).*(x(:,1)./100)).^2).*(rZ - Rs) +(( Cmax .* ( x(:,1)./100 ) ).^2).*( w.^2.*(rZ - Rs) ) - (((x(:,2)./10).*k).*(x(:,1)./100)); 

    F(:,2) = (x(:,2).*k).^2.*(iZ - w.*Ls) + (x(:,2).*k).^2.*x(:,1).*((w.*Ls)./200) + x(:,1).*((w.*Ls)/200).*(w.*Cmax).^2 + (w.*Cmax).^2 .*(iZ -(w.*Ls));

    if nargout > 1 % need Jacobian
        jacF = [- k - (k.^2.*x(:,2).*x(:,1).*(Rs - rZ))./50,                         - (k.^2.*x(:,2).^2.*(Rs - rZ))./100 - (Cmax.^2.*w.^2.*(Rs - rZ))./100;
                2.*k.^2.*x(:,2).*(iZ - Ls.*w) + (k.^2.*Ls.*x(:,2).*w.*x(:,1))./100,(Ls.*Cmax.^2.*w.^3)./200 + (Ls.*k.^2.*x(:,2).^2.*w)./200];
    end
end

那么我对 fsolve 的配置是这样的

    options = optimoptions('fsolve','Display','iter-detailed','PlotFcn',@optimplotfirstorderopt);
%     options.StepTolerance = 1e-13;
    options.OptimalityTolerance = 1e-12;
    options.FunctionTolerance = 6e-11;
    options.MaxIterations = 100000;
    options.MaxFunctionEvaluations = 400;%*400;
    options.Algorithm = 'trust-region-dogleg';%'trust-region'%'levenberg-marquardt';%
%     options.FiniteDifferenceType= 'central';
    options.SpecifyObjectiveGradient = true;
    fun= @objective;

    x0 = [x1',x2'];

    % Solve the function fun
    [gwc,fval,exitflag,output,jacobianEval] =fsolve(fun,x0,options);

作为方程的值

Rs =
    0.1640
Ls =
   1.1000e-07
Cmax =
   7.0000e-11
w =
   1.7040e+08
rZ =
   12.6518
iZ = 
   14.5273
K =
    0.1007
x0 = 
    70.56 0.0759

我的问题来了，因为我不明白为什么 fsolve 似乎没有像我预期的那样迭代 x(:,1)。我确实知道上述系统和参数的解决方案应该是x1=58.8和x2=0.0775。 为了测试方法的收敛性，我设置为初始猜测x0 = [x1*(1+20/100) 0.0759] = [70.56 0.0759]（x1 中的误差为 20%，x2 上的值更高），但是 fsolve 给出的解决方案是初始点，这是为什么呢？我在我的设置中做错了什么吗？

提前致谢

编辑：添加了具有通用系数的方程

Answer 1

使用您的代码，我试图找出问题所在，但似乎没有问题，除了“病态”行为。

我试图观察你的函数的行为，并注意到：

x = -50:50;
y = -50:50;

[X,Y]=meshgrid(x,y);
F1 = zeros(size(X));
F2 = zeros(size(X));
for i=1:size(X,1)
    for j=1:size(X,2)
        f = objective([X(i,j) Y(i,j)]);
        F1(i,j) = f(1);
        F2(i,j) = f(2);
    end
end

figure; 
subplot(1,2,1)
surf(X,Y,F1); shading interp; xlabel('x1');ylabel('x2');zlabel('F_1'); colorbar
subplot(1,2,2)
surf(X,Y,F2); shading interp; xlabel('x1');ylabel('x2');zlabel('F_2'); colorbar

F1 是向量函数的第一个元素，F2 是第二个元素。 请注意，F1 在整个范围内几乎是恒定的（从 0 到 1 变化很小）。另请注意，对于F2 的大部分表面，您有亮黄色部分，这意味着此功能也没有太大变化。对于任何给定的初始条件，F 函数的范数都足够小，因此该区域中的任何点都将被认为“对于函数的零足够好”。 我们还注意到，您等式中的某些值具有非常不同的数量级，

Rs =
    0.1640
Ls =
   1.1000e-07
Cmax =
   7.0000e-11
w =
   1.7040e+08
rZ =
   12.6518

这使得解决起来更加困难。最好的解决方案是尝试改变变量来标准化输入和输出。这应该缩放自变量和向量值函数，以改善模型的数值不准确性。

---

好的，你现在已经提供了你的方程，你的问题似乎比我最初想象的更糟糕（在数字意义上）。

你的方程式是：

鉴于x1 是第二个变量的函数，您实际上只有 1 个自变量，而不是之前看起来的 2 个。因此，您可以编写F1 = f1(x2) 和F2=f2(x2)，因为这两个函数都只是一个变量的函数。你有两种选择来解决这个由 2 个方程和 1 个变量组成的方程组（注意到这里的问题了吗？）

1. 第一个选项 - 简单的选项 - 分别求解每个方程。通过这样做，您将获得x2 的两个不同值，一个满足第一个方程，另一个满足第二个方程。这没有帮助。

2. 第二个选项 - 最难的 - 是同时求解两个方程。这里的难点在于，您必须在求解方程时定义一个要满足的精确标准。看，你有两个方程，但只有一个变量，那么你怎么知道你找到的解是“最优”的呢？一个常用的标准是求解least-squares sense 中的两个方程，即，您发现x2 的值使S=F1^2+F2^2 的平方和最小。 现在，具有通用系数的方程与您提供的代码不完全匹配，所以我不知道这些ai、bi 系数中的哪一个是Rs、Ls、Cmax 等，但是你通常可以以最小二乘的方式解决您的方程：

   a1 = 兰特； a2 = 兰特； a3 = 兰特； b1 = 兰特； b2 = 兰特； b3 = 兰特； b4 = 兰特； b5 = 兰特；

   fun = @(x2) myfun(x2,a1,a2,a3,b1,b2,b3,b4,b5); x0 = 1; [X,FVAL,EXITFLAG] = fminsearch(fun,x0)

   函数 S = myfun(x2,a1,a2,a3,b1,b2,b3,b4,b5) x1 = a2*x2./(a1*x2.^2+a3);

   f1 = a1*(x1.*x2).^2 - a2*x1.*x2 + a3.*x1.^2;
   f2 = b1*(x1.*x2).^2 + b2*x1.*(x1.*x2).^2 + b3*x1.^3 + b4*x1.^2 + b5*x1;
   
   S = f1^2+f2^2;
   

   结束

请注意，您必须正确定义常数和初始估计值。由于x1 是x2 的函数，因此您可以在函数体内计算它。然后，计算函数的每个分量，然后创建平方和。函数fminsearch 求平方和的最小值——方程的解。