单纯形法与线性规划入门

Posted 2020-12-17 chenxiaoran666

前言

说实话其实这个算法真的不难，只要你有耐心把它啃完。。。

线性规划

线性规划的标准型长这个样子：

[egin{align}max &sum_{j=1}^nc_jx_j&\s.t. &sum_{j=1}^na_{i,j}x_jle b_i&i=1,2,...,m\&x_jge0&j=1,2,...,nend{align} ]

它的松弛型则是长这个样子：

[egin{align}max &sum_{j=1}^nc_jx_j&\s.t. &x_{n+i}=b_i-sum_{j=1}^na_{i,j}x_j&i=1,2,...,m\&x_jge0&j=1,2,...,n+mend{align} ]

考虑和原本的标准型相比，松弛型有什么优势。

容易发现，原本限制条件中的不等号，被我们转化成了等号。

这一步转化是之后所有操作的基础。

(Pivot)（转轴）

考虑到我们需要最大化的式子(sum_{j=1}^nc_jx_j)。

如果我们能够通过某些手段把(x)的系数全部都变成负数，那么因为(x_jge0)，只要令(x_{1sim n}=0)即可得到最大值。

但我们该怎样把(x)的系数都变成负数呢？

其实，只要通过等量代换就可以了。

考虑(x_jge 0)并不是仅针对(x_{1sim n})，而同时要求(x_{n+1sim n+m})也满足这个条件。

而它们之间存在一定的等量关系，就可以试着通过等量代换来调整所需最大化的式子。

我们规定把(x_e)用(x_{n+l})代换的过程叫做(Pivot(l,e))。

已知：

[x_{n+l}=b_l-sum_{j=1}^na_{l,j}x_j ]

我们移出右式中的(x_e)，将它和(x_{n+l})调换位置得到：

[x_e=frac{b_l}{a_{l,e}}-frac1{a_{l,e}}x_{n+l}-sum_{j=1}^n[j ot=e]frac{a_{l,j}}{a_{l,e}}x_{j} ]

然后我们再把其余式子中的(x_e)全部用右式代换，就彻底抹消了(x_e)的存在，完成了一次转轴。

(Init)（初始化）

考虑如果存在某个(b_i<0)，我们就无法令(x_{1sim n}=0)，因为这样就不满足限制了。

仔细观察转轴操作，可发现这一操作过程中并不会改变(b_i)的正负性，因此只要我们在一开始把所有(b_i)都强制修改为正数即可。

为了达成这一目标，每次我们找出一个(b_l<0)，并找到一个(a_{l,e}<0)，然后只要执行(Pivot(l,e))即可。

如果找不到，显然无解。

(Simplex)（单纯形法）

单纯形法主过程的核心问题，就是每次应该对哪一对(l,e)执行转轴操作：

• 首先，需要满足(c_e>0)，且转轴后可以让新的(c_e<0)（这是我们执行转轴操作的初衷）。

• 其次，在所有可选的(l)中，我们要使(frac{b_l}{a_{l,e}})最小（要选择最紧的限制）。

只要不断进行(Pivot(l,e))，直到无法继续就可以了。

代码

【UOJ179】线性规划（据说因为出题人卡常卡精度，只要拿到(97)分就可以了？）

#include<bits/stdc++.h>
#define Tp template<typename Ty>
#define Ts template<typename Ty,typename... Ar>
#define Reg register
#define RI Reg int
#define Con const
#define CI Con int&
#define I inline
#define W while
#define N 20
#define DB double
#define eps 1e-8
using namespace std;
int n,m,op;DB a[N+5][N+5];
namespace SimplexMethod//单纯形法
{
	int id[2*N+5];DB s[N+5];
	I void Pivot(CI l,CI e)//转轴
	{
		RI i,j;DB t=a[l][e];swap(id[n+l],id[e]);//交换编号
		for(a[l][e]=1,i=0;i<=n;++i) a[l][i]/=t;//对这个式子变形
		for(i=0;i<=m;++i) if(i^l&&fabs(a[i][e])>eps)//枚举其他式子
			for(t=a[i][e],a[i][e]=j=0;j<=n;++j) a[i][j]-=t*a[l][j];//等量代换
	}
	I bool Init()//初始化，使所有b[i]为正
	{
		RI i,l,e;W(1)//不断操作
		{
			for(l=e=0,i=1;i<=m;++i) a[i][0]<-eps&&(!l||rand()&1)&&(l=i);if(!l) break;//找到b[l]<0
			for(i=1;i<=n;++i) a[l][i]<-eps&&(!e||rand()&1)&&(e=i);//找到a[l][e]<0
			if(!e) return puts("Infeasible"),0;Pivot(l,e);//找不到则无解，否则转轴
		}return 1;
	}
	I bool Simplex()//单纯形法主过程
	{
		RI i,l,e;DB Mn;W(1)//不断操作
		{
			for(l=e=0,Mn=1e9,i=1;i<=n;++i) if(a[0][i]>eps) {e=i;break;}if(!e) break;//找到c[e]>0
			for(i=1;i<=m;++i) a[i][e]>eps&&a[i][0]/a[i][e]<Mn&&(Mn=a[i][0]/a[i][e],l=i);//找到限制最紧的l
			if(!l) return puts("Unbounded"),0;Pivot(l,e);//如果所有a[l][e]<0无穷大，否则转轴
		}return 1;
	}
	I void Solve()//单纯形法
	{
		RI i;for(srand(20050521),i=1;i<=n;++i) id[i]=i;if(Init()&&Simplex())
		{
			if(printf("%.10lf
",-a[0][0]),!op) return;//输出答案
			for(i=1;i<=m;++i) s[id[n+i]]=a[i][0];for(i=1;i<=n;++i) printf("%.10lf ",s[i]);//输出方案
		}
	}
}
int main()
{
	RI i,j;for(scanf("%d%d%d",&n,&m,&op),i=1;i<=n;++i) scanf("%lf",a[0]+i);//读入
	for(i=1;i<=m;++i) {for(j=1;j<=n;++j) scanf("%lf",a[i]+j);scanf("%lf",a[i]);}//读入
	return SimplexMethod::Solve(),0;//单纯形法
}