参变量变分原理(二)
参变量变分原理(1)利用二次规划问题的最优性条件可以得到线性互补问题。实际上，利用线性规划的最优性条件也可以得到线性互补问题。线性互补问题可以利用基于单纯形法的方法解决，也可以利用内点法求解。接下来给出方法。考虑二次规划问题:式中，是的对称矩阵,是的矩阵,为维向量,为维向量。引入乘子和，定义拉格朗日函数再引入松弛变量，使这样，(1)的KKT条件可写成令代入(2)得其中，均为维列向量，则是阶矩阵。式(3)称为线性互补问题。首先，如果所有的，则可满足式(3)，这就是问题的解。但是,一般情况下，中至少有一个元素。Lemke引人了一个人工变量，构建一个新的方程组:其中，,为单位矩阵。利用上式，可构造一个初始的基本可行解和。这里是中最小的负值元素，如此可保证初始解中，是非负的。算法包括初始步和迭代步。在初始步中，利用枢轴操作将变换为基变量，记出基变量为。接下来，即将入基的非基变量就是上一次迭代中出基变量的互补变量。如果出基，则下一次迭代中入基，反之亦然。这样一来，可始终保持互补条件成立。求解线性互补问题的算法Step1入基，中最小元素对应的元素出基。如果中不存在负元素，则得到解为，。算法停止。Step2针对对应的列和第行，开展初等行变换(枢轴操作)。此时，的所有元素都已经转换为非负值。Step3由于第行对应的基变量出基，其对应的互补变量，即第列对应的元素入基(最开始迭代时，出基，入基)。Step4针对第列元素开展比值计算，计算方式为除以该列中的所有正数元素，确定枢轴变量，据此找到相应的出基变量。如果没有元素为正数，则线性互补问题无解。这种情况称为存在射线解，应停止计算。Step5开展初等行变换，使得枢轴变量为1，该列中的其他元素全部为0。如果上一次枢轴操作可使得基变量出基，则迭代完成，算法停止;否则，转至Step3。算例用方法求解下列凸二次规划问题转为线性互补问题引入人工变量，建立初始表如表1所示。表1初始表w1w2w3z1z2z3z0qw1100-21-3-1-1w20101-4-2[-1]-10w3001320-16入基，出基，以表中加的元素为主元进行消元运算，得换基后的表如表2所示。表2z0入基，w2出基w1w2w3z1z2z3z0qw11-10-3[5]-109w20-10-142110w30-11362016入基，按最小比值规则确定出基，经主元消去法运算，得换基后的表格如表3所示。表3z2入基，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所示。表4z1入基，w3出基w1w2w3z1z2z3z0qw11/14-5/283/28011/7033/14w2-1/2-1/4-1/400[2]13/2w3-3/141/285/28104/7013/14入基,按最小比值规则确定出基,经主元消去法运算,得换基后的表格如表5所示。表4z1入基，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代码importnumpyasnp##AA是上述表格1的矩阵defLemke(AA):number_of_rows=AA.shape[0]colm_of_poivt=AA.shape[1]-2NBS=np.arange(number_of_rows)iter=0max_iter=500row_of_poivt=0piv=0piv=np.min(AA[:,-1])row_of_poivt=np.argmin(AA[:,-1])#判断最后一列全为正数ifpiv==0:print("已经满足条件！")returnwhile1:iter+=1print(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_poivtprint(NBS)c1=1./AA[row_of_poivt,colm_of_poivt]AA[row_of_poivt,:]=c1*AA[row_of_poivt,:]foriinrange(number_of_rows):ifi!=row_of_poivt:AA[i,:]=AA[i,:]-AA[i,colm_of_poivt]*AA[row_of_poivt,:]ifout_base==AA.shape[1]-2:breakelse:colm_of_poivt=number_of_rows+out_baser1=list(filter(lambdax:x>0,AA[:,colm_of_poivt]))iflen(r1)==0:print("射线解!!!")return#以上是判断射线解(全小于0)dd=AA[:,-1]/AA[:,colm_of_poivt]ee=[3.4E38ifx<0elsexforxindd]row_of_poivt=np.argmin(ee)print(ee)print(AA)print(row_of_poivt)ifiter>max_iter:breakx=np.zeros(number_of_rows)foriinrange(number_of_rows):idx=np.argmin(NBS)x[i]=AA[idx,-1]NBS[idx]=888888888print(x)来源：数值分析与有限元编程
