单纯形法与线性规划入门

前言

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

线性规划

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

\[\begin{aligned}max\ \ &\sum_{j=1}^nc_jx_j&\\s.t.\ \ &\sum_{j=1}^na_{i,j}x_j\le b_i&i=1,2,...,m\\&x_j\ge0&j=1,2,...,n\end{aligned} \]

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

\[\begin{aligned}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_j\ge0&j=1,2,...,n+m\end{aligned} \]

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

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

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

\(Pivot\)（转轴）

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

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

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

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

考虑\(x_j\ge 0\)并不是仅针对\(x_{1\sim n}\)，而同时要求\(x_{n+1\sim 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\not=e]\frac{a_{l,j}}{a_{l,e}}x_{j} \]

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

\(Init\)（初始化）

考虑如果存在某个\(b_i<0\)，我们就无法令\(x_{1\sim 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\n",-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;//单纯形法
}