单纯形法c++

单纯形法（Simplex Algorithm）是一种常用的线性规划算法，它通过不断移动一个n维空间中的多面体（称为单纯形或棱锥）来找到最优解。在每个顶点处计算目标函数的值，然后朝着更优的方向移动。如果无法找到更优的解，则认为已经找到了最优解。

以下是一个简单的单纯形法的C++实现：

#include <iostream>
#include <vector>
#include <algorithm>

using namespace std;

const double eps = 1e-6; // 精度
const int MAXN = 1005; // 最大变量数和约束数

int n, m; // 变量数和约束数
double a[MAXN][MAXN]; // 约束矩阵
int basis[MAXN]; // 基变量的列编号
double ans; // 最优解

// 换入变量
int pivot(int r, int c) {
    double maxv = 0;
    int id = -1;
    for (int i = 1; i <= n; i++) {
        if (i != c && a[r][i] > maxv) {
            maxv = a[r][i];
            id = i;
        }
    }
    return id;
}

// 换出变量
int feasible() {
    while (true) {
        int r, c;
        double minv = 1e9;
        for (int i = 1; i <= m; i++) {
            if (a[i][0] < minv) {
                minv = a[i][0];
                r = i;
            }
        }
        if (minv > -eps) break; // 所有系数都大于等于0，已经找到最优解
        minv = 1e9;
        for (int i = 1; i <= n; i++) {
            if (a[r][i] < -eps && a[0][i] / a[r][i] < minv) {
                minv = a[0][i] / a[r][i];
                c = i;
            }
        }
        if (minv == 1e9) return -1; // 无法找到换出变量，问题无解
        int id = pivot(r, c);
        if (id == -1) return -2; // 无法找到换入变量，问题无界
        basis[r] = id;
        double tmp = a[r][c];
        for (int i = 0; i <= n; i++) a[r][i] /= tmp;
        for (int i = 1; i <= m; i++) {
            if (i != r) {
                tmp = a[i][c];
                for (int j = 0; j <= n; j++) a[i][j] -= a[r][j] * tmp;
            }
        }
    }
    return 0;
}

// 线性规划
double simplex() {
    while (true) {
        int idx = feasible();
        if (idx == -1) return 1e9; // 无法找到换出变量，问题无解
        if (idx == -2) return -1e9; // 无法找到换入变量，问题无界
        bool flag = true;
        for (int i = 1; i <= n; i++) {
            if (a[0][i] < -eps) {
                flag = false;
                int r, c;
                double minv = 1e9;
                for (int j = 1; j <= m; j++) {
                    if (a[j][i] < -eps && a[j][0] / a[j][i] < minv) {
                        minv = a[j][0] / a[j][i];
                        r = j;
                    }
                }
                int id = pivot(r, i);
                if (id == -1) return -2; // 无法找到换入变量，问题无界
                basis[r] = id;
                double tmp = a[r][i];
                for (int j = 0; j <= n; j++) a[r][j] /= tmp;
                for (int j = 1; j <= m; j++) {
                    if (j != r) {
                        tmp = a[j][i];
                        for (int k = 0; k <= n; k++) a[j][k] -= a[r][k] * tmp;
                    }
                }
                break;
            }
        }
        if (flag) break; // 所有系数都大于等于0，已经找到最优解
    }
    ans = a[0][0];
    for (int i = 1; i <= m; i++) {
        if (basis[i] <= n) {
            ans += a[i][0] * a[0][basis[i]];
        }
    }
    return ans;
}

int main() {
    cin >> n >> m;
    for (int i = 1; i <= n; i++) cin >> a[0][i];
    for (int i = 1; i <= m; i++) {
        for (int j = 1; j <= n + 1; j++) {
            cin >> a[i][j];
        }
        basis[i] = n + i;
    }
    ans = simplex();
    if (ans == 1e9) cout << "Infeasible" << endl; // 无法找到换出变量，问题无解
    else if (ans == -1e9) cout << "Unbounded" << endl; // 无法找到换入变量，问题无界
    else printf("%.6f\n", ans);
    return 0;
}

注意，上述代码只适用于标准形式的线性规划问题，即：

$$\max\limits_{x_1,x_2,\cdots,x_n}c_1x_1+c_2x_2+\cdots+c_nx_n$$

$$\text{s.t.}\ \begin{cases}a_{11}x_1+a_{12}x_2+\cdots+a_{1n}x_n\le b_1\\a_{21}x_1+a_{22}x_2+\cdots+a_{2n}x_n\le b_2\\\ \ \ \vdots\\a_{m1}x_1+a_{m2}x_2+\cdots+a_{mn}x_n\le b_m\end{cases}$$

其中，$c_1,c_2,\cdots,c_n,b_1,b_2,\cdots,b_m$均为常数，$a_{11},a_{12},\cdots,a_{mn}$为系数。如果需要解决非标准形式的线性规划问题，需要将其转化为标准形式。