单纯形算法
作用¶
单纯形法是解决线性规划问题的一个有效的算法。线性规划就是在一组线性约束条件下,求解目标函数最优解的问题。
线性规划的一般形式¶
在约束条件下,寻找目标函数
线性规划的可行域¶
满足线性规划问题约束条件的所有点组成的集合就是线性规划的可行域。若可行域有界(以下主要考虑有界可行域),线性规划问题的目标函数最优解必然在可行域的顶点上达到最优。
单纯形法就是通过设置不同的基向量,经过矩阵的线性变换,求得基可行解(可行域顶点),并判断该解是否最优,否则继续设置另一组基向量,重复执行以上步骤,直到找到最优解。所以,单纯形法的求解过程是一个循环迭代的过程。
线性规划的标准形式¶
在说明单纯形法的原理之前,需要明白线性规划的标准形式。因为单纯形算法是通过线性规划的标准形来求解的。一般,规定线性规划的标准形式为:
写成矩阵形式:
标准形的形式为:
-
1)目标函数要求
\max -
2)约束条件均为等式
-
3)决策变量为非负约束
普通线性规划化为标准形:
- 1)若目标函数为最小化,可以通过取负,求最大化
- 2)约束不等式为小于等于不等式,可以在左端加入非负变量,转变为等式,比如:
同理,约束不等式为大于等于不等式时,可以在左端减去一个非负松弛变量,变为等式。
- 3)若存在取值无约束的变量,可转变为两个非负变量的差,比如:
本文最开始的线性规划问题转化为标准形为:
单纯形法¶
几何意义¶
在标准形中,有
如果令非基变量等于
如果为可行解的话,
还是通过上述具体的线性规划问题来说明:
如果选择
所以,通过选择不同的基变量,可以获得不同的可行域的顶点。
如何判断最优¶
如前所述,基变量可由非基变量表示:
目标函数
当达到最优解时,所有的
当前的目标函数值为
如何选择新的基变量¶
如果存在多个
如何选择被替换的基变量¶
假如我们选择非基变量
继续通过上面的例子来说明:
从最后一行可以看到,
第一行是符号
第二行:若
第三行:若
尽管替换
终止条件¶
当目标函数用非基变量的线性组合表示时,所有的系数均不大于
如果,有一个非基变量的系数为
使用单纯形法来求解线性规划,输入单纯形法的松弛形式,是一个大矩阵,第一行为目标函数的系数,且最后一个数字为当前轴值下的
算法和使用单纯性表求解线性规划相同。
对于线性规划问题:
我们可以得到其松弛形式:
我们可以构造单纯形表,其中最后一行打星的列为轴值。
在单纯形表中,我们发现非轴值的
其实我们由于每个
我们可以发现,对于约束方程
继续计算,我们得到:
此时我们发现,所有非轴的
整个过程代码如下所示:
#include <bits/stdc++.h>
using namespace std;
vector<vector<double> > Matrix;
double Z;
set<int> P;
size_t cn, bn;
bool Pivot(pair<size_t, size_t> &p) { // 返回0表示所有的非轴元素都小于0
int x = 0, y = 0;
double cmax = -INT_MAX;
vector<double> C = Matrix[0];
vector<double> B;
for (size_t i = 0; i < bn; i++) {
B.push_back(Matrix[i][cn - 1]);
}
for (size_t i = 0; i < C.size(); i++) { // 在非轴元素中找最大的c
if (cmax < C[i] && P.find(i) == P.end()) {
cmax = C[i];
y = i;
}
}
if (cmax < 0) {
return 0;
}
double bmin = INT_MAX;
for (size_t i = 1; i < bn; i++) {
double tmp = B[i] / Matrix[i][y];
if (Matrix[i][y] != 0 && bmin > tmp) {
bmin = tmp;
x = i;
}
}
p = make_pair(x, y);
for (set<int>::iterator it = P.begin(); it != P.end(); it++) {
if (Matrix[x][*it] != 0) {
// cout<<"erase "<<*it<<endl;
P.erase(*it);
break;
}
}
P.insert(y);
// cout<<"add "<<y<<endl;
return true;
}
void pnt() {
for (size_t i = 0; i < Matrix.size(); i++) {
for (size_t j = 0; j < Matrix[0].size(); j++) {
cout << Matrix[i][j] << "\t";
}
cout << endl;
}
cout << "result z:" << -Matrix[0][cn - 1] << endl;
}
void Gaussian(pair<size_t, size_t> p) { // 行变换
size_t x = p.first;
size_t y = p.second;
double norm = Matrix[x][y];
for (size_t i = 0; i < cn; i++) { // 主行归一化
Matrix[x][i] /= norm;
}
for (size_t i = 0; i < bn; i++) {
if (i != x && Matrix[i][y] != 0) {
double tmpnorm = Matrix[i][y];
for (size_t j = 0; j < cn; j++) {
Matrix[i][j] = Matrix[i][j] - tmpnorm * Matrix[x][j];
}
}
}
}
void solve() {
pair<size_t, size_t> t;
while (1) {
pnt();
if (Pivot(t) == 0) {
return;
}
cout << t.first << " " << t.second << endl;
for (set<int>::iterator it = P.begin(); it != P.end(); it++) {
cout << *it << " ";
}
cout << endl;
Gaussian(t);
}
}
int main(int argc, char *argv[]) {
// ifstream fin;
// fin.open("./test");
cin >> cn >> bn;
for (size_t i = 0; i < bn; i++) {
vector<double> vectmp;
for (size_t j = 0; j < cn; j++) {
double tmp = 0;
cin >> tmp;
vectmp.push_back(tmp);
}
Matrix.push_back(vectmp);
}
for (size_t i = 0; i < bn - 1; i++) {
P.insert(cn - i - 2);
}
solve();
}
/////////////////////////////////////
// glpk input:
///* Variables */
// var x1 >= 0;
// var x2 >= 0;
// var x3 >= 0;
///* Object function */
// maximize z: x1 + 14*x2 + 6*x3;
///* Constrains */
// s.t. con1: x1 + x2 + x3 <= 4;
// s.t. con2: x1 <= 2;
// s.t. con3: x3 <= 3;
// s.t. con4: 3*x2 + x3 <= 6;
// end;
/////////////////////////////////////
// myinput:
/*
8 5
1 14 6 0 0 0 0 0
1 1 1 1 0 0 0 4
1 0 0 0 1 0 0 2
0 0 1 0 0 1 0 3
0 3 1 0 0 0 1 6
*/
/////////////////////////////////////
理论罗列¶
标准型¶
最大化
最大化
转换为标准型¶
若目标函数要求取最小值,那么可以对其取相反数变成取最大值。对于限制条件
松弛型¶
基本变量
非基变量
松弛变量
等式左侧为基本变量,右侧为非基本变量。
变量¶
- 替入变量
x_e - 替出变量
x_l
可行解¶
-
基本解:所有非基变量设为
0 -
基本可行解:所有
b_i \geq 0
注:单纯形法的过程中 B 和 N 不断交换,在 n 维空间中不断走,“相当于不等式上的高斯消元”。
转轴¶
选取一个非基本变量
Bland 规则 可以参看:最优化方法
初始化¶
在所有
算法实现¶
每个约束定义了
以下问题可以转换为单纯形:
- 最短路
- 最大流
- 最小费用最大流
- 多商品流
基本思想就是改写
转动时,
也就是说,约束条件只用
simplex
是主过程,基本思想是找到一个
例题「NOI2008」志愿者招募
题目大意:长度为
原始问题
对偶问题
把对应出的系数矩阵代入到单纯形算法就可以求出最优解了。
#include <algorithm>
#include <cmath>
#include <cstdio>
#include <cstring>
#include <iostream>
using namespace std;
typedef long long ll;
const int M = 10005, N = 1005, INF = 1e9;
const double eps = 1e-6;
inline 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] > eps) break;
if (e == n + 1) return v;
double mn = INF;
for (int i = 1; i <= m; i++)
if (a[i][e] > eps && mn > b[i] / a[i][e]) mn = b[i] / a[i][e], l = i;
if (mn == INF) return INF; // unbounded
pivot(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;
b[i] = read();
}
printf("%d", (int)(simplex() + 0.5));
}
对偶原理¶
最大化与最小化互换,常数与目标函数互换,改变不等号,变量与约束对应。
令
全幺模矩阵(Totally Unimodular Matrix)¶
充分条件:
-
仅有
-1,0,1 -
每列至多两个非零数
-
行可分为两个集合:
- 一列包含两个同号非零数,两行不在同一个集合
- 一列包含两个异号非零数,两行在同一个集合
线性规划中