参变量变分原理(二)

参变量变分原理(1)

利用二次规划问题的最优性条件可以得到线性互补问题。实际上，利用线性规划的最优性条件也可以得到线性互补问题。线性互补问题可以利用基于单纯形法的方法解决，也可以利用内点法求解。接下来给出    方法。考虑二次规划问题:

式中，    是    的对称矩阵,    是    的矩阵,    为    维向量,    为    维向量。

引入乘子    和    ，定义拉格朗日函数

再引入松弛变量    ，使

这样，(1)的KKT条件可写成

令

代入(2)得

其中，    均为    维列向量，    则是    阶矩阵。

式(3)称为线性互补问题。首先，如果所有的    ，则    可满足式(3)，这就是问题的解。但是,一般情况下，    中至少有一个元素    。Lemke引人了一个人工变量，构建一个新的方程组:

其中，    ,    为单位矩阵。利用上式，可构造一个初始的基本可行解    和    。这里    是    中最小的负值元素，如此可保证初始解中，    是非负的。

算法包括初始步和迭代步。在初始步中，利用枢轴操作将    变换为基变量，记出基变量为    。接下来，即将入基的非基变量就是上一次迭代中出基变量的互补变量。如果    出基，则下一次迭代中    入基，反之亦然。这样一来，可始终保持互补条件    成立。

 

求解线性互补问题的    算法

Step 1    入基，    中最小元素    对应的元素    出基。如果    中不存在负元素，则得到解为    ，    。算法停止。

Step 2 针对    对应的列和第    行，开展初等行变换(枢轴操作)。此时，    的所有元素都已经转换为非负值。

Step 3 由于第    行对应的基变量出基，其对应的互补变量，即第    列对应的元素入基(最开始迭代时，    出基，    入基)。

Step 4 针对第    列元素开展比值计算，计算方式为    除以该列中的所有正数元素，确定枢轴变量，据此找到相应的出基变量。如果没有元素为正数，则线性互补问题无解。这种情况称为存在射线解，应停止计算。

Step 5 开展初等行变换，使得枢轴变量为1，该列中的其他元素全部为0。如果上一次枢轴操作可使得基变量    出基，则迭代完成，算法停止;否则，转至Step3。

算例

用    方法求解下列凸二次规划问题

转为线性互补问题

引入人工变量    ，建立初始表如表1所示。

   入基，    出基，以表中加    的元素为主元进行消元运算，得换基后的表如表2所示。

   入基，按最小比值规则确定    出基，经主元消去法运算，得换基后的表格如表3所示。

   入基，按最小比值规则确定    出基，经主元消去法运算，得换基后的表格如表4所示。

   入基,按最小比值规则确定    出基,经主元消去法运算,得换基后的表格如表5所示。

   已出基，得到互补基本可行解

最优点

相应的目标函数最小值为    

Python代码

import numpy as np

## AA是上述表格1的矩阵

def Lemke(AA):

    number_of_rows = AA.shape[0]
    colm_of_poivt = AA.shape[1] - 2
    NBS= np.arange(number_of_rows )
    
    iter = 0
    max_iter = 500

    row_of_poivt = 0
    piv  = 0    
    piv = np.min( AA[:,-1] ) 
    row_of_poivt = np.argmin( AA[:,-1] )
    # 判断最后一列全为正数
    if piv == 0:
        print("已经满足条件！")
        return
    
    while 1:
        iter += 1
        print(f"第{iter}步:")
        print(f"{colm_of_poivt = }")
        print(f"{row_of_poivt = }")
        out_base = NBS[row_of_poivt]
        print(f"{out_base = }")
        NBS[row_of_poivt] = colm_of_poivt 
        print(NBS)
        c1 = 1. / AA[row_of_poivt, colm_of_poivt]
        AA[row_of_poivt,:] = c1 * AA[row_of_poivt,:]

        for i in range(number_of_rows):
            if i != row_of_poivt:
                AA[i,:] = AA[i,:] - AA[i,colm_of_poivt]  * AA[row_of_poivt, :] 

        
        if out_base == AA.shape[1] - 2 :
            break
        else:
            colm_of_poivt =   number_of_rows  + out_base
        


        r1 = list( filter(lambda x: x > 0,  AA[:,colm_of_poivt] ) )
        if len(r1) == 0:
            print("射线解!!!") 
            return
        # 以上是判断射线解(全小于0)

        dd =  AA[:, -1] / AA[:,colm_of_poivt] 
        ee = [3.4E38 if x < 0 else x for x in dd ]
        row_of_poivt = np.argmin(ee)

        print(ee)

        print(AA)
        print(row_of_poivt)


        if iter > max_iter:
            break

    x = np.zeros(number_of_rows)  
    for i in range(number_of_rows):
        idx  = np.argmin(NBS)
        x[i] = AA[idx, -1]
        NBS[idx] = 888888888
    print(x)