8.16.2 ACM-ICPC 线性规划 单纯形算法
简介
单纯形算法(Simplex Algorithm)是乔治·丹齐格(George Dantzig)于1947年提出的一种用于求解线性规划问题的算法。该算法在优化领域有着广泛的应用,并且由于其高效性和实际应用中的优秀表现,成为解决线性规划问题的经典方法。
线性规划问题的基本形式
一个标准的线性规划问题可以表述为:
最大化或最小化目标函数:
其中,Z是目标函数,cic_ici是目标函数的系数,aija_{ij}aij是约束条件的系数,bib_ibi是约束条件的右端项。
单纯形算法的基本概念
基本解与可行解
基本解是通过将一些变量设置为零并求解其他变量得到的解。可行解是满足所有约束条件的基本解。
基本可行解与极点
基本可行解是满足所有约束条件且非负的解,对应于几何意义上的一个极点。
目标函数的改进
单纯形算法通过从一个基本可行解移动到另一个基本可行解,逐步改进目标函数的值,直到找到最优解。
单纯形算法的步骤
1. 初始化
找到一个初始的基本可行解。通常通过引入松弛变量将不等式约束转换为等式约束来实现。
2. 检查最优性
检查当前解是否满足最优条件。如果所有目标函数系数均非负,则当前解为最优解。
3. 选择进入变量
选择一个使目标函数值增加最快的变量进入基。如果存在多个变量均能使目标函数值增加,则任意选择一个。
4. 选择离开变量
确定哪个变量将离开基。通过计算每个基本变量在新解中变为零时的约束,选择最小的非负比值所对应的变量作为离开变量。
5. 更新基
更新基变量和非基变量,重新计算新的基本可行解。
6. 重复步骤2至5
重复上述步骤,直到找到最优解或确定问题无界。
单纯形算法的实际应用
单纯形算法广泛应用于各种优化问题,如资源分配、生产计划、运输问题等。其高效性和鲁棒性使其成为求解线性规划问题的常用方法。
示例应用
假设我们有一个生产计划问题:
目标函数:
约束条件:
通过单纯形算法,我们可以逐步求解,最终找到最优解 x1x_1x1 和 x2x_2x2 的值,使目标函数 ZZZ 最大化。
结论
单纯形算法是一种强大且高效的线性规划求解方法。通过不断改进基本可行解,它能够在多种应用中找到最优解。理解单纯形算法的基本步骤和概念对于解决实际问题具有重要意义。
在实际应用中,单纯形算法虽然在最坏情况下可能表现出指数时间复杂度,但其在大多数实际问题中的表现非常优异。通过对算法的深入理解和实际应用,可以更好地解决线性规划问题 (ACM Digital Library) (UCI Bren School) (Codeforces)。
8.16.2 ACM-ICPC 线性规划 单纯形算法
定义
单纯形法是解决线性规划问题的一个有效算法。线性规划就是在一组线性约束条件下,求解目标函数最优解的问题。
线性规划的一般形式
在约束条件下,寻找目标函数 zzz 的最大值:
线性规划的可行域
满足线性规划问题约束条件的所有点组成的集合就是线性规划的可行域。若可行域有界(以下主要考虑有界可行域),线性规划问题的目标函数最优解必然在可行域的顶点上达到最优。
单纯形法的思想
单纯形法通过设置不同的基向量,经过矩阵的线性变换,求得基可行解(可行域顶点),并判断该解是否最优,否则继续设置另一组基向量,重复执行以上步骤,直到找到最优解。所以,单纯形法的求解过程是一个循环迭代的过程。
可行域示意图
线性规划的标准形式
在说明单纯形法的原理之前,需要明白线性规划的标准形式。因为单纯形算法是通过线性规划的标准形式来求解的。一般,规定线性规划的标准形式为:
写成矩阵形式:
标准形式的特征:
- 目标函数要求 max\maxmax
- 约束条件均为等式
- 决策变量为非负约束
普通线性规划化为标准形式
若目标函数为最小化,可以通过取负求最大化。约束不等式为小于等于不等式,可以在左端加入非负变量,转变为等式。
例如:
同理,约束不等式为大于等于不等式时,可以在左端减去一个非负松弛变量,变为等式。若存在取值无约束的变量,可转变为两个非负变量的差。
例如:
本文最开始的线性规划问题转化为标准形式为:
单纯形法的几何意义与例子
在标准形式中,有 mmm 个约束条件(不包括非负约束),nnn 个决策变量,且 n≥mn \geq mn≥m。首先选取 mmm 个基变量 xj′(j=1,2,…,m) x_j^{'}(j = 1, 2, \ldots, m )xj′(j=1,2,…,m),基变量对应约束系数矩阵的列向量线性无关。通过矩阵的线性变换,基变量可由非基变量表示:
如果为可行解的话,CiC_iCi 大于 0。几何意义就是在具体的线性规划问题中,通过选择不同的基变量,可以获得不同的可行域顶点。
例如:
如果选择 x2x_2x2、x3x_3x3 为基变量,那么令 x1x_1x1、x4x_4x4 等于 0,可以去求解基变量 x2x_2x2、x3x_3x3 的值。对系数矩阵做行变换如下所示,
X1=0X_1=0X1=0 表示可行解在 y 轴上;X4=0X_4=0X4=0 表示可行解在 x1+2x2=9x_1+2x_2=9x1+2x2=9 的直线上。那么,求得的可行解即表示这两条直线的交点,也是可行域的顶点。
如何判断最优解
基变量可由非基变量表示:
目标函数 zzz 也可以完全由非基变量表示:
当达到最优解时,所有的 σj\sigma_jσj 应小于等于 0,当存在 jjj,σj>0\sigma_j > 0σj>0 时,当前解不是最优解。
如何选择新的基变量
如果存在多个 σj>0\sigma_j > 0σj>0,选择最大的 σj>0\sigma_j > 0σj>0 对应的变量作为基变量,这表示目标函数随着 xj′x_j^{'}xj′ 的增加,增长最快。
如何选择被替换的基变量
假如选择非基变量 xs′x_s^{'}xs′ 作为下一轮的基变量,那么被替换基变量 xj′x_j^{'}xj′ 在下一轮中作为非基变量,等于 0。选择 xj′x_j^{'}xj′ 的原则是替换后尽量使 xs′x_s^{'}xs′ 值最大,但要保证替换基变量后的解仍是可行解,因此应该选择最紧的限制。
终止条件
当目标函数用非基变量的线性组合表示时,所有的系数均不大于 0,则表示目标函数达到最优。
如果有一个非基变量的系数为 0,其他的均小于 0,表示目标函数的最优解有无穷多个。这是因为目标函数的梯度与某一边界正交,在这个边界上,目标函数的取值均相等,且为最优。
单纯形法的具体实现
标准型:m+n 个约束 n 个变量用 x 向量表示,A 是一个 m×n 的矩阵,c 是一个 n 的向量,b 是一个 m 的向量,最大化 cx 满足约束 Ax ≤ b,x > 0。
最大化:
满足如下约束条件:
松弛型:基本变量 B,|B|=m,一个约束对应一个,表示松弛量,叫做松弛变量(基本变量)。
非基变量 N,|N|=n,x_n + i = b_i - \sum a_{ij}x_j \geq 0
松弛变量 xn+ix_{n+i}xn+i
等式左侧为基本变量,右侧为非基本变量。
变量的替换
替入变量 xex_exe(非基变量)
替出变量 xlx_lxl(基本变量)
可行解
基本解:所有非基变量设为 0,基本变量为右侧的常数。
基本可行解:所有 bi≥0b_i \geq 0bi≥0。
注:单纯形法的过程中 B 和 N 不断交换,在 n 维空间中不断走,「相当于不等式上的高斯消元」。
转轴
选取一个非基本变量 xex_exe 为替入变量,基本变量 xlx_lxl 为替出变量,将其互换,为了防止循环,根据 Bland 规则,选择下标最小的变量。
Bland 规则:可以参看最优化方法。
初始化
在所有 bi<0b_i < 0bi<0 的约束中随机选一个作为 xlx_lxl,再随机选一个 ale<0a_{le} < 0ale<0 作为 xex_exe,然后 pivot(l,e) 后 bib_ibi 就变正了。如果不存在这样的 xlx_lxl,那么说明原问题无解。
算法实现
每个约束定义了 n 维空间中的一个半空间(超平面),交集形成的可行域是一个凸区域称为单纯形。目标函数是一个超平面,最优解在凸区域定点处取得。通过不断的转轴操作,在 n 维凸区域的顶点上不断移动(转轴),使得基本解的目标值不断变大,最终达到最优解。如果基本解的目标值没有上限,那么说明原问题无界。
例题
例如,NOI2008志愿者招募问题。
代码实现
#include <algorithm>
#include <cmath>
#include <cstdio>
#include <cstring>
#include <iostream>
using namespace std;
const int M = 10005, N = 1005, INF = 1e9;
int read() { // 快读
char c = getchar();
int x = 0, f = 1;
while (c < '0' || c > '9') {
if (c == '-') f = -1;
c = getchar();
}
while (c >= '0' && c <= '9') {
x = x * 10 + c - '0';
c = getchar();
}
return x * f;
}
int n, m;
double a[M][N], b[M], c[N], v;
void pivot(int l, int e) { // 转轴操作函数
b[l] /= a[l][e];
for (int j = 1; j <= n; j++)
if (j != e) a[l][j] /= a[l][e];
a[l][e] = 1 / a[l][e];
for (int i = 1; i <= m; i++)
if (i != l && fabs(a[i][e]) > 0) {
b[i] -= a[i][e] * b[l];
for (int j = 1; j <= n; j++)
if (j != e) a[i][j] -= a[i][e] * a[l][j];
a[i][e] = -a[i][e] * a[l][e];
}
v += c[e] * b[l];
for (int j = 1; j <= n; j++)
if (j != e) c[j] -= c[e] * a[l][j];
c[e] = -c[e] * a[l][e];
// swap(B[l],N[e])
}
double simplex() {
while (true) {
int e = 0, l = 0;
for (e = 1; e <= n; e++)
if (c[e] > (double)0) break;
if (e == n + 1) return v; // 此时v即为最优解
double mn = INF;
for (int i = 1; i <= m; i++) {
if (a[i][e] > (double)0 && mn > b[i] / a[i][e]) {
mn = b[i] / a[i][e]; // 找对这个e限制最紧的l
l = i;
}
}
if (mn == INF) return INF; // unbounded
pivot(l, e); // 转动l,e
}
}
int main() {
n = read();
m = read();
for (int i = 1; i <= n; i++) c[i] = read();
for (int i = 1; i <= m; i++) {
int s = read(), t = read();
for (int j = s; j <= t; j++) a[i][j] = 1; // 表示第i种志愿者在j时间可以服务
b[i] = read();
}
printf("%d", (int)(simplex() + 0.5));
}
对偶原理
最大化与最小化互换,常数与目标函数互换,改变不等号,变量与约束对应。
全幺模矩阵(Totally Unimodular Matrix)
当一个矩阵的任意一个子方阵的行列式都为 ±1,0\pm1,0±1,0 时,我们称这个矩阵是全幺模的。如果单纯形矩阵是全幺模的,那么单纯形就具有整数解。
证明:在线性规划中,最优解通常位于由约束条件形成的可行域的边界上。具体来说,它位于这个高维多面体的顶点处。这些顶点可以通过将部分线性独立的不等式约束转化为等式,然后求解得到的线性方程组来确定。如果约束矩阵是全幺模的(即每个方阵子矩阵的行列式为 0、1 或 -1),那么根据克莱姆法则(Cramer's rule),这些线性方程组的解将是整数。
线性规划中 A 为全幺模矩阵,则单纯形法过程中所有系数 ∈−1,0,1\in -1,0,1∈−1,0,1,可以去除系数为 0 的项进行优化!
注:任何最大流、最小费用最大流的线性规划都是全幺模矩阵。
习题练习
UOJ#179. 线性规划
参考资料
- 线性规划之单纯形法【超详解 + 图解】
- 2016 国家集训队论文
- 算法导论