用c++语言编写程序实现单纯形法
时间: 2023-06-13 15:07:37 浏览: 63
单纯形法是一种线性规划求解方法,能够在多项式时间内求解一类特殊的线性规划问题。下面是一个简单的单纯形法实现,仅供参考:
```c++
#include <iostream>
#include <vector>
#include <algorithm>
#include <iomanip>
using namespace std;
const double EPS = 1e-6; // 精度
// 矩阵类
class Matrix {
public:
Matrix(int n, int m) : n_(n), m_(m), data_(n * m) {}
double& operator() (int i, int j) {
return data_[i * m_ + j];
}
const double& operator() (int i, int j) const {
return data_[i * m_ + j];
}
int rows() const {
return n_;
}
int cols() const {
return m_;
}
Matrix transpose() const {
Matrix res(m_, n_);
for (int i = 0; i < n_; ++i)
for (int j = 0; j < m_; ++j)
res(j, i) = (*this)(i, j);
return res;
}
private:
int n_, m_;
vector<double> data_;
};
// 单纯形法求解
bool simplex(Matrix& A, Matrix& b, Matrix& c, double& res) {
int n = A.rows(), m = A.cols();
Matrix B(n, n), N(n, m - n), cB(n, 1), cN(m - n, 1), xB(n, 1), xN(m - n, 1), yT(1, n);
vector<int> basis(n), nonbasis(m - n);
for (int i = 0; i < n; ++i) {
basis[i] = i;
for (int j = 0; j < n; ++j)
B(i, j) = A(i, j);
xB(i, 0) = b(i, 0);
}
for (int i = 0; i < m - n; ++i) {
nonbasis[i] = i + n;
for (int j = 0; j < n; ++j)
N(j, i) = A(j, i + n);
xN(i, 0) = 0;
}
cB = c.transpose().get_submatrix(basis);
cN = c.transpose().get_submatrix(nonbasis);
while (true) {
// 计算 yT
yT = cB.transpose() * B.inverse();
// 计算 delta
double delta = 0;
int pos = -1;
for (int i = 0; i < m - n; ++i) {
double d = cN(i, 0) - yT * N.get_column(i);
if (d > delta) {
delta = d;
pos = i;
}
}
if (pos == -1) { // 最优解
res = cB.transpose() * xB;
return true;
}
// 计算 xN
Matrix dN = N.get_column(pos);
Matrix dB = B.inverse() * dN;
double theta = 1e20;
pos = -1;
for (int i = 0; i < n; ++i) {
if (dB(i, 0) > EPS && xB(i, 0) / dB(i, 0) < theta) {
theta = xB(i, 0) / dB(i, 0);
pos = i;
}
}
if (pos == -1) // 无界
return false;
// 更新基变量
xB = xB - theta * dB;
xB(pos, 0) = theta;
basis[pos] = nonbasis[pos];
nonbasis[pos] = -1;
// 更新基矩阵
swap(B(pos, pos), N(pos, pos - n));
for (int i = 0; i < m - n; ++i) {
if (i == pos)
continue;
Matrix dn = N.get_column(i) - N.get_column(pos) * N(pos, i);
Matrix db = B.inverse() * dn;
db(pos, 0) = -1;
swap(xN(i, 0), xB(pos, 0));
swap(nonbasis[i], basis[pos]);
N.set_column(i, dn);
B.set_column(pos, db);
}
cB = c.transpose().get_submatrix(basis);
cN = c.transpose().get_submatrix(nonbasis);
}
}
int main() {
int n, m;
cin >> n >> m;
Matrix A(n, m), b(n, 1), c(m, 1);
for (int i = 0; i < n; ++i)
for (int j = 0; j < m; ++j)
cin >> A(i, j);
for (int i = 0; i < n; ++i)
cin >> b(i, 0);
for (int i = 0; i < m; ++i)
cin >> c(i, 0);
double res;
if (simplex(A, b, c, res))
cout << fixed << setprecision(6) << res << endl;
else
cout << "Infeasible or unbounded!" << endl;
return 0;
}
```
该程序的输入格式为:
```
n m
a11 a12 ... a1m
a21 a22 ... a2m
...
an1 an2 ... anm
b1
b2
...
bn
c1 c2 ... cm
```
其中,$n$ 和 $m$ 分别表示约束条件的个数和变量的个数,$a_{ij}$ 表示第 $i$ 个约束条件中第 $j$ 个变量的系数,$b_i$ 表示第 $i$ 个约束条件的右边的常数,$c_j$ 表示第 $j$ 个变量的系数。程序输出最优解或者提示无解或无界。
相关推荐
![](https://csdnimg.cn/download_wenku/file_type_ask_c1.png)
![](https://csdnimg.cn/download_wenku/file_type_ask_c1.png)
![](https://csdnimg.cn/download_wenku/file_type_ask_c1.png)
![](https://csdnimg.cn/download_wenku/file_type_ask_c1.png)
![](https://csdnimg.cn/download_wenku/file_type_ask_c1.png)
![](https://csdnimg.cn/download_wenku/file_type_ask_c1.png)
![](https://csdnimg.cn/download_wenku/file_type_ask_c1.png)
![](https://csdnimg.cn/download_wenku/file_type_ask_c1.png)
![](https://csdnimg.cn/download_wenku/file_type_ask_c1.png)
![](https://csdnimg.cn/download_wenku/file_type_ask_c1.png)
![](https://csdnimg.cn/download_wenku/file_type_ask_c1.png)
![](https://csdnimg.cn/download_wenku/file_type_ask_c1.png)