首页/文章/ 详情

参变量变分原理(二)

4月前浏览3939

参变量变分原理(1)

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

 

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

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

 

再引入松弛变量    ,使

 

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

 

 

代入(2)得

 

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

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

 

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

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

 

求解线性互补问题的    算法

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

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

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

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

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

算例

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

 
 

转为线性互补问题

 

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

表1 初始表

w1w2w3z1z2z3z0q
w1100-21-3-1-1
w20101-4-2[-1]-10
w3001320-16

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

表2 z0入基,w2出基

w1w2w3z1z2z3z0q
w11-10-3[5]-109
w20-10-142110
w30-11362016

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

表3 z2入基,w1出基

w1w2w3z1z2z3z0q
w11/5-1/50-3/51-1/505/9
w2-4/5-1/507/5014/5114/5
w3-6/51/51[28/5]016/5026/5

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

表4 z1入基,w3出基

w1w2w3z1z2z3z0q
w11/14-5/283/28011/7033/14
w2-1/2-1/4-1/400[2]13/2
w3-3/141/285/28104/7013/14

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

表4 z1入基,w3出基

w1w2w3z1z2z3z0q
w13/28-9/567/56010-1/149/4
w2-1/4-1/8-1/80011/23/4
w3-1/143/281/4100-2/71/2

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

 

最优点

 

相应的目标函数最小值为    

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)  


来源:数值分析与有限元编程
pythonUM
著作权归作者所有,欢迎分享,未经许可,不得转载
首次发布时间:2024-06-23
最近编辑:4月前
太白金星
本科 慢慢来
获赞 5粉丝 11文章 325课程 0
点赞
收藏
未登录
还没有评论
课程
培训
服务
行家
VIP会员 学习 福利任务 兑换礼品
下载APP
联系我们
帮助与反馈