simplex method

Posted on Edited on

单纯形法线性规划

1. 单纯形法求解线性规划步骤

废话不多说，从一个例子开始：

求解下列线性规划：

首先预处理，令： \[ M=25x_1+33x_2+18x_3\\ \] 然后在第一个约束条件中添加松弛变量(slack variables)\(x_4(>0)\)，得到： \[ 2x_1+3x_2+4x_3+x4=60 \] 第二个约束条件和第三个约束条件以此类推，然后把预处理完的问题写成下列的形式：

去掉变量，构造单纯形表，下面的表格被称为初始单纯形表（initial simplex tableau）：

在最后一行找到最小的数（-33）所在的列（第2列）。

令\(b_i\)表示最后一列第\(i\)行的数，令\(a_{i2}\)表示第2列第\(i\)行的数，找到\({b_i}/{a_{i2}}\)的值最小的行（第一行，60/3 < 50/2 < 46/1），如下图所示：

然后第1行进行归一化（除以\(a_{12}\))，并用对其他行高斯消元：

如果最后一行还有负数，找到最小的负数所在的列（\(j\))，再找到\({b_i}/{a_{ij}}\)的最小值（非负）所在的行（\(i\)），归一化后对其他行进行高斯消元。

最后一行的系数中没有负数，结束。

\(M\)的最大值即为\(\frac{4854}{7}\)，\(x_1,x_2,x_3\)的值分别为\(\frac{78}{7},\frac{88}{7},0\)。

2. 单纯形算法的几何解释

通过几何来理解更加直观，往往能建立更好的数学直觉，下面举一个例子。在变量数目较少的情况下我们可以使用“图解法”，如下面的问题：

对应的图示，其中蓝色的区域被称为可行区域，其中的每个点对应一组可行解：

可得最大值在\(B(20,30,0)\)取到，为130。

我们使用单纯形算法计算一遍：

\(x_1\) \(x_2\) \(x_3\) \(x_4\) \(x_5\) \(M\)
1 1 1 1 0 0 50
1 2 4 0 1 0 80
-2 -3 -4 0 0 1 0

初始情况下\(x_1\) \(x_2\) \(x_3\) \(x_4\) \(x_5\)的值分别为0,0,0,50,80。\([x_1,x_2,x_3]^T\)为\([0,0,0]^T\)正好对应图中的原点。

选择第3列，第2行，归一化，消元得:

\(x_1\) \(x_2\) \(x_3\) \(x_4\) \(x_5\) \(M\)
\(\frac{3}{4}\) \(\frac{1}{2}\) 0 1 \(-\frac{1}{4}\) 0 30
\(\frac{1}{4}\) \(\frac{1}{2}\) 1 0 \(\frac{1}{4}\) 0 20
-1 -1 0 0 1 1 80

此时\([x_1,x_2,x_3]^T\)为\([0,0,20]^T\)正好对应图中的\(D\)点。

选择第1列，第1行，进行同样的操作得：

\(x_1\) \(x_2\) \(x_3\) \(x_4\) \(x_5\) \(M\)
1 \(\frac{2}{3}\) 0 \(\frac{4}{3}\) \(-\frac{1}{3}\) 0 40
0 \(\frac{1}{3}\) 1 \(-\frac{1}{3}\) \(\frac{1}{3}\) 0 10
0 \(-\frac{1}{3}\) 0 \(\frac{4}{3}\) \(\frac{2}{3}\) 1 120

此时\([x_1,x_2,x_3]^T\)为\([40,0,10]^T\)正好对应图中的\(E\)点。

选择第2列，第2行，继续操作：

\(x_1\) \(x_2\) \(x_3\) \(x_4\) \(x_5\) \(M\)
1 0 -2 2 -1 0 20
0 1 3 -1 1 0 30
0 0 1 1 1 1 130

此时\([x_1,x_2,x_3]^T\)为\([20,30,0]^T\)正好对应图中的\(B\)点。

进行到这里，算法结束，正好就是最大值。

由单纯形法的计算过程可以看出\([x_1,x_2,x_3]^T\)沿着”单纯形“的边不断进行调整变化，直到达到最优点，即：

\([0,0,0]^T\Rightarrow[0,0,20]^T\Rightarrow[40,0,10]^T\Rightarrow[20,30,0]^T\)

3. 单纯形算法代码模板和例题

全部代码给出如下（代码来源[2])：

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
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
#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();
    system("pause");
}
/////////////////////////////////////
// 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
*/
/////////////////////////////////////

注：上述代码是从[2]找的，和上述步骤有些许不同，下面是代码实验:

求解下列问题：

输入：

1
2
3
4
6 3
2 3 4 0 0 0 
1 1 1 1 0 50
1 2 4 0 1 80

输出：

4. 线性规划问题的对偶问题（dual problem）

每一个标准的（canonical）线性规划问题（原问题，primal problem）都有一个相互关联的对偶问题（dual problem）。下面是原问题和对偶问题形式上的比较：

举一个简单的例子，原问题如下：

其对偶问题为：

原问题的最大值即为对偶问题的最小值。

5. 无界（unbounded）和不可行（infeasible）

5.1 无界

单纯形表出现下面的情况时说明无界，最优解为无穷。

\(x_1\) \(x_2\) \(x_3\) \(x_4\) \(x_5\) \(M\)
1 1 -1 1 0 0 50
1 2 -4 0 1 0 80
-2 -3 -4 0 0 1 0

按照算法应该选中第3列的某一行进行归一化，但是第三列全是负值，说明无界。

5.2 不可行

单纯形表出现下面的情况时说明可行解集为空，无解。

\(x_1\) \(x_2\) \(x_3\) \(x_4\) \(x_5\) \(M\)
1 1 1 1 0 0 -50
1 2 4 0 1 0 80
-2 -3 -4 0 0 1 0

如上图所示，第1行\(b_1\)为负数，但是\(a_{1j}(j=1,2,\cdots,5)\)均为正。

6. 例题

模板题

6.1 使用样例3来说明不可行情况下的单纯形表：

初始单纯形表（没有最后一行目标函数，为了方便省略了）：

\(x_1\) \(x_2\) \(x_3\) \(x_4\) \(x_5\) \(x_6\) \(M\)
-2 1 0 1 0 0 0 -4
1 1 0 0 1 0 0 4
1 -2 0 0 0 1 0 -4

最右列\(b\)含有小于0的情况，要先进行初始化。\(b_1=-4，a_{11}=-2\)，小于0，因此要把第1列对应的\(x_1\)引入。那么哪一行归一化呢？计算所有\(b_i/a_{i1}>0\)的情况（这里\(i=1,2\)），取最小的那一行归一化（第1行），然后高斯消元，得到：

\(x_1\) \(x_2\) \(x_3\) \(x_4\) \(x_5\) \(x_6\) \(M\)
1 \(-\frac{1}{2}\) 0 \(-\frac{1}{2}\) 0 0 0 2
0 \(\frac{3}{2}\) 0 \(\frac{1}{2}\) 1 0 0 2
0 \(-\frac{3}{2}\) 0 \(\frac{1}{2}\) 0 1 0 -6

重复上述操作，\(b_3=-6\)，并且\(a_{32}<0\)把这一行对应的\(x_2\)引入。这里应该选择第2行进行归一化，因为\(b_2/a_{22}<b_3/a_{32}\)，第2行归一化之后对其它行进行高斯消元，这步操作之后我们发现第3行为：

\(x_1\) \(x_2\) \(x_3\) \(x_4\) \(x_5\) \(x_6\) \(M\)
0 0 0 1 1 1 0 -4

\(b_3\)依然小于0，但是找不到一个\(a_{3j}\)小于0了，说明可行域为空，无解。

6.2 使用样例4来说明无界情况下的单纯形表：

初始单纯形表：

\(x_1\) \(x_2\) \(x_3\) \(M\)
1 0 1 0 1
0 -1 0 1 0

目标函数行（第2行）对应最小负数列（第2列）为-1，但是第2列上没有正数，说明无界，最优解为无穷。

代码：

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
123
124
125
126
127
//见[1]理论罗列
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <iostream>
#include <algorithm>
#include <ctime>
using namespace std;
#define Maxn 25
const double eps = 0.00000001, INF = 1e15;

int n, m;

int id[Maxn * 2];//id[n+1]~id[n+m]记录了m个约束对应的主元？？
double a[Maxn][Maxn];
//第一维是限制，B集合
//第二维是元素，N集合
//a[0][xx] -> c 目标函数系数
//a[xx][0] -> b 限制等式常数
//a[xx][yy] -> A 限制等式系数向量
//最大化 sigma(ci*xi),i属于N
//限制 xj=bj-sigma(aji*xi) ,j属于B

double myabs(double x) { return x > 0 ? x : -x; }

void Pivot(int l, int e)
{
    //转轴l和e
    swap(id[n + l], id[e]);
    double t = a[l][e];
    a[l][e] = 1;
    for (int j = 0; j <= n; j++)
        a[l][j] /= t;
    for (int i = 0; i <= m; i++)
        if (i != l && (myabs(a[i][e]) > eps))//行不等于l，并且a[i][e] ！= 0
        {
            t = a[i][e];
            a[i][e] = 0;
            for (int j = 0; j <= n; j++)
                a[i][j] -= a[l][j] * t;
        }
}

//初始化-辅助问题
bool init()
{
    while (1)
    {
        int e = 0, l = 0;
        for (int i = 1; i <= m; i++)
            if ((a[i][0] < -eps) && (!l || (rand() & 1)))//在所有bi<0中随机选择一个作为换出,没加括号之前，烂代码的代表！！小于号和与运算的优先级
                l = i;
        if (!l)//bi都大于0，初始化结束。
            break;
        for (int j = 1; j <= n; j++)
            if ((a[l][j] < -eps) && (!e || (rand() & 1)))//在所有a[i]中随机选择一个作为换出
                e = j;
        if (!e)//bl小于0但是这一行的alj都大于0，无解，可行域为空
        {
            printf("Infeasible\n");
            return 0;
        }
        Pivot(l, e);
    }
    return 1;
}

//最优化
bool simplex()
{
    while (1)
    {
        int l = 0, e = 0;
        double mn = INF;
        for (int j = 1; j <= n; j++)
            if (a[0][j] > eps)//大于0
            {
                e = j;
                break;
            }
        if (!e)
            break; //如果目标变量c都小于0 找到答案
        for (int i = 1; i <= m; i++)
            if ((a[i][e] > eps) && (a[i][0] / a[i][e] < mn))
                mn = a[i][0] / a[i][e], l = i; //找a[i][0]/a[i][e]最小的i进行转轴
        if (!l)
        {
            printf("Unbounded\n");
            return 0;
        }
        //如果所有的a[i][e]都小于0，说明最优值正无穷
        Pivot(l, e);
    }
    return 1;
}

double ans[Maxn];

int main()
{
    srand(time(0));
    int t;
    scanf("%d%d%d", &n, &m, &t);
    for (int i = 1; i <= n; i++)
        scanf("%lf", &a[0][i]);
    for (int i = 1; i <= m; i++)
    {
        for (int j = 1; j <= n; j++)
            scanf("%lf", &a[i][j]);
        scanf("%lf", &a[i][0]);
    }
    for (int i = 1; i <= n; i++)
        id[i] = i;
    if (init() && simplex())
    {
        printf("%.8lf\n", -a[0][0]);
        if (t)
        {
            for (int i = 1; i <= m; i++)
                ans[id[n + i]] = a[i][0];
            for (int i = 1; i <= n; i++)
                printf("%.8lf ", ans[i]);
        }
    }
    //system("pause");
    return 0;
}

参考资料

[1] linear algebra and its application chapter 9 Optimization

[2] 单纯形算法

[3] 线性规划之单纯形法【超详解+图解】

[4] 单纯形法