uoj 179 线性规划 单纯形法模板

缘起

了解到高中学习的线性规划问题可以有统一的解法（其实大学运筹学已经教过这种方法了, 只是当时觉得这种算法没什么美感~）. 于是来学习一下acm中如何将单纯形法的板子耍的酷一点~ uoj 179 线性规划

分析

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
这是一道模板题。

（这个题现在标程挂了。。哪位哥哥愿意提供一下靠谱的标程呀？） 现在的新标程有没有问题啊？

本题中你需要求解一个标准型线性规划：

有 n 个实数变量 x1,x2,…,xn 和 m 条约束，其中第 i 条约束形如 ∑_{j=1,..,n}aijxj≤bi。

此外这 n 个变量需要满足非负性限制，即 xj≥0。

在满足上述所有条件的情况下，你需要指定每个变量 xj的取值，使得目标函数 F=∑_{j=1,..,n}cjxj的值最大。

输入格式
第一行三个正整数 n,m,t。其中 t=0或者1。

第二行有 n 个整数 c1,c2,⋯,cn，整数间均用一个空格分隔。

接下来 m 行，每行代表一条约束，其中第 i 行有n+1个整数 ai1,ai2,⋯,ain,bi，整数间均用一个空格分隔。

输出格式
如果不存在满足所有约束的解，仅输出一行 "Infeasible"。

如果对于任意的 M，都存在一组解使得目标函数的值大于 M，仅输出一行"Unbounded"。

否则，第一行输出一个实数，表示目标函数的最大值F。当第一行与标准答案的相对误差或绝对误差不超过1e−6，你的答案被判为正确。

如果t=1，那么你还需要输出第二行，用空格隔开的n个非负实数，表示此时 x1,x2,…,xn 的取值，如有多组方案请任意输出其中一个。

如果 t=0，或者出现 Infeasible 或 Unbounded 时，不需要输出第二行。

样例一
input
2 2 1
1 1
2 1 6
-1 2 3

output
4.2
1.8 2.4 

explanation
两条约束分别为 2x1+x2≤6,−x1+2x2≤3

当 x1=1.8,x2=2.4时目标函数 x1+x2 取到最大值 4.24.2。

样例二
input
2 2 1
1 -1
1 1 4
-1 -2 -2

output
4.0
4.0 0.0

explanation
注意 xj≥0的限制。

样例三
input
3 3 1
0 0 1
-2 1 0 -4
1 1 0 4
1 -2 0 -4

output
Infeasible

样例四
input
2 1 1
0 1
1 0 1

output
Unbounded

限制与约定
对于所有数据，1≤n,m≤20，0≤|aij|,|bi|,|cj|≤100，t∈{0,1}。

本题包含 4 个子任务，每个 25 分。

子任务 1,3 满足 bi≥0。

子任务 2,4 没有特殊限制。

子任务 1,2 中 t=0。

子任务 3,4 中 t=1

时间限制：1s
空间限制：256MB

本文基于我对 【2】这篇好文章的内容的吸收写下的.

高中时代（初中?）就学过线性规划，但是那个时候学习的都是二维的线性规划. 之所以仅仅是二维, 是因为二维线性规划可以画坐标图解决. 但是一旦到了高维——画图就不好使了. 于是数学家丹齐格在1947年就发明了这种算法——单纯形法（simplex method）求解一般的线性规划问题.

所谓线性规划问题（下简称LP问题）就是求一个线性的目标函数在若干线性的约束条件下的某种最值（最大值或者最小值）

而所谓单纯形（simplex） 就是n维空间中n-1维超平面围成的凸区域（n-1维超平面围成的区域显然是凸的）.

所以LP（linear projection 线性规划）问题本质就是一个目标函数在一个凸区域上的极值问题. 不难知道, 一定在simplex的边界取得极值. 最优解可能是唯一存在，也可能无穷多个，甚至可能无界解（因此没有最优解，只有更优解）. 还有可能无解（例如一个LP问题既要求x1>=1又让x1<=0, 单纯形为空集, 其实只要单纯形不空, 就一定有可行解, 只是唯一最优解还是无穷多最优解还是无界解的区别而已）. 关于LP问题解的分类的上述叙述其实就是单纯形基本定理. 虽说是定理, 但是直觉上很容易理解.

首先, 关于单纯形的算法的严格数学理论参见算法导论第29章. 这里和最大流一样（【1】）, 仅仅采用意识流但是不失主线的叙事手法讲述单纯形算法这个故事.

标准型

首先我们将要处理的LP问题归为如下标准型.

其中A是m*n 矩阵, c和x是n维列向量. b是m维列向量. 注意, 只对x的符号做了限制，其余A、b、c都没有做任何限制. 我们完全可以使用 (A,b,c) 简写一个单纯形的标准型.

任意一个LP问题都可以从下面四个方面来转换为标准型.

  1. 目标函数是最小化, 而不是最大化.

    例如目标是极小化 -2x1+3x2, 则其实只需要求最大化的 2x1-3x2, 然后相反数就是答案

  2. 可能有变量不具备>=0约束

    例如xi这个变量不具备xi>=0的约束, 则可以令 xi=xi1-xi2, 其中xi1>=0, xi2>=0 即可, 然后将用到xi的地方都替换为xi1-xi2即可

  3. 可能有等式约束

    f(x1,…,xn)==c 等价于 f(x1,…,xn)>=c && f(x1,…,xn)<=c 2个条件

  4. 可能有>=bi的约束

    f(x1,…,xn)>=b 等价于 -f(x1,…,xn)<=-b

综上4点, 我们就可以将任何一个LP问题转换为单纯形的标准型.

松弛型

上一节我们将任意的LP问题转换为了单纯形的标准型. 下面就要求解单纯形的标准型. 因为单纯形的标准型是<=, 所以求解其实并不方便. 遂牛人们想出了继续将单纯形的标准型转换为松弛型. 转换方法如下

不是 Ax<=b 吗? 遂存在一个m维的列向量X=(xn+1,…,xn+m) 满足 Ax+X=b. X的m个分量就叫做松弛变量.（图论中的最短距离算法中也出现过松弛这一术语. 其实算法导论的29.2节以及【3】都告诉我们单源最短路问题其实可以转换为LP问题）

按照这种思想，我们可以如下定义出一个单纯形的松弛型

其中z是目标函数, {cj, j属于N} 就是标准型中的c, {aij}就是标准型中的A, {bi}就是标准型中的b（一个m维列向量）. N={1,….,n}, B={n+1,…,n+m}. 并且x1,…,xn+m 都是非负变量. 这里, {xi, $x\in N$} 我们称之为非基变量（N的含义就是non-basis）, 而 {xi, $i\in B$} 称之为基变量（B的含义就是basis） 目标函数z在整个单纯形算法的过程中都仅仅和非基变量有关, 和基变量无关. 而且基变量的特征是在(29.43)中的每个约束中仅仅出现一次,而且出现在最左侧.

同LP问题的单纯形的标准型，我们以 (N,B,A,b,c,v) 这个元组表示一个单纯形的松弛型.

ps: 比较单纯形的标准型和松弛型的表达元组之间的区别，前者是(A,b, c) 后者是(N,B,A,b,c,v), 后者多出了N、B、v三个东西. 为什么会多出这三个东西? 不急, 因为我们马上就会看到(N,B,A,b,c,v) 这5个东西都会在松弛型求解的过程中发生变化，这5个东西完整的刻画了当前松弛型的求解状态. 所以我们需要这5个东西共同来刻画一个松弛型的问题.

求解

不难看出标准型和松弛型是等价的. 下面来求解松弛型. 我们看看(29.42)这个松弛型的目标函数. 显然，因为变量都是非负的. 所以如果z中的cj(j在N中)都是<=0的话, 则已经求出了z的最大值——就是v.

下面假设存在某个cj系数>0, 那么为了增大v我们应该做什么事情呢? 当然是增大xj啊~ 但是它能无限增大吗? 如果能的话, 则本问题没有最优解只有更优解（其实一句话——它的解无界——unbounded）. 那么什么情况下会unbounded? 就是在增大xj的同时, xj越大，所有约束条件越能满足的时候. 那什么时候会这样呢? 就是所有的aij(i属于B)<=0的时候. 例如下面的松弛型

1
2
3
4
z=100+3x1-4x2
x3=46-(-3x1)+2x2
x4=1-3x1+x2
x1,x2,x3,x4>=0

增大x1显然是能让目标函数z增大的, 而且x1越大, x3越能>=0（因为x1前面的系数为-3<=0）. 所以如果出现了这种情况的话, 则算法直接返回unbounded.

如果不出现这种情况呢? 则我们可以根据bi(i属于B)们以及>0的aij们得到关于xj一个最紧的上界限制lim. 这个最紧的上界估计 xj<=lim 就刻画了xj最多能增长多少. 怎么理解这句话? 我们修改一下上面的例子

1
2
3
4
z=100+3x1-4x2
x3=46-3x1+2x2
x4=1-3x1+x2
x1,x2,x3,x4>=0

则x1前的系数3>0, 所以考虑增大x1, 但是x1不能无限增大——因为会导致x3、x4<0.

  • 为了不让x3<0, 我们有46-3x1>=0 即 x1<=46/3
  • 为了不让x4<0, 我们有1-3x1>=0 即 x1<=1/3

因为46/3>1/3, 所以我们得到关于x1一个最紧的上界估计就是 x1<=1/3, 即lim=1/3 即x1最多能增长到1/3.

然后发现, x4是限制x1要满足 x1<=lim的基变量. 所以考虑做pivot操作（OI圈将其称为转轴操作）. 所谓转轴操作，就是用x4根据松弛型的等式约束（这里是x4=1-3x1+x2）求解出x1来, 得到

1
x1=1/3+1/3*x2-1/3*x4

然后将z中的x1以及其他松弛型的等式约束中的x1都替换成上面的表达式（即x1用x4表达出来）。这就是pivot操作的全部. 所以pivot操作就是交换非基变量（x1）和基变量（x4）——让原本的基变量（x4）变为非基变量（x1）, 而基变量（x4）成为非基变量（x1）. 为什么需要这样做? 这需要抽象的数学表达. 我们在（29.42）、（29.43）构建的松弛型中做这种替（转）换（轴）操作之后松弛型变为（假设是非基变量xe和基变量xl之间的交换, $e\in N$, $l\in B$）

令 $N’=N\backslash\cup {l}, B’=B\backslash\cup {e} $ , 这就是新的非基变量和基变量. 不难知道pivot操作的复杂度是O(n*m)的. 上述公式称之为pivot公式.

这就是完成了一次pivot(N,B,A,b,c,l,e)操作. 值得注意的是新的v, 变成了 v+ce*bl/a_{le} 而ce>0, bl（我们先假设所有的bi{i属于B}都是>=0的, 即最初的松弛型的基本解就是可行解，所谓基本解就是让非基变量都变成0的解，则基变量就是bi们, 注意, b和基变量都是m维的列向量, 至于基变量不是可行解的情形我们后面再讨论）非负, a_{le}>0（这里a_{le}是使得上界限制的最紧的那个, 如果a_{le}<=0的话, 则unbounded），所以新的v>=老的v

即我们通过转轴操作增大了v. 再来看看我们的bi们是否满足约束（即仍然使得基本解是可行解），是的, 因为新的bi=(bi-a_{ie}/a_{le}*bl), 如果a_{ie}<=0, 则新的bi一定>=0, 否则, a_{ie}>0, 但是我们对a_{le}的选择策略是最紧的上界, 所以 bi/a_{ie} >= bl/a_{le}, 所以依旧能保证新的bi>=0, 即依旧基本解是可行解.

也就是说——如果我们假定初始的松弛型的基本解就是可行解的话, 则通过转轴操作, 我们能将松弛型转换为一个等价的松弛型（即得到的目标函数的解是一样的，下同），而且新的v还能增大. 而算法导论证明这个过程必将在有限步数内结束——要么发现所有的cj<=0了, 要么发现unbounded了.

下面来讨论一下——如果初始的松弛型基本解不是可行解的话会怎么样? 例如下面的松弛型

1
2
3
4
5
最大化		2x1-x2
约束
		  x3=2-2x1+x2
		  x4=-4-x1+5x2
		  xj>=0, j=1,2,3,4

上面的松弛型的基本解显然不是可行解——因为b2=-4<0. 我们该如何处理呢? 换言之, 我们该如何将其转换为一个基本解可行的松弛型呢? 答案是引入辅助松弛型（下简称aux）

1
2
3
4
5
最大化		-x0
约束
		  x3=2-2x1+x2+x0
		  x4=-4-x1+5x2+x0
		  xj>=0, j=0,1,2,3,4

因为bi最小值为-4, 所以对aux做pivot(4,0)（这里简写了N、B、A、b、c这5个参数, 仅仅保留了l和e，即基变量和非基变量）操作将anx转换为一个等价的松弛型（下简称aux1）

1
2
3
4
5
最大化		-4-x1+5x2-x4
约束
		  x3=6-x1-4x2+x4
		  x0=4+x1-5x2+x4
		  xj>=0, j=0,1,2,3,4

注意，变化后的aux满足基本解是可行解，为什么pivot(4,0) 就可以将aux变成基本解可行的等价松弛型呢? 这还是要从一般的数学公式出发. 首先 bl<0（不然最初的松弛型就是满足基本解可行的, 我们就不必大费周章的引入啥aux了），其次bl是所有bi们中最小的那一个. 则xl的约束变成

ps: 上面 x_{l+n} 是第l个基变量, 加上n是为了防止和非基变量 x1,..,xn区分开来.

则 -bl>0, 而至于其他的约束变成什么样子了呢?

因为bl是bi中最小的, 所以其他约束也满足新的bi>=0, 所以转换后的aux1必将满足基本解可行. 则下一步就是找出该aux1的最优解, 令aux2是反复对aux1调用pivot最后求出最优解的单纯形状态（即cj系数都<=0, 注意, 因为aux的最优解<=0是肯定的, 所以一定不会出现无界解的情况）. 因为aux和aux1等价, 而aux的目标函数<=0是显然的(因为aux的目标函数是-x0, 而x0>=0)， 如果不是0的话（即严格<0）, 则算法导论的引理29.11表明原松弛型不可行. 即无解.

在aux2中令 x0=0, 就可以将原松弛型转换为基本解可行的松弛型.

好了, 说了这么多, 来看看本题如何解决. 顺便搓一个单纯形法的板子。下面的板子是根据网上的大神改的. 实现上和上面的说法未必完全一致，但是思想是一致的.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
//#include "stdafx.h"
#include <stdio.h>
#include <algorithm>
#include <math.h>
using namespace std;
//#define LOCAL

const int maxn = 50; // n+m的上界
int n,m,type, id[maxn], tp[maxn];
double a[maxn][maxn];
const double eps = 1e-9, inf = 1e9;

inline int kk(bool x) // true返回1, false返回-1
{
	return x?1:-1;
}

void pivot(int r, int c) // 转轴r行c列
{
	swap(id[r+n], id[c]); // 将非基变量id[c]（即伊始的id[c]）变动到基变量id[r+n]去. 因为最初的id[i](i>n)是0，所以id[c]就变成0了,注意, 为什么基变量的编号不重要? 因为时时刻刻目标函数的值仅仅和非基变量有关, 和基变量无关
	double t = a[r][c]; // pivot中要除以的系数
	a[r][c] = 1;
	for (int i = 0;i<=n; i++)
	{
		a[r][i] /= t;
	} // pivot公式第三行, 系数都除以了t（即原本的a[r][c]）
	for (int i = 0; i<=m; i++)
	{
		if (abs(a[i][c])>eps && i!=r) // 如果a[i][c]不为0（如果是0的话, 即pivot公式中的aie为0, 则可以看到,诸多系数都是不变的）并且不是已经处理掉的第r行的话
		{ // 注意, 这里必须要引入math.h 不然abs运行的巨慢
			t = a[i][c];
			a[i][c] = 0; // 注意pivot公式第一行和第二行最后那个系数就知道这里赋值为0的作用了
			for (int j = 0; j<=n; j++)
			{
				a[i][j]+=kk(!i&&!j)*t*a[r][j]; // 注意, a[r][j]是已经除以过 t = a[r][c] 的(pivot公式第三行),对于i=0,j=0的情形,即pivot公式第一行的最左边,它是+而不是-,所以这里要乘以kk
			}
		}
	} // pivot公式的第一行和第二行, 即变换a[i][j]的除了第r行之外的行的系数
} // 这段代码对照文章中的pivot公式来看

void simplex() // 单纯形法, 主体就是2个while, 第一个while是使得单纯形满足基本解可行（这个过程可能发现不可行）, 第二个while就是使得cj们<=0或者发现无界解
{
	int i,j;
	double w,t;
	for (int k = 1; k<=n; k++)
	{
		id[k] = k;
	} // id[i]（1<=i<=n）是非基变量的编号, 例如 id[1]=2, 表示当前处于x1位置的其实是伊始的x2变量
	while(1)
	{
		i=0,j =0; // i是行, j是列
		w = -eps;
		for (int k = 1; k<=m; k++)
		{
			if (a[k][0]<w) // 如果真的<0
			{
				w = a[k][0];
				i = k;
			} // 获取最小的bi
		}
		if (!i)
		{
			break;
		} // 说明bi都非负, 则已经满足基本解可行, 就不需要转换了
		for (int k = 1;k<=n; k++)
		{
			if (a[i][k]<-eps)
			{
				j = k;
				break;
			}
		} // 找到bi所在行严格负的aik
		if (!j) // 如果bi<0并且它所在的一行a系数都>=0,则显然不可行, 因为bi所在的那一行约束是不可能成立的(bi<0, x们都>=0, 则x_{i+n}<0, 而不可能>=0了)
		{
			puts("Infeasible");
			return;
		}
		pivot(i,j); // 转轴第i行和第j列, 目的是使得bi变成正的（从这里也不难想象需要aik<0）
	} // 通过转轴操作使得基本型可行——即bi们都>=0
	while(1)
	{
		i =0, j =0;
		w = eps;
		for (int k = 1; k<=n; k++)
		{
			if (a[0][k]>w)
			{
				w = a[0][k];
				j = k;
			}
		} // 找到cj>0的那一列,而且是cj最大的那一列（这倒无所谓,只是提供一个选择的标准而已）
		if (!j)
		{
			break;
		} // 如果所有cj都<=0了, 则说明已经达到了最优解,单纯形算法就结束了
		w = inf; // 准备找最紧的约束了
		for (int k = 1;k<=m; k++) // 找到cj>0, 则准备增大xj, 但是能不能无限增大就看有没有最紧的约束了
		{
			if (a[k][j]>eps && (t=a[k][0]/a[k][j])<w) // 如果系数真的>0, 则a[k][j]和a[k][0](即bk们)就可以给出约束, 注意, 因为上面的while(1)已经保证了bk们>=0, 所以这里的约束一定是>=0的,不用担心不可行的情况
			{
				w = t;
				i = k;
			}
		}// 找到能给出最紧约束的行i
		if (!i)
		{
			puts("Unbounded");
			return;
		} // 如果a[k][j]都<=0, 则增大xj, 越大越能满足, 则存在无界解
		pivot(i,j); // 转轴操作, 增大目标函数中的v
	}
	printf("%.9f\n", a[0][0]); // 打印最优解
	for (int i = n+1; i<=n+m; i++)
	{
		tp[id[i]] = i-n; // tp[i]=j表示伊始的非基变量xi变成基变量,而这个基变量对应的是bj(从而最后要达到最优解需要xi取值为bj)
	} // id[i]可能是0, 也可能是1,...,n（表示伊始的非基变量变成了基变量）
	if (type) // 如果要打印最优解啥时候取到
	{
		for (int i = 1; i<=n; i++)
		{
			printf("%.9f ", tp[i]?a[tp[i]][0]:0); // 如果tp[i]不是0的话, 则表明xi是映射成了基变量, 否则说明它仍然处于非基变量的位置,注意, 最后伊始的基变量一定都变到非基变量的位置上. 因为最后非基变量都是0（这样,m个约束才都取等号,这样最后取最值的点才在单纯形的角点上）
		}
	}
}

int main()
{
#ifdef LOCAL
	freopen("d:\\data.in", "r", stdin);
	//	freopen("d:\\my.out", "w", stdout);
#endif
	scanf("%d%d%d", &n, &m ,&type);
	for (int i = 1; i<=n; i++)
	{
		scanf("%lf", a[0]+i); // 读入cj系数
	}
	for (int i = 1; i<=m; i++)
	{
		for (int j = 1; j<=n; j++)
		{
			scanf("%lf", a[i]+j); // 读入aij
		}
		scanf("%lf", a[i]); // 读入 bi
	}
	simplex();
	return 0;
}

ac情况

#365822 #179. 线性规划 yfs 97 0ms 600kb C++11 5.1kb 2019-10-02 21:14:39

97分,在一个extra test上被T掉了~

其实从上述单纯形算法不难看出, 单纯形法就是不断的从一个单纯形的角点移动到另一个角点的过程.

参考

【1】https://yfsyfs.github.io/2019/09/20/%E6%B4%9B%E8%B0%B7-3376-%E3%80%90%E6%A8%A1%E6%9D%BF%E3%80%91%E7%BD%91%E7%BB%9C%E6%9C%80%E5%A4%A7%E6%B5%81-Ford-Fulkerson%E7%AE%97%E6%B3%95/

【2】https://www.cnblogs.com/zzqsblog/p/5457091.html

【3】https://yfsyfs.github.io/2019/08/23/poj-3169-Layout-%E6%9C%80%E7%9F%AD%E8%B7%AF/