【BZOJ3265】志愿者招募加强版 线性规划 单纯形法 对偶原理

链接：

#include <stdio.h>int main(){    puts("转载请注明出处[vmurder]谢谢");    puts("网址：blog.csdn.net/vmurder/article/details/45673079");}

题解：

这道题是线性规划求目标函数最小值 ，对偶原理转一下就成了单纯形算法求线性规划最大值。

单纯形法：

首先这篇博客缺失了很多证明，只能讲述单纯形法的实现。

 // N个变量 M个限制double a[M][N]; // 这m个限制里变量的系数double b[M]; // 第i个限制加和 <=bidouble c[N]; // 目标函数里变量系数double ans; // 目标函数自带的系数void pivot(int x,int y); // 转轴操作void solve();

首先对于如下数据

3 3 2 3 41 1 2 21 2 3 51 3 3 2

我们可以得到以下的限制：
目标函数: 2x1+5x2+2x3+0
限制函数:
x1≤2
x1+x2≤3
x2+x3≤4

然后若 x4x5x6 均为任意值。
那么我们可以将式子表示为
状态 1: **************
目标函数: 2x1+5x2+2x3+0
限制函数:
x4+1x1=2
x5+1x1+1x2=3
x6+1x2+1x3=4

然后在单纯形的过程中函数们的各系数得到改变的过程如下：

状态 1: **************
目标函数: 2x1+5x2+2x3+0
限制函数:
x4+1x1=2
x5+1x1+1x2=3
x6+1x2+1x3=4
状态 2: **************
目标函数: −2x4+5x2+2x3+4
限制函数:
x1+1x4=2
x5+−1x4+1x2=1
x6+1x2+1x3=4
状态 3: **************
目标函数: 3x4+−5x5+2x3+9
限制函数:
x1+1x4=2
x2+−1x4+1x5=1
x6+1x4+−1x5+1x3=3
状态 4: **************
目标函数: −3x1+−5x5+2x3+15
限制函数:
x4+1x1=2
x2+1x1+1x5=3
x6+−1x1+−1x5+1x3=1
状态 5: **************
目标函数: −1x1+−3x5+−2x6+17
限制函数:
x4+1x1=2
x2+1x1+1x5=3
x3+−1x1+−1x5+1x6=1

线性规划最终答案：17

单纯形思想：

其实就是每次找一个在目标函数中存在的变量，跟不在其中的某个变量交换一下，然后换来换去最终目标函数中所有变量系数都为负了，那么目标函数最后加的常数就是答案辣。

单纯形实现：

具体的实现过程呢？

struct Simplex{    int n,m; // n个变量、m个限制    double a[M][N],b[M],c[M],ans;    void pivot(int x,int y)    {        int i,j;        double t;        b[x]/=a[x][y];        for(i=1;i<=n;i++)if(i!=y)a[x][i]/=a[x][y];        a[x][y]=1.0/a[x][y];        for(i=1;i<=m;i++)if(i!=x&&fabs(a[i][y])>eps)        {            b[i]-=a[i][y]*b[x],t=a[i][y];            for(j=1;j<=n;j++)a[i][j]-=t*a[x][j];            a[i][y]=-t*a[x][y];        }        ans+=c[y]*b[x],t=c[y];        for(i=1;i<=n;i++)c[i]-=t*a[x][i];        c[y]=-t*a[x][y];    }    double solve()    {        read();        double t;        for(int i,x,y;;)        {            for(i=1;i<=m;i++)if(c[i]>eps){y=i;break;}            if(i>m)return ans;            for(t=inf,i=1;i<=m;i++)                if(a[i][y]>eps&&t>b[i]/a[i][y])                    x=i,t=b[i]/a[i][y];            if(t==inf)return t;            else pivot(x,y);        }    }}simplex;

下面的实现里存在两个坑，请不要深究它们怎么证。

double solve() :

首先我们每次找一个标号最靠前的【坑1】目标函数中系数为正的变量 y (如第30行)，然后如果找不到，就表示已经找到了最优解【坑2】。

然后我们找这个变量在哪个限制里是”最紧”的【坑3】，即33、34、35这三行。
如果找不到任何一个”紧”的行，即 ∀a(i,y)<=eps 那么这个变量想多大就能多大，反正有辅助变量顶着，可以直接返回无界【坑4】，即答案是无穷大。

然后我们对此限制中的辅助变量和找到的 y 进行转轴操作。

void pivot(int x,int y) :

先把第 x 个限制的系数改一改，然后用它去消其它的限制，再消一下目标函数。结束。类似高斯消元？

对偶原理：

如果要求最小值，那么我们把模型对偶一下就可以辣。【坑5】

即限制矩阵 【转置】 一下，然后目标函数的系数和每个限制的最终 ≤ 的那个值也互换一下。

【转置】 ：矩阵的元素 a(i,j) 变为 a(j,i) 。

对于证明坑的总结：

【坑1】为什么要找标号最小的？
根据Bland法则，这么做可以避免被卡死循环。
【坑2】为什么标号都是负的就表示找到了最优解？
表示无法理解为什么不能给某个存在负系数的变量转一转。
【坑3】为什么要找”最紧”的？
不知道不知道。
【坑4】为什么直接返回无界了？
不可以别的某个变量再来转一转，然后就回归有界？
【坑5】对偶原理是毛线？
天然坑，深坑。

代码：

#include <cmath>#include <cstdio>#include <cstring>#include <iostream>#include <algorithm>#define N 1010 // 变量个数#define M 10100 // 限制个数#define inf 1e9#define eps 1e-5using namespace std;struct Simplex{    int n,m; // n个变量、m个限制    double a[M][N],b[M],c[M],ans;    void read()    {        int i,j,k,l,r;        scanf("%d%d",&m,&n);        for(i=1;i<=m;i++)scanf("%lf",&c[i]);        for(i=1;i<=n;i++)        {            scanf("%d",&k);            while(k--)            {                scanf("%d%d",&l,&r);                for(j=l;j<=r;j++)a[i][j]=1.0;            }            scanf("%lf",&b[i]);        }        swap(n,m);    }    void pivot(int x,int y)    {        int i,j;        double t;        b[x]/=a[x][y];        for(i=1;i<=n;i++)if(i!=y)a[x][i]/=a[x][y];        a[x][y]=1.0/a[x][y];        for(i=1;i<=m;i++)if(i!=x&&fabs(a[i][y])>eps)        {            b[i]-=a[i][y]*b[x],t=a[i][y];            for(j=1;j<=n;j++)a[i][j]-=t*a[x][j];            a[i][y]=-t*a[x][y];        }        ans+=c[y]*b[x],t=c[y];        for(i=1;i<=n;i++)c[i]-=t*a[x][i];        c[y]=-t*a[x][y];    }    double solve()    {        read();        double t;        for(int i,x,y;;)        {            for(i=1;i<=m;i++)if(c[i]>eps){y=i;break;}            if(i>m)return ans;            for(t=inf,i=1;i<=m;i++)                if(a[i][y]>eps&&t>b[i]/a[i][y])                    x=i,t=b[i]/a[i][y];            if(t==inf)return t;            else pivot(x,y);        }    }}simplex;int main(){//  freopen("test.in","r",stdin);    printf("%d\n",(int)(simplex.solve()+eps));    return 0;}