线性规划与单纯形
线性规划与单纯形
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
- 事实上，把普通形式转化为标准型是不难的，有几个问题，我们分别讨论怎么解决
- $ x_{i} $ 可能没有限制
把 $ x_{i} $ 全部拆成 $ x_{i} = x_{i+1}-x_{i+2} $，然后 $ x_{i+1},x_{i+2} \geqslant 0 $ 稍微思考一下就发现这很正确 - 有 $ x_{i} < 0 $ 这种坑爹限制
把 $ x_{i} $ 全部换成 $ -x_{i} $ - 可能限制是 $ \sum a_{i,j}x_{j} \geqslant b_{i} $ 的形式
不等式整个乘-1取反 - 可能要最小化C
取反取反！ - 限制可能是小于或者大于，没有等于
分情况讨论，如果是浮点数，两者没有区别
如果是整数，直接+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问题的矩阵好了
\[\begin{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
\]
- 写出系数矩阵：
\[\begin{bmatrix}
2 & 1 & 1 & 0 \\
-1 & 1 & 0 & 1
\end{bmatrix}
\]
- 把初始解搞出来
\[\left\{\begin{matrix}
x_{1}=0
\\ x_{2}=0
\\ x_{3}=2
\\ x_{4}=3
\end{matrix}\right.
\]
- 然后我们发现，既然现在的基变量的系数全是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\{\begin{matrix}
x_{1}=1
\\ x_{2}=0
\\ x_{3}=0
\\ x_{4}=4
\end{matrix}\right.
\]
- 可以发现，现在的状况和刚才是一样的：所有非基变量都是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\n");
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\n");
return false;
}
pivot(out, in);
}
return true;
}
void work() {
if (init_solution() && simplex()) {
printf("%.8lf\n", -a[0][0]);//注意一下算出来的解是真正答案的相反数
if (t == 1) {
rep1(i, m) ans[id[n + i]] = a[i][0];
rep1(i, n) printf("%.8lf ", ans[i]);
putchar('\n');
}
}
}
}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} $ 是整数的题，除非系数矩阵是全幺模矩阵（网络流题的矩阵都是全幺模矩阵）
许愿弥生改二
