线性规划初探

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

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

附赠uoj179的模板一个

  1 #include<iostream>
  2 #include<cstdio>
  3 #include<algorithm>
  4 #include<cmath>
  5 #include<cstring>
  6 #include<stdlib.h>
  7
  8 using namespace std;
  9 const double eps=1e-9;
 10 double a[50][50];
 11 int b[50],u[50],n,m,ty;
 12
 13 int read()
 14 {
 15     int x=0,f=1;char ch=getchar();
 16     while(ch<‘0‘||ch>‘9‘){if(ch==‘-‘)f=-1;ch=getchar();}
 17     while(ch>=‘0‘&&ch<=‘9‘){x=x*10+ch-‘0‘;ch=getchar();}
 18     return x*f;
 19 }
 20 int dcmp(double x)
 21 {
 22     if (fabs(x)<=eps) return 0;
 23     if (x>0) return 1; else return -1;
 24 }
 25
 26 void pivot(int x,int y)
 27 {
 28      swap(b[x],u[y]);
 29      double k=a[x][y]; a[x][y]=1;
 30      for (int i=0; i<=n; i++) a[x][i]/=k;
 31      for (int i=0; i<=m; i++)
 32        if (i!=x&&dcmp(a[i][y])!=0)
 33        {
 34           k=a[i][y]; a[i][y]=0;
 35           a[i][0]+=(i?-1:1)*k*a[x][0];
 36           for (int j=1; j<=n; j++)
 37             a[i][j]-=k*a[x][j];
 38        }
 39 }
 40 bool initial()
 41 {
 42      for (int i=1; i<=n; i++) u[i]=i;
 43      for (int i=1; i<=m; i++) b[i]=n+i;
 44      while (1)
 45      {
 46            int x=0,y=0;
 47            for (int i=1; i<=m; i++)
 48              if (dcmp(a[i][0])<0) x=i; //加break会TLE？
 49            if (!x) return 1;
 50            for (int i=1; i<=n; i++)
 51              if (dcmp(a[x][i])<0) y=i;//加break会TLE？
 52            if (!y) return 0;
 53            pivot(x,y);
 54      }
 55 }
 56
 57 int simplex()
 58 {
 59     if (!initial()) return 0;
 60     while (1)
 61     {
 62           int x=0,y=0;
 63           for (int i=1; i<=n; i++)
 64             if (dcmp(a[0][i])>0) {y=i; break;}
 65           if (!y) return 1;
 66           double mi=1e15;
 67           for (int i=1; i<=m; i++)
 68             if (dcmp(a[i][y])>0&&(!x||a[i][0]/a[i][y]<mi)) {mi=a[i][0]/a[i][y];x=i;}
 69           if (!x) return -1;
 70           pivot(x,y);
 71     }
 72 }
 73 int main()
 74 {
 75     n=read(); m=read(); ty=read();
 76     for (int i=1; i<=n; i++) a[0][i]=read();
 77     for (int i=1; i<=m; i++)
 78     {
 79         for (int j=1; j<=n; j++) a[i][j]=read();
 80         a[i][0]=read();
 81     }
 82     switch (simplex())
 83     {
 84            case 1:{
 85                 printf("%.8lf\n",a[0][0]);
 86                 if (ty)
 87                 {
 88                    for (int i=1; i<=n; i++) u[i]=0;
 89                    for (int i=1; i<=m; i++) if (b[i]<=n) u[b[i]]=i;
 90                    for (int i=1; i<=n; i++) printf("%.8lf ",u[i]?a[u[i]][0]:double(0));
 91                 }
 92                 break;
 93            }
 94            case 0:puts("Infeasible");break;
 95            case -1:puts("Unbounded");break;
 96     }
 97     system("pause");
 98     return 0;
 99 }
100     

我唯一不解的是时间复杂度问题，据说会被卡但实际上基本跑得飞起（下面会用实例证明）

还有模板的initialization上有一处我很不解就是为什么在48行和51行处的循环加break和不加break在时间会带来明显的差距……

（uoj上这两处加上break会TLE……）求神犇指教……

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

比如这个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，

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

 1 #include<iostream>
 2 #include<cstdio>
 3 #include<cstring>
 4 #include<cmath>
 5 #include<stdlib.h>
 6
 7 using namespace std;
 8 const double eps=1e-9;
 9 int b[11100],u[11100],n,m;
10 double a[10010][1010];
11
12 int dcmp(double x)
13 {
14     if (fabs(x)<=eps) return 0;
15     if (x>0) return 1; else return -1;
16 }
17
18 void pivot(int x,int y)
19 {
20      swap(b[x],u[y]);
21      double k=a[x][y];a[x][y]=1;
22      for (int i=0; i<=n; i++) a[x][i]/=k;
23      for (int i=0; i<=m; i++)
24          if (i!=x&&dcmp(a[i][y])!=0)
25          {
26             k=a[i][y]; a[i][y]=0;
27             a[i][0]+=(i?-1:1)*k*a[x][0];
28             for (int j=1; j<=n; j++) a[i][j]-=k*a[x][j];
29          }
30 }
31
32 void simplex()
33 {
34      for (int i=1; i<=n; i++) u[i]=i;
35      for (int i=1; i<=m; i++) b[i]=i+n;
36      while (1)
37      {
38            int x=0,y=0;
39            for (int i=1; i<=n; i++)
40              if (dcmp(a[0][i])>0) {y=i;break;}
41            if (!y) break;
42            double mi=1e20;
43            for (int i=1; i<=m; i++)
44              if (dcmp(a[i][y])>0&&(!x||mi>a[i][0]/a[i][y])) {x=i; mi=a[i][0]/a[i][y];}
45            if (!x) break;
46            pivot(x,y);
47      }
48 }
49
50 int main()
51 {
52     scanf("%d%d",&n,&m);
53     for (int i=1; i<=n; i++) scanf("%lf",&a[0][i]);
54     for (int i=1; i<=m; i++)
55     {
56         int s,t,c;
57         scanf("%d%d%d",&s,&t,&c);
58         for (int j=s; j<=t; j++) a[i][j]=1;
59         a[i][0]=c;
60     }
61     simplex();
62     printf("%d\n",(int)a[0][0]);
63     system("pause");
64     return 0;
65 }

1061