Skip to main content

__debug's Home Keep it simple, stupid

对单纯形算法的一些提示

教程

推荐阅读算法导论上的第 29 章, 讲得十分详细.

如果觉得内容太多, 第一遍阅读可以跳过 29.2 和 29.4 节.

前面的部分应该都很好理解, 主要是 \(\text{INITIALIZE-SIMPLEX}\) 稍微困难一些.

实现

即便看懂了算法导论, 也比较难打出简洁高效的代码…

如果是第一次, 还是建议阅读别人的代码. 顺便给出几个实现上的提示:

  1. 不用 \(A, b, c\) 三个数组来表示系数, 而只用一个 \(A\) 来表示所有的系数. (比如说, 我的代码里 \(A_{1..m}\) 是限制, \(A_{m + 1}\) 则是目标函数; 同时 \(A_{i, 1..n}\) 是变量前的系数, \(A_{i, n + 1}\) 则是常数项) 这样在 pivot 时会方便很多
  2. 接上面, 如果你只用一个 \(A\) 来表示所有的系数, 那么所有 \(A_{i, j}\) 前的符号最好一样 (比如说都是 \(+\)), 而不是像算导里的 \(x_i = b_i \color{red}{-} \sum_{j = 1}^{n} A_{i, j} x_j\)
  3. 算导中的 \(\text{INITIALIZE-SIMPLEX}\) 中最后一部分十分复杂, 即消除 \(x_0\) 的影响的部分; 在代码实现中, 不妨加入两条 \(x_0 > 0\) 和 \(x_0 < 0\) 的限制 (在程序中这里有一个很 trick 的写法, 详见 init 的最后一部分)

UOJ 179 是模板题, 贴下代码.

const int MAXN = 31, MAXM = 31;
const double INF = 1e18, EPS = 1e-8;

namespace Simplex
{

int n, m;
double a[MAXM][MAXN];
int idx[MAXN];                  // index of non-basic variable
int idy[MAXM];                  // index of basic variable

void pivot(int x, int y)
{
    double t = -1 / a[y][x];
    a[y][x] = -1;
    for (int i = 0; i <= n + 1; ++i) {
	a[y][i] *= t;
    }

    for (int i = 1; i <= m + 3; ++i) { // m must + 3
	if (i != y && fabs(a[i][x]) > EPS) {
	    double t0 = a[i][x];
	    a[i][x] = 0;
	    for (int j = 0; j <= n + 1; ++j) {
		a[i][j] += a[y][j] * t0;
	    }
	}
    }

    std::swap(idx[x], idy[y]);
}

double simplex()
{
    for (;;) {
	int x = -1;
	for (int i = 0; i <= n; ++i) {
	    if (a[m + 1][i] > EPS) {
		x = i;
		break;
	    }
	}
	if (x == -1)
	    break;

	int y = -1;
	double mindelta = INF;
	for (int i = 1; i <= m; ++i) {
	    if (a[i][x] < -EPS && chkmin(mindelta, -a[i][n + 1] / a[i][x]))
		y = i;
	}
	if (y == -1)            // unbounded
	    return +INF;

	pivot(x, y);
    }

    return a[m + 1][n + 1];
}

bool init(int n_, int m_)
{
    n = n_;
    m = m_;
    for (int i = 0; i <= n; ++i) {
	idx[i] = i;
    }
    for (int i = 1; i <= m + 2; ++i) {
	idy[i] = i + n;
    }

    // setting the objective funtion to -x_0
    a[m + 1][0] = -1;
    // adding +x_0 to each constraints
    for (int i = 1; i <= m; ++i) {
	a[i][0] = +1;
    }

    // let k be the index of minimum b_i
    int k = 1;
    for (int i = 2; i <= m; ++i) {
	if (a[i][n + 1] < a[k][n + 1])
	    k = i;
    }

    pivot(0, k);
    if (simplex() < -EPS)     // infeasible
	return false;

    // adding a constraint x_0 = 0
    for (int i = 0; i <= n + 1; ++i) {
	a[m + 2][i] = -a[m + 1][i];
    }
    // setting the objective funtion to the original
    m += 2;

    return true;
}

void printAnswer(int n_)
{
    static double ans[MAXN * 2 + 5];

    for (int i = 1; i <= m; ++i) {
	ans[idy[i]] = a[i][n + 1];
    }

    for (int i = 1; i <= n_; ++i) {
	printf("%.10f ", ans[i]);
    }
    putchar('\n');
}

}

int N, M, T;

void input()
{
    read(N); read(M); read(T);
    for (int i = 1; i <= N; ++i) {
	// we don't need the objective when initializing
	scanf("%lf", &Simplex::a[M + 3][i]);
    }
    for (int i = 1; i <= M; ++i) {
	for (int j = 1; j <= N; ++j) {
	    scanf("%lf", &Simplex::a[i][j]);
	    // convert standard form to augmented form
	    // attention that the sign of coefficient is different compared with CLRS
	    Simplex::a[i][j] = -Simplex::a[i][j];
	}
	scanf("%lf", &Simplex::a[i][N + 1]);
    }
}

void solve()
{
    if (Simplex::init(N, M) == false) {
	puts("Infeasible");
	return;
    }

    double ans = Simplex::simplex();

    if (ans == INF)
	puts("Unbounded");
    else {
	printf("%.10f\n", ans);
	if (T)
	    Simplex::printAnswer(N);
    }
}

完整代码见 179.cpp.

Comments

Comments powered by Disqus