运筹系列31：内点法python代码

1. 简介

用内点法求解线性规划问题理论上的计算复杂度为$O(n^{3.5}{L}^2)$，其中n是变量的维数，L是输入长度。而单纯形法本质上还是个搜索问题，其计算复杂度是$O(2^n)$。
内点法总结起来有两大类，如下：
（1）使用拉格朗日法将不等式去除，然后使用KKT条件将原问题转为方程组，然后用牛顿法求解。
（2）类似信赖域方法，每次在一定范围（比如使用尺度变换生成一个球）内移动到最优位置，迭代进行。
基本概念可以参见：https://blog.csdn.net/WASEFADG/article/details/90376051

本文依次介绍Logarithmic barrier method，Affine scaling method，Projective method。

2. 经典的Barrier法

详细的介绍可以参考 https://zhuanlan.zhihu.com/p/32685234，对应的代码https://github.com/PrimerLi/linear-programming/tree/master/src
下面两张图是原理：


为了能够找到一个初始可行解，添加一个barrier函数：

一般选用Logarithmic barrier method，如下。其中$\varphi (x)=-{\Sigma }_i\ln (x_i)$称为barrier函数。至于为什么用这个函数做barrier，可以参考这里

使用KKT条件进行转化，将最优化问题转换为求解方程。

使用牛顿法求解，不同$\mu$下的最优解构成的集合称为central path。
这里我们假设问题为线性：
$\min c^{T}x$，s.t. $Ax\le b$
barrier函数设置为$I_t(x)=-log(-u)/t$，目标函数为$tcx-\Sigma \log (-Ax+b)$，并且有：


内点法是选取一个比较小的初始参数$t_0$，求出这个参数对应的解$x_{t_0}$，逐步迭代从小$t_0$ 到大 $t$可以得到一系列的解，最后得到的解$x_t$称作是淬火解。这个解应该可以满足我们的精度需求。

def newton(t, A, b, c, x0):
    maxIterationNumber = 50
    eps = 1.0e-8
    for count in range(maxIterationNumber):
        slack = 1.0/(A.dot(x0) - b)
        gradient = t*c - A.transpose().dot(slack) # 梯度
        H = A.transpose().dot(np.diag(np.square(slack))).dot(A) 
        delta = np.linalg.inv(H).dot(gradient)
        x0 = x0 - delta
        error = np.linalg.norm(delta)
        if (-error < eps):
            break
    return x0,error

def solver(t0, A, b, c, x0, factor):
    upperBound = 1000*t0
    t = [t0]
    solutions = [x0]
    eps = 1.0e-10
    while(t0 < upperBound):
        x0,error = newton(t0, A, b, c, x0)
        solutions.append(x0)
        print (x0)
        t0 = factor*t0
        t.append(t0)
        if (error < eps):
            break
    return x0

下面是个例子：
$\min x+y$，
s.t.
$0.81x$