【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

我们可以得到以下的限制:
目标函数: 2x1+5x2+2x3+0
限制函数:
x12
x1+x23
x2+x34

然后若 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-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;
}
评论 3
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值