C++实现单纯性算法

这篇文章是针对我前段时间提出的问题——(4条消息) 实现单纯性算法时遇到的精度问题-编程语言-CSDN问答

的一个最终回答，实现了大M法和单纯性算法

#include <iostream>
#include <iomanip>
#include <vector>
#include <cmath>

using std::vector;
using std::cin;
using std::cout;
using std::endl;

#define M 1e+9

/**
 * 1. func scale: to scale a row, make the pivot variable 1;
 * 2. func subtract: to subtract one row from another in order to make the main variable on other lines 0;
**/
void scale(vector<vector<double>>& matrix, int row_index, double scale_num)
{
    // scale_num is  the value of the pivot variable
    for (int i = 0; i < matrix[row_index].size(); i++)
    {
        matrix[row_index].at(i) /= scale_num;
    }
}

void substract(vector<vector<double>>& matrix, int row_index_1, int pivot_row, double muti)
{
    /**
     * row_index_1 is some non-pivot row;
     * row_index_2 is the pivot row;
    **/
    
    for (int i = 0; i < matrix[row_index_1].size(); i++)
    {
        matrix[row_index_1].at(i) -= matrix[pivot_row].at(i) * muti;
    }
    
}

void gaussian(vector<vector<double>>& simplex_matrix, int pivot_column, int pivot_row)
{
    /*this function should make the pivot column value in other rows 0*/
    scale(simplex_matrix, pivot_row, simplex_matrix[pivot_row][pivot_column]);
    
    for (int pres_row = 0; pres_row < simplex_matrix.size(); pres_row++)
    {
        if(pres_row != pivot_row) 
        {
            double muti = simplex_matrix[pres_row][pivot_column] / simplex_matrix[pivot_row][pivot_column];
            substract(simplex_matrix, pres_row, pivot_row, muti);
        }
    }
}

/**
 * 1. input: input and store the constraints
 * 2. augment: augment the previous matrix, store the augmented form of the matrix(simplex matrix), initliaze the basis and variable values
**/
void input(vector<vector<double>>& simplex_matrix, vector<bool>& need_artificial_variable, int variable_num, int constraint_num)
{
    int RHS_index = simplex_matrix[0].size() - 1;
    int objective_function_index = simplex_matrix.size() - 1;

    /*input the constraint*/
    //cout<<"input constraint"<<endl;
    for (size_t i = 0; i < constraint_num; i++)
    {
        for (size_t j = 0; j < variable_num; j++)
        {
            cin>>simplex_matrix[i][j];
        }
    }
    
    /*input the RHS values*/
    //cout<<"input RHS values"<<endl;
    for (size_t i = 0; i < constraint_num; i++)
    {
        cin>>simplex_matrix[i][RHS_index];
        if(simplex_matrix[i][RHS_index] < 0) //If the RHS_value is negative, we need to transform it into a positive value.
        {
            for (size_t j = 0; j < simplex_matrix[i].size(); j++)
            {
                simplex_matrix[i][j] *= -1;
            }
            need_artificial_variable[i] = true;
        }
    }
    
    /*input the objective function*/
    //cout<<"input objective function"<<endl;
    for (size_t i = 0; i < variable_num; i++)
    {
        int temp;
        cin>>temp;
        simplex_matrix[objective_function_index][i] = -temp;
    }
}

void augment(vector<vector<double>>& simplex_matrix, vector<bool>& need_artificial_variable, vector<int>& basis, vector<double>& variable_values, int variable_num, int constratint_num)
{
    int RHS_index = simplex_matrix[0].size() - 1;
    int objective_function_index = simplex_matrix.size() - 1;

    /*add the slack variables*/
    int slack_index = variable_num;
    for (size_t i = 0; i < constratint_num; i++)
    {
        if(!need_artificial_variable[i])
        {
            simplex_matrix[i][slack_index] = 1;
        }
        else
        {
            simplex_matrix[i][slack_index] = -1;
        }
        basis[i] = slack_index;//The slack variable is going to be the basis, if there is no artificial variable in the constraint.
        slack_index++;
    }

    /*add the artificial variables*/
    int artificial_index = slack_index;
    for (size_t i = 0; i < constratint_num; i++)
    {
        if(need_artificial_variable[i])
        {
            simplex_matrix[i][artificial_index] = 1;
            simplex_matrix[objective_function_index][artificial_index] = M;
            basis[i] = artificial_index; //update the basis, if there is an artificial variable in the constraint, it is always the initial basis.
            /*Eliminate the coefficents in the objective function*/
            substract(simplex_matrix, objective_function_index, i, M);
            artificial_index++;
        }
    }

    /*initliaze the basis variable values*/
    for (size_t i = 0; i < basis.size(); i++)
    {
        variable_values[basis[i]] = simplex_matrix[i][RHS_index];
    }
}

/**
 * 1. test_optimal: test whether the solution is optimal
 * 2. test_unbounded: test whether the problem is unbounded
 * 3. test_feasible: test whether the problem has no feasiable solution
**/
bool test_optimal(vector<vector<double>>& simplex_matrix)
{
    /*checks the objective function*/
    int positive_counter = 0;
    int RHS_index = simplex_matrix[0].size() - 1;
    int objective_function_index = simplex_matrix.size() - 1;

    for (size_t i = 0; i < simplex_matrix[objective_function_index].size() - 1; i++)
    {
        if (simplex_matrix[objective_function_index][i] >= -1e-5)
        {
            positive_counter++;
        }
    }
    
    if (positive_counter == simplex_matrix[objective_function_index].size() - 1)
    {
        return true;
    }
    
    return false;
}

bool test_unbounded(vector<vector<double>>& simplex_matrix, int pivot_column)
{
    int negative_counter = 0;
    int RHS_index = simplex_matrix[0].size() - 1;
    int objective_function_index = simplex_matrix.size() - 1;

    for (size_t i = 0; i < simplex_matrix.size() - 1; i++)
    {
        if (simplex_matrix[i][pivot_column] <= 0)
        {
            negative_counter++;
        }
    }
    
    if (negative_counter == simplex_matrix.size() - 1)
    {
        return true;
    }
    return false;

}

bool test_fasiable(vector<int>& basis, vector<double>& variable_values, int variable_num, int constraint_num)
{
    int slack_num = variable_num + constraint_num;
    for (size_t i = 0; i < basis.size(); i++)
    {
        if (basis[i] >= slack_num && std::abs(variable_values[basis[i]]) >= 1e-5)
        {
            return true;
        }
    }
    
    return false;
}

/**
 * 1. find_in_variable: find the in-basis variable(which way to go)
 * 2. find_out_variable: find the out-basis variable(how far to go)
**/
int find_in_variable(vector<vector<double>>& simplex_matrix, vector<int>& basis, vector<double>& variable_values, int variable_num, int constraint_num)
{
    /*find the non-basis variable with the minimum coefficent*/
    int RHS_index = simplex_matrix[0].size() - 1;
    int objective_function_index = simplex_matrix.size() - 1;

    double min_value = 1e+9;
    int min_index = -1;
    for (size_t i = 0; i < simplex_matrix[objective_function_index].size() - 1; i++)
    {
        if(simplex_matrix[objective_function_index][i] < min_value)
        {
            bool in_basis = false;
            for (size_t j = 0; j < basis.size(); j++)
            {
                if (basis[j] == i)
                {
                    in_basis = true;
                    break;
                }
            }

            if(!in_basis && simplex_matrix[objective_function_index][i] < 0)
            {
                min_value = simplex_matrix[objective_function_index][i];
                min_index = i;
            }
        }
    }
    
    return min_index;
}

int find_out_variable(vector<vector<double>>& simplex_matrix, vector<int>& basis, vector<double>& variable_values, int pivot_column, int variable_num, int constraint_num)
{
    /*find the minium ratio in the pivot column*/
    int RHS_index = simplex_matrix[0].size() - 1;
    int objective_function_index = simplex_matrix.size() - 1;

    double min_ratio = 1e+9;
    int min_index = -1;
    for (size_t i = 0; i < simplex_matrix.size() - 1; i++)
    {
        if(simplex_matrix[i][pivot_column] <= 0)
        {
            continue;
        }
        double ratio = simplex_matrix[i][RHS_index] / simplex_matrix[i][pivot_column];
        //cout<<"ratio:"<<ratio<<endl;
        if (ratio < min_ratio)
        {
            //cout<<"min!"<<endl;
            min_ratio = ratio;
            min_index = i;
        }
    }
    
    return min_index;
}

int simplex(vector<vector<double>>& simplex_matrix, vector<int>& basis, vector<double>& variable_values, int variable_num, int constraint_num)
{
    int RHS_index = simplex_matrix[0].size() - 1;
    int objective_function_index = simplex_matrix.size() - 1;
    int iter_counter = 0;

    while (true)
    {
        /*cout<<"----------------------------------------------------------------------------------------------------------"<<endl;
        for (size_t i = 0; i < simplex_matrix.size(); i++)
        {
            for (size_t j = 0; j < simplex_matrix[i].size(); j++)
            {
                cout<<std::setw(15)<<simplex_matrix[i][j];
            }
            cout<<endl;
        }
        cout<<"basis: ";
        for (size_t i = 0; i < basis.size(); i++)
        {
            cout<<basis[i]<<' ';
        }
        cout<<endl;

        cout<<"variable_values: ";
        for (size_t i = 0; i < variable_values.size(); i++)
        {
            cout<<variable_values[i]<<' ';
        }
        cout<<endl;*/
        
        if (test_optimal(simplex_matrix))
        {
            return 1;
        }
        int in_variable = find_in_variable(simplex_matrix, basis, variable_values, variable_num, constraint_num);
        //cout<<"in variable: "<<in_variable<<endl;

        if (test_unbounded(simplex_matrix, in_variable))
        {
            return 2;
        }
        int out_variable = find_out_variable(simplex_matrix, basis, variable_values, in_variable, variable_num, constraint_num);
        //cout<<"out variable: "<<out_variable<<endl;

        gaussian(simplex_matrix, in_variable, out_variable);

        variable_values[basis[out_variable]] = 0;
        basis[out_variable] = in_variable;
        for (size_t i = 0; i < basis.size(); i++)
        {
            variable_values[basis[i]] = simplex_matrix[i][RHS_index] / simplex_matrix[i][basis[i]];
        }
    }

}

int main()
{
    int variable_num, constraint_num;
    cin>>constraint_num>>variable_num;

    vector<vector<double>> simplex_matrix(constraint_num + 1, vector<double>(variable_num + constraint_num*2 + 1));
    vector<int> basis(constraint_num); //The current basis
    vector<bool> need_artificial_variable(constraint_num, false);
    vector<double> variable_values(variable_num + constraint_num*2);
    int RHS_index = simplex_matrix[0].size() - 1;
    int objective_function_index = simplex_matrix.size() - 1;

    input(simplex_matrix, need_artificial_variable, variable_num, constraint_num);
    augment(simplex_matrix, need_artificial_variable, basis, variable_values, variable_num, constraint_num);
    
    int return_code = simplex(simplex_matrix, basis, variable_values, variable_num, constraint_num);
    for (size_t i = 0; i < basis.size(); i++)
    {
        variable_values[basis[i]] = simplex_matrix[i][RHS_index];
    }
    /*for (int i = 0; i < variable_num; i++)
    {
        cout<<std::fixed<<std::showpoint<<std::setprecision(15)<<variable_values[i]<<' ';
    }
    cout<<endl;*/

    if (test_fasiable(basis, variable_values, variable_num, constraint_num))
    {
        return_code = 3;
    }

    switch (return_code)
    {
    case 1:
        cout<<"Bounded solution"<<endl;
        
        for (int i = 0; i < variable_num; i++)
        {
            cout<<std::fixed<<std::showpoint<<std::setprecision(15)<<variable_values[i]<<' ';
        }
        cout<<endl;
        break;    
    case 2:
        cout<<"Infinity"<<endl;
        break;
    case 3:
        cout<<"No solution"<<endl;
        break;
    default:
        break;
    }

    return 0;   
}

代码如上