线性规划之单纯形算法

Posted 2022-08-12 clover_hxy

说到单纯形算法，首先就先从线性规划开始介绍。

什么是线性规划？  在给定有限的资源和竞争约束情况下，很多问题都可以表述为最大化或最小化某个目标。如果可以把目标指定为某些变量的线性函数，而且如果可以将资源约束指定为这些变量的等式或不等式，则得到了一个线性规划问题。

求解线性规划的两种常用格式：标准型和松弛型。标准型中所有的约束都上不等式，而在松弛型中约束都是等式。

标准型：    最大化    sigma(ci *xi)

                   约束条件   ai *xi<=bi   

                                     xi>=0

因为标准型的约束都是不等式，那么如果出现了等式怎么办呢？  

例如  aixi==bi   我们可以把这个式子转化成两个不等式  ai*xi<=bi  ai*xi>=bi    这是等价的

那如果出现ai*xi>=bi   的形式又该怎么办呢，其实根据不等式的性质取反即可    即  -ai*xi<=-bi   

因为单纯形算法是根据松弛型来写的所有这里对于标准型的问题不在赘述。

再来看看松弛型：

松弛型就是除了非负约束外，其他条件都是等式

sigma（j=1...n)  aij *xj<=bj

可以转写成    s=bj-sigma(j=1...n)aij *xj   s>=0 的形式  其中s是一个松弛变量。

松弛型的基本形式：

z=v+sigma(j=1...n) cj*xj

for (i=1;i<=m;i++)  xi+n=bi-sigma(j=1....n)aij*xj    其中xi+n表示的就是松弛变量,也叫辅助变量。

对了这里插入一些定义帮助理解

有了以上这些铺垫，我们就可以正式介绍单纯型算法了。

单纯型算法是解决线性规划问题的经典算法，它是建立在矩阵操作的基础上。使用单纯形算法首先要把所有的约束条件转化为松弛型，然后通过迭代不断更新最优解。从线性规划的松弛型中得到基本解：将每个非基本变量设为0，并从等式约束中计算基本变量的值。每轮迭代把一个松弛型转换成一个等价的松弛型。关联的基本可行解的目标值不会小于上一轮的迭代，为了增大目标值，我们选择一个非基本变量，使得从0开始增加变量值，目标值也会增加。变量的能够增加的幅度受限于其他的约束条件。特别的，我们增加它，直到某个基本变量变为0.然后重写松弛型，交换此基本变量与选定的非基本变量。

说太多都不如举个例子来的实在。

最大化  3x1+x2+2x3

满足约束   x1 +x2+2x3<=30

           2x1+2x2+5x3<=24

           4x1+x2+2x3<=36

           x1,x2,x3>=0

转换成松弛型

z=3x1+x2+2x3

x4=30-x1-x2-2x3

x5=24-2x1-2x2-5x3

x6=36-4x1-x2-2x3

注意：一个很关键的地方，我们虽然将最初的形式改成了一个松弛型，但是我们的矩阵A,B,C存储的都是最初的值，而不是移项后的值

每次迭代，我们的目标就是重新整理线性规划，使一个基本解有更大的目标值，我们选择一个目标函数中系数为正的非基本变量xe,而且尽可能增大xe的值而不违反任何约束。变量xe成为基本变量，而某个其他变量xl成为非基本变量。

回到上面的例子，让我们考虑增加x1的值，当x1增加的时候，x4,x5,x6的值都减小。因为每个变量都有一个非负约束，所以我们不允许它们中的任何一个值变成负值。尝试增加x1的值，发现最紧约束是第三个式子，因此我们互换x1和x6

于是得到x1=9-x2/4-x3/2-x6/4,然后把约束中的x1替换，得到一个等价的线性规划。

这个操作过程成为转轴操作：

转轴操作实际上就是把约束条件L变形，表示为x_e=…的形式，然后再逐个等式进行代换。首先对选定的等式L进行代换，因为把L变形的时候需要把x_e的系数化为1，也就是把等式里的每一个元素都除以了A_(L,e)，于是我们第一步就要进行这项操作。注意因为修改别的项的时候要用到A_(L,e)，所以A_(L,e)应该最后被修改，而因为我们把x_e代换成了y_{L（y就是基本变量）}，而y_L的初始系数为1，所以我们应该把A_(L,e)的系数修改为1/A_(L,e)。接下来是修改其它等式，因为我们把x_e代进去了，而所有用来表示x_e的元素都和以前是异号的（因为移到了x_e的异侧），所以要把每一个等式中对应的那个量减去等式L中对应的量乘以在这个等式中x_e的系数，举个例子，要修改等式K中的常数B_K，那么我们应该把B_K修改为B_K-B_L*A_(K,e)。仍然是最后修改A_(K,e)，并且把它修改成y_L应该对应的系数，即-A_(L,e)*A_(K,e)。这里为什么要加个负号呢其实刚开始我也是百思不得其解，最终明白是这样的，我们虽然是交换变量，但实际上不是真的交换，而是通过改变系数的方式来达到同下标，不同线性规划的效果。一最上面的式子为例

x6=36-4x1-x2-2x3  我们要交换的是x1,x6,换完后变成了x1=9-x2/4-x3/2-x6/4  也就是原来的这里一行的A由4,1,2变成了1/4，1/4,1/2,但是我们在往别的式子中代换的时候用的实际的系数是x1=9-x2/4-x3/2-x6/4 中的系数而不是，最初我们读入的不等式的系数，这就是为什么会出现负号，B_K-B_L*A_{(K,e)这里是减也是同理，那里直接用-A(L,e)*A(K,e)是因为x6前面的系数一定为负，而其他的因为带入后式子中不一定只有一项xi，所以是减。这样解释不知道是否能理解。}

最后一步操作是修改目标函数中对应的系数，用v存储那个常数，增量是原先x_e的系数乘上等式L中修改后的常数，即C_e*B_L。然后把其他系数也修改，对于系数C_i仍然是减去等式L中对应的量即A_(L,i)乘以在目标函数中X_e的系数即C_e，最后修改C_e即可。最后当找不到大于0的系数C时，常数v就是答案。

 说道这里，再说都是徒劳了

伪代码如下：

单纯形伪代码
void pivot()

修改选中式子的约束（除以选中项的系数）；
修改选中式子中除选中项的项的系数（除以选中项的系数）；
修改替入变量的系数（替入替出变量互换后其实就是修改了选中项的系数）；


循环每一个式子

除了选中的这个式子&&当前循环到的式子中存在选中项；
将当前式子的常数项修改；
循环当前式子中除选中项的每一项，修改每一项的系数；
修改选中项的系数；



计算目标函数中常数的增量；
修改目标函数中每一项的系数；
修改目标函数中选中项的系数；


void simplex（）

while (1)

在目标函数中找出系数不为0的一项，如果没有结束程序；
在所有的式子中找出包含当前选中项（系数不为0）且最紧的一项；
pivot()；

（感谢 zyf大神的帮助）

c++代码如下：

#include<iostream>  
#include<cstdio>  
#include<cstring>  
#include<algorithm>  
#include<cmath>  
#define inf 1e10  
#define esp 1e-7  
using namespace std;  
int n,m;  
double a[10003][1003],b[10003],c[10010],ans,v;  
void priov(int l,int e)  
  
   b[l]/=a[l][e];  
   for (int i=1;i<=m;i++)  
    if (i!=e)  a[l][i]/=a[l][e];  
   a[l][e]=1/a[l][e];  
   for (int i=1;i<=n;i++)  
    if (i!=l&&fabs(a[i][e])>esp)  
       
        b[i]-=b[l]*a[i][e];  
        for (int j=1;j<=m;j++)  
         if (j!=e)  a[i][j]-=a[i][e]*a[l][j];  
        a[i][e]=-a[i][e]*a[l][e];  
       
    v+=c[e]*b[l];   
    for (int i=1;i<=m;i++)  
     if (i!=e)  c[i]-=c[e]*a[l][i];  
    c[e]=-c[e]*a[l][e];  
  
double simple()  
  
  int l,i,e;  
  double t;  
  while(true)  
     
     for (i=1;i<=m;i++)  
      if (c[i]>esp)  break;  
     e=i;  
     if (e==m+1) return v;  
     t=inf;  
     for (i=1;i<=n;i++)  
      if ((a[i][e]>esp)&&t>b[i]/a[i][e])  
       t=b[i]/a[i][e],l=i;  
     if (t==inf)  return inf;  
     priov(l,e);  
     //cout<<l<<" "<<e<<endl;  
     //cout<<v<<endl;   
     
  
int main()  
  
  scanf("%d%d",&m,&n);  
  for (int i=1;i<=m;i++)  
   scanf("%lf",&c[i]);  
  for (int i=1;i<=n;i++)  
     
     int l,r;  
     scanf("%d%d%lf",&l,&r,&b[i]);  
     for (int j=l;j<=r;j++)  
      a[i][j]++;  
     
  ans=simple();  
  printf("%.0lf",ans);  
  return 0;  
  

另外在说明一个问题：我们在很多时候都会碰到所有基变量为0并不是一组可行解的状况，这时候就会有一个强大的定理排上用场了——对偶定理。

如果目标函数的系数构成的矩阵为C——它是一个行向量，约束条件的常数构成的矩阵为B——它是一个列向量，约束条件中所有的系数构成的矩阵为A，那么我们把A的I,J转置，再把B和C交换，把求最大值变成求最小值（或把求最小值变成求最大值），就得到了一个新的线性规划。它和原来的线性规划的基变量个数可能是不同的，但最优解是相同的。举个例子，对于线性规划：Maximize Cx(S.T. Ax<=B,x>=0)，我们可以把它转化成等价的线性规划：Minmize B^Tx(S.T. A^Ty<=C^T,y>=0)。