線性規劃專題——SIMPLEX 單純形演算法(四)——實現
阿新 • • 發佈:2018-12-17
實現的原理見:
線性規劃專題——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; }