【规划】单纯形法
求解线性规划的单纯形法(这玩意不是都能网络流写的吗(谁跟你这样说的
首先是线性规划的定义
利用松弛变量将约束中的不等式换成等式,则得到
在约束等式左边的 $x[n+1]…x[n+m]$ 叫做基变量,等式右边的 $x[1]…x[n]$ 叫做非基变量
线性规划有三种情况:有解、无解(不存在一个满足限制条件的解)、无界(目标函数可以无穷大)
显然当 $b_i$ 非负时令所有基变量为 $b_i$,非基变量为0,即可得到一个满足条件的初始解
有个重要的操作叫转轴(Pivot)
我们可以把一个基变量 $x[b]$ 和非基变量 $x[n]$ 互换,用 $x[b]$ 和其他非基变量代换这个 $x[n]$,这样 $x[b]$ 就成了非基变量,而 $x[n]$ 成了基变量
为了最大化目标函数,我们考虑不停地找一个系数为正的非基变量,然后增大这个量
最优化的具体过程:
举个栗子:
找到 $x1$,然后限制一为 $x1 \le 30$,限制二为 $x1 \le 12$,限制三为 $x1 \le 9$,所以我们用 $x1$ 代换下界最紧的 $x6$
得到下面的新线性规划:
然后找到 $x3$,得到:
增大 $x2$ 得到
最后得到最大值为28,此时 $x1=8,~x2=4,~x3=0$
但是,在 $bi$ 为负的时候,把所有基变量设为0不是一组合法的初始解
我们选择一个 $bi$ 为负的基变量 $x[i+n]$,然后我们在该约束右边找一个系数为正的非基变量,然后进行转轴操作,如果没有系数为正显然就无解了
时间复杂度的话 Pivot 操作的复杂度为 $O(nm)$,但是最优化操作中 Pivot 操作的调用次数可能会成为指数级
其他多项式时间算法:
椭球算法 $O(n^6 m^2)$
内点算法 $O(n^{3.5}m^2)$
改进的内点算法 $O(n^{3.5}m)$
例题:UOJ #179. 线性规划
#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; }
其余科学阅读资料: