线性规划

单纯形

标准型与松弛型

标准型是

$$
\sum_i a_{i, j}x_j\le b_i\
\max \sum_i c_ix_i
$$

松弛型是

$$
\sum_i a_{i, j}x_j=b_i\
\max \sum_i c_ix_i
$$

标准型转化为松弛型, 只要在所有限制左侧一个松弛变量即可. $\sum_i a_{i, j}x_j + x\le b_i$ .

一般写的单纯型求解的应该是标准型转化成的松弛型, (当然, 也有可能因为UOJ上模板题是输入标准型)

核心思想

每一条不等式限制限制的是多维空间的一个”半空间”(类比二维空间的半平面), 若干个这个的交应该是一个凸多维几何体.

而因为目标函数是线性的, 所以必然有一个顶点可以取到最优解. (可以理解一下, 线性函数等值线全部平行, 也就是朝着一个固定的方向函数值最大, 可以沿着等值线平移到一个顶点上). 因为是凸的所以只要贪心的走就能走到最后的最优解.

单纯形说的是, 选择一组线性无关的变量, 为基变量, 其他变量都可以由基变量表示.

单纯形实现中, 我们让基变量形式上都是松弛变量(对于每一个基变量, 其只在一个限制中系数为1的出现, 其他限制中系数都为0), 于是任意时刻每个基变量对应一个限制.

把这样就可以把它们都设成0, 此时基变量的取值也是确定的–它们对应的限制中的常数, 因为其他变量都成0了. 此时对应了凸函数的一个顶点. 于是不断切换基变量就可以到不同的顶点.

因为最后目标函数最大, 所以我们希望最后所有非基变量在目标函数中的系数都是负的, 这样就可以都设成0了.

于是算法过程是我们通过不断切换基变量在多维几何体的顶点上游走.

如何切换基变量-Pivot

我们每次选择一个非基变量替换已有的一个基变量.

考虑选择什么非基变量最优, 也就是让哪个变大最优(非基变量取值都是0), 我们要求 $\max$ , 所以贪心选择目标函数中系数最大的.

然后要替换掉一个基变量, 考虑非基变量最终取值由至少一个限制限制住, 那么就让这个限制的基变量出来换成它(因为最后它的取值就是这个限制的常数)

那么因为形式上基变量系数只在对应限制上是1其他为0, 于是先给限制整个除一下系数变成1, 然后再用这一限制消去其他所有对应限制中的这个基变量即可.

具体流程

初始化

要先找到一个可行解, 再进行.

考虑如果所有变量都是0不是一个可行解, 那么一定因为至少一个 $b_i<0$ .

对于一个 $b_i<0$ , 只要在对应限制里找到一个系数为负的变量pivot变成基变量就能让 $b_i>0$ . 如果没有系数为负的, 那么一定无解.

最优化

每次按照上面写的寻找非基变量和对应限制循环pivot即可. 终止条件是目标函数系数全不正(非基变量当系数不正, 基变量系数按照定义全0)

输出

发现, 基变量的取值是对应限制的常数.

非基变量取值是0.

答案是目标函数的常数.

实现(UOJ板子)

输入一个标准形, 求最大值和变量取值.

看了uoj上的实现, 前 $k$ 快是一个模子的, 然而很费解. 比如读入进来的明明是标准形式却给所有系数取负等.

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

#include <iomanip>
#include <iostream>
using namespace std;

//uoj卡精度部分
#define float __float128
const float eps = 1e-16, inf = 1e9;
int sgn(float x) {
    return x < -eps ? -1 : (x > eps ? 1 : 0);
}

const int N = 50;   //N是n+m
float a[N][N];      //系数矩阵,a[i][j]表示变量id[j]在第i个限制中的系数. a[0][i]表示目标函数, a[0][0]是目标函数常数,a[i][0]是第i个限制的b_i常数
int id[N];          //id[1...n]是系数矩阵第i列对应的变量编号,id[n+1...m]是第i个限制的基变量编号
int n, m;           //变量个数(不包含转化过去的松弛变量)

//代码理解困难可能主要在于a[...][i]表示的变量是在变的.

//把变量id[idx]设为行l当基变量,(用行l消掉其他行对应id[idx]的系数)
void pivot(int l, int idx) {
    swap(id[n + l], id[idx]);
    
    float c = a[l][idx];
    for (int i = 0; i <= n; i++)
        a[l][i] /= c;

    //a[...][idx]列要从变量id[idx]的系数变成对应限制l基变量的系数,所以后面还特判j=idx这一列.
    //注意id[n+l]这个基变量本来都系数是00000100000这样的
    a[l][idx] = 1 / c;
    for (int i = 0; i <= m; i++) {
        if (i == l  \vert  \vert  !sgn(a[i][idx]))
            continue;
        float x = a[i][idx];
        a[i][idx] = -a[l][idx] * x;
        for (int j = 0; j <= n; j++)
            if (j != idx)
                a[i][j] -= a[l][j] * x;
    }
}

//找可行解
bool init() {
    for (int i = 1; i <= n + m; i++)
        id[i] = i;
    while (true) {
        int line = 1, var = 0;
        //找常数最小的行
        for (int i = 1; i <= m; i++)
            if (a[i][0] < a[line][0])
                line = i;
        //已经可行
        if (sgn(a[line][0]) >= 0)
            return true;
        for (int i = 1; i <= n; i++) {
            //id[var]<id[i]是和循环i的顺序对应的,不能变,否则会死循环
            //记住是找的是最大的id[var]
            if (sgn(a[line][i]) < 0 && (var == 0  \vert  \vert  (id[var] < id[i]))) {
                var = i;
            }
        }
        if (!var)
            return false;
        pivot(line, var);
    }
}
bool simplex() {
    while (true) {
        int var = 1, line = 0, temp;
        float lim = inf;
        //找目标函数中系数最大的非基变量
        for (int i = 1; i <= n; i++)
            if (a[0][i] > a[0][var]) {
                var = i;
            }


        if (sgn(a[0][var]) <= 0)
            return true;
        
        //找变量id[var]最紧的限制
        for (int i = 1; i <= m; i++) {
            //这里id[i]>id[line]同样是找最大的,避免死循环,且这个和上面那个保持一致(要么都最大正循环,要么都找最小反着循环)
            if (sgn(a[i][var]) > 0 && (temp = sgn(a[i][0] / a[i][var] - lim), (temp < 0  \vert  \vert  temp == 0 && (!line \vert  \vert id[i] > id[line])))) {
                line = i;
                lim=a[i][0]/a[i][var];
            }
        }
        if (!line)
            return false;
        pivot(line, var);
    }
}
float ans[N];
int main() {
    int t;
    long double in;
    cin >> n >> m >> t;
    for (int i = 1; i <= n; i++)
        cin >> in, a[0][i] = in;
    for (int i = 1; i <= m; i++) {
        for (int j = 1; j <= n; j++)
            cin >> in, a[i][j] = in;
        cin >> in, a[i][0] = in;
    }
    if (!init()) {
        cout << "Infeasible" << endl;
        return 0;
    }
    if (!simplex()) {
        cout << "Unbounded" << endl;
        return 0;
    }
    cout << fixed << setprecision(10) << (long double)-a[0][0] << endl; //答案是-a[0][0],因为我们是把a[0][0]放在等式右边了(消元的时候用的...=b_i的b_i消的a[0][0]).
    if(t==0)return 0;
    for (int i = 1; i <= m; i++) {
        ans[id[i + n]] = a[i][0];
    }
    for (int i = 1; i <= n; i++)
        cout << fixed << setprecision(10) << (long double)ans[i] << " ";
}

速度测试

习题

志愿者招募

线性规划裸题, 变量是每一类志愿者的个数. 限制是每一天的志愿者都到达所需.

因为限制的常数(每一天的志愿者需求量)都非负, 可以跳过寻找可行解的初始化过程.