【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个限制加和 <=bi
double c[N]; // 目标函数里变量系数
double ans; // 目标函数自带的系数
void pivot(int x,int y); // 转轴操作
void solve();

首先对于如下数据

3 3 

2 3 4

1 1 2 2

1 2 3 5

1 3 3 2

我们可以得到以下的限制：
目标函数: $2x_1+5x_2+2x_3+0$
限制函数:
$x_1\leq2$
$x_1+x_2\leq3$
$x_2+x_3\leq4$

然后若 $x_4 x_5 x_6$ 均为任意值。
那么我们可以将式子表示为
状态 1: **************
目标函数: $2x_1+5x_2+2x_3+0$
限制函数:
$x_4+1x_1=2$
$x_5+1x_1+1x_2=3$
$x_6+1x_2+1x_3=4$

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

状态 1: **************
目标函数: $2x_1+5x_2+2x_3+0$
限制函数:
$x_4+1x_1=2$
$x_5+1x_1+1x_2=3$
$x_6+1x_2+1x_3=4$
状态 2: **************
目标函数: $-2x_4+5x_2+2x_3+4$
限制函数:
$x_1+1x_4=2$
$x_5+-1x_4+1x_2=1$
$x_6+1x_2+1x_3=4$
状态 3: **************
目标函数: $3x_4+-5x_5+2x_3+9$
限制函数:
$x_1+1x_4=2$
$x_2+-1x_4+1x_5=1$
$x_6+1x_4+-1x_5+1x_3=3$
状态 4: **************
目标函数: $-3x_1+-5x_5+2x_3+15$
限制函数:
$x_4+1x_1=2$
$x_2+1x_1+1x_5=3$
$x_6+-1x_1+-1x_5+1x_3=1$
状态 5: **************
目标函数: $-1x_1+-3x_5+-2x_6+17$
限制函数:
$x_4+1x_1=2$
$x_2+1x_1+1x_5=3$
$x_3+-1x_1+-1x_5+1x_6=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这三行。
如果找不到任何一个”紧”的行，即 $\forall a(i,y)<=eps$ 那么这个变量想多大就能多大，反正有辅助变量顶着，可以直接返回无界【坑4】，即答案是无穷大。

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

$void~pivot(int~x,int~y)$ :

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

对偶原理：

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

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

【转置】 ：矩阵的元素 $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-5
using 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;
}