1. 程式人生 > >線性規劃專題——SIMPLEX 單純形演算法(四)——實現

線性規劃專題——SIMPLEX 單純形演算法(四)——實現

實現的原理見:
線性規劃專題——SIMPLEX 單純形演算法(一)
線性規劃專題——SIMPLEX 單純形演算法(二)
線性規劃專題——SIMPLEX 單純形演算法(三)——示例、注意點

實現過程

重點是實現PIVOT 指定主元的高斯約旦消元,以及INITICIALSIMPLEX初始化SIMPLEX,注意 在實現中,本人並沒有構造輔助線性規劃,因為看起來比較麻煩。

#include <stdio.h>
#include <stdlib.h>
#include <malloc.h>
#include <vector>
#include <algorithm>
#include <set>
#include <iostream>
using namespace std;
#define N (202)
#define INFTY (1e30)
#define EXP (1e-10)
//INFIY 正無窮
//EXP	無窮小量,主要為了對付double 在判斷時可能出現誤差
double co_matrix[N][N] = { 0 };
double b_matrix[N]= { 0 };
double c_matrix[N] = { 0 };
set<int> BI;//儲存基本列的下標
double z;//儲存最小值
double theta[N] = { 0 };
double lambda[N] = { 0 };
int INITLIZESIMPLEX(int m,int n);
vector<double > CALCULATEX(int m,int n);
void PIVOT(int i,int j,int m,int n);//以co_matrix的aij為主元,執行一次高斯約旦消元
int find_min_index(double * vec,int len)
{
	double tmp = vec[0];
	int index = 0;
	for (int i = 0; i < len; i++)
	{
		if (tmp>vec[i])
		{
			tmp = vec[i], index = i;//選第1個最小的
		}
	}
	return index;
}
int find_max_index(double * vec, int len)
{
	double tmp = vec[0];
	int index = 0;
	for (int i = 0; i < len; i++)
	{
		if (tmp<vec[i])
		{
			tmp = vec[i], index = i;//選第1個最大的
		}
	}
	return index;
}
pair < vector<double>, double > SIMPLEX(int m, int n, int init_flag = 0)
{
	vector<double> x;
	int feasible = INITLIZESIMPLEX(m, n);//可能有解的時候
	while (feasible)
	{
		int min_ci = find_min_index(c_matrix, n);
		if (c_matrix[min_ci] >= 0)
			//所有的ci都大於0,說明目標函式已經收斂了 找到了最優值
		{
			x = CALCULATEX(m,n);
			break;
		}
		//針對第min_ci列 計算最小的非負theta
		for (int i = 0; i < m; i++)
		{
			if (co_matrix[i][min_ci]>0)
				//這樣計算的theta才是非負的
			{
				theta[i] = b_matrix[i] / co_matrix[i][min_ci];
			}
			else
				//當co_matrix[i][min_ci]為0,或者為負的時候,都不用去考慮
			{
				theta[i] = INFINITY;
			}
		}
		int max_thetai = find_min_index(theta, m);
		if (abs(theta[max_thetai] - INFTY) < EXP)
		{
			x.clear();
			z = INFTY;
			printf("unbounded");
			break;
		}
		PIVOT(max_thetai, min_ci,m,n);//做一次高斯消元
	}
	return pair<vector<double> ,double> (x, -z);
}
int INITLIZESIMPLEX(int m,int n)
{
	int min_bi = find_min_index(b_matrix, m);//找到bi最小的那個

	BI.clear();
	for (int i = m; i < n; i++)
	{
		BI.insert(i);
	}

	if (b_matrix[min_bi] >= 0)
		//左右的b都是非負的,那麼這個問題就可以自然的獲取一組頂點和解
	{
		return 1;
	}
	//使用黑科技,直接從bi最小的那一行 任選一個 小於0 的co_matrix[i]來操作
	//並沒有構建輔助線性規劃,因為太麻煩了
	int col_index = -1;
	for (int i = 0; i < n; i++)
	{
		if (co_matrix[min_bi][i] < 0)
		{
			col_index = i;
			break;
		}
	}
	if (col_index == -1)
	{
		BI.clear();
		printf("Infeasible");
		return 0;//無解
	}
	//否則執行一次轉軸操作
	PIVOT(min_bi, col_index,m,n);
	return 1;
}
vector<double > CALCULATEX(int m,int n)
{
	vector<double> x;
	for (int i = 0; i < n; i++)
	{
		x.push_back(0);
	}
	for (auto it = BI.begin(); it != BI.end();it++)
	{
		int i = *it;
		{
			for (int j = 0; j < m; j++)
			{
				if (abs(co_matrix[j][i] -1 ) < EXP)
				{
					x[i] = b_matrix[j];
					if (x[i] < 0)
					{
						x.clear();
						printf("Infeasible");
						break;
						//無解
					}
				}
			}
		}
	}
	return x;
}
void PIVOT(int row, int col,int m,int n)
{
	//基的變換
	for (auto it = BI.begin(); it != BI.end(); it++)
		//從已有的基裡面找到那個需要被替換出去的基
	{
		if (abs(co_matrix[row][*it] - 1) < EXP)
		{
			BI.erase(*it);
			break;
		}
	}
	BI.insert(col);

	//係數矩陣的高斯約旦消元
	for (int i = 0; i < n; i++)
		//把當前行除以co_matrix[row][col]
	{
		if (i != col)
		{
			co_matrix[row][i] /= co_matrix[row][col];
		}
	}
	b_matrix[row] /= co_matrix[row][col];
	co_matrix[row][col] = 1;

	//把其他行 減去第row行,校區第row行的第col個元素,以及b_matrix的高斯約旦
	for (int i = 0; i < m; i++)
	{
		if (i != row)
		{
			double factor = co_matrix[i][col] / co_matrix[row][col];
			for (int j = 0; j < n; j++)
			{
				co_matrix[i][j] -= co_matrix[row][j] * factor;
			}
			b_matrix[i] -= b_matrix[row] * factor;
		}
	}
	//c_matrix的高斯約旦
	double factor = c_matrix[col] / co_matrix[row][col];
	for (int j = 0; j < n; j++)
	{
		c_matrix[j] -= co_matrix[row][j] * factor;
	}
	//z的變換
	z -= b_matrix[row] * factor;



}
int main()
{
	//TEST 1.
	/*int m = 4, n = 7;
	co_matrix[0][0] = 1, co_matrix[0][1] = 1, co_matrix[0][2] = 1, co_matrix[0][3] = 1, co_matrix[0][4] = 0, co_matrix[0][5] = 0, co_matrix[0][6] = 0;
	co_matrix[1][0] = 1, co_matrix[1][1] = 0, co_matrix[1][2] = 0, co_matrix[1][3] = 0, co_matrix[1][4] = 1, co_matrix[1][5] = 0, co_matrix[1][6] = 0;
	co_matrix[2][0] = 0, co_matrix[2][1] = 0, co_matrix[2][2] = 1, co_matrix[2][3] = 0, co_matrix[2][4] = 0, co_matrix[2][5] = 1, co_matrix[2][6] = 0;
	co_matrix[3][0] = 0, co_matrix[3][1] = 3, co_matrix[3][2] = 1, co_matrix[3][3] = 0, co_matrix[3][4] = 0, co_matrix[3][5] = 0, co_matrix[3][6] = 1;
	b_matrix[0] = 4, b_matrix[1] = 2, b_matrix[2] = 3, b_matrix[3] = 6;
	c_matrix[0] = -1, c_matrix[1] = -14, c_matrix[2] = -6, c_matrix[3]= 0, c_matrix[4]= 0, c_matrix[5] = 0, c_matrix[6] = 0;*/

	//TEST 2.
	/*int m = 2, n = 4;
	co_matrix[0][0] = -1, co_matrix[0][1] = -1, co_matrix[0][2] = 1, co_matrix[0][3] = 0;
	co_matrix[1][0] = 1, co_matrix[1][1] = 1, co_matrix[1][2] = 0, co_matrix[1][3] = 1;
	b_matrix[0] = -1, b_matrix[1] = 2;
	c_matrix[0] = 1, c_matrix[1]= 2, c_matrix[2] = 0, c_matrix[3]= 0;
	*/
	//TEST 3.無解
	int m = 2, n = 4;
	co_matrix[0][0] = -1, co_matrix[0][1] = -1, co_matrix[0][2] = 1, co_matrix[0][3] = 0;
	co_matrix[1][0] = 1, co_matrix[1][1] = 1, co_matrix[1][2] = 0, co_matrix[1][3] = 1;
	b_matrix[0] = -2, b_matrix[1] = -1;
	c_matrix[0] = 1, c_matrix[1] = 2, c_matrix[2] = 0, c_matrix[3] = 0;

	auto rst = SIMPLEX(m, n);
	if (rst.first.size())
		printf("min:%0.5f\n", rst.second);

	system("pause");
	return 0;
}