java单纯形法_单纯形法 - fjzzq2002 - 博客园

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

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

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

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

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

8c2f8174003f2160596d8d047abc7d6c.png

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

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

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

867ce901e85ef27cf9986f9393db2c33.png

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

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

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

aa593e5a3d558f5ae871cc98bfc2a791.png

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

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

7ead3b25e7222d676e1eb926cb888bd2.png

额那不等号怎么办呢?

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

affe9335cec28f96dd251649c8e7a5e4.png

在松弛型等式左边的那些变量我们把它们叫做基变量,在上面的松弛型表示中就是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]成了基变量。

一开始我们知道0f9829182dd002cf6638a57c417d62c7.png

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

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

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

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

fc4ff82946d442b72fe5293ae8cf85dd.png

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

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

这个线性规划是这样的:

7cf04098c03f18fa67ad5618dd6b02be.png

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

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

86231ef694167e8a7e6001375200c0cc.png

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

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

6e0e1abc28f4b187f758585c74969e70.png

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

efe3778d64e034d54a0a51969fe85518.png

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

1ca7f17809a9de16b1cbf1574e48f95a.png

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

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

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

7d70248c03f5586580716342fd63bed4.png

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

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

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

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

例如:

e04507b206b9b2cd544688537671d854.png

我们引入辅助线性规划:

04790200e8b5ba662fed0a69a1ec4768.png

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

b408a89528f585246ca7e384fac1df91.png

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

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

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

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

85092f61e7f2919b4d071587fbcb6279.png

而-bi显然为正的。

而其他约束会变为:

548986973320892b789e3c5284fdf5c7.png

由于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个基变量,把系数取负就可以了。

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

#include #include#include#include#include#include#include#include#include#include

using namespacestd;

typedefdoubleld;#define SZ 233

intn,m,t,id[SZ];

ld a[SZ][SZ],vv[SZ];const ld eps=1e-7;intdcmp(ld x)

{if(xeps) 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];

}

}voidsolve()

{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

}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]);

}intmain()

{

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();

}

评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

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

抵扣说明:

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

余额充值