二阶锥规划算法

二阶锥规划的定义

二阶锥规划问题的形式为

$\min_{x} f^{T} x$

需满足以下约束

$\begin{matrix} ‖ A_{sc} (i) \cdot x - b_{sc} (i) ‖ \leq d_{sc}^{T} (i) \cdot x - γ (i) \\ A \cdot x \leq b \\ Aeq \cdot x = beq \\ lb \leq x \leq ub . \end{matrix}$

f、x、b、beq、lb 和 ub 是向量，A 和 Aeq 是矩阵。对于每个 i，矩阵 A_sc(i)、向量 b_sc(i) 和 d_sc(i) 以及标量 γ(i) 都处于使用 secondordercone 创建的二阶锥约束中。

换句话说，该问题具有线性目标函数和线性约束，以及一组形式为 $‖ A_{sc} (i) \cdot x - b_{sc} (i) ‖ \leq d_{sc}^{T} (i) \cdot x - γ (i)$ 的二阶锥约束。

`coneprog` 算法

coneprog 求解器使用 Andersen、Roos 和 Terlaky [1] 中描述的算法。该方法是一种类似于内点 linprog 算法的内点算法。

标准格式

该算法首先将问题置于标准形式中。该算法添加了非负松弛变量，因此问题具有以下形式

$\min_{x} f^{T} x$

需满足以下约束

$\begin{array}{l} A \cdot x = b \\ x \in K . \end{array}$

求解器扩大线性系数向量 f 和线性约束矩阵 A 的大小，以考虑松弛变量。

区域 K 是洛伦兹锥 公式 1 与非负正交的叉积。转换每个凸锥

$‖ A_{sc} (i) \cdot x - b_{sc} (i) ‖ \leq d_{sc}^{T} (i) \cdot x - γ (i)$

到洛伦兹锥公式 1，创建变量向量 t₁, t₂, …, t_n+1：

$\begin{matrix} t_{1} = d^{T} x - γ \\ t_{2 : (n + 1)} = A_{sc} x - b_{sc} . \end{matrix}$

这里，每个锥 i 的变量数量 n 是 A_sc(i) 中的行数。根据其定义，变量向量 t 满足不等式

‖ t_{2 : (n + 1)} ‖ \leq t_{1} .

(1)

公式 1 是 (n+1) 个变量中的洛伦兹锥的定义。变量 t 出现在问题中，代替了凸区域 K 中的变量 x。

在内部，该算法还在锥约束的重构中使用了旋转的洛伦兹锥，但本主题不讨论这种情况。有关详细信息，请参阅 Andersen、Roos 和 Terlaky [1]。

添加松弛变量时，算法会根据需要否定变量，并添加适当的常量，以便：

只有一个边界的变量的下界为零。
具有两个边界的变量的下界为零，而使用松弛变量则没有上界。
将没有边界的变量放置在洛伦兹锥中，并以松弛变量作为约束变量。这个松弛变量不是任何其他表达式、目标或约束的一部分。

对偶问题

对偶锥是

$K_{*} = {s : s^{T} x \geq 0 \forall x \in K} .$

其对偶问题是

$\max_{y} b^{T} y$

使得

$A^{T} y + s = f$

其中

$s \in K_{*} .$

对偶最优解是满足对偶约束并最大化对偶目标的点 (y,s)。

齐次自对偶表示

为了处理可能不可行或无界的问题，该算法添加了两个变量 τ 和 κ，并将问题构造为同质（等于零）和自对偶。

\begin{matrix} A x - b τ = 0 \\ A^{T} y + s - f τ = 0 \\ - f^{T} x + b^{T} y - κ = 0 \end{matrix}

(2)

以及约束

(x; τ) \in \bar{K}, (s; κ) \in {\bar{K}}_{*} .

(3)

这里， $\bar{K}$ 是与非负实线相连的锥 K，它是 (x;τ) 的空间。类似地， ${\bar{K}}_{*}$ 是与非负实线相连的锥 $K_{*}$ ，它是 (s;κ) 的空间。在这个表示中，以下引理表明 τ 是可行解的缩放，而 κ 是不可行问题的指标。

引理（[1] 引理 2.1）

让 (x, τ, y, s, κ) 成为公式 2 的可行解，同时满足公式 3 中的约束。

x^Ts + τκ = 0.
如果 τ > 0，则 (x, y, s)/τ 是标准形式二阶锥问题的原始对偶最优解。
如果 κ > 0，则以下严格不等式中至少有一个成立：
b^Ty > 0
f^Tx < 0.
如果第一个不等式成立，则标准形式的原始二阶锥问题不可行。如果第二个不等式成立，则标准形式的对偶二阶锥问题不可行。

总之，对于可行问题，变量 τ 会缩放原始标准形式问题和齐次自对偶问题之间的解。对于不可行的问题，最后的迭代（x、y、s、τ、κ）为原始标准形式问题提供了不可行性的证明。

起点

迭代的起点是可行点：

对于每个非负变量，x = 1；每个洛伦兹锥中的第一个变量为 1；否则为 0。
y = 0。
每个锥为 s = (1,0,…,0)，每个非负变量为 1。
τ = 1。
κ = 1。

中心路径

该算法试图遵循中心路径，即以下方程中 γ 的参数解，从 1 减小至 0。

\begin{matrix} A x - b τ = γ (A x_{0} - b τ_{0}) \\ A^{T} y + s - c τ = γ (A^{T} y_{0} + s_{0} - f τ_{0}) \\ - f^{T} x + b^{T} y - κ = γ (- f^{T} x_{0} + b^{T} y_{0} - κ_{0}) \\ X S e = γ μ_{0} e \\ τ κ = γ μ_{0} . \end{matrix}

(4)

每个变量带有 0 下标，表示该变量的起点。
变量 X 和 S 分别是由 x 和 s 向量形成的箭头矩阵。对于向量 x = [x₁,x₂,…,x_n]，箭头矩阵 X 具有定义
$X = mat (x) = [\begin{matrix} x_{1} & x_{2 : n}^{T} \\ x_{2 : n} & x_{1} I \end{matrix}] .$
根据其定义，X 是对称的。
变量 e 是每个锥坐标中为 1 的向量，对应于 x₁ 洛伦兹锥坐标。
变量 μ₀ 的定义
$μ_{0} = \frac{x_{0}^{T} s_{0} + τ_{0} κ_{0}}{k + 1},$
其中 k 是 x₀ 中非零元素的数量。

中心路径始于起点，终于齐次自对偶问题的最优解。

Andersen、Roos 和 Terlaky [1] 在引理 3.1 中表明，互补条件 x^Ts = 0（其中 x 和 s 位于 Lorentz 锥的乘积 L 中）等价于条件

$X_{i} S_{i} e_{i} = S_{i} X_{i} e_{i} = 0$

对于每个锥 i。这里 X_i = mat(x_i)，x_i 是与洛伦兹锥 i 相关的变量，S_i = mat(s_i)，e_i 是适当维度的单位向量[1,0,0,…,0]。该讨论表明，中心路径在其终点处满足互补条件。

搜索方向

为了当参数 γ 从 1 减小到 0 时获得靠近中心路径的点，该算法使用牛顿法。要查找的变量已标记为 (x、τ、y、s、κ)。让 d_x 表示 x 变量的搜索方向，依此类推。然后，牛顿步骤求解由公式 4 得出的以下线性系统。

$\begin{matrix} A d_{x} - b d_{τ} = (γ - 1) (A x_{0} - b τ_{0}) \\ A^{T} d_{y} + d_{s} - f d_{τ} = (γ - 1) (A^{T} y_{0} + s_{0} - f τ_{0}) \\ - f^{T} d_{x} + b^{T} d_{y} - d_{κ} = (γ - 1) (- f^{T} x_{0} + b^{T} y_{0} - κ) \\ X_{0} d_{s} + S_{0} d_{x} = - X_{0} S_{0} e + γ μ_{0} e \\ τ_{0} d_{κ} + κ_{0} d_{t} a u = - τ_{0} κ_{0} + γ μ_{0} . \end{matrix}$

该算法通过朝 d 方向迈出一步来获得下一个点。

$[\begin{matrix} x_{1} \\ τ_{1} \\ y_{1} \\ \begin{array}{l} s_{1} \\ κ_{1} \end{array} \end{matrix}] = [\begin{matrix} x_{0} \\ τ_{0} \\ y_{0} \\ \begin{array}{l} s_{0} \\ κ_{0} \end{array} \end{matrix}] + α [\begin{matrix} d_{x} \\ d_{τ} \\ d_{y} \\ \begin{array}{l} d_{s} \\ d_{κ} \end{array} \end{matrix}]$

对于某个步骤 $α \in [0, 1]$ 。

为了实现数值稳定性和加速收敛，该算法根据 Nesterov 和 Todd [8] 中的建议来缩放步骤。此外，该算法根据 Mehrotra 的预测校正器 [7] 的变体来校正步骤。（有关更多详细信息，请参阅 Andersen、Roos 和 Terlaky 的 [1]。）

步骤求解器变体

前面的讨论涉及指定了值 'augmented' 的 LinearSolver 选项。求解器具有其他值，可以改变步骤计算以适应不同类型的问题。

'auto' (默认) - coneprog 选择步骤求解器：
- 如果问题是稀疏的，步长求解器为 'prodchol'。
- 否则，步长求解器为 'augmented'。
'normal' - 求解器使用 'augmented' 步骤的变体，适用于问题稀疏的情况。请参阅 Andersen、Roos 和 Terlaky [1]。
'schur' - 求解器使用改进的 Schur 补方法来处理具有少量密集列的稀疏问题。此方法也适用于大型锥。请参阅 Andersen [2]。
'prodchol' - 求解器使用 Goldfarb 和 Scheinberg ([4] 和 [5]) 中描述的方法处理具有少量密集列的稀疏问题。此方法也适用于大型锥。

迭代显示和停止条件

在每次迭代 k 中，该算法计算三个相对收敛测度：

原始不可行性
${Infeas}_{Primal}^{k} = \frac{‖ A x_{k} - b τ_{k} ‖}{\max (1, ‖ A x_{0} - b τ_{0} ‖)} .$
对偶不可行性
${Infeas}_{Dual}^{k} = \frac{‖ A^{T} y_{k} + s_{k} - f τ_{k} ‖}{\max (1, ‖ A^{T} y_{0} + s_{0} - f τ_{0} ‖)} .$
间隙不可行性
${Infeas}_{Gap}^{k} = \frac{| - f^{T} x_{k} + b^{T} y_{k} - κ_{k} |}{\max (1, | - f^{T} x_{0} + b^{T} y_{0} - κ_{0} |)} .$

您可以在命令行中通过指定迭代显示来查看这三个统计数据。

options = optimoptions('coneprog','Display','iter');

当问题可行且求解器收敛时，所有三个都应该趋近于零。对于可行问题，变量 κ_k 趋近于零，变量 τ_k 趋近于正常数。

一个停止条件与间隙不可行性有一定关系。停止条件是当以下最优性测度降低到最优容差差以下时。

${Optimality}^{k} = \frac{| f^{T} x_{k} - b^{T} y_{k} |}{τ_{k} + | b^{T} y_{k} |} = \frac{| f^{T} x_{k} / τ_{k} - b^{T} y_{k} / τ_{k} |}{1 + | b^{T} y_{k} / τ_{k} |} .$

该统计数据衡量了客观值的精确度。

在下列条件下，求解器也会停止并声明问题不可行。三个相对不可行性测度小于 c = ConstraintTolerance，并且

$τ_{k} \leq c \max (1, κ_{k}) .$

如果是 b^Ty_k > 0，则求解器声明原始问题不可行。如果是 f^Tx_k < 0，则求解器声明对偶问题不可行。

当

$μ_{k} \leq c μ_{0}$

和

$τ_{k} \leq c \max (1, κ_{k}) .$

在这种情况下，coneprog 报告问题在数值上不稳定（退出标志 -10）。

当至少一个不可行性测度大于 ConstraintTolerance 并且计算的步长太小时，就会出现剩余的停止条件。在这种情况下，coneprog 报告搜索方向太小，无法取得进一步进展（退出标志 -7）。

参考

[1] Andersen, E. D., C. Roos, and T. Terlaky. On implementing a primal-dual interior-point method for conic quadratic optimization. Math. Program., Ser. B 95, pp. 249–277 (2003). https://doi.org/10.1007/s10107-002-0349-3

[2] Andersen, K. D. A modified schur-complement method for handling dense columns in interior-point methods for linear programming. ACM Transactions on Mathematical Software (TOMS), 22(3):348–356, 1996.

[3] Ben-Tal, Aharon, and Arkadi Nemirovski. Convex Optimization in Engineering: Modeling, Analysis, Algorithms. (1998).

[4] Goldfarb, D. and K. Scheinberg. A product-form cholesky factorization method for handling dense columns in interior point methods for linear programming. Mathematical Programming, 99(1):1–34, 2004.

[5] Goldfarb, D. and K. Scheinberg. Product-form cholesky factorization in interior point methods for second-order cone programming. Mathematical Programming, 103(1):153–179, 2005.

[6] Luo, Zhi-Quan, Jos F. Sturm, and Shuzhong Zhang. Duality and Self-Duality for Conic Convex Programming. (1996).

[7] Mehrotra, Sanjay. “On the Implementation of a Primal-Dual Interior Point Method.” SIAM Journal on Optimization 2, no. 4 (November 1992): 575–601. https://doi.org/10.1137/0802028.

[8] Nesterov, Yu. E., and M. J. Todd. “Self-Scaled Barriers and Interior-Point Methods for Convex Programming.” Mathematics of Operations Research 22, no. 1 (February 1997): 1–42. https://doi.org/10.1287/moor.22.1.1.