单纯形法

看了集训队答辩,感觉要学习的有杜教筛高级版、线性规划、FFT、仙人掌、高级版线段树

不出意外的话一个月内博客内都不会有别的东西了QAQ

首先是喜闻乐见的单纯形法解线性规划。

今年(2016年)和线性规划有关的集训队论文有两篇,大家可以自行翻一下集训队论文(当然如果你没有拿到你可以去UOJ群下载啊),下面的大部分内容都是参阅akf那篇

线性规划的标准型一般长得像这样:

image

一般我们拿到的都是像标准型这样的问题,例如网络流问题,我们就是要最大化流量并且让每个点流量守恒。

假设我们把汇点到源点连一条容量为+∞的边那么所有点就都流量守恒了。我们用c(u,v)表示(u,v)这条边的容量,f(u,v)表示这条边的流量。

那么我们不难得到这样的一个标准型线性规划

image

不难把它转化成标准型这样的模型。

而最小费用流则也是一个标准型,我们还是把汇点到源点加一条容量为+∞,费用为0的边。

那么也可以得到一个比较标准的线性规划:

image

(较论文中的配文有修改)

如果你要最小费用最大流的话只要加一个约束条件$\sum_{(u,v)}{f(u,v)}=maxflow$image

额那不等号怎么办呢?

对于不等约束条件${\sum_{j=1}^n{a_{ij}x_j}}\leq{b_i}$image

在松弛型等式左边的那些变量我们把它们叫做基变量,在上面的松弛型表示中就是x[n+1]…x[n+m],非基变量就是在等式右边的那些,在上面的表示中是x[1]…x[n]。

显然当bi非负时令所有基变量为bi,非基变量为0,即可得到一个满足条件的初始解。(至于bi可能为负的情况下面在单独说)

单纯形法有一个重要的操作,叫转轴(pivot)操作。就是说我们可以把一个基变量x[b]和非基变量x[n]互换,用x[b]和其他非基变量代换这个x[n],这样x[b]就成了非基变量,而x[n]成了基变量。

一开始我们知道image

那么简单代换一下会发现image,然后我们也可以把其他约束条件中的x[n]这样代换掉,于是就完成了这个转轴操作。

然后有了这个转轴过程的帮助,我们就可以实现线性规划算法了。

首先为了最大化目标函数,我们考虑不停地找一个系数为正的非基变量,然后增大这个x。

具体的这个最优化过程如下:

image

(注意邹逍遥原论文在这里的做法是选择满足image的,ysy告诉我这样是不科学的,akf的论文里写的也是ce>0,于是做了一些修改)

这是为什么呢?akf在论文中给出了一个例子,从中我们可以大概理解一下。

这个线性规划是这样的:

image

我们考虑第一步我们选择换正系数的x1,因为增大x1必然会增大目标函数。

我们分别考虑三个限制,显然第一个限制下x1<=30,第二个限制下x1<=12,第三个限制x1<=9,然后x1<=9的下界最紧,所以我们用x1代换x6,得到下面的新线性规划:

image

哇目标函数增大到了27,可喜可贺。

接下来假设我们来增大x3,类似地可以得到:

image

然后可以增大x2,可以得到:

image

由于所有非基变量系数均为负的,所以我们不可能再得到更小的解了。

image

咦,有没有注意到上面漏了什么东西说要下面讲啊…

对了,在bi为负的时候,如果你把所有基变量带0,而非基变量就不一定是一组合法的初始解。

这时候我们就需要一个初始化操作,初始化操作的基本思想是引入一个辅助线性规划:

image

如果我们求得了这个辅助线性规划的最优解,那么如果最优解中x0>0显然原线性规划无解。

如果最优解中x0=0,那么x[1]…x[n+m]就是原线性规划的一个可行解。

而我们容易构造出辅助线性规划的一个可行解。

我们只要把x0作为换入变量,找到bi的最小值bl,把x[i+l]用x0来替代。

例如:

image

我们引入辅助线性规划:

image

bi最小值是-4,那么我们把x4用x0来替代:

image

然后我们就可以得到一组x的初始解用来进行最优化操作。

讲道理:为什么bi这样之后一定为正?

首先bi的最小值bl一定为负(废话,否则就不用进行最优化操作了)

然后我们考虑第l个约束会变为:

image

而-bi显然为正的。

而其他约束会变为:

image

由于bl是最小值,所以bi>=bl,所以-bl+bi>=0。

而我们发现UOJ上似乎有一种更加神奇的初始化方法?

我们的目标显然是让所有系数bi都为非负的。

我们选择一个bi为负的基变量x[i+n],然后我们在该约束右边找一个系数为正的非基变量,然后进行转轴操作,如果没有系数为正显然就无解了。

额其实这和这个初始化操作本质上是一样的。

所以这样就可以完成整个线性规划的过程。

我们来分析一下时间复杂度?

pivot操作的复杂度显然是O(NM)的,但是最优化操作中pivot操作的调用次数可能会成为指数级。

但是我们可以发现要达到这个指数级的调用次数,边权也必须是指数级的。所以在OI中往往跑得比谁都快。所以“能在1s内跑出范围为几百的数据”。

然而一般线性规划是可以在多项式范围内求解的,只不过…

椭球算法 O(n^6*m^2)

内点算法 O(n^3.5*m^2)

改进的内点算法 O(n^3.5*m)

在oi中一般并不能有什么卵用

所以单纯形法就这么讲完啦。

实现的时候可以发现,写单纯形并不要真的把松弛型建出来,只要假装第i个约束条件就对应第i个基变量,把系数取负就可以了。

科学的代码 http://uoj.ac/submission/61872 (worse)

我的代码(有一些细节写的比较丑

#include <iostream>
#include <stdio.h>
#include <stdlib.h>
#include <algorithm>
#include <string.h>
#include <vector>
#include <math.h>
#include <limits>
#include <set>
#include <map>
using namespace std;
typedef double ld;
#define SZ 233
int n,m,t,id[SZ];
ld a[SZ][SZ],vv[SZ];
const ld eps=1e-7;
int dcmp(ld x)
{
    if(x<-eps) return -1;
    if(x>eps) return 1;
    return 0;
}
void pivot(int r,int c) //基变量r和非基变量c 
{
    swap(id[r+n],id[c]);
    ld x=-a[r][c]; a[r][c]=-1;
    for(int i=0;i<=n;i++) a[r][i]/=x;
    for(int i=0;i<=m;i++)
    {
        if(dcmp(a[i][c])&&i!=r);else continue;
        x=a[i][c]; a[i][c]=0;
        for(int j=0;j<=n;j++) a[i][j]+=x*a[r][j];
    }
}
void solve()
{
    for(int i=1;i<=n;i++) id[i]=i;
    while(1) //init
    {
        int x=0,y=0;
        for(int i=1;i<=m;i++)
        {
            if(dcmp(a[i][0])<0&&(!x||(rand()&1))) x=i; 
        }
        if(!x) break;
        for(int i=1;i<=n;i++)
        {
            if(dcmp(a[x][i])>0&&(!y||(rand()&1))) {y=i; break;}
        }
        if(!y) {puts("Infeasible"); return;}
        pivot(x,y);
    }
    while(1) //simplex
    {
        int x=0,y=0;
        for(int i=1;i<=n;i++)
        {
            if(dcmp(a[0][i])>0&&(!x||(rand()&1))) x=i; 
        }
        if(!x) break;
        double w,t; bool f=1;
        for(int i=1;i<=m;i++)
        {
            if(dcmp(a[i][x])<0&&((t=-a[i][0]/a[i][x]),f||t<w)) w=t, y=i, f=0;
        }
        if(!y) {puts("Unbounded"); return;}
        pivot(y,x);
    }
    printf("%.9lf\n",a[0][0]);
    if(!t) return;
    for(int i=1;i<=n;i++) vv[i]=0;
    for(int i=n+1;i<=n+m;i++) vv[id[i]]=a[i-n][0];
    for(int i=1;i<=n;i++) printf("%.9lf ",vv[i]);
}
int main()
{
    scanf("%d%d%d",&n,&m,&t);
    for(int i=1;i<=n;i++) scanf("%lf",&a[0][i]);
    for(int i=1;i<=m;i++)
    {
        for(int j=1;j<=n;j++) scanf("%lf",&a[i][j]), a[i][j]*=-1;
        scanf("%lf",&a[i][0]);
    }
    solve();
}
posted @ 2016-05-04 06:33 fjzzq2002 阅读(...) 评论(...) 编辑 收藏