simplex method

Posted on 2022-07-10 Edited on 2024-07-21

单纯形法线性规划

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

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

求解下列线性规划：

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

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

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

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

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

如果最后一行还有负数，找到最小的负数所在的列（ $j$ )，再找到 $b_{i} / a_{i j}$ 的最小值（非负）所在的行（ $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])：

#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]找的，和上述步骤有些许不同，下面是代码实验:

求解下列问题：

输入：

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_{1 j} (j = 1, 2, \dots, 5)$ 均为正。

6. 例题

模板题

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

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

$x_{1}$	$x_{2}$	$x_{4}$	$x_{5}$	$x_{6}$
-2	1	1	0	0	-4
1	1	0	1	0	4
1	-2	0	0	1	-4

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

$x_{1}$	$x_{2}$	$x_{4}$	$x_{5}$	$x_{6}$
1	$- \frac{1}{2}$	$- \frac{1}{2}$	0	0	2
0	$\frac{3}{2}$	$\frac{1}{2}$	1	0	2
0	$- \frac{3}{2}$	$\frac{1}{2}$	0	1	-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_{3 j}$ 小于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]理论罗列
#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] 单纯形法