内点法（interior point method）求解二次规划，附python代码

河北一帆

已于 2022-07-21 08:07:10 修改

阅读量4.5k

点赞数 4

文章标签：最优化二次规划

于 2022-07-20 11:41:22 首次发布

本文链接：https://blog.csdn.net/qq_41816368/article/details/125888172

版权

内点法介绍这篇博文写的很好https://blog.csdn.net/dymodi/article/details/46441783，在这篇文章的基础上，本文给出障碍函数法代码和一个算例。

障碍函数内点法的主要思想是：把不等式约束放进目标函数里。以下面的问题为例

不等式约束放进目标函数，变成如下形式：

当 $f_{i}(x)>0$ 时，函数 $I_{-}(u)$ 特别大，我们希望 $I_{-}(u)$ 这个函数像下图红虚线的样子，

但那样不可导，所以我们用如下形式的函数去近似。

图中， $u=-3$ 时由下往上的的3根蓝线，对应的 $t$ 值是依次增大。障碍函数内点法在迭代的时候，t初值给一个小值，t逐渐增大，使可行域越来越接近红虚线。

在 t 的每一步迭代中，原问题转化为如下子问题：

求解子问题，得到一个最优解 $x^{*}$ ，这个 $x^{*}$ 是 t 下一步迭代的x初值。

子问题的求解可以用拉格朗日乘子法。

拉格朗日乘子法中， $F(x, \lambda)=0$ 通过牛顿迭代法（Newton-Raphson method）求解。

牛顿迭代法解非线性方程组

$x_{n+1}=x_{n} - \alpha * (H(f))^{-1}\bigtriangledown f$

H(f)是Hessian矩阵。

得到子问题的最优解 $x^{*}$ ，然后继续迭代 t

障碍函数内点法的算法步骤如下：（节选自Stephen Boyd的《Convex Optimization》）

迭代终止条件中的m就是不等式约束的个数那个m，为什么以这个条件为终止的原因，参见文章开头提到的博客，Boyd的书里也有详细的证明，是通过原问题的对偶问题证明得到的。 $\mu$ 是一个大于1的参数，使t每次迭代增大。

算例：

求解如下二次规划问题（解为 x1=0.5，x2=2.25）：

转化为如下子问题：

$t(x_{1}^{2}-x_{1}x_{2}+2x_{2}^{2}-x_{1}-10x^{2})-log(6-3x_{1}-2x_{2})-log(x_{1})-log(x_{2})$

完整python代码如下：

import numpy as np
import time


def grad_f(t, x1, x2):
    return np.array([[t * (2 * x1 - x2 - 1) + 3 / (6 - 3 * x1 - 2 * x2) - 1 / x1],
                    [t * (-1 * x1 + 4 * x2 - 10) + 2 / (6 - 3 * x1 - 2 * x2) - 1 / x2]])


def Hessian_f(t, x1, x2):
    return np.array([[2*t - 9 / pow((6-3*x1-2*x2), 2) + 1 / pow(x1, 2), -t - 6 / pow((6-3*x1-2*x2), 2)],
                    [-t - 6 / pow((6-3*x1-2*x2), 2), 4*t - 4 / pow((6-3*x1-2*x2), 2) + 1 / pow(x2, 2)]])


def NewtonRaphson(t, x1, x2):
    gf = grad_f(t, x1, x2)

    Hf = Hessian_f(t, x1, x2)
    Hf_inv = np.linalg.inv(Hf)

    deltaX = 0.1 * np.matmul(Hf_inv, gf)
    res = np.linalg.norm(deltaX, 2)

    return x1 - deltaX[0, 0], x2 - deltaX[1, 0], res


if __name__ == "__main__":
    time_start = time.time()
    t = 2
    x1 = 1
    x2 = 2
    while True:
        while True:
            x1, x2, res = NewtonRaphson(t, x1, x2)
            if res < 0.0001:
                break

            # print(x1, x2, res)
            # print("------")

        if 3.0 / t < 0.0001:
            time_end = time.time()
            print('consume time:', time_end - time_start)
            print("t:{}, x1:{}, x2:{}".format(t, x1, x2))
            break
        t = 2 * t