Octave中带约束函数最小化：求解最短过河路径问题

需找到河流边界曲线上的过渡点，使起点x0（曲线一侧陆地）到终点x1（曲线另一侧水域）的人员通行时间最短。已知参数包括陆地速度Vb、水中速度Vr、水流速度Vw，以及定义河流边界的函数句柄fshore。

初始用Octave的fminbnd实现求解，但存在未满足的约束：从过渡点到x1的路径必须位于曲线下方。尝试过sqp但不确定约束定义，且希望不依赖额外工具包。

初始实现代码 #

function res = maliboo(fshore,Vr,x0,x1,Vb,Vw)

x = linspace(x0(1)-2, x1(1)+2);
y = fshore(x);
plot(x,y)
ylim([min(min(y), x1(2))-1 max(max(y), x0(1))+1])
hold on
plot(x0(1), x0(2), 'o')
plot(x1(1), x1(2), '*')
v2 = sqrt(Vr.^2+Vw.^2);

t = @(xx) (sqrt((x0(1) - xx).^2 + (fshore(xx)-x0(2)).^2))./ Vb + (sqrt((x1(1) - xx).^2 + (fshore(xx)-x1(2)).^2))./ v2;
xxx = fminbnd(t, x0(1), x1(1));
##  h = @(xx) fshore(xx + (x1(1)-xx)/2) - (x1(2)-fshore(xx))/(x1(1)-xx) * (xx + (x1(1)-xx)/2) - fshore(xx)-((x1(2)-fshore(xx))/(x1(1)-xx))*xx;
##  xxx = sqp((x1(1)-x0(1))/2, t);
plot(xxx, fshore(xxx), 'x')
plot([x0(1) xxx], [x0(2) fshore(xxx)], '--', 'linewidth', 2)
plot([xxx x1(1)], [fshore(xxx) x1(2)], '--', 'linewidth', 2)

约束实现方案提示 #

1. 约束的数学本质 #

要保证过渡点(xxx, fshore(xxx))到x1的线段全程在曲线下方，等价于：线段上任意横坐标x对应的纵坐标 ≤ fshore(x)。

线段的纵坐标插值公式为：

y_line(x) = fshore(xxx) + (x1(2) - fshore(xxx))/(x1(1) - xxx) * (x - xxx)

约束即对所有x ∈ [min(xxx, x1(1)), max(xxx, x1(1))]，满足y_line(x) ≤ fshore(x)。

2. 无额外工具包的解法：带惩罚项的目标函数 #

无需额外工具包，可给原时间函数添加惩罚项，强制算法选择满足约束的点：

• 定义惩罚函数，通过检查线段上的采样点（比如中点）判断约束是否违反，违反则添加大额惩罚：

  penalty = @(xx) begin
      x_mid = xx + (x1(1) - xx)/2;
      y_mid_line = fshore(xx) + (x1(2) - fshore(xx))/(x1(1) - xx) * (x_mid - xx);
      y_mid_curve = fshore(x_mid);
      % 若线段中点在曲线上方，返回大惩罚值，否则返回0
      max(0, y_mid_line - y_mid_curve) * 1e6;  % 1e6可根据实际调整
  end;
  

• 改造目标函数：

  t_penalized = @(xx) t(xx) + penalty(xx);
  

• 用fminbnd求解改造后的函数：

  xxx = fminbnd(t_penalized, x0(1), x1(1));
  

  注：可增加采样点数量（比如三分点、四分点）提升约束检查的准确性。

3. 使用sqp的约束定义（若允许使用） #

若用sqp，可定义非线性不等式约束（约束函数返回值需≤0）：

% 以线段中点为检查点，可扩展多个点强化约束
h = @(xx) begin
    x_mid = xx + (x1(1) - xx)/2;
    y_line_mid = fshore(xx) + (x1(2)-fshore(xx))/(x1(1)-xx)*(x_mid - xx);
    y_curve_mid = fshore(x_mid);
    y_line_mid - y_curve_mid;  % 返回值≤0则满足约束
end;

调用sqp时传入约束：

xxx = sqp((x1(1)-x0(1))/2, t, [], [], h, []);

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