线性规划与单纯形法学习笔记

单纯形法这方面的资料网上不是很多（也许是我没找到qwq），蒟蒻觉得这个算法也不是特别好理解，现在也只是敢说大致了解了它的原理会套一下板子做题（证明过程什么的完全不会），感觉也没有了解到这个算法的本质，总之还是学得很浅显。决定先记录下来，以后再回来深入一点地学习，希望不要鸽了（qwq）。

首先是推荐学习博客：

原理：https://www.hrwhisper.me/introduction-to-simplex-algorithm/

代码及习题：http://www.cnblogs.com/ECJTUACM-873284962/p/7097864.html#autoid-1-1-0

总结：https://www.cnblogs.com/candy99/p/lp.html 

单纯形法是用来解决线性规划问题的一个算法，简单总结起来就是先找到一个可行解，然后通过不断旋转(pivot)找到更优解，直到找不到更优解时结束。

如果是囫囵吞枣版做题的话就是：根据题目求出线性方程--->看是否需要对偶性(单纯形的标准型是目标函数最大化，题目若是目标最小的话就要对偶)--->套板子  。

这里要注意理解对偶性(借用Candy大佬的原话)：

 

其实就是常数和目标函数互换，把系数矩阵转置。

UOJ-179 线性规划

线性规划输出方案的模板题，抄的大佬代码当板子用。qwq

#include<iostream>
#include<cstdio>
#include<algorithm>
#include<cmath>
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
", -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

#include<bits/stdc++.h>
using 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

简单模板题，每个区间限制列一条方程，每个数再列一条方程，套板子求解即可。

#include<bits/stdc++.h>
using 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<i+t1;j++) a[i][j]=1;
        b[i]=t2;
    }
    for (int i=1;i<=n;i++) a[++m][i]=1,b[m]=1;  //限制每个数只能选一次 
    
    printf("%d",(int)(simplex()+0.5));
}

BZOJ-3112

这题题目要求最小值，所以列出方程用对偶性转化，套个板子就OK了。

#include<bits/stdc++.h>
using 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该题评论区有大佬说单纯形法是错误的（蒟蒻一脸懵逼）。

#include<bits/stdc++.h>
using 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<iostream>
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cmath>
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
",(int)(simplex()+0.5));
    }
}