线性规划与单纯形法

标准型 最大化 ∑j=1ncjxj 满足约束 ∑j=1naijxj≤bi i=1,2,...,m ​ xj≥0 j=1,2,...,m 松弛型 最大化 ∑j=1ncjxj 满足约束 xi+n=bi−∑j=1naijxj ​ xj≥0 j=1,2,...,n+m 线性规划的解空间是一个凸性区域,局部最优解就是全局最优解。 单纯形 基变量:在松弛型等式左侧的所有变量。 非基变量:在松弛型等式右侧的所有变量。 如果常数项都非负,我们可以拿出一个基本解,即所有基变量的值为右侧常数项,非基变量的值为0。 转轴(pivot)操作 选择一个基变量xB和一个非基变量xN,将其互换(称这个非基变量为换入变量,基变量为换出变量)。 xB=bi−∑j=1naijxj xN=(bi−∑j≠Naijxj−xB)/aiN 将目标函数以及其余所有限制中xN用这个等式替代。 simplex 操作 从基本解出发,每次找一个目标函数系数为正的非基变量和一个对他约束最强的基变量进行转轴操作。 保证每次转轴操作都会使得目标函数增大。 initialization操作 初始化,即找到一组基本解。 做法:引入一个辅助线性规划: 最大化 −x0 满足约束 xi+n=bi−∑j=1naijxj+x0 ​ xj≥0 如果原线性规划存在一个解为(x1,...,xn+m),那么(0,x1,..,xn+m)为辅助线性规划的一个可行解。 我们最大化−x0那么最优解就是这个可行解。 构造辅助线性规划的初始解,我们把x0当做换入变量,找到最小的bl,把xl+n当做换出变量执行一次转轴操作。 x0=−bl+∑j=1naljxj+xl+n 对于其他的约束,变为 xi+n=−bl+bi+∑j=1(alj−aij)xj+xl+n 操作之后辅助线性规划有一个可行的初始基本解了,执行simplex操作即可。 时间复杂度 O(simplex)×O(nm)=O(???)=O(可过) 可以证明单纯形法一定能够在有限步之后终止,但是最坏情况算法的时间复杂度为指数级别的。 对偶问题 给定一个原始线性规划: 最小化 cTx 满足约束 Ax≥b x≥0 其对偶线性规划为: 最大化 bTy 满足约束 ATy≤c y≥0 线性规划弱对偶性:若x=(x1,...xn)是原问题的一个可行解,y=(y1,...ym)是对偶问题的可行解。那么: ∑j=1ncjxj≥∑i=1mbiyi 线性规划对偶性:若x=(x1,...xn)是原问题的一个最优解,y=(y1,...ym)是对偶问题的最优解。那么: ∑j=1ncjxj=∑i=1mbiyi 根据这个性质可以将一系列问题的初始化操作去掉,转化成对偶后的线性规划进行求解。 网络流解线性规划问题 https://www.cnblogs.com/suika/p/9079315.html 对于一些线性规划问题,我们通常能够转化成 每个变量的都出现两次,且系数分别为+1和−1。 就是这样的模型,我们可以用网络流的方法巧妙的解决。 首先网络流有性质:对于除了源点和汇点的其他点,有 流入的总量等于流出的总量。 这让我们有一个思路,把点当成每一个方程,把边看成给予限制的变量,对于每个变量,把+1看成流出−1看成流入,找到对应方程的位置连边。 常见的问题转化成线性规划 最短路: 最大化 dt 满足约束 ds=0 dvi≤dui+li 最大流 设fu,v表示u到v的实际流量。 最大化 ∑ufs,u 满足约束 fu,v=fv,u ​ ∑i≠s,tfu,i=0 全幺模矩阵 全幺模矩阵是指任何一个行数和列数相等的满秩子矩阵行列式的值都是1或−1的矩阵。 对于一个线性规划中的A为全幺模矩阵,那么所有约束系数都是1,0,−1中的一个。 也就是说,直接执行线性规划就能得到相同的整数规划问题的正确答案。 模板 http://uoj.ac/problem/179 #include #include #include using namespace std; #define N 105 #define eps 1e-8 typedef double f2; f2 a[N][N]; int n,m,type,id[N],tp[N]; void pivot(int r,int c) { swap(id[r+n],id[c]); f2 tmp=-a[r][c]; a[r][c]=-1; int i,j; for(i=0;i<=n;i++) a[r][i]/=tmp; for(i=0;i<=m;i++) if(a[i][c]&&i!=r) { tmp=a[i][c]; a[i][c]=0; for(j=0;j<=n;j++) a[i][j]+=tmp*a[r][j]; } } void solve() { int i,j,k; for(i=1;i<=n;i++) id[i]=i; while(1) { i=j=0; f2 tmp=-eps; for(k=1;k<=m;k++) if(a[k][0]eps) {j=k; break;} if(!j) {puts("Infeasible"); return ;} pivot(i,j); } while(1) { i=j=0; f2 tmp=eps; for(k=1;k<=n;k++) if(a[0][k]>tmp) tmp=a[0][j=k]; if(!j) break; tmp=1e9; for(k=1;k<=m;k++) if(a[k][j]<-eps&&-a[k][0]/a[k][j] #include #include using namespace std; #define N 10050 #define M 100050 #define S (10049) #define T (10048) #define inf (1<<30) int head[N],to[M],nxt[M],flow[M],val[M],cnt=1,path[N],Q[N],l,r,dis[N],inq[N],n,m,a[N],c[M],s[N],t[N]; inline void add(int u,int v,int f,int w) { to[++cnt]=v; nxt[cnt]=head[u]; head[u]=cnt; flow[cnt]=f; val[cnt]=w; to[++cnt]=u; nxt[cnt]=head[v]; head[v]=cnt; flow[cnt]=0; val[cnt]=-w; } bool spfa() { memset(dis,0x3f,sizeof(dis)); memset(path,0,sizeof(path)); l=r=0; Q[r++]=S; dis[S]=0; inq[S]=1; while(l!=r) { int x=Q[l++],i; if(l==n+10) l=0; inq[x]=0; for(i=head[x];i;i=nxt[i]) { if(flow[i]&&dis[to[i]]>dis[x]+val[i]) { dis[to[i]]=dis[x]+val[i]; path[to[i]]=i^1; if(!inq[to[i]]) { Q[r++]=to[i]; if(r==n+10) r=0; inq[to[i]]=1; } } } } return path[T]; } void mcmf() { int minc=0,maxf=0; while(spfa()) { int nf=1<<30,i; for(i=T;i!=S;i=to[path[i]]) { nf=min(nf,flow[path[i]^1]); } for(i=T;i!=S;i=to[path[i]]) { flow[path[i]]+=nf; flow[path[i]^1]-=nf; minc+=nf*val[path[i]^1]; } maxf+=nf; } printf("%d\n",minc); } int main() { scanf("%d%d",&n,&m); int i; for(i=1;i<=n;i++) { scanf("%d",&c[i]); add(i+1,i,inf,0); } for(i=1;i<=m;i++) { scanf("%d%d%d",&s[i],&t[i],&a[i]); add(s[i],t[i]+1,inf,a[i]); } for(i=1;i<=n+1;i++) { if(c[i]-c[i-1]>0) { add(S,i,c[i]-c[i-1],0); }else { add(i,T,c[i-1]-c[i],0); } } mcmf(); } BZOJ3550[ONTAK2010]Vacation&&BZOJ1283:序列 https://lydsy.com/JudgeOnline/problem.php?id=3550 题解:https://www.cnblogs.com/suika/p/9080922.html 代码: #include #include #include using namespace std; #define N 1050 #define M 100050 #define S (n-m+3) #define T (n-m+4) #define inf (1<<30) int head[N],to[M],nxt[M],cnt=1,path[N],dis[N],Q[N],l,r,inq[N],flow[M],val[M]; int n,m,K,a[N]; inline void add(int u,int v,int f,int w) { to[++cnt]=v; nxt[cnt]=head[u]; head[u]=cnt; flow[cnt]=f; val[cnt]=w; to[++cnt]=u; nxt[cnt]=head[v]; head[v]=cnt; flow[cnt]=0; val[cnt]=-w; } bool spfa() { memset(dis,0x80,sizeof(dis)); memset(path,0,sizeof(path)); l=r=0; Q[r++]=S; dis[S]=0; inq[S]=1; while(l!=r) { int x=Q[l++],i; if(l==T+1) l=0; inq[x]=0; for(i=head[x];i;i=nxt[i]) { if(flow[i]&&dis[to[i]]n-m) add(i-m+1,n-m+2,1,a[i]); else add(i-m+1,i+1,1,a[i]); } mcmf(); } BZOJ1937: [Shoi2004]Mst 最小生成树 https://lydsy.com/JudgeOnline/problem.php?id=1937 同BZOJ3118。 空间开成[9900][770]可过。 需要判负0的情况。 代码: #include #include #include #include using namespace std; typedef double f2; #define N 55 #define M 770 //maxm=766 #define eps 1e-8 int head[N],to[N<<1],nxt[N<<1],cnt,n,m,fa[N],val[N<<1]; int idx[N],dep[N],tot; f2 a[9950][M]; struct A { int x,y,k,z,ist; void read() { scanf("%d%d%d",&x,&y,&z); } }e[M]; inline void add(int u,int v,int w) { to[++cnt]=v; nxt[cnt]=head[u]; head[u]=cnt; val[cnt]=w; } void dfs(int x,int y) { int i; fa[x]=y; dep[x]=dep[y]+1; for(i=head[x];i;i=nxt[i]) if(to[i]!=y) { dfs(to[i],x); idx[to[i]]=val[i]; } } int lca(int x,int y) { while(x!=y) { if(dep[x]eps) {j=k; break;} // if(!j) {puts("NO"); return ;} pivot(i,j,n,m); } while(1) { i=j=0; t=eps; for(k=1;k<=n;k++) if(a[0][k]>t) t=a[0][j=k]; if(!j) break; t=1e9; for(k=1;k<=m;k++) if(a[k][j]<-eps&&-a[k][0]/a[k][j]e[i].y) swap(e[i].x,e[i].y); ID[e[i].x][e[i].y]=i; } int x,y; for(i=1;iy) swap(x,y); e[ID[x][y]].ist=1; } for(i=1;i<=m;i++) { if(e[i].ist) { add(e[i].x, e[i].y, i); add(e[i].y, e[i].x, i); } } dfs(1,0); for(i=1;i<=m;i++) { if(e[i].ist==0) { build(i); a[0][i]=-1; }else { a[0][i]=-1; } } simplex(m,tot); } https://www.codechef.com/FEB12/problems/FLYDIST http://poj.org/problem?id=3689 https://www.codechef.com/JUNE15/problems/CHEFBOOK BZOJ3118: Orz the MST https://lydsy.com/JudgeOnline/problem.php?id=3118 分析: 考虑一条非树边会约束哪些树边。 然后就可以做了,注意求最小要把c都取负数求最大。 代码: #include #include #include using namespace std; typedef double f2; #define N 305 #define M 1005 #define eps 1e-8 int head[N],to[N<<1],nxt[N<<1],cnt,n,m,fa[N],val[N<<1]; int idx[N],dep[N],tot; f2 a[4050][M]; struct A { int x,y,k,z,a,b,ist; void read() { scanf("%d%d%d%d%d%d",&x,&y,&z,&ist,&a,&b); } }e[M]; inline void add(int u,int v,int w) { to[++cnt]=v; nxt[cnt]=head[u]; head[u]=cnt; val[cnt]=w; } void dfs(int x,int y) { int i; fa[x]=y; dep[x]=dep[y]+1; for(i=head[x];i;i=nxt[i]) if(to[i]!=y) { dfs(to[i],x); idx[to[i]]=val[i]; } } int lca(int x,int y) { while(x!=y) { if(dep[x]eps) {j=k; break;} // if(!j) {puts("NO"); return ;} pivot(i,j,n,m); } while(1) { i=j=0; t=eps; for(k=1;k<=n;k++) if(a[0][k]>t) t=a[0][j=k]; if(!j) break; t=1e9; for(k=1;k<=m;k++) if(a[k][j]<-eps&&-a[k][0]/a[k][j] #include #include #include using namespace std; typedef double f2; typedef long double f3; #define N 1050 #define M 10050 #define eps 1e-9 #define inf 1e10 f2 a[M][N]; int n,m,x; void pivot(int r,int c,int n,int m) { f2 t=-a[r][c]; a[r][c]=-1; int i,j; for(i=0;i<=n;i++) a[r][i]/=t; for(i=0;i<=m;i++) if(abs(a[i][c])>eps&&i!=r) { t=a[i][c]; a[i][c]=0; for(j=0;j<=n;j++) a[i][j]+=t*a[r][j]; } } void simplex(int n,int m) { int i,j,k; f2 t; while(1) { i=j=0; t=-eps; for(k=1;k<=m;k++) if(a[k][0]eps) {j
50000+
5万行代码练就真实本领
17年
创办于2008年老牌培训机构
1000+
合作企业
98%
就业率

联系我们

电话咨询

0532-85025005

扫码添加微信