线性规划初探

看完《算法导论》肯定会写单纯形

因为单纯形不仅好写而且《算法导论》里讲的很清楚

附赠uoj179的模板一个

#include<iostream>
#include<cstdio>
#include<algorithm>
#include<cmath>
#include<cstring>
#include<stdlib.h>

using namespace std;
const double eps=1e-9;
double a[50][50];
int b[50],u[50],n,m,ty;

int read()
{
    int x=0,f=1;char ch=getchar();
    while(ch<'0'||ch>'9'){if(ch=='-')f=-1;ch=getchar();}
    while(ch>='0'&&ch<='9'){x=x*10+ch-'0';ch=getchar();}
    return x*f;
}
int dcmp(double x)
{
    if (fabs(x)<=eps) return 0;
    if (x>0) return 1; else return -1;
}

void pivot(int x,int y)
{
     swap(b[x],u[y]);
     double k=a[x][y]; a[x][y]=1;
     for (int i=0; i<=n; i++) a[x][i]/=k;
     for (int i=0; i<=m; i++)
       if (i!=x&&dcmp(a[i][y])!=0)
       {
          k=a[i][y]; a[i][y]=0;
          a[i][0]+=(i?-1:1)*k*a[x][0];
          for (int j=1; j<=n; j++)
            a[i][j]-=k*a[x][j];
       }
}
bool initial()
{
     for (int i=1; i<=n; i++) u[i]=i;
     for (int i=1; i<=m; i++) b[i]=n+i;
     while (1)
     {
           int x=0,y=0;
           for (int i=1; i<=m; i++)
             if (dcmp(a[i][0])<0) x=i;
           if (!x) return 1;
           for (int i=1; i<=n; i++)
             if (dcmp(a[x][i])<0) y=i;
           if (!y) return 0;
           pivot(x,y);
     }  
}
               
int simplex()
{
    if (!initial()) return 0;
    while (1)
    {
          int x=0,y=0; 
          for (int i=1; i<=n; i++)
            if (dcmp(a[0][i])>0) {y=i; break;}
          if (!y) return 1;
          double mi=1e15;
          for (int i=1; i<=m; i++)
            if (dcmp(a[i][y])>0&&(!x||a[i][0]/a[i][y]<mi)) {mi=a[i][0]/a[i][y];x=i;}
          if (!x) return -1;
          pivot(x,y);
    }  
}          
int main()
{
    n=read(); m=read(); ty=read(); 
    for (int i=1; i<=n; i++) a[0][i]=read();
    for (int i=1; i<=m; i++)
    {
        for (int j=1; j<=n; j++) a[i][j]=read(); 
        a[i][0]=read();
    }
    switch (simplex())
    {
           case 1:{
                printf("%.8lf\n",a[0][0]);
                if (ty)
                {
                   for (int i=1; i<=n; i++) u[i]=0;
                   for (int i=1; i<=m; i++) if (b[i]<=n) u[b[i]]=i;
                   for (int i=1; i<=n; i++) printf("%.8lf ",u[i]?a[u[i]][0]:double(0));
                }
                break;
           }
           case 0:puts("Infeasible");break;
           case -1:puts("Unbounded");break;
    }
    system("pause");
    return 0;
}

据说会被卡但实际上基本跑得飞起（下面会用实例证明）

下面是实际应用，凡是最大流，费用流的问题大概都能用线性规划解决，而且会很快变裸题……

比如这个1061，很明显把每天的要求人数bi作为约束，每种志愿者雇佣数量做变量xi

也就是求最小化∑ci*xi i=1..m

并且x要满足约数条件 ∑aij*xi>=bi i=1..n, j=1..m, ai,j=(sj<=i<=ti)?1:0，且xj非负

但是这与标准线性规划刚好相反，标准应该是全是<=且最大化

当然如果你看过《算法导论》关于单纯形最优解证明的话，你就会知道在对偶线性规划下这很简单

对于标准型线性规划：最大化∑ci*xi i=1..n ∑aij*xi<=bi i=1..m, j=1..n，x非负

我们很容易转化成对偶情况：最小化∑bi*xi i=1..m, ∑aij*xi>=bi i=1..n, j=1..m，x非负

这两种线性规划最优值相等（不禁联想到最大流和最小割的关系）

于是我们直接上单纯形即可

但是也许有疑虑，n<=1000,m<=10000能过吗&……

然而答案是c++版本的单纯形跑了1212ms，而我以前pascal版本的费用流跑了1764ms，

所以单纯形还是非常非常厉害的，大胆的使用吧……

#include<iostream>
#include<cstdio>
#include<cstring>
#include<cmath> 
#include<stdlib.h> 

using namespace std;
const double eps=1e-9;
int b[11100],u[11100],n,m;
double a[10010][1010];

int dcmp(double x)
{
    if (fabs(x)<=eps) return 0;
    if (x>0) return 1; else return -1;
}

void pivot(int x,int y)
{
     swap(b[x],u[y]);
     double k=a[x][y];a[x][y]=1;
     for (int i=0; i<=n; i++) a[x][i]/=k;
     for (int i=0; i<=m; i++)
         if (i!=x&&dcmp(a[i][y])!=0)
         {
            k=a[i][y]; a[i][y]=0;
            a[i][0]+=(i?-1:1)*k*a[x][0];
            for (int j=1; j<=n; j++) a[i][j]-=k*a[x][j];
         }
}
                
void simplex()
{
     for (int i=1; i<=n; i++) u[i]=i;
     for (int i=1; i<=m; i++) b[i]=i+n;
     while (1)
     {
           int x=0,y=0;
           for (int i=1; i<=n; i++)
             if (dcmp(a[0][i])>0) {y=i;break;}
           if (!y) break;
           double mi=1e20;
           for (int i=1; i<=m; i++)
             if (dcmp(a[i][y])>0&&(!x||mi>a[i][0]/a[i][y])) {x=i; mi=a[i][0]/a[i][y];}
           if (!x) break;
           pivot(x,y);
     }  
}           
     
int main()
{
    scanf("%d%d",&n,&m);
    for (int i=1; i<=n; i++) scanf("%lf",&a[0][i]);
    for (int i=1; i<=m; i++)
    {
        int s,t,c;
        scanf("%d%d%d",&s,&t,&c);
        for (int j=s; j<=t; j++) a[i][j]=1;
        a[i][0]=c;
    }
    simplex();
    printf("%d\n",(int)a[0][0]);
    system("pause");
    return 0;
}

 bzoj3112 题目在这里http://blog.csdn.net/popoqqq/article/details/44312949

很显然和1061大同小异

#include<iostream>
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cmath>
#include<stdlib.h>

using namespace std;
const double eps=1e-9;
double a[1010][10010];
int n,m;

int dcmp(double x)
{
     if (fabs(x)<=eps) return 0;
     if (x>0) return 1; else return -1;
}

void pivot(int x,int y)
{
     double k=a[x][y]; a[x][y]=1;
     for (int i=0; i<=n; i++) a[x][i]/=k;
     for (int i=0; i<=m; i++)
       if (i!=x&&dcmp(a[i][y])!=0)
       {
          k=a[i][y]; a[i][y]=0;
          a[i][0]+=(i?-1:1)*k*a[x][0];
          for (int j=1; j<=n; j++) a[i][j]-=k*a[x][j];
       }
}
void simplex()
{
     while (1)
     {
           int x=0,y=0;
           for (int i=1; i<=n; i++)
             if (dcmp(a[0][i])>0) {y=i;break;}
           if (!y) break;
           double mi=1e20;
           for (int i=1; i<=m; i++)
             if (dcmp(a[i][y])>0&&(!x||mi>a[i][0]/a[i][y])) {x=i; mi=a[i][0]/a[i][y];}
           if (!x) break;
           pivot(x,y);
    }
}

int main()
{
    scanf("%d%d",&m,&n);
    for (int i=1; i<=m; i++) scanf("%lf",&a[i][0]);
    for (int i=1; i<=n; i++)
    {
        int l,r,d;
        scanf("%d%d%d",&l,&r,&d);
        for (int j=l; j<=r; j++) a[j][i]=1;
        a[0][i]=d;
    }
    simplex();
    printf("%d\n",(int)a[0][0]);
    return 0;
}