simplex.hpp¶
Maximize cx subject to Ax <= b and x >= 0 with a two-phase simplex tableau; reports optimal, infeasible, or unbounded.
求线性规划 max c·x，约束 Ax≤b、x≥0，并区分最优、无解和无界。
Implementation¶
#ifndef NOYA_SIMPLEX_HPP
#define NOYA_SIMPLEX_HPP 1
/// @complexity Time: Exponential worst case; typically polynomial pivot work per iteration.
/// Space: O(rows * columns).
#include <algorithm>
#include <cassert>
#include <cmath>
#include <limits>
#include <utility>
#include <vector>
namespace noya {
enum class linear_program_status { optimal, infeasible, unbounded };
struct linear_program_result {
linear_program_status status = linear_program_status::infeasible;
long double objective = 0;
std::vector<long double> solution;
};
namespace simplex_internal {
struct tableau {
static constexpr long double epsilon = 1e-12L;
int constraint_count;
int variable_count;
std::vector<int> basic;
std::vector<int> nonbasic;
std::vector<std::vector<long double>> value;
tableau(const std::vector<std::vector<long double>> &a,
const std::vector<long double> &b, const std::vector<long double> &c)
: constraint_count(int(b.size())), variable_count(int(c.size())),
basic(constraint_count), nonbasic(variable_count + 1),
value(constraint_count + 2,
std::vector<long double>(variable_count + 2)) {
for (int row = 0; row < constraint_count; row++) {
for (int column = 0; column < variable_count; column++) {
value[row][column] = a[row][column];
}
basic[row] = variable_count + row;
value[row][variable_count] = -1;
value[row][variable_count + 1] = b[row];
}
for (int column = 0; column < variable_count; column++) {
nonbasic[column] = column;
value[constraint_count][column] = -c[column];
}
nonbasic[variable_count] = -1;
value[constraint_count + 1][variable_count] = 1;
}
void pivot(int leaving_row, int entering_column) {
long double inverse = 1 / value[leaving_row][entering_column];
for (int row = 0; row < constraint_count + 2; row++) {
if (row == leaving_row) {
continue;
}
for (int column = 0; column < variable_count + 2; column++) {
if (column != entering_column) {
value[row][column] -= value[leaving_row][column] *
value[row][entering_column] * inverse;
}
}
}
for (int column = 0; column < variable_count + 2; column++) {
if (column != entering_column) {
value[leaving_row][column] *= inverse;
}
}
for (int row = 0; row < constraint_count + 2; row++) {
if (row != leaving_row) {
value[row][entering_column] *= -inverse;
}
}
value[leaving_row][entering_column] = inverse;
std::swap(basic[leaving_row], nonbasic[entering_column]);
}
bool optimize(int phase) {
int objective_row = phase == 1 ? constraint_count + 1 : constraint_count;
while (true) {
int entering = -1;
for (int column = 0; column <= variable_count; column++) {
if (phase == 2 && nonbasic[column] == -1) {
continue;
}
if (entering == -1 ||
value[objective_row][column] <
value[objective_row][entering] - epsilon ||
(std::abs(value[objective_row][column] -
value[objective_row][entering]) <= epsilon &&
nonbasic[column] < nonbasic[entering])) {
entering = column;
}
}
if (entering == -1) {
return true;
}
if (value[objective_row][entering] >= -epsilon) {
return true;
}
int leaving = -1;
for (int row = 0; row < constraint_count; row++) {
if (value[row][entering] <= epsilon) {
continue;
}
if (leaving == -1) {
leaving = row;
continue;
}
long double first =
value[row][variable_count + 1] / value[row][entering];
long double second =
value[leaving][variable_count + 1] / value[leaving][entering];
if (first < second - epsilon || (std::abs(first - second) <= epsilon &&
basic[row] < basic[leaving])) {
leaving = row;
}
}
if (leaving == -1) {
return false;
}
pivot(leaving, entering);
}
}
linear_program_result solve() {
int lowest_row = 0;
for (int row = 1; row < constraint_count; row++) {
if (value[row][variable_count + 1] <
value[lowest_row][variable_count + 1]) {
lowest_row = row;
}
}
if (constraint_count > 0 &&
value[lowest_row][variable_count + 1] < -epsilon) {
pivot(lowest_row, variable_count);
if (!optimize(1) ||
value[constraint_count + 1][variable_count + 1] < -epsilon) {
return {linear_program_status::infeasible, 0, {}};
}
if (std::abs(value[constraint_count + 1][variable_count + 1]) > epsilon) {
return {linear_program_status::infeasible, 0, {}};
}
for (int row = 0; row < constraint_count; row++) {
if (basic[row] != -1) {
continue;
}
int entering = 0;
for (int column = 1; column <= variable_count; column++) {
if (std::abs(value[row][column]) >
std::abs(value[row][entering]) + epsilon ||
(std::abs(std::abs(value[row][column]) -
std::abs(value[row][entering])) <= epsilon &&
nonbasic[column] < nonbasic[entering])) {
entering = column;
}
}
if (std::abs(value[row][entering]) > epsilon) {
pivot(row, entering);
}
}
}
if (!optimize(2)) {
return {linear_program_status::unbounded,
std::numeric_limits<long double>::infinity(),
{}};
}
std::vector<long double> solution(variable_count);
for (int row = 0; row < constraint_count; row++) {
if (basic[row] < variable_count) {
solution[basic[row]] = value[row][variable_count + 1];
}
}
return {linear_program_status::optimal,
value[constraint_count][variable_count + 1], std::move(solution)};
}
};
} // namespace simplex_internal
/// @brief Maximize c*x subject to A*x <= b and x >= 0 with a two-phase
/// simplex tableau; reports optimal, infeasible, or unbounded.
inline linear_program_result
simplex(const std::vector<std::vector<long double>> &a,
const std::vector<long double> &b, const std::vector<long double> &c) {
assert(a.size() == b.size());
for (const auto &row : a) {
assert(row.size() == c.size());
}
return simplex_internal::tableau(a, b, c).solve();
}
} // namespace noya
#endif // NOYA_SIMPLEX_HPP
#include <algorithm>
#include <cassert>
#include <cmath>
#include <limits>
#include <utility>
#include <vector>
/// @complexity Time: Exponential worst case; typically polynomial pivot work per iteration.
/// Space: O(rows * columns).
namespace noya {
enum class linear_program_status { optimal, infeasible, unbounded };
struct linear_program_result {
linear_program_status status = linear_program_status::infeasible;
long double objective = 0;
std::vector<long double> solution;
};
namespace simplex_internal {
struct tableau {
static constexpr long double epsilon = 1e-12L;
int constraint_count;
int variable_count;
std::vector<int> basic;
std::vector<int> nonbasic;
std::vector<std::vector<long double>> value;
tableau(const std::vector<std::vector<long double>> &a,
const std::vector<long double> &b, const std::vector<long double> &c)
: constraint_count(int(b.size())), variable_count(int(c.size())),
basic(constraint_count), nonbasic(variable_count + 1),
value(constraint_count + 2,
std::vector<long double>(variable_count + 2)) {
for (int row = 0; row < constraint_count; row++) {
for (int column = 0; column < variable_count; column++) {
value[row][column] = a[row][column];
}
basic[row] = variable_count + row;
value[row][variable_count] = -1;
value[row][variable_count + 1] = b[row];
}
for (int column = 0; column < variable_count; column++) {
nonbasic[column] = column;
value[constraint_count][column] = -c[column];
}
nonbasic[variable_count] = -1;
value[constraint_count + 1][variable_count] = 1;
}
void pivot(int leaving_row, int entering_column) {
long double inverse = 1 / value[leaving_row][entering_column];
for (int row = 0; row < constraint_count + 2; row++) {
if (row == leaving_row) {
continue;
}
for (int column = 0; column < variable_count + 2; column++) {
if (column != entering_column) {
value[row][column] -= value[leaving_row][column] *
value[row][entering_column] * inverse;
}
}
}
for (int column = 0; column < variable_count + 2; column++) {
if (column != entering_column) {
value[leaving_row][column] *= inverse;
}
}
for (int row = 0; row < constraint_count + 2; row++) {
if (row != leaving_row) {
value[row][entering_column] *= -inverse;
}
}
value[leaving_row][entering_column] = inverse;
std::swap(basic[leaving_row], nonbasic[entering_column]);
}
bool optimize(int phase) {
int objective_row = phase == 1 ? constraint_count + 1 : constraint_count;
while (true) {
int entering = -1;
for (int column = 0; column <= variable_count; column++) {
if (phase == 2 && nonbasic[column] == -1) {
continue;
}
if (entering == -1 ||
value[objective_row][column] <
value[objective_row][entering] - epsilon ||
(std::abs(value[objective_row][column] -
value[objective_row][entering]) <= epsilon &&
nonbasic[column] < nonbasic[entering])) {
entering = column;
}
}
if (entering == -1) {
return true;
}
if (value[objective_row][entering] >= -epsilon) {
return true;
}
int leaving = -1;
for (int row = 0; row < constraint_count; row++) {
if (value[row][entering] <= epsilon) {
continue;
}
if (leaving == -1) {
leaving = row;
continue;
}
long double first =
value[row][variable_count + 1] / value[row][entering];
long double second =
value[leaving][variable_count + 1] / value[leaving][entering];
if (first < second - epsilon || (std::abs(first - second) <= epsilon &&
basic[row] < basic[leaving])) {
leaving = row;
}
}
if (leaving == -1) {
return false;
}
pivot(leaving, entering);
}
}
linear_program_result solve() {
int lowest_row = 0;
for (int row = 1; row < constraint_count; row++) {
if (value[row][variable_count + 1] <
value[lowest_row][variable_count + 1]) {
lowest_row = row;
}
}
if (constraint_count > 0 &&
value[lowest_row][variable_count + 1] < -epsilon) {
pivot(lowest_row, variable_count);
if (!optimize(1) ||
value[constraint_count + 1][variable_count + 1] < -epsilon) {
return {linear_program_status::infeasible, 0, {}};
}
if (std::abs(value[constraint_count + 1][variable_count + 1]) > epsilon) {
return {linear_program_status::infeasible, 0, {}};
}
for (int row = 0; row < constraint_count; row++) {
if (basic[row] != -1) {
continue;
}
int entering = 0;
for (int column = 1; column <= variable_count; column++) {
if (std::abs(value[row][column]) >
std::abs(value[row][entering]) + epsilon ||
(std::abs(std::abs(value[row][column]) -
std::abs(value[row][entering])) <= epsilon &&
nonbasic[column] < nonbasic[entering])) {
entering = column;
}
}
if (std::abs(value[row][entering]) > epsilon) {
pivot(row, entering);
}
}
}
if (!optimize(2)) {
return {linear_program_status::unbounded,
std::numeric_limits<long double>::infinity(),
{}};
}
std::vector<long double> solution(variable_count);
for (int row = 0; row < constraint_count; row++) {
if (basic[row] < variable_count) {
solution[basic[row]] = value[row][variable_count + 1];
}
}
return {linear_program_status::optimal,
value[constraint_count][variable_count + 1], std::move(solution)};
}
};
} // namespace simplex_internal
/// @brief Maximize c*x subject to A*x <= b and x >= 0 with a two-phase
/// simplex tableau; reports optimal, infeasible, or unbounded.
inline linear_program_result
simplex(const std::vector<std::vector<long double>> &a,
const std::vector<long double> &b, const std::vector<long double> &c) {
assert(a.size() == b.size());
for (const auto &row : a) {
assert(row.size() == c.size());
}
return simplex_internal::tableau(a, b, c).solve();
}
} // namespace noya