Processing math: 100%

NanoApe's Blog

既是咸鱼又是辣鸡

【规划】单纯形法

NanoApe posted @ 2016年6月10日 18:00 in 蒟蒻不撕烤智熵何来 , 3028 阅读

求解线性规划的单纯形法(这玩意不是都能网络流写的吗(谁跟你这样说的

首先是线性规划的定义

利用松弛变量将约束中的不等式换成等式,则得到

在约束等式左边的 x[n+1]x[n+m] 叫做基变量,等式右边的 x[1]x[n] 叫做非基变量

线性规划有三种情况:有解、无解(不存在一个满足限制条件的解)、无界(目标函数可以无穷大)

显然当 bi 非负时令所有基变量为 bi,非基变量为0,即可得到一个满足条件的初始解

有个重要的操作叫转轴(Pivot)

我们可以把一个基变量 x[b] 和非基变量 x[n] 互换,用 x[b] 和其他非基变量代换这个 x[n],这样 x[b] 就成了非基变量,而 x[n] 成了基变量

为了最大化目标函数,我们考虑不停地找一个系数为正的非基变量,然后增大这个量

最优化的具体过程:

举个栗子:

找到 x1,然后限制一为 x130,限制二为 x112,限制三为 x19,所以我们用 x1 代换下界最紧的 x6

得到下面的新线性规划:

然后找到 x3,得到:

增大 x2 得到

最后得到最大值为28,此时 x1=8, x2=4, x3=0

但是,在 bi 为负的时候,把所有基变量设为0不是一组合法的初始解

我们选择一个 bi 为负的基变量 x[i+n],然后我们在该约束右边找一个系数为正的非基变量,然后进行转轴操作,如果没有系数为正显然就无解了

时间复杂度的话 Pivot 操作的复杂度为 O(nm),但是最优化操作中 Pivot 操作的调用次数可能会成为指数级

其他多项式时间算法:

椭球算法 O(n6m2)

内点算法 O(n3.5m2)

改进的内点算法 O(n3.5m)

 

例题:UOJ #179. 线性规划

Simplex
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
#include <cstdio>
#include <algorithm>
 
#define rep(i, l, r) for(int i=l; i<=r; i++)
     
using namespace std;
 
inline int read()
{
    int x=0; bool f=0; char ch=getchar();
    while (ch<'0' || '9'<ch) f|=ch=='-', ch=getchar();
    while ('0'<=ch && ch<='9') x=x*10+ch-'0', ch=getchar();
    return f?-x:x;
}
 
#define maxn 21
#define eps 1e-8
 
int n, m, op, tot, q[maxn], idx[maxn], idy[maxn]; double a[maxn][maxn], A[maxn];
 
inline void Pivot(int x, int y)
{
    swap(idy[x],idx[y]);
    double tmp=a[x][y]; a[x][y]=1/a[x][y];
    rep(i, 0, n) if (y!=i) a[x][i]/=tmp;
    tot=0; rep(i, 0, n) if (i!=y && (a[x][i]>eps || a[x][i]<-eps)) q[++tot]=i;
    rep(i, 0, m)
    {
        if ((x==i) || (a[i][y]<eps && a[i][y]>-eps)) continue;
        rep(j, 1, tot) a[i][q[j]]-=a[x][q[j]]*a[i][y];
        a[i][y]=-a[i][y]/tmp;
    }
}
 
int main()
{
    n=read(), m=read(), op=read();
    rep(i, 1, n) a[0][i]=read();
    rep(i, 1, m)
    {
        rep(j, 1, n) a[i][j]=read();
        a[i][0]=read();
    }
     
    rep(i, 1, n) idx[i]=i;
    rep(i, 1, m) idy[i]=i+n;
    while (true)
    {
        int x=0, y=0;
        rep(i, 1, m) if (a[i][0]<-eps && ((!x)||(rand()&1))) x=i; if (!x) break;
        rep(i, 1, n) if (a[x][i]<-eps && ((!y)||(rand()&1))) y=i; if (!y) return puts("Infeasible"),0;
        Pivot(x,y);
    }
     
    while (true)
    {
        int x=0, y=0; double mn=1e15;
        rep(i, 1, n) if (a[0][i]>eps) {y=i; break;} if (!y) break;
        rep(i, 1, m) if (a[i][y]>eps && a[i][0]/a[i][y]<mn) mn=a[i][0]/a[i][y], x=i; if (!x) return puts("Unbounded"),0;
        Pivot(x,y);
    }
     
    printf("%.8lf\n", -a[0][0]); if (!op) return 0;
    rep(i, 1, m) if (idy[i]<=n) A[idy[i]]=a[i][0];
    rep(i, 1, n) printf("%.8lf%c", A[i], i==n?'\n':' ');
     
    return 0;
}

 

其余科学阅读资料:

http://www.cnblogs.com/zzqsblog/p/5457091.html


登录 *


loading captcha image...
(输入验证码)
or Ctrl+Enter