单纯形法是求解线性规划问题最常用的算法，在实际应用中，特别是解决中小规模的线性规划问题时，其有相当优异的表现．但Klee-Minty举出了单纯形法最差的一个方面，即导致迭代次数较问题的规模呈指数倍增长[1]．从那时起，研究者们开始不断尝试找到更好、更快的方法来解决线性规划问题．随后的研究推动了内点法的发展，内点法在一些大规模的线性规划问题求解中显著优于单纯形算法．尽管如此，文献[2-3]证明了单纯形法求解线性规划问题的迭代次数的期望值是多项式形式．此外，内点法在重新求解带扰动数据的线性规划问题中未表现出如单纯形法一样的效率[4-5]．因此，算法的核心部分，即核心思想的改进仍然非常重要．当前，关于单纯形法核心法则的研究越来越受到线性规划领域研究者们的关注[6]，更不乏一些将单纯形法和内点法相结合的成功尝试[7-8]，关于原始对偶算法收敛性的一些研究，其计算结果亦令人鼓舞[9-10]．为获得更高的计算效率并实现进一步的突破，有必要改变单纯形法和内点法的基本思路．基于此，文献[11]提出了一种改进后的通用算法用于线性规划问题，该对偶原始算法(DPA)通过正交三角(QR)分解处理一系列约束最小二乘问题来确定搜索方向，此搜索方向与原始单纯形法、对偶单纯形法、原始对偶算法和内点法完全不同．该方法从任意原始可行解(不一定从基本解或严格的内解)开始，然后沿着严格的下降方向(不一定沿边缘或内部)移动到目标函数．DPA不受原始简并性影响，因此对经典单纯形法的退化问题更有效．此外，DPA可以同时得到一对精确的原始和对偶最优解，而无需内点法所要求的繁琐过程．理论上，DPA结合了单纯形法和内点法的优点．文献[11]中给出了一组随机构造的问题，其初步数值结果优异．然而，虽然该对偶原始算法在理论上非常成功，并且比传统的单纯形法需要迭代次数更少，但是须要进一步考虑计算实现的问题．该算法的计算瓶颈在于如何获得相关约束最小二乘问题的残差，若利用正常的非线性优化方法，则在每次迭代中耗时很多，失去了算法的优越性．本研究给出该对偶原始算法的一种有效实现方法，即提出了一种主要处理一系列无约束最小二乘问题的子算法，而不是通过经典非线性优化算法来处理约束最小二乘问题．此外，该子算法中使用的数据结构与文献[11]中提出的对偶原始算法的主体部分相同，可以很容易地将其嵌入到主程序中，因此大大减少了消耗的计算时间，降低计算成本．本研究先介绍解决线性规划问题的对偶原始算法；然后描述了一种通过求解一系列无约束最小二乘问题来计算搜索方向的子算法，并将子算法结合到主过程中，从而获得新版本的对偶原始算法；最后，通过给出一组随机构造的问题结果的说明性示例和数值展示所提出新方法的优越性和有效性．1 对偶原始算法考虑线性规划问题的标准形式mincTx;s.t. Ax=b,x≥0，(1)及其对偶问题maxbTy;s.t. ATy+s=c,s≥0，(2)其中c,x,s∈Rn，b,y∈Rm，A∈Rm×n，mn，rank(A)=m．下面简要介绍对偶原始算法．设x˙是问题(1)的原始可行解，M和N分别为两个指标集：M={i|x¯i0,i∈{1,2,⋯,n}}; (3)N={i|x¯i=0,i∈{1,2,⋯,n}}．(4)系数矩阵A和s的单位系数矩阵可以相应地分别划分为[AM,AN]和[IM,IN]两个部分．向量x，s和c也被一致地划分．现在考虑约束最小二乘问题minsN≥0||[AT,IN][y,sN]T-c||2．(5)对偶原始算法主要步骤如下：步骤1 设x˙是问题(1)的原始可行解；步骤2 计算上述问题(5)的解的残差r；步骤3 当残差r=[AT,IN][y,sN]T-c=0时停止计算，此时的解x˙和s¯=0    (j∈M),s¯j   (j∈N)分别是原始最优解和对偶最优解；步骤4 当r≥0时停止计算，问题(1)无下界；步骤5 通过α=min{-x˙j/rj|rj0,j∈M}确定步长α；步骤6 令x˙:=x˙+αr，然后从步骤1开始循环．从问题(1)的初始可行解开始，DPA经过有限步骤后在以下任一处终止：a. 步骤3，达到原始最优解和对偶最优解；或者b. 步骤4，检测问题(1)的下无界性．有关上述结果的详细证明，参见文献[11]．2 对偶原始算法的有效实现2.1　对偶原始算法的局限性用于线性规划的对偶原始算法可有效避免简并性问题且收敛于有限步数．与传统的单纯形法相比，它在退化问题上效率更高．对偶原始算法可衍生出许多其他线性规划算法，包括两种最成功的方法——单纯形法和内点法的变体，以及与二者完全不同的其他方法．理论上讲，大多数线性规划算法都可以通过选择不同情况变换后得到．例如，求解线性规划问题的Criss-cross算法就是一个交替进行原始和对偶迭代的算法，它也是对偶原始算法的一种特殊形式[12]．然而，尽管如上所述和文献[11]中的讨论，对偶原始算法拥有诸多优点，但仍然只能说其理论上有助于线性规划求解．该算法所需的迭代步骤通常比单纯形法及其他形式少得多，但就CPU时间而言，它有时却比后二者更耗时．对偶原始算法的核心是计算下降方向，即求约束最小二乘问题(5)的残差，这也是在每次迭代中最耗时的部分．如果能够较好地处理这个计算瓶颈，那么对偶原始算法的效率将大大提高．受此启发，对文献[11]中提出的对偶原始算法做了进一步改进，提出了一种新的对偶原始算法(NDPA)，该算法只须要处理一系列无约束的最小二乘问题，并且核心程序的数据结构不变．2.2　子算法令x˙是问题(1)的原始可行解．这里首先展示其对应问题(5)的一些基本结果．假定(y˜T,s˙NT)T是问题(5)的一个可行解，r˜是对应的残差．现在考虑一个辅助的无约束最小二乘问题min||[AT,IN][y,zN]T+r˜||2，(6)式中：y∈Rm；zN∈RN;N为集合N的势，集合zN的元素须满足当(s˜N)j=0时，(zN)j=0．现在假定(y˙T,z˙NT)T是问题(6)的解并将其记作d˙=(y˙T,z˙NT)T，(7)即可得出以下结论．引理1 令d˙如式(7)中定义，若d˙≠0，则d˙是问题(5)在(y˙T,s˙NT)T的可行下降方向．证明 显然，根据问题(6)的定义，若d˙≠0，则d˙是问题(5)在(y˙T,s˙NT)T的一个可行方向，因为若(s˜N)j=0，则(zN)j=0．另一方面，通过最小二乘问题的基本理论结果，已知(y˙T,z˙NT)T是问题(6)的解，这意味着[AT,IN]T[AT,IN](y˙T,z˙NT)T=-[AT，IN]Tr˜．因此得到(y˙T,z˙NT)T(-[AT,IN]Tr˜)=(y˙T,z˙NT)T∙[AT,IN]T[AT,IN](y˙T,z˙NT)T≥0. (8)注意[AT,IN]Tr˜恰好是问题(5)的目标函数在点(y˜T,s˜NT)T的梯度．因此d˙=(y˙T,z˙NT)T是问题(5)在(y˜T,s˜NT)T的下降方向．若d˙=0，无法通过上述方法得到可行的下降方向，于是可利用目标函数在当前点的梯度向量构造一个可行的下降方向．令g=[AT,IN]Tr˜表示目标函数在当前点的梯度向量．对g作分割，得到g=(g1T,g2T)T，其中g1∈Rm，g2∈RN．根据w=(w1T,w2T)T∈Rm+N定义w，其中w1=g1∈Rm,w2∈RN，(w2)j=0    ((g2)j≥0);1    ((g2)j0). (9)引理2 若d˙=0，w=0，则当前点(y˜T,s˜NT)T是问题(5)的最优解．证明 考虑问题(5)的拉格朗日算子L(z,λ)=||A˜Tz-c||2-λTz，式中：A˜T=[AT,IN]；z=[y,sN]T．选取λ*=2g作为对应的拉格朗日乘子．下面说明当前点满足KKT条件．首先，∇z(||A˜Tz-c||2-λTz)=2A˜(A˜Tz-c)-λ*=2A˜r˜-λ*=2g-2g=0；其次，根据以上w的定义知w=0的条件表明g1=0，g2≥0．因此，λ*=2g≥0．接下来检验最后一个条件λ*Tz*=0，式中z*=[y*,sN*]T．定义指标集J={j(s˜N)j}=0，记问题(6)中的矩阵为A˜(J)T，即从矩阵A˜T中删去对应指标集J的列项．所以，问题(6)可写作一个无约束最小二乘问题min||A˜(J)Tu+r˜||2，(10)其中u∈Rm+N-J．由最小二乘问题的基本理论，有-A˜(J)r˜=A˜(J)A˜(J)Tu*，(11)式中u*为问题(10)的最优解．进一步，条件d˙=(yT,zNT)T=0意味着u*=0．因此，由式(11)有-A˜(J)r˜=0．注意到λ*=2g=2A˜r˜且-A˜(J)r˜与使得(s˜N)j0成立的j相关，可得结果λ*Tz*=0．至此证明完成．若d˙=0，w≠0，表明g≠0．不失一般性，假设gi≠0，定义d˙=-giei，(12)式中ei∈Rm+N为第i个分量为1、其他分量为0的单位向量．此时显然上述定义的向量d˙是当前点(y˜T,s˜NT)T的可行下降方向．由以上讨论，可构造如下子算法．令(y˜T,s˜NT)T为问题(5)的初始可行解．步骤1 求对应的辅助无约束最小二乘问题(6)的解(y˙T,z˙NT)T；步骤2 若(y˙T,z˙NT)T≠0，定义搜索方向d˙=(y˙T,z˙NT)T，转到步骤5；步骤3 若式(9)中定义的向量w为零，运行停止，当前点(y˜T,s˜NT)T即问题(5)的最优解；步骤4 由式(12)定义搜索方向d˙；步骤5 由α=min{-(s˜N)j/(d˙)j|(d˙)j0，j∈{m+1,m+2,⋯,m+|N|}}定义步长α；步骤6 更新(y˜T,s˜NT)T,(y˜T,s˜NT)T:=(y˜T，s˜NT)T+αd˙，回到步骤1继续进行．上述子算法迭代步骤是有限的，因为目标函数在每次迭代和每次主动约束(即有效约束)数量增加时都严格减少．此外，有效约束集(满足所有有效约束的点的集合)的数量亦有限．2.3　新对偶原始算法下面完整地阐述本文提出的经改进后的新对偶原始算法．设x˙是问题(1)的原始可行解．步骤1 分别由式(3)和(4)定义M和N；步骤2 对任意给定的点(yT,zNT)T，zN≥0(问题(5)的一个初始可行解)，调用子算法计算问题(5)的最优解(y*T,z*NT)T；步骤3 若(y*T,z*NT)T=0，则停止计算，当前解x˙最优；步骤4 定义搜索方向r=[AT,IN]y*Tz*NT-c；步骤5 当r≥0时停止计算，此时问题(1)无下界；步骤6 由α=min{-x¯j/r˙j|r˙j0,j∈M}确定步长α；步骤7 更新x¯，令x¯:=x¯+αr，回到步骤1开始循环．NDP算法只须要处理一系列无约束最小二乘问题，而对于无约束最小二乘问题，目前已经有解决它们的有效方法，在此不再赘述．正因如此，NDP算法能以更低的时间成本找到搜索方向．3 计算结果在此，将修正后的单纯形法和NDP算法进行比较，简述部分数值测试并报告计算结果，进而验证NDP算法在随机生成的线性规划问题上具有较高效率．用于实验的随机生成的线性规划问题均为标准形式．所执行并比较的两个算法说明如下：代码RSA为修正版单纯形法，代码NDPA为第一阶段的NDP算法．以上程序均用PC系统的MATLAB 9.4编译和运行，所有的计算都在一个处理节点上执行而不进行任何切换．试验所用的计算机为Windows7操作系统，处理器为2.3 GHz的Intel处理器，内存4 GB．为观察和进一步了解改进后的NDP算法在数值结果上的表现，将其与修正版的单纯形法比较．首先，分别采用这两种算法解决小规模问题；然后，采用一组随机生成的线性规划问题进行更进一步的数值实验．样例1考虑以标准形式(1)随机生成的线性规划问题A=[A7×41,A7×42,A7×43]7×12；b=[b1×41，b1×32]1×7T； c=c1×41,c1×42,c1×431×12T，式中：A7×41=-2.756 1-0.459 31.175 71.712 5-0.870 3-0.526 7-1.546 1-1.331 7-2.422 5-1.160 02.133 01.861 4-0.158 9-1.073 41.859 81.778 72.512 52.155 8-0.604 2-0.977 2-0.042 70.849 7-1.168 4-0.223 54.250 70.785 1-0.971 1-2.587 0；A7×42=-7.263 60.661 9-1.995 92.720 0-0.476 0-1.233 9-0.440 50.746 3-2.708 90.132 5-0.605 42.099 9-0.693 30.809 00.519 50.511 63.665 90.286 30.873 3-2.860 10.378 30.143 60.957 3-1.012 68.926 2-0.077 71.745 9-2.857 1；A7×43=-4.213 4-3.271 63.548 73.156 7-1.035 50.496 80.416 51.236 4-1.470 9-3.020 31.974 11.406 5-0.550 5 -2.051 60.502 81.125 53.378 12.283 7-1.828 4-2.238 90.352 81.154 7-0.180 7-1.766 15.003 54.771 3-4.890 2-3.403 0；b1×41=[0.083 1,0.470 6, 0.215 6, 1.094 9]；b1×32=[0.718 9,0.250 5,0.030 8]；c1×41=[-0.543 9,-0.100 7,-0.655 6,0.937 6]；c1×42=[-0.288 6,0.000 0,0.000 0,0.000 0]；c1×43=[0.000 0,0.000 0,0.000 0,0.000 0]．NDPA的求解过程如下．在第一阶段的主要迭代后，可得到一个可行解x0=[x1×411,x1×412,x1×413]1×12T，式中：x1×411=[0.000 0,0.537 9,0.197 5,0.000 0]；x1×412=[0.253 7,0.516 3,1.052 9,0.226 2]；x1×413=[0.121 4,0.319 2,0.698 3,0.685 1]．然后再做带残差的4次迭代r1=[r1×411,r1×412,r1×413]1×12T；r2=[r1×421,r1×422,r1×423]1×12T；r3=[r1×431,r1×432,r1×433]1×12T;r4=1.0e-14×[r1×441,r1×442,r1×443]1×12T，式中：r1×411=[0.086 5,0.009 4,0.017 1,0.000 0]；r1×412=[0.142 8,0.000 0,-0.061 2,-0.013 0]；r1×413=[-0.133 8,-0.004 7,0.217 2,-0.056 1]；r1×421=[0.106 9,-0.056 5,-0.005 9,0.000 0];r1×422=[0.044 0,0.000 0,-0.014 8,0.036 8]；r1×423=[0.000 0,0.030 1,0.200 1,-0.046 3];r1×431=[0.112 5,0.000 0,-0.040 5,0.000 0];r1×432=[0.043 8,0.000 0,0.008 8,0.089 2];r1×433=[0.000 0,-0.024 6,0.143 2,-0.043 7];r1×441=[0.000 0,-0.023 6,0.022 2,-0.111 0]；r1×442=[0.016 7,-0.022 2,0.000 0,0.000 0]；r1×443=[0.127 7,0.000 0,0.000 0,-0.044 4]，其中r4≈0．因此，可得最优解如下xopt=[x1×421,x1×422,x1×423]1×12T，式中：x1×421=[1.713 2,-0.000 0,0.411 0,0.000 0]；x1×422=[0.697 8,-0.000 0,1.190 3,0.938 4]；x1×423=[0.000 0,0.829 5,3.361 0,0.000 0]．采用NDPA可得最优值为-1.402 6，CPU占用时间为0.21 s．相对地，最终获得相同的最优解和相同的最优目标值，采用修正版单纯形法解决上述问题需要10次迭代，Matlab编码下CPU占用时间为0.24 s．本研究分别采用RSA和NDPA进行了10次试验，获得全部计算结果的统计数据见表1(样例1为该试验组的问题2)，表中：P为试验问题序号；Sm+n为规模，即矩阵A的行和列之和；IRSA和INDPA分别为完全求解每个试验问题时RSA和NDPA所需的主要迭代次数；TRSA和TNDPA分别为运行时间．10.13245/j.hust.210603.T001表1RSA和NDPA计算结果统计PSm+nIRSAINDPATRSA/sTNDPA/s114520.120.132191050.240.213353180.670.584402130.380.505453160.530.4268047110.730.4979052110.330.31814073351.060.65916092372.541.711020067132.211.13从表1可看出：与RSA相比，NDPA需要的平均迭代次数和运行时间均明显减少；换言之，对于大多数问题而言NDPA明显优于RSA，除了问题(1)和(4)在使用前者解决时需要稍微多一些运行时间．进一步地，观察表1可看出：NDPA在运行时间上相对于RSA的优势随测试问题规模的增大而明显增加，即随问题规模的增大，NDPA的迭代次数和运行时间较RSA有明显减少，当问题规模增加到200时，运行时间可减少48.87%，迭代次数仅为RSA的1/3，这一趋势完全符合对NDP算法的期望．4 结语本研究提出了一种计算性能优越的新对偶原始算法(NDPA)，规避了以往线性规划算法的受原始退化影响、迭代次数随规模增大而大幅增长、占用CPU时间较长导致效率低等不足，并通过随机生成的线性规划问题试验，证明了该算法不仅可减少迭代次数，亦可大幅降低运行时间，且随问题规模的增大其计算效率明显提升．此外，在文献[9，13-15]等较为复杂的算法的实现以及文献[16-18]的一些实际应用中，也可尝试加入NDPA的思路，对复杂算法中的一部分做相应改进，进而提升整体计算效率．