线性规划与单纯形

线性规划与单纯形

0x00 前言

• 二轮前心有点乱，做不动题，来写写算法笔记调整一下心情
• 你可能需要：
• 矩阵基础
• 高斯消元基础
• 高中线性规划基础
• 来理解这个东西

0x01 简介

• 线性规划描述这样一类问题：有一些变量 $ x_{i} $ ，有一个目标函数 $ C = sum c_{i}x_{i} $， 有一些限制 $ sum a_{i,j}x_{j} leqslant or geqslant b_{i} $，我们要最大化/最小化C
• 比方说，你要找一个npy，npy有一个心情值和对你的喜爱程度，有2种物品，你给ta买一个物品i会花掉你$ c_{i} $的钱，增加npy $ a_{i,2} $ 的心情值，$ a_{i,2} $ 的喜爱程度，问你最少要花多少钱才能让ta的心情值>=b1，对你的喜爱程度>=b2，为了便于理解，这里的变量和上面的定义是对应的
• 这个问题就相当于，设你买了x1个物品1，x2个物品2，要求在

[a_{1,1}x_{1}+a_{2,1}x_{2} geqslant b_{1} ]

[a_{1,2}x_{1}+a_{2,2}x_{2} geqslant b_{2} ]

的前提下，最小化

[C = c_{1}x_{1}+c_{2}x_{2} ]

• 这就是一个线性规划的典型形式

0x02 简化版问题

• 我知道你们都学过简化版线性规划（平面版的），所以我就不讲那个愚蠢的用一条直线去切可行解区域的方法了。

0x03 松弛型与标准型

• 线性规划有两种规范的形式，即松弛型与标准型
• 这一条讲给有一定基础的人：我这里的标准型和松弛型可能和大学课本里的定义不一样，我以我的为准
• 先讲一下标准型，所谓的标准型就是在刚才的定义的基础上

有一些变量 $ x_{i} $ ，有一个目标函数 $ C = sum c_{i}x_{i} $， 有一些限制 $ sum a_{i,j}x_{j} leqslant or geqslant b_{i} $，我们要最大化/最小化C

• 加一些更严格的限制条件，变成这样

有一些变量 $ x_{i} geqslant 0 $ ，有一个目标函数 $ C = sum c_{i}x_{i} $， 有一些限制 $ sum a_{i,j}x_{j} leqslant b_{i} $，我们要最大化C

• 事实上，把普通形式转化为标准型是不难的，有几个问题，我们分别讨论怎么解决

1. $ x_{i} $ 可能没有限制
   把 $ x_{i} $ 全部拆成 $ x_{i} = x_{i+1}-x_{i+2} $，然后 $ x_{i+1},x_{i+2} geqslant 0 $ 稍微思考一下就发现这很正确
2. 有 $ x_{i} < 0 $ 这种坑爹限制
   把 $ x_{i} $ 全部换成 $ -x_{i} $
3. 可能限制是 $ sum a_{i,j}x_{j} geqslant b_{i} $ 的形式
   不等式整个乘-1取反
4. 可能要最小化C
   取反取反！
5. 限制可能是小于或者大于，没有等于
   分情况讨论，如果是浮点数，两者没有区别
   如果是整数，直接+1/-1即可

• 然后我们就得到了标准形式，然并卵
• 我们还要讨论一下松弛形式
• 如果你学过高斯消元，你会知道，n个线性无关的方程可以解出n个变量，但是n个不等式啥都干不了
• 但是我们可以给它加上一些松弛变量，然后让不等式变成等式
• 显然

[a_{1,1}x_{1}+a_{2,1}x_{2} leqslant b_{1} ]

[a_{1,1}x_{1}+a_{2,1}x_{2}+x_{3} = b_{1} ]

• 是等价的，这里就可以看出为啥要统一成小于等于了，因为还有 $ x_{i} geqslant 0 $ 的性质，在小于等于的基础上 $ +x_{3} $ 才是不违反限制的
• 然后就都变成等式啦！

0x04 基解与可行解

• 但是这里有问题了，假如一共有n个变量，m个限制，因为松弛的需要，又额外加了m个变量，于是现在又n+m个变量，m个方程，根据克拉玛定理（或者说，显然），这个方程不止一组解，我们不妨先弄一组解出来

• 我们来考虑一下系数 $ a_{i,j} $ 的矩阵，就用刚才那个npy问题的矩阵好了

[egin{bmatrix} a_{1,1} & a_{2,1} & 1 & 0 \ a_{1,2} & a_{2,2} & 0 & 1 end{bmatrix} ]

• 注意这个是加了松弛变量以后的矩阵
• 发现前面那一堆a不太好弄，指不准行列式就是0，但是后面有一个I，好解&行列式不会为0，那么$ x1 = 0,x2 = 0,x3 = b1,x4 = b2 $ 是一组很不错的解
• 然后呢？然后这个解没什么卵用（C是0），而且可能不可行（b1，b2可能小于0）
• 这样的一个解里有m个变量不是0（叫做基变量），n个变量是0（叫做非基变量）
• 我们把任何一组这样的解叫做基解，如果它可行，则称作基可行解
• 不妨来考虑一下这个东西的几何意义
• 你已经知道了二维平面上的最优解总是在凸包顶点上（或者会有无穷多最优解）
• 事实上，这相当于在一个n+m维空间里切n+m维的超多面体，理解一下这个，这很重要
• 脑洞一下就能发现，最优解也类似的，肯定在超多面体的顶点上
• 而我们解出来的基解一定是在顶点上的
• 通过一些步骤我们可以让解在多面体上移动，最终跑到最优解
• 好，几何到此为止，因为再下去就不可言传了，你得自己体会，而我越讲你越乱

0x05 单纯形法

• 假设有这样一个线性规划
• 在这里为了简便，我们假设所有 $ b_{i} geqslant 0 $，小于0的后面会讲到
• 为了易于理解与书写，我们把变量设成具体的数字

[C = x_{1}+2x_{2} ]

[2x_{1}+x_{2} leqslant 2 ]

[-x_{1}+x_{2} leqslant 3 ]

• 转化成松弛型：

[C = x_{1}+2x_{2} ]

[2x_{1}+x_{2} + x_{3} = 2 ]

[-x_{1}+x_{2} + x_{4} = 3 ]

• 写出系数矩阵：

[egin{bmatrix} 2 & 1 & 1 & 0 \ -1 & 1 & 0 & 1 end{bmatrix} ]

• 把初始解搞出来

[left{egin{matrix} x_{1}=0 \ x_{2}=0 \ x_{3}=2 \ x_{4}=3 end{matrix} ight. ]

• 然后我们发现，既然现在的基变量的系数全是0，我们可以变相的认为改变基变量不会对C产生影响，为什么我们不从C的表达式中选一个系数大于0的变量增大呢？
• 现在有这样一些等式

[x_{3} = 2 - 2x_{1} - x_{2} ]

[x_{4} = 3 + x_{1} + x_{2} ]

• 如果我们选择 $ x_{1} $ 并增大，那么 $ x_{3} $ 与 $ x_{4} $ 也会相应的变化，对 $ x_{1} $ 变化的幅度产生限制，如果 $ x_{1} $ 在允许的范围内尽可能变大，那么最终会导致 $ x_{3} $ 与 $ x_{4} $ 中的一个变成0
• 那么如果 $ x_{1} $ 的变化没有限制呢？那答案就是正无穷了
• 现在，我们发现由于 $ x_{3} $ 的限制，$ x_{1} $只能变化到1就会卡住（不然 $ x_{3} $ 就小于0了）
• 然后我们把所有式子中的 $ x_{1} $ 通过 $ x_{3} = 2 - 2x_{1} - x_{2} $ 给取代掉
• 首先把 $ x_{1} $ 的系数消成1，变形成

[x_{1} = 1 - frac{1}{2}x_{3} - frac{1}{2}x_{2} ]

• 然后把所有的 $ x_{1} $ 换掉

[C = 1 - frac{1}{2}x_{3} - frac{3}{2} x_{2} ]

[x_{1} = 1 - frac{1}{2}x_{3} - frac{1}{2}x_{2} ]

[x_{4} = 4 - frac{1}{2}x_{3} + frac{1}{2}x_{2} ]

• 现在的解是：

[left{egin{matrix} x_{1}=1 \ x_{2}=0 \ x_{3}=0 \ x_{4}=4 end{matrix} ight. ]

• 可以发现，现在的状况和刚才是一样的：所有非基变量都是0，所有基变量在C中的系数都是0，基变量的系数构成一个I，但是C的常数项变大了，我们可以一直这么做下去，直到所有的C中变量系数都小于0为止
• 这个就是单纯形法

0x06 上代码！##

#include <iostream>
#include <cstdlib>
#include <cstdio>
#include <algorithm>
#include <string>
#include <cstring>
#include <cmath>
#include <map>
#include <stack>
#include <set>
#include <vector>
#include <queue>
#include <time.h>
#define eps 1e-10
#define INF 1e15
#define MOD 1000000007
#define rep0(j,n) for(int j=0;j<n;++j)
#define rep1(j,n) for(int j=1;j<=n;++j)
#define pb push_back
#define set0(n) memset(n,0,sizeof(n))
#define ll long long
#define OJ UOJ
using namespace std;
int read() {
	char c = getchar(); int f = 1, x = 0;
	while (!isdigit(c)) { if (c == '-') f = -1; c = getchar(); }
	while (isdigit(c)) { x = x * 10 + c - '0'; c = getchar(); }
	return f*x;
}
char get_ch() {
	char c = getchar();
	while (!isalpha(c)) c = getchar();
	return c;
}
const int MAXINT = 25;
int n, m, t;
struct Linear_Programming {
	double a[MAXINT][MAXINT]; //m*n的矩阵，[1,1]到[m,n]是系数矩阵，[0,1]到[0,n]是c，[1,0]到[m,0]是b，[0,0]是目标函数值
	int id[MAXINT * 2];
	double ans[MAXINT];//1~n是非基变量，n+1~n+m是基变量
	void pivot(int out, int in) {
		swap(id[n + out], id[in]);//把非基变量丢进去，基变量丢出来
		double tmp = a[out][in];
		a[out][in] = 1;//把基变量那一列换出来
		for (int i = 0; i <= n; i++) a[out][i] /= tmp;//把需要换进去的变量系数消成1
		for (int i = 0; i <= m; i++) if (i != out && fabs(a[i][in]) > eps) {//处理系数矩阵，消元，事实上把c和目标函数一起处理了
			tmp = a[i][in]; a[i][in] = 0;
			for (int j = 0; j <= n; j++) a[i][j] -= tmp*a[out][j];
		}
	}
	bool init_solution() {
		rep1(i, n) id[i] = i;
		while (1) {
			int in = 0, out = 0;
			for (int i = 1; i <= m; i++) if (a[i][0] < -eps && (!out || (rand() & 1)))out = i; //初始解不合法，随机选一个不合法的基变量换出来
			if (!out) break; //合法
			for (int i = 1; i <= n; i++) if (a[out][i] < -eps && (!in || (rand() & 1))) in = i;//把一个负系数的变量换进去来尝试使初始解合法
			if (!in) {//没有办法使初始解合法了
				printf("Infeasible
");
				return false;
			}
			pivot(out, in);
		}
		return true;
	}
	bool simplex() {//求解单纯形
		while (1) {
			int in = 0, out = 0;
			double mn = INF;
			for (int i = 1; i <= n; i++) if (a[0][i] > eps) { in = i; break; }//选一个系数大于0的丢进去
			if (!in) break;//没有系数大于0的，已经达到最优解
			for (int i = 1; i <= m; i++) 
				if (a[i][in] > eps && a[i][0] / a[i][in] < mn) {
					mn = a[i][0] / a[i][in]; out = i;
				}
			if (!out) {//没有任何限制，无界解
				printf("Unbounded
");
				return false;
			}
			pivot(out, in);
		}
		return true;
	}
	void work() {
		if (init_solution() && simplex()) {
			printf("%.8lf
", -a[0][0]);//注意一下算出来的解是真正答案的相反数
			if (t == 1) {
				rep1(i, m) ans[id[n + i]] = a[i][0];
				rep1(i, n) printf("%.8lf ", ans[i]);
				putchar('
');
			}
		}
	}
}lp;
int main() {
	srand(559320414);
	n = read(); m = read(); t = read();
	rep1(i, n) lp.a[0][i] = read();
	rep1(i, m) {
		rep1(j, n) {
			lp.a[i][j] = read();
		}
		lp.a[i][0] = read();
	}
	lp.work();
	return 0;
}
//uoj的模板题
//事实上只需要记录非基变量，反正编号为i的基变量的系数向量除了i这个位置是1其他都是0，c[i]肯定为0，根本不用记录
//求解的时候要统一成sigma axi <= b的形式，而且是最大化目标函数，可以用取反或者对偶原理解决
//对偶原理：把c向量和b向量互换，a矩阵转置（如果像本程序用一个大矩阵记录的话就是把这整个矩阵转置），最大化变成最小化（最小化变成最大化），最终的答案不变

• 这里使用了一些trick，a[1][1]开始的矩阵记录的才是系数急诊，a[0][i]记录的是C中的系数，a[i][0]记录的是b，a[0][0]记录的是C的相反数

0x07 大杀器

• 随机初始解：如果初始有一些b小于0，随机选个负系数变量交换进去
• 对偶原理：把整个系数矩阵转置，约束和C中的系数互换，最大化C与最小化C互换，答案不变
• 具体原因不知

0x08 什么时候能用这玩意？

• 一般只能做模板题
• 能做所有网络流题
• 不能做要求 $ x_{i} $ 是整数的题，除非系数矩阵是全幺模矩阵（网络流题的矩阵都是全幺模矩阵）