最优化问题 - 内点法

以下是一种循序推理的方式，来帮助你从基础概念出发，理解 内点法（Interior-Point Method, IPM） 是什么、为什么要用它，以及它是如何工作的。

1. 问题起点：带不等式约束的优化

假设你有一个带不等式约束的优化问题：
$$\begin{array}{rl}& \underset{x}{\min }f(x)\\ & \text{subject to }{g}_i(x)\le 0,i=1,\dots ,m.\end{array}$$

• $f(x)$：目标函数
• $g_i(x)\le 0$：不等式约束（例如，控制输入不能超范围）

我们想要找到一个满足所有约束且使目标函数最小的点 $x^{\ast }$。

2. 直接在边界上“走”的困难

一种思路是：最优点往往在可行域边界上（很多经典优化问题里，最优解会出现在约束生效的边界）。

• 但是，如果我们只在边界上“走”，常常要先找到并跟随所有活跃约束的交线，这在高维情况下非常复杂。
• 并且，如果问题是非线性的，直接处理“在边界上走”会更难，容易出现数值不稳定或无法收敛的问题。

3. 想法：能不能“在里面”搜索，最后再逼近边界？

为了避免“一上来就卡在边界”，有人提出：

• 先在可行域内部（所有约束 $g_i(x)<0$ 都“宽裕”）附近做搜索，那里没有太多束缚，比较容易用连续的梯度法去迭代寻优；
• 当我们逐渐找到一个更优解，需要靠近边界时，再“慢慢地”逼近它。

这样做，就不需要每一步都严格跟在边界上，也能保证可行性（不跑到约束外面去）。

4. 内点法的关键：障碍函数（Barrier Function）

为在域内部搜索，需要一种数学手段，让解“自动”不跑出界。这时就引入了障碍函数。

4.1 形式：对约束加一个“负对数”项

对每个不等式约束 $g_i(x)\le 0$，定义一个项：
$$-\log (-g_i(x)).$$

• 如果 $g_i(x)<0$，那么 $-g_i(x)>0$，对数是可以求的。
• 一旦 $g_i(x)$ 接近 0（也就是接近边界），$-g_i(x)$ 趋向 0，这个对数会趋向 $-\infty$，但我们在目标函数中是带一个负号“-”，变成了一个正无穷大的惩罚。
• 这样一来，就相当于在可行域的边缘设置了一道**“墙”，越靠近边界，代价就越大，解被迫**留在可行域内部。

4.2 组合成“障碍型”目标函数

把原本的目标函数 $f(x)$ 和障碍项结合：
$$B(x,\mu )=f(x)-\mu \sum _{i=1}^m\log (-g_i(x)).$$

• $\mu >0$ 是一个系数，称为“障碍参数（barrier parameter）”。
• 当 $\mu$ 较大时，障碍惩罚也大，解就离约束边界更远；当我们逐渐减小 $\mu$，系统会允许解更靠近边界。

5. 迭代逼近最优解

内点法并不是一次就把问题解决，而是：

1. 从一个“较宽松”的问题开始（$\mu$ 比较大，边界惩罚很强）。
2. 求解 ${\min }_{x}B(x,\mu )$，得到一个可行域内部的解。
3. 降低 $\mu$ 的值，重新求解，这时可以允许解更加靠近（或到达）真正的最优边界。
4. 多次迭代后，$\mu \to 0$，解会逼近真实最优解。

这一过程里，算法会用到类似于 牛顿法、KKT 条件 等工具来求解各轮的子问题。

6. 内点法 vs. 其他方法

• 单纯形法（Simplex）：在可行域的“顶点—边界—面”上走，比较适合线性规划；但在高维非线性问题上，效率可能较低。
• 内点法：通过“障碍函数”把解留在域内部，逐渐往边界逼近，常在非线性规划（NLP）中表现更好，也更适合大规模问题。

7. 在 MPC 中的应用

MPC 求解器（如 IPOPT）正是利用内点法来处理：

• 多步预测的系统动力学约束
• 输入/状态上下限
• 非线性目标函数

让无人机等系统在高维空间里高效地搜索最优控制输入。

🔑 小结：内点法的一句话总结

内点法是一种在可行域“内部”迭代搜索的求解策略，通过“障碍函数”阻挡解跑出可行域，并逐步放松障碍参数，最终逼近最优解和约束边界。

这就是内点法的核心“推理过程”：与其从边界开始，不如在“里层”走，让数值算法更稳定，再慢慢让解贴近约束边界，从而找到最优。

下面我将针对这段代码的逻辑和实现细节，结合“内点法(障碍函数) + 固定步长梯度下降”这一思路，做一个比较细致的解析。

1. 整体思路与算法框架

目标问题（从注释和注解上看）是：

$$\begin{array}{rl}& \min f(x)=x_1+x_2+x_3\\ & \text{s.t.}{x}_1^2+x_2^2+x_3^2\le 1.\end{array}$$

• 可行域是单位球 { x ∈ R 3 : ∥ x ∥ ≤ 1 } \{x \in \mathbb{R}^3 : \|x\|\le 1 \} {x∈R3:∥x∥≤1}；

• 希望用内点法来解决不等式约束 $x_1^2+x_2^2+x_3^2\le 1$。
  通常在内点法中，会把不等式 (g(x) \le 0)（这里 $g(x)=x^2+y^2+z^2-1$ ）转写成形如 (h(x) > 0)，再将 (\log h(x)) 当作障碍项（barrier）。本例里定义
  $$h(x)=1-(x_1^2+x_2^2+x_3^2),$$
  于是$h(x)>0$等价于 $x_1^2+x_2^2+x_3^2<1$。

• 内点法的障碍目标函数定义为
  $$B(x,\mu )=f(x)-\mu \ln (h(x)),$$
  其中 $\mu >0$ 是障碍系数(barrier parameter)。$\mu$ 越小，$-\mu \log h(x)$的惩罚作用越强，最终会将可行解“推”到最接近边界的最优位置上。

• 在固定一个$\mu$ 的情况下，常用牛顿法或梯度法去最小化 $B(x,\mu )$ 。在算法外层循环中，则逐步减小 $\mu$（比如乘个衰减因子），让解逐步逼近约束边界并最终得到较准确的可行最优解。

在这份代码中：

1. 外层循环 (共 max_outer 次)：

   • 每次先用当前 $\mu$ 在球内做若干步梯度下降，尝试求解 $\min B(x,\mu )$；
   • 之后衰减 $\mu$ $\leftarrow$ $\mu \cdot$ $\text{(mu\_decay)}$，进入下一轮；
   • 重复，直到$\mu$ 足够小或者达到外层迭代上限。

2. 内层循环 (共 max_inner 次)：

   • 计算当前 $B$ 的梯度 $\nabla B(x,\mu )$；
   • 做一步固定步长的梯度下降：$x\leftarrow x-\alpha \nabla B(x,\mu )$；
   • 如果发现新点 $x_{\text{new}}$不可行（或者离边界过近），就缩小步长再试；
   • 如果两次迭代 $\Vert x_{\text{new}}-x\Vert$ 非常小(< tol)，就视为收敛并 break。

这样就得到一条渐进接近最优解的迭代轨迹，存放在 x_history 中。

2. 代码主干解析

从最外层开始看起，核心部分在函数

function x_history = interior_point_3d_solve()
    % 参数
    mu_init = 1.0;      % 初始障碍参数
    mu_decay = 0.2;     % 每轮迭代后 mu 的衰减因子
    alpha    = 0.001;   % 固定梯度下降步长
    tol      = 1e-6;    % 收敛判据
    max_outer= 10;      % 外层循环次数
    max_inner= 50;      % 每次 mu 下最大内层迭代次数

    % 初始化可行解(球内)
    x = [0;0;0];  
    mu = mu_init;

    x_history = [];

    for outer = 1:max_outer
        for inner = 1:max_inner
            g = grad_B(x, mu);  
            x_new = x - alpha*g;

            % 若越过球边界 => h(x_new) <= 0
            if h_3d(x_new) <= 1e-9
                x_new = x - 0.1*alpha*g;  % 缩小步长再试
            end

            % 收敛判断
            if norm(x_new - x) < tol
                x = x_new;
                break;  % 跳出内层循环
            end

            x = x_new;
            x_history = [x_history; x']; 
        end

        % 减小 mu
        mu = mu * mu_decay;
        if mu < 1e-12
            break;    % mu 已经很小，不必再迭代
        end
    end

    % 加入最终点
    x_history = [x_history; x'];
end

1. mu_init 设为 1.0，之后每次外层循环会 mu = mu * mu_decay，即乘以 0.2。这样大概迭代几次后就会让 mu 变得很小。
2. alpha=0.001 是固定的梯度下降步长，相对比较小，所以我们用到 max_inner=50 步来让它收敛到一个合适精度。
3. 收敛阈值 tol=1e-6：若相邻两步$\Vert x_{\text{new}}-x\Vert$ 小于这个值，就认为内层已经收敛。
4. 每次更新 x 后都会把它记录到 x_history 里，用于可视化迭代轨迹。

这里值得注意的是：如果单纯用固定步长，可能碰到越过可行域边界（即 (x^2+y2+z^2>1)）的风险。为此，代码做了一个简单检查：

if h_3d(x_new) <= 1e-9
    x_new = x - 0.1*alpha*g;  % 步长缩小10倍
end

当然，这只是一个很“简单粗暴”的处理，工业级内点法通常要做更精细的线搜索或牛顿校正，这里是为了示例演示。

3. 障碍函数与梯度的计算

3.1 h_3d(x)

function val = h_3d(x)
    val = 1.0 - (x(1)^2 + x(2)^2 + x(3)^2);
end

即$h(x)=1-r^2$，其中 $r^2=x_1^2+x_2^2+x_3^2$。球内保证 $h(x)>0$。

3.2 B_3d(x, mu)

function val = B_3d(x, mu)
    val = f_3d(x) - mu*log(h_3d(x));
end

对应障碍目标$B(x,\mu )=f(x)-\mu \ln (h(x))$。

3.3 grad_B(x, mu) 的推导

从数学上看，如果
$$B(x,\mu )=f(x)-\mu \ln (h(x)),$$
则其梯度为
$$\nabla B(x,\mu )=\nabla f(x)-\mu \nabla \ln (h(x)).$$
其中
$$\nabla \ln (h(x))=\frac{1}{h(x)}\nabla h(x).$$
而 (h(x) = 1 - r^2)，(\nabla h(x) = -2x)。所以
$$\nabla \ln (h(x))=\frac{-2x}{1-r^2}.$$
带上负号“$-\mu$”一起，就得到对障碍项的贡献为$+\frac{2\mu x}{1-r^2}$。

如果我们要最小化 $f(x)=x_1+x_2+x_3$，其梯度就是 $(1,1,1)$。
于是
$$\nabla B(x,\mu )=(1,1,1)+\frac{2\mu }{1-r^2}(x_1,x_2,x_3).$$

在代码中，可以看到：

function g = grad_B(x, mu)
    hx = h_3d(x);  % 1 - r^2
    dB_dx0 = 1.0 + (2.0 * mu * x(1) / hx);
    dB_dx1 = 1.0 + (2.0 * mu * x(2) / hx);
    dB_dx2 = 1.0 + (2.0 * mu * x(3) / hx);

    g = [dB_dx0; dB_dx1; dB_dx2];
end

正好对应上面的公式：1.0 就是 $\nabla f(x)$ 的那一部分，$(2.0\ast mu\ast x(i)/hx)$ 对应障碍项梯度那一部分。

4. 运行与结果

如果修正了 f_3d，那么在球内最小化 (x+y+z) 的最优解，理论上会落在球面上与 ((1,1,1)) 方向相反的地方，也就是球面朝着 ((-1,-1,-1)) 的方向。其最优解应当是

$$\left(-\frac{1}{\sqrt{3}},-\frac{1}{\sqrt{3}},-\frac{1}{\sqrt{3}}\right)$$
因为在约束 $x^2+y^2+z^2=1$ 上，(x+y+z) 的最小值就是 $-\sqrt{3}$。迭代跑起来后，你应该能看到解最终趋近这个点，轨迹会从球心出发，一路在球内前进，并在迭代后期逐渐贴近球面。

如果把可视化部分加上（即下面这段示例）：

x_history = interior_point_3d_solve();

figure('Color','w','Name','Interior-Point 3D Demo');
hold on; grid on; axis equal;

% 绘制单位球面
[Xs, Ys, Zs] = sphere(50);
surf(Xs, Ys, Zs, ...
     'FaceAlpha',0.1, 'EdgeColor','none', 'FaceColor','c');

xlabel('x_0'); ylabel('x_1'); zlabel('x_2');
title('Minimize x_0 + x_1 + x_2 subject to x_0^2 + x_1^2 + x_2^2 \le 1');

% 画迭代轨迹
x0_hist = x_history(:,1);
x1_hist = x_history(:,2);
x2_hist = x_history(:,3);
plot3(x0_hist, x1_hist, x2_hist, 'b-o','LineWidth',1.5);

% 最后一点标红
plot3(x0_hist(end), x1_hist(end), x2_hist(end), ...
      'ro', 'MarkerSize',8, 'MarkerFaceColor','r');

legend('Unit Sphere','Iter Process','Final Solution');
view(35,25);

就能看到一个球面和从原点出发，沿着负对角线方向慢慢收敛到球面那一点的轨迹。

5. 小结

1. 内点法原理：通过把约束$x_1^2+x_2^2+x_3^2\le 1$ 转化为对数障碍$-\mu \ln (1-r^2)$，并在外层迭代不断减小 $\mu$。这能保证迭代点始终在球内 ($h(x)>0$)，同时在 $\mu \to 0$ 时逐渐收敛到可行域边界上的最优点。

2. 实现细节：

   • 用固定步长 + 简单可行性校正(越界时缩步长)；
   • 用 max_outer 与 max_inner 控制多重循环；
   • 用 tol 判断迭代收敛；
   • 在每次迭代都记录 x 到 x_history 中，用于可视化。

3. 需要修正的地方：
   代码中的 f_3d(x) 与注释/梯度公式不一致，应改回

   function val = f_3d(x)
       val = x(1) + x(2) + x(3);
   end
   

   这样才与“最小化 (x+y+z)”的需求相吻合。

修正后运行，即可得到一个从球心(原点)出发，最终靠近 $(-\frac{1}{\sqrt{3}},-\frac{1}{\sqrt{3}},-\frac{1}{\sqrt{3}})$ 的迭代过程。

参考：为什么最优解是 $(-1/\sqrt{3},-1/\sqrt{3},-1/\sqrt{3})$？

因为在单位球约束下，若要最小化$x_1+x_2+x_3$，相当于“在球面上找与(1,1,1)方向夹角最大的点”——也就是与 $(1,1,1)$ 反方向的单位向量。它正好是 $\frac{-1}{\sqrt{3}}(1,1,1)$。目标值是
$$x_1+x_2+x_3=-\sqrt{3}.$$
这与几何直觉、拉格朗日乘子法都能得到一样的结果。

6. 结语

• 该代码很好地演示了使用对数障碍(log-barrier)的内点法思路：用一系列的无约束子问题（带障碍项）来逼近有约束优化，并在外层迭代中逐渐减小障碍系数 (\mu)，使解贴近约束边界的最优解。
• 不过，在正式应用时，往往不会用固定步长的简单梯度下降，而会用牛顿法或线搜索，以获得更好的数值稳定性和收敛速度。
• 代码本身最大的问题是 f_3d 与注释/公式不匹配，只要改回 f_3d(x) = x(1)+x(2)+x(3) 即可与文档保持一致，也能正确体现“最小化$\sum x_i$”的意图。

在使用对数障碍（log‐barrier）形式的内点法时，$\mu$（有时也记作 $t$ 或 $\nu$）通常被称为障碍参数（barrier parameter）。它的核心作用是：

1. 控制对数障碍项的“强度”
   在障碍型目标函数
   $$B(x,\mu )=f(x)-\mu \ln (h(x))$$
   中，(\mu) 决定了 (-\mu \ln\bigl(h(x)\bigr)) 这部分惩罚力度的大小。

   • 当 (\mu) 较大时，对(\ln(h(x)))的惩罚力度相对较小，算法对靠近约束边界的“敏感度”不高，所以在收敛初期可以更自由地在可行域里移动。
   • 当 (\mu) 变小时，(-\mu \ln(h(x)))会变得更尖锐，迫使解更加贴近且“贴合”可行域的边界（若这是最优解所在的位置）。

2. 帮助逐步逼近最优解并保持可行
   内点法的思路是：一开始用较大的 (\mu)（障碍作用较弱）来保证算法稳步地在可行域里进行搜索；之后在外层循环中逐步减小 (\mu)，使对数障碍项逐渐变得陡峭，从而把解“推”到真正需要的边界附近并逼近最优解。

3. 数值上的平衡

   • 如果 (\mu) 过小，一开始 (-\mu \ln(h(x))) 的势垒就会非常强，导致很难在可行域里移动，且容易产生数值不稳定（如(\ln(h(x)))趋向负无穷）。
   • 如果 (\mu) 过大，到收敛后期也无法精确地在边界附近找到最优解。所以通常会有一条**“(\mu)衰减”路径**（比如 (\mu \leftarrow \beta \mu)，(\beta<1)）来让解逐步逼近最优值。

简而言之：(\mu) 是控制“障碍强度”的调节器。随着 (\mu) 从大到小的不断衰减，解会逐渐向可行域边界靠拢并最终获得精确的约束最优解。

完整代码

function x_history = interior_point_3d_solve()
%{
 使用内点法(障碍函数 + 固定步长梯度下降)求解:
 
   min  f(x) = x(1) + x(2) + x(3)
   s.t. x(1)^2 + x(2)^2 + x(3)^2 <= 1.
 
 内点法: 定义障碍型目标
   B(x, mu) = f(x) - mu * log( h(x) ), 其中
   h(x) = 1 - (x0^2 + x1^2 + x2^2).
 
 输出:
   x_history: 每个迭代得到的 (x0, x1, x2).
%}

    % 参数
    mu_init = 1.0;      % 初始障碍参数
    mu_decay = 0.2;     % 每轮迭代后 mu 的衰减因子
    alpha = 0.001;       % 梯度下降步长
    tol = 1e-6;         % 收敛阈值
    max_outer = 10;     % 外层循环(更新 mu)次数
    max_inner = 50;     % 每次 mu 下最大迭代次数

    % 初始化可行解 (x0, x1, x2)，球内，例如原点
    x = [0; 0; 0];  
    mu = mu_init;

    x_history = [];

    for outer = 1:max_outer
        for inner = 1:max_inner
            g = grad_B(x, mu);  % 计算梯度
            x_new = x - alpha*g;

            % 如果越过球边界 => h(x_new) <= 0 => x_new^2+y^2+z^2 >=1
            % 简单地缩小步长再试
            if h_3d(x_new) <= 1e-9
                x_new = x - 0.1*alpha*g;
            end

            % 收敛判断
            if norm(x_new - x) < tol
                x = x_new;
                break;
            end

            x = x_new;
            x_history = [x_history; x']; % 记录轨迹(行向量)
        end

        % 降低 mu 使解更逼近约束边界
        mu = mu * mu_decay;
        if mu < 1e-12
            break;
        end
    end

    % 最后再将最终点加入
    x_history = [x_history; x'];
end

%% 目标函数 f(x)
function val = f_3d(x)
    % f(x0, x1, x2) = x0 + x1 + x2
    val = x(1) + x(2) + x(3);
end

%% 障碍函数项 h(x) = 1 - (x0^2 + x1^2 + x2^2)
function val = h_3d(x)
    val = 1.0 - (x(1)^2 + x(2)^2 + x(3)^2);
end

%% 障碍型目标 B(x, mu) = f(x) - mu*ln(h(x))
function val = B_3d(x, mu)
    val = f_3d(x) - mu*log(h_3d(x));
end

%% B(x, mu) 的梯度
function g = grad_B(x, mu)
%{
 B(x) = (x0 + x1 + x2) - mu * ln(1 - r^2),
  其中 r^2 = x0^2 + x1^2 + x2^2
 => dB/dx0 = 1 + [2 mu x0 / (1 - r^2)]
 => dB/dx1 = 1 + [2 mu x1 / (1 - r^2)]
 => dB/dx2 = 1 + [2 mu x2 / (1 - r^2)]
%}
    hx = h_3d(x);
    r2 = x(1)^2 + x(2)^2 + x(3)^2;
    % hx = 1 - r2 > 0  (只要在球内)

    dB_dx0 = 1.0 + (2.0 * mu * x(1) / hx);
    dB_dx1 = 1.0 + (2.0 * mu * x(2) / hx);
    dB_dx2 = 1.0 + (2.0 * mu * x(3) / hx);

    g = [dB_dx0; dB_dx1; dB_dx2];
end



%{
 演示如何在 3D 中用“障碍函数 + 简单梯度下降”的内点法
 来最小化 f(x) = x0 + x1 + x2
 subject to x0^2 + x1^2 + x2^2 <= 1.
 
 可行域是单位球 (x0^2 + x1^2 + x2^2 <= 1)。
 我们会在图中绘制球面，并用散点绘制迭代轨迹。
%}

% 1) 调用求解函数，得到每步迭代的解 x(k)
x_history = interior_point_3d_solve();

% 2) 可视化
figure('Color','w','Name','Interior-Point 3D Demo');
hold on; grid on; axis equal;  % 3D 坐标中，最好设 axis equal

% 2.1 绘制单位球面（x0^2 + x1^2 + x2^2 = 1）
[Xs, Ys, Zs] = sphere(50);  
% sphere() 生成一个半径为 1 的球面网格
surf(Xs, Ys, Zs, 'FaceAlpha',0.1, 'EdgeColor','none', 'FaceColor','c');
% 给球面一个半透明的青色

xlabel('x_0'); ylabel('x_1'); zlabel('x_2');
title('Minimize x_0 + x_1 + x_2 subject to x_0^2 + x_1^2 + x_2^2 \le 1');

% 2.2 绘制迭代轨迹
x0_hist = x_history(:,1);
x1_hist = x_history(:,2);
x2_hist = x_history(:,3);

nPoints = size(x_history,1);
if nPoints > 1
    % 中间过程点用蓝色散点
    plot3(x0_hist(1:end-1), x1_hist(1:end-1), x2_hist(1:end-1), ...
        'bo-', 'LineWidth',1.5, 'MarkerSize',4, 'MarkerFaceColor','b');
end

% 最后一点用红色标记
plot3(x0_hist(end), x1_hist(end), x2_hist(end), ...
      'ro', 'MarkerSize',8, 'MarkerFaceColor','r');

legend('Unit Sphere (Constraint)', 'Iter Process', 'Final Solution');
view(35, 25);  % 调整3D视角