参变量变分原理(二)
参变量变分原理(1)利用二次规划问题的最优性条件可以得到线性互补问题。实际上，利用线性规划的最优性条件也可以得到线性互补问题。线性互补问题可以利用基于单纯形法的方法解决，也可以利用内点法求解。接下来给出 方法。考虑二次规划问题: 式中， 是 的对称矩阵, 是 的矩阵, 为 维向量, 为 维向量。引入乘子 和 ，定义拉格朗日函数 再引入松弛变量 ，使 这样，(1)的KKT条件可写成 令 代入(2)得 其中， 均为 维列向量， 则是 阶矩阵。式(3)称为线性互补问题。首先，如果所有的 ，则 可满足式(3)，这就是问题的解。但是,一般情况下， 中至少有一个元素 。Lemke引人了一个人工变量，构建一个新的方程组: 其中， , 为单位矩阵。利用上式，可构造一个初始的基本可行解 和 。这里 是 中最小的负值元素，如此可保证初始解中， 是非负的。算法包括初始步和迭代步。在初始步中，利用枢轴操作将 变换为基变量，记出基变量为 。接下来，即将入基的非基变量就是上一次迭代中出基变量的互补变量。如果 出基，则下一次迭代中 入基，反之亦然。这样一来，可始终保持互补条件 成立。 求解线性互补问题的 算法Step 1 入基， 中最小元素 对应的元素 出基。如果 中不存在负元素，则得到解为 ， 。算法停止。Step 2 针对 对应的列和第 行，开展初等行变换(枢轴操作)。此时， 的所有元素都已经转换为非负值。Step 3 由于第 行对应的基变量出基，其对应的互补变量，即第 列对应的元素入基(最开始迭代时， 出基， 入基)。Step 4 针对第 列元素开展比值计算，计算方式为 除以该列中的所有正数元素，确定枢轴变量，据此找到相应的出基变量。如果没有元素为正数，则线性互补问题无解。这种情况称为存在射线解，应停止计算。Step 5 开展初等行变换，使得枢轴变量为1，该列中的其他元素全部为0。如果上一次枢轴操作可使得基变量 出基，则迭代完成，算法停止;否则，转至Step3。算例用 方法求解下列凸二次规划问题 转为线性互补问题 引入人工变量 ，建立初始表如表1所示。表1 初始表w1w2w3z1z2z3z0qw1100-21-3-1-1w20101-4-2[-1]-10w3001320-16 入基， 出基，以表中加 的元素为主元进行消元运算，得换基后的表如表2所示。表2 z0入基，w2出基w1w2w3z1z2z3z0qw11-10-3[5]-109w20-10-142110w30-11362016 入基，按最小比值规则确定 出基，经主元消去法运算，得换基后的表格如表3所示。表3 z2入基，w1出基w1w2w3z1z2z3z0qw11/5-1/50-3/51-1/505/9w2-4/5-1/507/5014/5114/5w3-6/51/51[28/5]016/5026/5 入基，按最小比值规则确定 出基，经主元消去法运算，得换基后的表格如表4所示。表4 z1入基，w3出基w1w2w3z1z2z3z0qw11/14-5/283/28011/7033/14w2-1/2-1/4-1/400[2]13/2w3-3/141/285/28104/7013/14 入基,按最小比值规则确定 出基,经主元消去法运算,得换基后的表格如表5所示。表4 z1入基，w3出基w1w2w3z1z2z3z0qw13/28-9/567/56010-1/149/4w2-1/4-1/8-1/80011/23/4w3-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) 来源：数值分析与有限元编程
