线性规划-单纯型

不依靠矩阵，能讲明白线性方程组吗

1.线性规划简介

线性规划用于处理以下问题：当线性方程组的方程个数小于变量个数，甚至线性方程组中包含不等号的时候，如何在所有可能的解中选出最优的。这里的“最优”形如 \(c_1x_1+c_2x_2+\dots +c_nx_n\) 取到最大值。

为了方便，线性规划还要求给出的解中 \(\forall x_i\ge 0\)，这虽然是个限制，但是下文可以看到这个限制反而会让运算变简单。

形式化地，线性规划对变量 \(x_1\dots x_n\) 求最大/最小的 \(c_1x_1+c_1x_2+\dots +c_nx_n\)，满足

\[\begin{cases} a_{1,1}x_1+a_{1,2}x_2+\dots+a_{1,n}x_n\le b_1\\ a_{2,1}x_1+a_{2,2}x_2+\dots+a_{2,n}x_n\le b_2\\ \vdots\\ a_{t_1,1}x_1+a_{t_1,2}x_2+\dots+a_{t_1,n}x_n\ge b_{t_1}\\ a_{t_1+1,1}x_1+a_{t_1+1,2}x_2+\dots+a_{t_1+1,n}x_n\ge b_{t_1+1}\\ \vdots\\ a_{t_1+t_2,1}x_1+a_{t_1+t_2,2}x_2+\dots+a_{t_1+t_2,n}x_n\ge b_{t_1+t_2}\\ a_{t_1+t_2+1,1}x_1+a_{t_1+t_2+1,2}x_2+\dots+a_{t_1+t_2+1,n}x_n= b_{t_1+t_2+1}\\ a_{t_1+t_2+2,1}x_1+a_{t_1+t_2+2,2}x_2+\dots+a_{t_1+t_2+2,n}x_n= b_{t_1+t_2+2}\\ \vdots\\ a_{m,1}x_1+a_{m,2}x_2+\dots+a_{m,n}x_n= b_{m}\\ x_1,x_2,\dots,x_n\ge 0 \end{cases}\]

看起来一堆东西，其实就是三个二维数组，所以写成矩阵会比较简单。但我想试试不用矩阵去讲线性规划。

题外话，这里面的线性可以理解为“一次”，就是说全部方程的 \(x_i\) 都是一次项。

2.标准线性规划

发现有两种最值，三种方程比较繁琐，可以减少情况。首先是求最值，如果想求最小值，将所有 \(c_i\) 变为相反数，问题就变为求最大值。注意这里并没有改变 \(x_i\ge 0\) 的限制。类似的， \(t_1\) 到 \(t_1+t_2\) 的大于等于号部分，也可以把系数全都取反，大于号变成小于号，所以现在只需要考虑小于等于号和等号。

因为小于等于其实就是差大于等于零，可以把差新设成一个变量，这样把方程中的小于等于变成了线性规划的固有性质 \(x\ge 0\)。

具体地，如果要求 \(a_1x_1+a_2x_2\dots a_nx_n\le b_1\)，就是要求 \(a_1x_1+a_2x_2\dots a_nx_n+x_{n+1}=b_1\) 且 \(x_{n+1}\ge 0\)。

这样，我们可以得到只有等号的线性方程组。添加的这个变量 \(x_{n+1}\) 称作松弛变量。

这样有 \(n’\) 个变量，\(m\) 条等式方程的线性规划，就是标准线性规划。

3.解的形态

这部分先假设 \(m<n\)。

然后考虑怎么解标准线性规划，这是一个线性方程组的问题。在 \(n\times n\) 的情况下，我们可以通过高斯消元使第 \(i\) 行只保留 \(x_i\) 系数不为 \(0\)。

那么 \(n\times m\) 的情况下，我们也可以进行 \(m\) 轮对 \(x_1,x_2\dots x_m\) 的高斯消元，使得第 \(i\) 行中只有 \(x_i\) 和 \(x_{m+1},x_{m+2}\dots x_{n}\) 系数不为 \(0\)。比如说洛谷的样例：

\[\begin{cases} 2x_1+x_2+x_3+0x_4=6\\ -x_1+2x_2+0x_3+x_4=3 \end{cases}\]

我们可以高斯消元变成

\[\begin{cases} 2x_1+0x_2+\frac{4}{5}x_3-\frac{1}{3}x_4=\frac{18}{5}\\ 0x_1+\frac{5}{2}x_2+\frac{1}{2}x_3+x_4=6 \end{cases}\]

这样如果 \(x_3,x_4\) 已知的话，就是

\[\begin{cases} 2x_1+0x_2=\frac{18}{5}-\frac{4}{5}x_3+\frac{1}{3}x_4\\ 0x_1+\frac{5}{2}x_2=6-\frac{1}{2}x_3-x_4 \end{cases}\]

变为

\[\begin{cases} x_1=\frac{9}{5}-\frac{2}{5}x_3+\frac{1}{6}x_4\\ x_2=\frac{12}{5}-\frac{1}{5}x_3-\frac{2}{5}x_4 \end{cases}\]

其中右侧是常数。也就是说不考虑 \(xi\ge 0\) 的限制，任意给定 \(x_3,x_4\)，都存在唯一 \(x1,x2\) 的取值。

我们发现每一行只要选择一个非零元去消其它行即可，第一行不一定要选 \(x_1\)，第二行也不一定要选 \(x_2\)，所以不考虑 \(0\) 的情况，我们可以用 \(m\) 个方程对 \(m\) 个变量 \(x_i\) 进行消元，也就是说，只要任意指定 \(n-m\) 个变量 \(x_{t_1},x_{t_2}\dots x_{t_n-m}\) 确定，剩下被消过的 \(m\) 个变量也可以确定。

考虑零的情况，其实我们也可以像高斯消元一样，假设第 \(i\) 个我们要消 \(x_{t_i}\)，那么在剩下的 \(i\) 到 \(m\) 行中任意选择一个 \(a_{j,t_i}\not=0\) 的行与第 \(i\) 行交换即可。

所以除非这一列全都为 \(0\)，或者像 \(x_1+x_2+x_3=4,2x_1+2x_2+x_4=5\) 一样，如果 \(x_1\) 和 \(x_2\) 中一个被消元，由于系数呈倍数关系，另一个对应的列会全变为 \(0\)，所以无法同时消元。这是因为 \(x_1+x_2\) 总是同时出现，故而只能解出 \(x_1+x_2\) 而不知道两者具体值。

否则对任意指定的 \(n-m\) 变量，都可以在随便填的情况下唯一地解出剩下 \(m\) 个变量的值。

所以在规避这两种情况后，我们有 \(n-m\) 个可以任意指定的变量，如果仍然没有 \(x_i\ge 0\) 的限制，完全可以直接由上面的第四个等式组算出来最终取值。

比如样例中要求 \(x_1+x_2\) 最大，就是 \(\frac{21}{5}-\frac{3}{5}x_3-\frac{7}{30}x_4\) 最大。

那么直接把 \(x_3\) 和 \(x_4\) 设为无穷小就行了。有限制的情况，不仅要考虑 \(x_3,x_4\)，还要考虑 \(x_1,x_2\) 的限制，所以要用到单纯型做法。

4.瞎说算法思想

主要思想是，\(x_3\) 和 \(x_4\) 的设定其实是独立的，如果自由设置的话，完全可以设成 \(0\) 考虑一个怎么设置，比如你设置 \(x_3\)，那么 \(x_1,x_2\) 都是关于 \(x_3\) 的一次函数，能求出来 \(x_3\) 的范围，再把 \(x_3\) 调到范围左侧/右侧。由于一次函数结构简单，范围要么没端点，要么端点是由 \(x_1,x_2\) 中的一个（或多个）解出来的，所以调完范围后 \(x_1\) 或者 \(x_2\) 会变成 \(0\)。

没端点的情况显然就是最终的解可以无限大/无限小，遇到就可以程序结束了。

这样做一次不一定最优，可能还要再找 \(n-m\) 个自由变量。为了方便，也为了尽可能选那些没到上界的，可以每次选 \(n-m\) 个为 \(0\) 的变量，之后每次更改一个变量，再把另一个变量变成 \(0\)，仍然存在这样的 \(n-m\) 个变量。

这样做到这 \(n-m\) 个变量变大都没收益得到的其实就是答案。

为什么呢？

首先答案一定有 \(n-m\) 个是为 \(0\) 的，因为想要达到最优，你选择的自由变量，不是把它们变成 \(0\) 就是把它们变成无穷大。所以但凡有超过 \(m\) 个没取到 \(0\) 的变量，把它当做自由变量，变大变小总有一个好的，故而肯定变成 \(0\) 更优。

总结一下过程，这相当于动态维护一个大小为 \(n-m\) 的集合，每次操作后的答案，就是集合内的数都为 \(0\) 时的答案。这样贪心操作是集合中删一个数，加一个数，并没有后效性，不会因为你把谁排除出了集合你就没法达到最终答案对应的集合了。

也不可能会出现集合想变成最终答案，一定会有一次操作将答案变小的情况，因为假设你让 \(x_i\) 增加并移除集合，那最终它一定又变成零。不然你最终随便指定自由变量，再让 \(x_i\) 减小，一次函数的 \(k\) 也没有变，仍可以让它减小并增大答案。

假设它变成零那一次被排除集合的是 \(x_j\)，那么这次完全可以改成让 \(x_j\) 移除集合。虽然 \(x_j\) 也可能会让答案变小，但这样证明下去总会矛盾。

P.S.注意这里的证明是“一定存在操作路径”，如果你的路径是随便找的，则不能假设这样的路径压缩。

于是大概就得到了答案，这个过程我们需要维护自由变量集合，还要维护剩下的 \(m\) 个数的一次函数。

5.单纯型

这部分再做一个小小的归化，初始的等式可以变成大于等于号和小于等于号，再变成两个添加了松弛变量的等式方程。

假设归化后等式方程共有 \(m\) 个，变量共 \(n\) 个，这样首先保证 \(n>m\)，其次这样做后全部的 \(m\) 个松弛变量系数都为 \(1\)，且一个等式只包含一个松弛变量。比如

\[\begin{cases} 2x_1+x_2+x_3+0x_4=6\\ -x_1+2x_2+0x_3+x_4=3\\ \end{cases}\]

其中 \(x_3,x_4\) 可以直接解出来。

这样我们初始选择 \(n-m\) 个非松弛变量作为自由变量，\(m\) 个松弛变量在 \(m\) 个方程中已经被消元好了，故而不会出现 3 里面说的无法消元的情况。

那么来模拟一下样例。初始时候集合是 \(\{x_1,x_2\}\)，有

\[\begin{cases} x_3=6-2x_1-x_2\\ x_4=3+x_1-2x_2 \end{cases}\]

这样改变 \(x_1\) 和 \(x_2\) 的贡献都是 \(1\)，因为就是 \(c_1\) 和 \(c_2\)，选择改变 \(x_1\)，由 \(x_3\) 的限制解出来 \(x_1\) 最多增加 \(3\)。

于是集合变成 \(x_2\) 和 \(x_3\)，有

\[\begin{cases} x_1=3-\frac{1}{2}x_2-\frac{1}{2}x_3\\ x_4=3-2x_2+x_1=6-\frac{5}{2}x_2-\frac{1}{2}x_3 \end{cases} \]

这样 \(x_2\) 的贡献是 \(c_2-\frac{1}{2}c_1=\frac{1}{2}\)，\(x_3\) 的贡献是 \(-\frac{1}{2}=-1\)，所以选择改变 \(x_2\)，由 \(x_4\) 的限制解出来最多增加 \(\frac{12}{5}\)。

于是集合变成 \(x_3\) 和 \(x_4\)，有

\[\begin{cases} x_1=\frac{9}{5}-\frac{2}{5}x_3+\frac{1}{6}x_4\\ x_2=\frac{12}{5}-\frac{1}{5}x_3-\frac{2}{5}x_4 \end{cases}\]

发现增加哪一个贡献都是负的，所以得到答案。

这个结果和我们 2 中做的高斯消元是一样的，可见每个时刻的答案只和集合有关，和过程无关。

求一次函数系数可以先把在集合里的移项到右边，再对左边高斯消元。

于是有代码如下，B[i][j]表示左边高斯消元后第 \(i\) 个式子的右边对左侧为 \(x_j\) 式子的贡献。

int vis[205];
struct jz{double a[205][205];}A,B;
double c[205],b[205],sb[205];
//sb 是 k 这一列中c[j]*B[k][j]的和
int sim(double &ans,int m,int n)
{
	for(int i=1;i<=m;i++) sb[i]=0;
	for(int i=1;i<=n;i++) if(!vis[i])
	for(int j=1;j<=m;j++)
	{
		sb[j]+=c[i]*B.a[j][i];
	}
	double mx=0;int pos=-1;
	for(int i=1;i<=n;i++)
	if(vis[i])
	{
		double det=0;
		for(int j=1;j<=m;j++) det-=A.a[j][i]*sb[j];
		//因为移到右边所以系数变负数，注意右边只有常数和xi两项
		//这里只需要斜率所以没记录常数项
		det+=c[i];
		if(det>mx) mx=det,pos=i;
	}
	double slo=mx;int pos0=pos;
	if(pos==-1) return 0;
	mx=inf;pos=-1;
	for(int i=1;i<=n;i++)
	if(!vis[i])
	{
		double fm=0,fz=0;
		for(int j=1;j<=m;j++)
		fz+=B.a[j][i]*b[j],
		fm+=A.a[j][pos0]*B.a[j][i];
		//这里分子是右边的常数项和，分母是挪过去的x_pos系数和
		if(fz<-eps) fz=-fz,fm=-fm;
		if(fm<eps) continue;
		//这种情况是x_pos越大贡献越大，所以xi无法限制x_pos
		if(fz<eps) {pos=i;mx=0;break;}
		fz/=fm;
		if(fz<mx) pos=i,mx=fz;
	}
	//如果没有能限制xi的x_pos则解无限
	if(pos==-1) puts("Unbounded"),exit(0);
	ans+=slo*mx;
	pivot(pos0,pos,m,n);
	return 1;
}

这里我维护的不在自由元内的变量一定是可解的，不会出现 2 中的反例。因为既然从集合删除 \(x_i\) 的时候 \(x_j\) 是瓶颈，在能解出其它值的情况下，\(x_i\) 和 \(x_j\) 一定仍然是一次函数关系，这样新增需要表示的 \(x_i\) 一定可以被表示，而原来需要 \(x_i\) 来表示的也都可以用 \(x_j\) 表示。

6.优化高斯消元

其实可以不通过高斯消元维护 \(B\) 的值，因为每次只拿走一个元素拿进来一个元素。

比如说向集合添加 \(x_a\)，并将 \(x_b\) 从集合内删除，现在不在集合内，需要解方程的有 \(x_1,x_1\dots x_t,x_b\)。仍然按照原先的 \(B\) 算 \(x_i\) 的话，因为 \(x_b\) 回到了等式左边，算出的值和正确的值应该查了 \(x_b\) 乘以一个常数。

而按照原先的 \(B\) 算 \(x_a\) 的方法，则能直接计算出 \(x_b\)，因为 \(x_a\) 被挪到等式右边，左边只留下了 \(x_b\)。

于是可以先知道怎么计算 \(x_b\)，再把计算 \(x_b\) 的这一列按照系数乘到其它列，消除其它列多余的 \(x_b\) 即可。

这样能在 \(O(m^2)\) 的复杂度内完成一次操作，我们把操作称作转轴。

代码实现如下：

int vis[205];
struct jz{double a[205][205];}A,B;
void pivot(int x,int y,int m,int n)
{//转轴操作维护逆矩阵保证复杂度
	vis[x]=0;
	vis[y]=1;
	double tmp=0;
	for(int i=1;i<=m;i++)
	B.a[i][x]=B.a[i][y],B.a[i][y]=0,
	tmp+=B.a[i][x]*A.a[i][x];
	for(int i=1;i<=m;i++) B.a[i][x]/=tmp;
	for(int i=1;i<=n;i++)
	if(!vis[i]&&i!=x)
	{
		tmp=0;
		for(int j=1;j<=m;j++)
		tmp+=B.a[j][i]*A.a[j][x];
		for(int j=1;j<=m;j++)
		B.a[j][i]-=B.a[j][x]*tmp;
	}
}

7.初始解

上面是 \(b_i\ge 0\) 的情况，初始的集合选择所有非松弛变量是可行的，但是如果 \(b_i<0\)，这么选会导致 \(x_i\) 是负数而无解。

技巧是新增变量 \(x_{n+1}\)，并在每个式子左侧减去 \(x_{n+1}\)，比如样例3（去掉了没用的 \(x_3\)，下标顺延了）

\[\begin{cases} -2x_1+x_2+x_3=-4\\ x_1+x_2+x_5=4\\ x_1-2x_2+x_6=-4\\ \end{cases} \]

变成

\[\begin{cases} -2x_1+x_2+x_3-x_7=-4\\ x_1+x_2+x_5-x_7=4\\ x_1-2x_2+x_6-x_7=-4\\ \end{cases} \]

你选 \(b_i<0\) 中 \(b_i\) 最小的那一行，并让那一行的 \(x_i=0\) 的话，就存在 \(x_7=4\) 的一组解，因为此时每个方程都只有松弛变量 \(x_j\) 以及 \(x_7\)，然后 \(x_7\) 移项到右边，右边的数值肯定会变成正数，所以存在一组初始解。

然后先运行一个让 \(-x_7\) 最大的单纯型，如果 \(-x_7\) 能达到 \(0\) 则得到了没有 \(x_7\) 的一组解，可以重新运行单纯型；否则就是原线性规划无解。

因为如果线性规划有解，\(x_7=0\) 也肯定是加入 \(x_7\) 后线性规划的解嘛。

于是就得到了初始解，接着就能完美地通过洛谷模板题辣！

代码：

#include<bits/stdc++.h>
using namespace std;
int vis[205];
struct jz{double a[205][205];}A,B;
void pivot(int x,int y,int m,int n)
{//转轴操作维护逆矩阵保证复杂度
	vis[x]=0;
	vis[y]=1;
	double tmp=0;
	for(int i=1;i<=m;i++)
	B.a[i][x]=B.a[i][y],B.a[i][y]=0,
	tmp+=B.a[i][x]*A.a[i][x];
	for(int i=1;i<=m;i++) B.a[i][x]/=tmp;
	for(int i=1;i<=n;i++)
	if(!vis[i]&&i!=x)
	{
		tmp=0;
		for(int j=1;j<=m;j++)
		tmp+=B.a[j][i]*A.a[j][x];
		for(int j=1;j<=m;j++)
		B.a[j][i]-=B.a[j][x]*tmp;
	}
}
const double eps=1e-7,inf=1e18;
double c[205],b[205],sb[205];
//sb 是 k 这一列中c[j]*B[k][j]的和
int sim(double &ans,int m,int n)
{
	for(int i=1;i<=m;i++) sb[i]=0;
	for(int i=1;i<=n;i++) if(!vis[i])
	for(int j=1;j<=m;j++)
	{
		sb[j]+=c[i]*B.a[j][i];
	}
	double mx=0;int pos=-1;
	for(int i=1;i<=n;i++)
	if(vis[i])
	{
		double det=0;
		for(int j=1;j<=m;j++)
		det-=A.a[j][i]*sb[j];
		det+=c[i];
		if(det>mx) {mx=det,pos=i;break;}
		//使用 bland 原则
	}
	double slo=mx;int pos0=pos;
	if(pos==-1) return 0;
	mx=inf;pos=-1;
	for(int i=1;i<=n;i++)
	if(!vis[i])
	{
		double fm=0,fz=0;
		for(int j=1;j<=m;j++)
		fz+=B.a[j][i]*b[j],
		fm+=A.a[j][pos0]*B.a[j][i];
		if(fz<-eps) fz=-fz,fm=-fm;
		if(fm<eps) continue;
		if(fz<eps) {pos=i;mx=0;break;}
		fz/=fm;
		if(fz<mx) pos=i,mx=fz;
	}
	if(pos==-1) puts("Unbounded"),exit(0);
	ans+=slo*mx;pivot(pos0,pos,m,n);
	return 1;
}
double tmp[202];
int main()
{
	int n,m,pos=-1;cin>>n>>m;
	double mi=inf;c[n+m+1]=-1;
	for(int i=1;i<=n;i++) scanf("%lf",&tmp[i]),vis[i]=1;
	for(int i=1;i<=m;i++)
	{
		for(int j=1;j<=n;j++)
		scanf("%lf",&A.a[i][j]);
		A.a[i][n+i]=1;
		scanf("%lf",&b[i]);
		if(b[i]<mi) mi=b[i],pos=i;
		A.a[i][n+m+1]=-1;
	}
	if(mi<-eps)
	{//初始解是 n+m+1 取最小的那个-b，这样松弛变量都是正的
		vis[n+pos]=1;vis[n+m+1]=0;
		for(int i=1;i<=m;i++)
		{
			if(i==pos) B.a[i][n+m+1]=-1;
			else B.a[i][n+i]=1,
			B.a[pos][n+i]=-1;
		}
	}
	else
	{
		vis[n+m+1]=1;
		for(int i=1;i<=m;i++)
		B.a[i][n+i]=1;
	}
	double ans=mi;
    //这里单纯型的ans是-x_{n+m+1}，所以它的初始值就是最小的-b
	while(1){if(!sim(ans,m,n+m+1))break;}
	if(ans<-eps) puts("Infeasible"),exit(0);
	if(!vis[n+m+1])
	{
		for(int i=1;i<=n+m;i++)
		if(vis[i])
		{pivot(i,n+m+1,m,n+m+1);
		break;}
	}
	ans=0;
	for(int i=1;i<=n;i++) c[i]=tmp[i];
	for(int i=n+1;i<=n+m;i++) c[i]=0;
	for(int i=1;i<=n;i++)
	if(!vis[i])
	{
		double now=0;
		for(int j=1;j<=m;j++)
		now+=B.a[j][i]*b[j];
		ans+=c[i]*now;
        //这里先求一下初始解对应的ans
	}
	while(1){if(!sim(ans,m,n+m))break;}
	printf("%.12lf\n",ans);
	for(int i=1;i<=n;i++)
	if(!vis[i])
	{
		double now=0;
		for(int j=1;j<=m;j++)
		now+=B.a[j][i]*b[j];
		printf("%.12lf ",now);
	}
	else printf("0 ");
	return 0;
}

8.总结

似乎没有什么题用到单纯型呢，只有 oiwiki 上的模板题 志愿者招募，而且也可以用网络流做。

但是我觉得单纯型算法对于理解线性规划是非常有意义的，对于没学过线性代数的高中生来说，能够加深对于向量的理解。

并且线性规划本身很有意义啊。所以就这样（？