单纯形法这方面的资料网上不是很多(也许是我没找到qwq),蒟蒻觉得这个算法也不是特别好理解,现在也只是敢说大致了解了它的原理会套一下板子做题(证明过程什么的完全不会),感觉也没有了解到这个算法的本质,总之还是学得很浅显。决定先记录下来,以后再回来深入一点地学习,希望不要鸽了(qwq)。
首先是推荐学习博客:
原理:
代码及习题:
总结:
单纯形法是用来解决线性规划问题的一个算法,简单总结起来就是先找到一个可行解,然后通过不断旋转(pivot)找到更优解,直到找不到更优解时结束。
如果是囫囵吞枣版做题的话就是:根据题目求出线性方程--->看是否需要对偶性(单纯形的标准型是目标函数最大化,题目若是目标最小的话就要对偶)--->套板子 。
这里要注意理解对偶性(借用Candy大佬的原话):
其实就是常数和目标函数互换,把系数矩阵转置。
UOJ-179 线性规划
线性规划输出方案的模板题,抄的大佬代码当板子用。qwq
#include#include #include #include using namespace std;const int MAXN = 51, INF = 1e9 + 10;const double eps = 1e-8;int N, M, opt;int id[MAXN];double ans[MAXN], a[MAXN][MAXN];void Pivot(int l, int e) { swap(id[N + l], id[e]); double t = a[l][e]; a[l][e] = 1; for(int i = 0; i <= N; i++) a[l][i] /= t;//交换基变量与非基变量 for(int i = 0; i <= M; i++) { if(i != l && abs(a[i][e]) > eps) { t = a[i][e]; a[i][e] = 0;//t表示系数 for(int j = 0; j <= N; j++) a[i][j] -= a[l][j] * t;//消去需要消去的非基变量 } }//带入的过程就是消去非基变量 }bool init() { //寻找初始化解,若bi < 0,找到一个ai < 0,转换它们,不断重复直到所有的bi都 > 0 //此时x1 = 0, x2 = 0···就是一组解 while(1) { int l = 0, e = 0; for(int i = 1; i <= M; i++) if(a[i][0] < -eps && (!l || (rand() & 1))) l = i; if(!l) break; for(int i = 1; i <= N; i++) if(a[l][i] < -eps && (!e || (rand() & 1))) e = i; if(!e) {puts("Infeasible"); return 0;} //这里所有的Xi都是正的,而bi是负的,这与松弛型相矛盾 Pivot(l, e); } return 1;}bool simplex() { while(1) { int l = 0, e = 0; double mn = INF; for(int i = 1; i <= N; i++) if(a[0][i] > eps) {e = i; break;} if(!e) break; for(int i = 1; i <= M; i++) if(a[i][e] > eps && a[i][0] / a[i][e] < mn) mn = a[i][0] / a[i][e], l = i;//找到下界最紧的松弛 if(!l) {puts("Unbounded"); return 0;} Pivot(l, e); } return 1;}int main() { srand(19260817); scanf("%d",&N); scanf("%d",&M); scanf("%d",&opt); for(int i = 1; i <= N; i++) scanf("%lf",&a[0][i]);//最大化C1X1 + C2X2··· for(int i = 1; i <= M; i++) { for(int j = 1; j <= N; j++) scanf("%lf",&a[i][j]); scanf("%lf",&a[i][0]);// <= b } for(int i = 1; i <= N; i++) id[i] = i; if(init() && simplex()) { printf("%.8lf\n", -a[0][0]); if(opt) { for(int i = 1; i <= M; i++) ans[id[i + N]] = a[i][0];//tag for(int i = 1; i <= N; i++) printf("%.8lf ", ans[i]); } } return 0;}
BZOJ-1061 志愿者招募
列出方程式后,用对偶性转换,套板子求解。依旧抄的大佬的代码,当作简易版本的线性规划板子用。qwq
#includeusing namespace std;typedef long long ll;const int M=10005,N=1005,INF=1e9;const double eps=1e-6;int n,m; //变量个数,方程个数 double a[M][N],b[M],c[N],v; //对偶后:a[i]系数矩阵,b[i]常数,c[i]目标 void pivot(int l,int e){ b[l]/=a[l][e]; for(int j=1;j<=n;j++) if(j!=e) a[l][j]/=a[l][e]; a[l][e]=1/a[l][e]; for(int i=1;i<=m;i++) if(i!=l&&fabs(a[i][e])>0){ b[i]-=a[i][e]*b[l]; for(int j=1;j<=n;j++) if(j!=e) a[i][j]-=a[i][e]*a[l][j]; a[i][e]=-a[i][e]*a[l][e]; } v+=c[e]*b[l]; for(int j=1;j<=n;j++) if(j!=e) c[j]-=c[e]*a[l][j]; c[e]=-c[e]*a[l][e]; //swap(B[l],N[e])}double simplex(){ while(true){ int e=0,l=0; for(e=1;e<=n;e++) if(c[e]>eps) break; if(e==n+1) return v; //所有c都<=0,得到最优解 double mn=INF; for(int i=1;i<=m;i++) if(a[i][e]>eps&&mn>b[i]/a[i][e]) mn=b[i]/a[i][e],l=i; if(mn==INF) return INF; //unbounded pivot(l,e); }}int main(){ scanf("%d%d",&n,&m); for(int i=1;i<=n;i++) scanf("%lf",&c[i]); for(int i=1;i<=m;i++){ int s,t; scanf("%d%d",&s,&t); for(int j=s;j<=t;j++) a[i][j]=1; scanf("%lf",&b[i]); } printf("%d",(int)(simplex()+0.5));}
题目练习:
BZOJ-3550
简单模板题,每个区间限制列一条方程,每个数再列一条方程,套板子求解即可。
#includeusing namespace std;typedef long long ll;const int N=1e3+10,INF=0x3f3f3f3f;const double eps=1e-6;int n,m; //变量个数,方程个数 double a[N][N],b[N],c[N],v; //对偶后:a[i]系数矩阵,b[i]常数,c[i]目标 void pivot(int l,int e){ b[l]/=a[l][e]; for(int j=1;j<=n;j++) if(j!=e) a[l][j]/=a[l][e]; a[l][e]=1/a[l][e]; for(int i=1;i<=m;i++) if(i!=l&&fabs(a[i][e])>0){ b[i]-=a[i][e]*b[l]; for(int j=1;j<=n;j++) if(j!=e) a[i][j]-=a[i][e]*a[l][j]; a[i][e]=-a[i][e]*a[l][e]; } v+=c[e]*b[l]; for(int j=1;j<=n;j++) if(j!=e) c[j]-=c[e]*a[l][j]; c[e]=-c[e]*a[l][e]; //swap(B[l],N[e])}double simplex(){ while(true){ int e=0,l=0; for(e=1;e<=n;e++) if(c[e]>eps) break; if(e==n+1) return v; //所有c都<=0,得到最优解 double mn=INF; for(int i=1;i<=m;i++) if(a[i][e]>eps&&mn>b[i]/a[i][e]) mn=b[i]/a[i][e],l=i; if(mn==INF) return INF; //unbounded pivot(l,e); }}int main(){ int t1,t2; scanf("%d%d",&t1,&t2); n=3*t1; m=n-t1+1; for (int i=1;i<=n;i++) scanf("%lf",&c[i]); for (int i=1;i<=m;i++) { for (int j=i;j
BZOJ-3112
这题题目要求最小值,所以列出方程用对偶性转化,套个板子就OK了。
#includeusing namespace std;typedef long long ll;const int M=1005,N=10005,INF=1e9;const double eps=1e-6;int n,m; //变量个数,方程个数 double a[M][N],b[M],c[N],v; //对偶后:a[i]系数矩阵,b[i]常数,c[i]目标 void pivot(int l,int e){ b[l]/=a[l][e]; for(int j=1;j<=n;j++) if(j!=e) a[l][j]/=a[l][e]; a[l][e]=1/a[l][e]; for(int i=1;i<=m;i++) if(i!=l&&fabs(a[i][e])>0){ b[i]-=a[i][e]*b[l]; for(int j=1;j<=n;j++) if(j!=e) a[i][j]-=a[i][e]*a[l][j]; a[i][e]=-a[i][e]*a[l][e]; } v+=c[e]*b[l]; for(int j=1;j<=n;j++) if(j!=e) c[j]-=c[e]*a[l][j]; c[e]=-c[e]*a[l][e]; //swap(B[l],N[e])}double simplex(){ while(true){ int e=0,l=0; for(e=1;e<=n;e++) if(c[e]>eps) break; if(e==n+1) return v; //所有c都<=0,得到最优解 double mn=INF; for(int i=1;i<=m;i++) if(a[i][e]>eps&&mn>b[i]/a[i][e]) mn=b[i]/a[i][e],l=i; if(mn==INF) return INF; //unbounded pivot(l,e); }}int main(){ scanf("%d%d",&m,&n); for(int i=1;i<=m;i++) scanf("%lf",&b[i]); for(int j=1;j<=n;j++){ int s,t; scanf("%d%d",&s,&t); for(int i=s;i<=t;i++) a[i][j]=1; scanf("%lf",&c[j]); } printf("%d",(int)(simplex()+0.5));}
BZOJ 3265
这题是BZOJ1061的增强版,其实也没差多少。但是BZOJ该题评论区有大佬说单纯形法是错误的(蒟蒻一脸懵逼)。
#includeusing namespace std;typedef long long ll;const int M=10005,N=1005,INF=1e9;const double eps=1e-6;int n,m; //变量个数,方程个数 double a[M][N],b[M],c[N],v; //对偶后:a[i]系数矩阵,b[i]常数,c[i]目标 void pivot(int l,int e){ b[l]/=a[l][e]; for(int j=1;j<=n;j++) if(j!=e) a[l][j]/=a[l][e]; a[l][e]=1/a[l][e]; for(int i=1;i<=m;i++) if(i!=l&&fabs(a[i][e])>0){ b[i]-=a[i][e]*b[l]; for(int j=1;j<=n;j++) if(j!=e) a[i][j]-=a[i][e]*a[l][j]; a[i][e]=-a[i][e]*a[l][e]; } v+=c[e]*b[l]; for(int j=1;j<=n;j++) if(j!=e) c[j]-=c[e]*a[l][j]; c[e]=-c[e]*a[l][e]; //swap(B[l],N[e])}double simplex(){ while(true){ int e=0,l=0; for(e=1;e<=n;e++) if(c[e]>eps) break; if(e==n+1) return v; //所有c都<=0,得到最优解 double mn=INF; for(int i=1;i<=m;i++) if(a[i][e]>eps&&mn>b[i]/a[i][e]) mn=b[i]/a[i][e],l=i; if(mn==INF) return INF; //unbounded pivot(l,e); }}int main(){ scanf("%d%d",&n,&m); for(int i=1;i<=n;i++) scanf("%lf",&c[i]); for(int i=1;i<=m;i++){ int t; scanf("%d",&t); while (t--) { int s,t; scanf("%d%d",&s,&t); for(int j=s;j<=t;j++) a[i][j]=1; } scanf("%lf",&b[i]); } printf("%d",(int)(simplex()+0.5));}
POJ 1275
这题跟BZOJ1061很像,线性方程应该不难想。注意一些细节:小时会循环,代价就是人数,矩阵记得清零。还有这个板子注意v也要清零。
#include#include #include #include #include using namespace std;typedef long long ll;const int M=1005,N=30,INF=1e9;const double eps=1e-6;int n,m,sum[30]; //变量个数,方程个数 double a[M][N],b[M],c[N],v; //对偶后:a[i]系数矩阵,b[i]常数,c[i]目标 void pivot(int l,int e){ b[l]/=a[l][e]; for(int j=1;j<=n;j++) if(j!=e) a[l][j]/=a[l][e]; a[l][e]=1/a[l][e]; for(int i=1;i<=m;i++) if(i!=l&&fabs(a[i][e])>0){ b[i]-=a[i][e]*b[l]; for(int j=1;j<=n;j++) if(j!=e) a[i][j]-=a[i][e]*a[l][j]; a[i][e]=-a[i][e]*a[l][e]; } v+=c[e]*b[l]; for(int j=1;j<=n;j++) if(j!=e) c[j]-=c[e]*a[l][j]; c[e]=-c[e]*a[l][e]; //swap(B[l],N[e])}double simplex(){ while(true){ int e=0,l=0; for(e=1;e<=n;e++) if(c[e]>eps) break; if(e==n+1) return v; //所有c都<=0,得到最优解 double mn=INF; for(int i=1;i<=m;i++) if(a[i][e]>eps&&mn>b[i]/a[i][e]) mn=b[i]/a[i][e],l=i; if(mn==INF) return INF; //unbounded pivot(l,e); }}int main(){ int T; cin>>T; while (T--) { memset(a,0,sizeof(a)); memset(b,0,sizeof(b)); memset(c,0,sizeof(c)); memset(sum,0,sizeof(sum)); n=24; for (int i=1;i<=24;i++) scanf("%lf",&c[i]); scanf("%d",&m); for (int i=1;i<=m;i++) { int s; scanf("%d",&s); s++; for (int j=0;j<8;j++) { int t=(s+j-1)%24+1; a[i][t]=1; sum[t]++; } b[i]=1; } bool ok=0; for (int i=1;i<=24;i++) if (c[i]>sum[i]) { ok=1; break; } if (ok) { puts("No Solution"); continue; } v=0; printf("%d\n",(int)(simplex()+0.5)); }}