一、内点法简介 

二、仿射内点法

ppt

KARMARKAR标准型转换与迭代算法实现

• karmarkar标准型转换

• karmarkar迭代过程

• main函数

karmarkar标准型转换

标准型转换的原理就不说了，参见高等数学规划教材。

def trans_to_karmarkar(a, b, c, m, n):

    # 约束矩阵A 等式右端b(列向量) 目标参数c(列向量) 约束数量m 变量数量n

    at = np.transpose(a)
    ct = np.transpose(c)
    bt = np.transpose(b)

    fu_b = -1 * b
    fu_c = -1 * c

    fu_at = -1*at
    fu_bt = -1*bt

    i_n = np.identity(n)
    zero_1 = np.zeros((m, 2*m+n))  # 第一行的0
    zero_2 = np.zeros((n, n))  # 第二行的0
    zero_3 = np.zeros((1, n))  # 第三行的0
    # print(zero_3)
    zero_4 = np.array([[0]])

    e_total = np.ones(2*m+2*n)
    e_last = np.ones((1, 2*m+2*n+2))

    # print(b.shape,c.shape,zero_3.shape)

    b_gang = np.concatenate((b, c, zero_4), axis=0)
    # print(b_gang.shape)
    fu_b_bang = -1*b_gang

    # 按行拼接再按列拼接
    A_plus1 = np.concatenate((a, zero_1), axis=1)
    # print(A_plus1.shape)
    A_plus2 = np.concatenate((zero_2, at, fu_at, i_n), axis=1)
    # print(A_plus2.shape)
    A_plus3 = np.concatenate((ct, fu_bt, bt, zero_3), axis=1)
    # print(A_plus3.shape)
    A_plus = np.concatenate((A_plus1, A_plus2, A_plus3), axis=0)
    # print(A_plus.shape)
    B_ae = b_gang - np.reshape(np.dot(A_plus, e_total), (m+n+1, 1))

    # print(B_ae.shape)

    B_1 = np.concatenate((A_plus, B_ae, fu_b_bang), axis=1)
    # print(B_1.shape)
    B = np.concatenate((B_1, e_last), axis=0)

    c_aim_1 = np.zeros((2*m+2*n, 1))
    temp = np.array([[1], [0]])
    c_aim = np.concatenate((c_aim_1, temp), axis=0)

    return B, c_aim

karmarkar迭代过程

迭代过程的算法也不说了，重要的要提一下迭代的收敛条件：

• 有些教材里给的是，其中是足够大正整数，

• 这里倾向于目标值变动比率足够小。
  在迭代过程中会出现步长越界的情况，需要调整参数alpha和上述两个收敛条件的关键参数，用来权衡越界情况和目标值优化情况。

def karmarkar(B, c_aim):
    # 近似解的效果严重依赖于参数L 和 步长alpha
    # 收敛条件可采用L参数法 或比较目标函数值的变动比率
    # 调节两个参数的时候 alpha相当于learning rate 尽量小
    # L或者变动比率 L越小即绝对值越大收敛条件越苛刻 容易出现步长越界
    # 同理变动比率越小越容易越界 但越容易逼近目标最优值

    m = len(B)
    n = len(B[1])
    A = B[:m-1]
    # print("输入A",A)
    r = 1/math.sqrt(n*(n-1))
    # print('r',r)
    # print(A)

    L = -22
    alpha = 1/32
    X_0 = np.ones((n, 1))/n
    # print("x0",X_0)
    # X_0=np.array([[0.6220],[1/3],[0.0447]])
    Y_ORIGINAL = np.dot(np.transpose(c_aim), X_0)
    # print("y00",Y_ORIGINAL)
    FLAG = Y_ORIGINAL*math.pow(2, L)  # 这个是教材给的收敛条件用的
    # print("FLAG",FLAG)
    Y_each=[]
    X_each=[]
    print("--------------迭代开始---------------\n")
    while 1:
        X_1 = karmarkar_one_step(X_0, n, A, c_aim, alpha, r)
        Y_0 = np.dot(np.transpose(c_aim), X_0)
        Y_1 = np.dot(np.transpose(c_aim), X_1)

        assert Y_1 < Y_0  # 步长越界时该处报错

        X_0 = X_1
        Y_each.append(Y_1)
        X_each.append(X_1)
        print("当前目标函数值：", Y_1)

        bias = (Y_0-Y_1)/Y_1

        if bias < 0.0014:
            print("--------------迭代结束---------------")
            print("最优近似变量值\n", X_1, "\n最优近似 karmarkar目标值\n", Y_1)
            break
        else:
            continue
    X_LAST = X_1
    #Y_LAST = np.dot(np.transpose(c_aim),X_LAST)
    # print('original X AND Y',X_LAST,Y_LAST)

    return X_LAST,Y_each,X_each

返回的是最优karmarkar标准型的变量向量，每一步迭代的karmarkar函数值和变量值（为了后边输出迭代过程的图表）。
每一步的迭代过程

def karmarkar_one_step(X_0, n, A, c_aim, alpha, r):
    # 可行点X0 m约束数量 n变量数量 A约束阵 c_aim目标函数参数 alpha

    e_n = np.ones((1, n))
    # print("e_n",e_n)
    X_0 = np.reshape(X_0, (n,))
    # print("X_0",X_0)
    D = np.diag(np.transpose(X_0))

    # print("d:",D)
    D_ni = np.linalg.inv(D)
    # print("D_ni",D_ni)
    A_1 = np.dot(A, D)
    # print("A",A,"A_1",A_1)
    B = np.concatenate((A_1, e_n), axis=0)
    # print("B",B)
    B_t = np.transpose(B)
    temp_B = np.linalg.inv(np.dot(B, B_t))
    temp_pb = np.dot(B_t, temp_B)
    pb = np.dot(temp_pb, B)
    # print("pb",pb)
    n_1 = len(pb)
    n_2 = len(pb[1])
    # print("长度",n_1,n_2)

    I = np.identity(n_1)
    temp_dy = I-pb
    # print("temp_dy",temp_dy)
    dy = -1*np.dot(np.dot(temp_dy, D), c_aim)
    # print("dy",dy)
    dy_norm = np.linalg.norm(dy)

    e = np.ones((n_1, 1))/n_1
    # print('e',e)
    et = np.transpose(np.ones((n_1, 1)))
    Y = e+alpha*r*dy/dy_norm
    # print('Y',Y)
    X = np.dot(D, Y)/np.dot(np.dot(et, D), Y)

    return X

main函数

def main():
    # m是约束条件数量 n是变量数量
    
    ##第六题参数集
    a=np.array([[3,8,-1,0],
                [5,2,0,-1],
                ])
    b=np.array([[4],
        [7]])
    c=np.array([[9],[6],[0],[0]])
    m=2
    n=4
    ##第七题参数集
    '''a = np.array([[1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0],
                  [0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 1],
                  [0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0],
                  [0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0],
                  [0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0],
                  [1, 0, 1, 0, -1, -1, -1, 0, 0, 0, 0, 0],
                  [0, 1, 0, 1, 0, 0, 0, -1, -1, -1, 0, 0]])
    b = np.array([[200000],
                  [60000],
                  [50000],
                  [100000],
                  [50000],
                  [0],
                  [0]])
    c = np.array([[0], [5], [4], [2], [3], [4], [5], [2], [1], [2], [0], [0]])
    m = 7
    n = 12'''

    '''
    ##测试参数集
    a=np.array([[1,1,1,1,0],
                [-1,1,2,0,1],
                ])
    b=np.array([[5],
        [6]])
    c=np.array([[1],[-1],[0],[0],[0]])
    m=2
    n=5'''

    B, c_aim = trans_to_karmarkar(a, b, c, m, n)

    # print("karmarkar_约束矩阵：\n",B)
    # print("karmarkar_目标系数：\n",np.transpose(c_aim))
    '''B = np.array([[1,-2,1],
        [1,1,1]])
    c_aim=np.array([[0],[0],[1]])'''

    X,karmarkar_Y_each_step,karmarlar_X_each_step = karmarkar(B, c_aim)
    X_LAST = X[:n]/X[-1]
    X_each_step = list(map(lambda x: x[:n]/x[-1],karmarlar_X_each_step))
    Y_each_step = list(map(lambda x: np.dot(np.transpose(c), x),X_each_step))
    Y = np.dot(np.transpose(c), X_LAST)
    print("原问题变量:\n", X_LAST, "\n原问题目标值:\n", Y)

    path='./save_karmarkar_original.csv'
    with open(path,'w',newline='') as f:
        writer=csv.writer(f)
        for i in range(len(Y_each_step)):
            writer.writerow([karmarkar_Y_each_step[i][0][0],Y_each_step[i][0][0]])


    '''path='./save.csv'
    with open(path,'w',newline='') as f:
        writer=csv.writer(f)
        for i in range(len(B)):
            writer.writerow(B[i])'''
    # print(len(B),len(B[1]))
    fig = plt.figure(constrained_layout=True)
    gs = fig.add_gridspec(ncols=1,nrows=2)

    x_data = range(len(Y_each_step))

    fig_ax1 = fig.add_subplot(gs[0,0])
    fig_ax1.plot(np.reshape(x_data,901),np.reshape(Y_each_step,901))
    fig_ax2 = fig.add_subplot(gs[1,0])
    fig_ax2.plot(np.reshape(x_data,901),np.reshape(karmarkar_Y_each_step,901))
    plt.show()

顺便连图一起画了，画得粗糙，有时间再优化。

内点法 (Interior Point Method)是一种求解线性规划或非线性凸优化问题的算法。它是由John von Neumann发明的，他利用戈尔丹的线性齐次系统提出了这种新的求解线性规划的方法。后被Narendra Karmarkar于1984年推广应用到线性规划，即Karmarkar算法。

1972年，V.Klee和G.Minty给出一个例子，他们构造一个线性规划问题，用单纯形方法求解，需要的计算时间为，其中是变量的维数。这个例子表明，单纯形算法虽然在实际应用中非常有效，至今占有绝对优势，但在理论上它还不是多项式算法。于是产生这样的问题：对于线性规划，能否找到多项式算法？

1979年，前苏联数学家第一个给出了求解线性规划的多项式算法，这就是所谓的椭球算法。它的计算复杂性为，其中是变量的维数，是输入长度。这个算法在理论上是重要的，但是计算结果很不理想，远不及单纯形方法有效。

算法上突破性的进展和当代科学技术发展的需要，又给人们提出进一步的问题：能否找到实用上也确实有效的多项式时间算法？正是在这样的背景下，产生了内点法，它的计算复杂性是。

内点法中有一个惩罚函数，用于描述凸集。与单纯形法不同，它通过遍历内部可行区域来搜索最优解。我们将在约束非线性规划内容讲到。

1984年，印度数学家N.Karmarkar针对线性规划问题提出了一种新的多项式时间算法，在实际计算效率方面，Karmarkar算法显示出可与单纯形法竞争的巨大潜力，Karmarkar算法的提出是线性规划理论研究的突破，而且对于处理非线性优化问题也显示出强大的生命力和广阔的应用前景。

Karmarkar算法是解如下形式的LP：

其中为行满秩矩阵，.设(1)有内点可行解且其最优值为零.