Matlab带边界约束的非线性方程组fsolve变换求解问题咨询

问题描述 #

我需要在Matlab中求解带边界约束（即lb<=x<=ub）的非线性方程组。Matlab官方建议将问题转换为最小化问题，使用lsqnonlin或fmincon求解，但我希望尝试结合变量变换的fsolve方法。

我选取了一个含2个未知量的非线性方程组示例，该方程组共有4个解，施加x>=0的约束后，仅x*=(10,20)满足条件。为求解约束条件下的方程组F(x)=0 s.t. x>=0，我尝试调用fsolve求解F(z.^2)=0，得到解z*后再转换回x=z.^2（平方变换确保x始终非负），初始点选择了可行值并避开(0,0)。

但该变换方法并未生效，尝试带边界的lsqnonlin也未成功——仅当初始点非常接近真实解(10,20)时才有效，但在更复杂的问题中我无法预知解的位置。

我的疑问：是否应该使用其他变换（比如x=exp(z)）？变换是否必须为一一映射？毕竟平方变换显然不是。

最小工作示例 #

%% fsolve with bound constraints
% 方程组有4个解，仅(10,20)满足x(1)>=0, x(2)>=0的约束
% (-1,-2) 不可行
% (10,-2) 不可行
% (-1,20) 不可行
% (10,20) 可行
% 变换思路：将原问题F(x)=0（x无界）转换为求解F(z.^2)=0（z无界，z.^2>=0）
clear
clc
close 
x0 = [5,10];
%opts = optimoptions(@fsolve,'Display','off');
%% 变换方法
x0 = sqrt(x0);
[x_star,~,flag] = fsolve(@fbnd_trans,x0); 
x_star = x_star.^2;
if flag<=0
    warning('方程组求解失败！')
end
disp('解x* = ')
disp(x_star)
disp('残差 = ')
Residuals = 
disp(fbnd(x_star))
 
%% lsqnonlin求解
lb = [0,0]';
ub = [inf,inf]';
[x_star,~,~,flag] = lsqnonlin(@fbnd,x0,lb,ub);
if flag<=0
    warning('方程组求解失败！')
end
disp('解x* = ')
disp(x_star)
disp('残差 = ')
disp(fbnd(x_star))

%------------------------- 子函数 ----------------------------------% 
function F = fbnd(x)
F(1) = (x(1)+1)*(10-x(1))*(1+x(2)^2)/(1+x(2)^2+x(2));
F(2) = (x(2)+2)*(20-x(2))*(1+x(1)^2)/(1+x(1)^2+x(1));
end
function F = fbnd_trans(x)
% 施加x>=0约束的变换
x = x.^2;
F(1) = (x(1)+1)*(10-x(1))*(1+x(2)^2)/(1+x(2)^2+x(2));
F(2) = (x(2)+2)*(20-x(2))*(1+x(1)^2)/(1+x(1)^2+x(1));
end

问题分析与解决方案 #

1. 平方变换失效的原因 #

平方变换不是全局一一映射，会导致变换后的方程组F(z²)=0存在多个z对应同一个x的情况，fsolve可能收敛到对应非可行x的z解（比如z=(±1,±2)对应x=(1,4)，并非目标解）。此外，平方变换在z=0附近导数为0，会导致数值求解时的雅可比矩阵奇异，严重影响收敛稳定性。

2. 变换的选择原则 #

优先选择严格单调的一一映射，保证x与变换后的变量（如z）一一对应，避免多解混淆求解器，同时保证导数非零，避免雅可比奇异：

• 对于x>=0的约束：x=exp(z)是最优选择之一，指数变换严格单调递增，且导数始终非零，不会出现数值稳定性问题；也可以用x = z² + ε（ε为极小正数），但稳定性不如指数变换。
• 对于lb<=x<=ub的有限区间约束：可以用x = lb + (ub-lb)*(1+tan(z))/(1-tan(z))或sigmoid类变换，将有限区间映射到整个实数域。

3. lsqnonlin失效的解决思路 #

lsqnonlin是局部优化器，原问题存在多个局部残差极小点，当初始点远离目标解时容易收敛到非可行解的局部极小。可以尝试：

• 多设置几个不同的初始点进行尝试；
• 结合全局优化工具（如globalSearch或MultiStart），配合lsqnonlin或fmincon遍历解空间，找到满足约束的全局最优解。

4. 修正后的指数变换示例 #

将变换改为x=exp(z)，修改后的变换函数及调用代码如下：

%% 指数变换方法
x0 = [5,10];
z0 = log(x0); % 初始点转换
[x_z_star,~,flag] = fsolve(@fbnd_exp_trans,z0); 
x_star = exp(x_z_star); % 转换回原变量
if flag<=0
    warning('方程组求解失败！')
end
disp('指数变换法解x* = ')
disp(x_star)
disp('残差 = ')
disp(fbnd(x_star))

% 指数变换的子函数
function F = fbnd_exp_trans(z)
x = exp(z);
F(1) = (x(1)+1)*(10-x(1))*(1+x(2)^2)/(1+x(2)^2+x(2));
F(2) = (x(2)+2)*(20-x(2))*(1+x(1)^2)/(1+x(1)^2+x(1));
end

内容的提问来源于stack exchange，提问作者Alessandro