修正单纯型法C++实现

修正单纯型法C++实现

修正单纯型法C++实现
/*****************************************************************************

* @brief    : Simplex algorithm

* @author   : mine

* @date     : 2020/10/18 10:35

* @version  : ver 1.0

*****************************************************************************/

#include <iostream>
#include <Eigen/Eigen/Dense>
#include <vector>
#include <ctime>
#include <windows.h>
#include <iomanip>

using namespace std;
using namespace Eigen;

int main()
{
	SetConsoleTitle(L"Simplex algorithm");

	/************************************************************************
	 * m:约束不等式个数 
	 * n:结构向量的个数                                                                     
	 ************************************************************************/

	int m = 0, n = 0; 
	cout << "－－－－－－－－－－单纯形法的参数输入－－－－－－－－－－－" << endl;
	cout << "请按提示输入下列参数：" << endl << endl;
	cout << "  变量个数n:      " << "   n =  ";
	cin >> n;
	cout << endl << "  约束不等式个数m:" << "   m =  ";
	cin >> m;

	/************************************************************************
	* A:系数矩阵   A = [B|N]
	* b:右端项矩阵
	* c:价格系数矩阵  c = [c_B|c_N]
	* B:基，非奇异 m*m 矩阵  --> 基变量x^*
	* N:非基，m*(n-m) 矩阵   --> 非基变量 x
	* I_B:基变量下标集合
	* I_N:非基变量下标集合
	* r:既约费用向量 r_j = c_j - (c_B)^T * B^(-1) * A_j  
	* p:单纯形乘子
	* d:单纯形法改进方向
	************************************************************************/

	MatrixXf A(m, n);
	VectorXf b(m), c(n), c_B(m), c_N(n - m);
	VectorXi I_B(m), I_N(n - m);

	//规范输出格式
	//IOFormat CommaInitFmt(StreamPrecision, DontAlignCols, ", ", ", ", "", "", " << ", ";");


	cout << endl << "  约束方程矩阵的系数:" << endl << endl << "          ";
	
	for (int i = 0; i <n; i++)
	{
		cout << "   x" << i + 1;
	}

	cout << "  = " << " 常数项" << endl;

	//输入约束方程系数
	for (int i = 0; i < m; i++)
	{
		cout << "    不等式" << i + 1 << "  ";

		for (int j = 0; j < n ; j++)
		{
			cin >> A(i, j);
		}
		cin >> b(i);
	}
	//cout << "系数矩阵 A = " << endl << A;
	//cout << endl;
	//cout << "右端项矩阵 b = " << endl << b;

	//初始化I_B
	cout << endl << "  选择基变量为：\t" << endl;
	for (int i = 0; i < m; i++)
	{
		cout << "\t\tx";
		int index = 0;
		cin >> index;
		I_B(i) = index;
	}
	//cout << "基变量下标矩阵 I_B = " << endl  << I_B;

	int index = 0;
	//初始化I_N
	for (int i = 0; i < n; i++)
	{

		if (i + 1 != I_B(index) )
		{
			I_N(i - index) = i + 1;
		}
		else
		{
			index += 1;
		}
	}
	//cout << endl << "非基变量下标矩阵 I_N = " << endl  << I_N;

	////初始化I_N
	//for (int i = 0; i < n; i++)
	//{
	//	int num = 0;
	//	if (i != I_B(i))
	//	{
	//		I_N(i) = i;
	//	}
	//	else
	//	{
	//		num += 1;
	//	}
	//}
	VectorXf p(m);
	VectorXf r(n - m);
	MatrixXf d(m, n - m);
	VectorXf d_q(m);
	VectorXf step(m);


	//输入目标函数系数 
	cout << endl << endl << "  目标函数的系数:" << endl << endl << "                ";
	for (int i = 1; i <= n; i++)
	{
		cout << "x" << i << "   ";
	}

	cout << endl << "  目标函数:   ";
	for (int i = 0; i < n; i++)
	{
		cin >> c(i);
	}
	//cout << endl << "价格系数矩阵 c = " << endl << c << endl;

	clock_t tStart = clock();

	//初始化B,c_B
	MatrixXf B(m, m);
	for (int i = 0; i < m; i++)
	{
		B.col(i) = A.col(I_B(i)-1);
		c_B(i) = c(I_B(i)-1);
	}
	//cout << endl << "基 B = " << endl << B << endl;
	//cout << endl << "基变量系数 c_B = " << endl << c_B;

	//初始化N,c_N
	MatrixXf N(m, n - m);
	for (int i = 0; i < n - m; i++)
	{
		N.col(i) = A.col(I_N(i)-1);
		c_N(i) = c(I_N(i)-1);
	}
	//cout << endl << "非基 N = " << endl << N << endl;
	//cout << endl << "非基变量系数 c_N = " << endl << c_N;


	while (true)
	{
		//初始化p   p = p^T
		p.transpose() = c_B.transpose()*B.inverse();
		cout << endl << " 单纯形算子 p = " << endl << p << endl;

		//求既约费用r

		for (int i = 0; i < n - m; i++)
		{
			r(i) = c_N(i) - p.transpose() * N.col(i);
		}

		ptrdiff_t max, min;
		float maxOfr = r.maxCoeff(&max);
		float minOfr = r.minCoeff(&min);
		float z = 0;
		if (minOfr >= 0)
		{
			cout << "\t现行解为最优解,最优解为：" << endl;
			for (int i = 0; i < m; i++)
			{
				cout << "\t\t\tx" << I_B(i) << "  = " << b(i) << endl;
				
				z += c(I_B(i) - 1)*b(i);
			}
			for (int i = 0; i < n - m; i++)
			{
				cout << "\t\t\tx" << I_N(i) << "  = " << "0" << endl;
			}
			cout << "  目标函数值为 " << z << endl;
			break;
		}

		for (int i = 0; i < n - m; i++)
		{
			if (r(i) < 0)
			{
				cout << endl << "进基变量为 x" << I_N(i) << endl;
				d = -1 * B.inverse()*N;
				d_q = d.col(i);
				float minOfd_q = d_q.minCoeff(&min);
				if (minOfd_q >= 0)
				{
					cout << "\t\t\t计算完毕，该问题无界" << endl;
				}

				for (int j = 0; j < m; j++)
				{
					if (d_q(j) <= 0)
					{
						step(j) = -1 * b(j) / d_q(j);
					}
				}

				float minOfstep = step.minCoeff(&min);
				float maxOfstep = step.minCoeff(&max);

				cout << endl << "退基变量为 x" << I_B(min) << endl;
				int index = I_B(min);
				I_B(min) = I_N(i);
				I_N(i) = index;

				//cout << I_B << endl;
				//cout << I_N << endl;

				VectorXf e(m);
				float t = 0;
				e = B.col(min);
				B.col(min) = N.col(i);
				N.col(i) = e;
				t = c_B(min);
				c_B(min) = c_N(i);
				c_N(i) = t;
				for (int k = 0; k < m; k++)
				{
					if (k != min)
					{
						b(k) = b(k) + minOfstep*d_q(k);
					}
					else
					{
						b(k) = minOfstep;
					}
				}
				//b(min) = r(i);

				//cout << B << endl;
				//cout << N << endl
				break;


				//break;
			}
		}
	}
	

	clock_t tEnd = clock();        
	cout << "程序总运行时间(s)：" << float(tEnd - tStart) / CLOCKS_PER_SEC << endl;

	system("pause");
	return 0;
}