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