线性规划求解器

.gitignore

@@ -0,0 +1 @@
+build

CMakeLists.txt

@@ -0,0 +1,36 @@
+cmake_minimum_required(VERSION 3.16.0)
+
+project(Solver)
+
+set(CMAKE_CXX_STANDARD 17)
+set(CMAKE_CXX_STANDARD_REQUIRED ON)
+
+# include(./FindGUROBI.cmake)
+# include_directories(${GUROBI_INCLUDE_DIRS})
+# include_directories(${CMAKE_CURRENT_SOURCE_DIR}/inc)
+
+set(SOURCE_DIR "${PROJECT_SOURCE_DIR}/src")
+set(INCLUDE_DIR "${PROJECT_SOURCE_DIR}/inc")
+
+file(GLOB SOURCES "${SOURCE_DIR}/*.cpp$")
+file(GLOB HEADERS "${SOURCE_INCLUDE}/*.hpp$")
+
+#include_directories(
+# ${GUROBI_INCLUDE_DIRS}
+#)
+
+message(${SOURCE_DIR})
+
+aux_source_directory(${SOURCE_DIR} src)
+aux_source_directory(${INCLUDE_DIR} inc)
+
+
+add_executable(Solver
+   src/main.cpp src/solver.cpp inc/solver.hpp
+  )
+
+#target_link_libraries(Solver 
+# ${GUROBI_LIBRARIES}
+#)
+
+target_include_directories(Solver PRIVATE ${SOURCE_DIR} ${INCLUDE_DIR})

FindGUROBI.cmake

@@ -0,0 +1,73 @@
+# ------------------------------------------------------------------------------
+# This file sets up Gurobi for CMake. Once done this will define
+#
+#   GUROBI_FOUND           - system has GUROBI
+#   GUROBI_INCLUDE_DIRS    - the GUROBI include directories
+#   GUROBI_LIBRARIES       - Link these to use GUROBI
+#
+#  In your CMakeLists file, you need to add, e.g. (modify it if necessary):
+#        if (GUROBI_FOUND)
+#            message(STATUS "Gurobi include dir: " ${GUROBI_INCLUDE_DIRS})
+#            message(STATUS "Gurobi libraries: " ${GUROBI_LIBRARIES})
+#            target_compile_definitions(${PROJECT_NAME} PUBLIC HAS_GUROBI)
+#            target_include_directories(${PROJECT_NAME} PRIVATE ${GUROBI_INCLUDE_DIRS})
+#            target_link_libraries(${PROJECT_NAME} PRIVATE ${GUROBI_LIBRARIES})
+#        endif()
+# ------------------------------------------------------------------------------
+
+
+# Is it already configured?
+if (NOT GUROBI_FOUND)
+
+    # Hardcoded search paths
+    set(SEARCH_PATHS_FOR_HEADERS
+        "$ENV{GUROBI_HOME}/include"
+        "/Library/gurobi1001/mac64/include"
+        "C:\\gurobi1001\\win64\\include"
+    )
+
+    set(SEARCH_PATHS_FOR_LIBRARIES
+        "$ENV{GUROBI_HOME}/lib"
+        "/Library/gurobi1001/mac64/lib"
+        "C:\\gurobi1001\\win64\\lib"
+    )
+
+    find_path(GUROBI_INCLUDE_DIR gurobi_c++.h
+      PATHS ${SEARCH_PATHS_FOR_HEADERS}
+    )
+
+
+    find_library( GUROBI_C_LIBRARY
+                  NAMES gurobi100 gurobi105 libgurobi
+                  PATHS ${SEARCH_PATHS_FOR_LIBRARIES}
+                  )
+
+    find_library( GUROBI_CXX_LIBRARY_DEBUG
+                NAMES gurobi_c++ gurobi_c++mdd2017
+                PATHS ${SEARCH_PATHS_FOR_LIBRARIES}
+                )
+
+    find_library( GUROBI_CXX_LIBRARY_RELEASE
+                NAMES gurobi_c++ gurobi_c++md2017
+                PATHS ${SEARCH_PATHS_FOR_LIBRARIES}
+                )
+
+    # setup header file directories
+    set(GUROBI_INCLUDE_DIRS ${GUROBI_INCLUDE_DIR})
+
+    # setup libraries files
+    set(GUROBI_LIBRARIES
+            debug ${GUROBI_CXX_LIBRARY_DEBUG}
+            optimized ${GUROBI_CXX_LIBRARY_RELEASE}
+            ${GUROBI_C_LIBRARY}
+            )
+
+endif()
+
+# Check that Gurobi was successfully found
+include(FindPackageHandleStandardArgs)
+find_package_handle_standard_args(GUROBI DEFAULT_MSG GUROBI_INCLUDE_DIRS)
+find_package_handle_standard_args(GUROBI DEFAULT_MSG GUROBI_LIBRARIES)
+
+# Hide variables from CMake-Gui options
+mark_as_advanced(GUROBI_LIBRARIES GUROBI_INCLUDE_DIRS GUROBI_INCLUDE_DIR)

inc/solver.hpp

@@ -0,0 +1,123 @@
+#pragma once
+#include <vector>
+using std::vector;
+
+class Expr;
+class Var;
+
+Expr operator+(const Expr& x, const Expr& y);
+Expr operator-(const Expr& x, const Expr& y);
+Expr operator+(const Expr& x);
+Expr operator+(Var x, Var y);
+Expr operator+(Var x, double a);
+Expr operator+(double a, Var x);
+Expr operator-(const Expr& x);
+Expr operator-(Var x);
+Expr operator-(Var x, Var y);
+Expr operator-(Var x, double a);
+Expr operator-(double a, Var x);
+Expr operator*(double a, Var x);
+Expr operator*(Var x, double a);
+Expr operator*(const Expr& x, double a);
+Expr operator*(double a, const Expr& x);
+Expr operator/(Var x, double a);
+Expr operator/(const Expr& x, double a);
+
+typedef vector<vector<double>> Matrix;
+
+
+enum class Rtn {
+	LOADED,
+	OPTIMAL,
+	INFEASIBLE,
+	INF_OR_UNBD,
+	UNBOUNDED,
+	CUTOFF,
+	ITERATION_LIMIT,
+	NODE_LIMIT,
+	TIME_LIMIT,
+	SOLUTION_LIMIT,
+	INTERRUPTED,
+	NUMERIC,
+	SUBOPTIMAL,
+	INPROGRESS,
+	USER_OBJ_LIMIT,
+	WORK_LIMIT,
+	MEM_LIMIT,
+	PIVOT,
+};
+
+class Var {
+public:
+	Var(double coef = 1) :col(0), coeffs(coef) {};
+	double coeffs;
+	int col;
+};
+
+class Expr
+{
+private:
+	double constant;
+	std::vector<double> coeffs;
+	std::vector<Var> vars;
+
+public:
+	Expr(const Expr& expr) = default;
+	Expr(double constant = 0.0);
+	Expr(Var var, double coeff = 1.0);
+
+	friend class Model;
+
+	friend Expr operator+(const Expr& x, const Expr& y);
+	friend Expr operator+(const Expr& x);
+	friend Expr operator+(Var x, Var y);
+	friend Expr operator+(Var x, double a);
+	friend Expr operator+(double a, Var x);
+	friend Expr operator-(const Expr& x, const Expr& y);
+	friend Expr operator-(const Expr& x);
+	friend Expr operator-(Var x);
+	friend Expr operator-(Var x, Var y);
+	friend Expr operator-(Var x, double a);
+	friend Expr operator-(double a, Var x);
+	friend Expr operator*(double a, Var x);
+	friend Expr operator*(Var x, double a);
+	friend Expr operator*(const Expr& x, double a);
+	friend Expr operator*(double a, const Expr& x);
+	friend Expr operator/(Var x, double a);
+	friend Expr operator/(const Expr& x, double a);
+
+
+	Expr operator=(const Expr& rhs);
+	void operator+=(const Expr& expr);
+	void operator-=(const Expr& expr);
+	void operator*=(double mult);
+	void operator/=(double a);
+	Expr operator+(const Expr& rhs);
+	Expr operator-(const Expr& rhs);
+};
+
+class Model
+{
+public:
+	Model();
+	Var* addVars(int col);
+	void addConstr(const Expr& expr, char sense, double rhs);
+	void setObjective(Expr obje, int sense = 0);
+	void print();
+	Rtn optimize(bool isMax = true);
+	double get();
+
+private:
+	double _simplex();
+	Rtn _pivot(std::pair<size_t, size_t>& p);
+	void _gaussian(std::pair<size_t, size_t> p);
+
+
+	Var* vars;
+
+	Matrix matrix;
+	size_t cn, bn;
+	vector<int> basic;
+	Rtn rtn;
+	double res;
+};

src/main.cpp

@@ -0,0 +1,35 @@
+#include <iostream>
+#include <cassert>
+
+#include "solver.hpp"
+
+using namespace std;
+
+int main(int argc, char* argv[])
+{
+	Model mdl;
+
+	Var* vars = mdl.addVars(3);
+
+	mdl.addConstr(vars[0] + vars[1] + vars[2], '<', 4);
+	mdl.addConstr(vars[0], '<', 2);
+	mdl.addConstr(vars[2], '<', 3);
+	mdl.addConstr(3 * vars[1] + vars[2], '<', 6);
+
+	mdl.setObjective(vars[0] + 14 * vars[1] + 6 * vars[2]);
+
+	switch(mdl.optimize(false)) {
+	case Rtn::OPTIMAL:			
+		cout << "find solution: " << mdl.get() << endl; 
+		break;
+	case Rtn::INFEASIBLE:	
+		cout << "no solution" << endl;					
+		break;
+	case Rtn::UNBOUNDED:	
+		cout << "unbounded solution" << endl;			
+		break;
+	default:
+		assert(false);
+	}
+	return 0;
+}

src/solver.cpp

@@ -0,0 +1,424 @@
+#include "solver.hpp"
+#include <iostream>
+#include <cassert>
+#include <algorithm>
+
+using std::make_pair;
+using std::pair;
+using std::cout;
+using std::endl;
+
+Expr operator+(const Expr& x, const Expr& y)
+{
+	Expr exp;
+	exp.coeffs.resize(std::max(x.coeffs.size(), y.coeffs.size()), 0);
+	for (int c = 0; c < exp.coeffs.size(); c++) {
+		if (c < x.coeffs.size() && c < y.coeffs.size()) {
+			exp.coeffs[c] = x.coeffs[c] + y.coeffs[c];
+		}
+		else {
+			exp.coeffs[c]= c < x.coeffs.size()? x.coeffs[c] : y.coeffs[c];
+		}
+	}
+	exp.constant = x.constant + y.constant;
+	return exp;
+}
+
+Expr operator+(const Expr& x)
+{
+	return x;
+}
+
+Expr operator+(Var x, Var y)
+{
+	Expr exp;
+	exp.coeffs.resize(std::max(x.col + 1, y.col + 1), 0);
+	exp.coeffs[x.col] = x.coeffs;
+	exp.coeffs[y.col] = y.coeffs;
+	return exp;
+}
+
+Expr operator+(Var x, double a)
+{
+	Expr exp;
+	exp.coeffs.resize(x.col + 1);
+	exp.constant = a;
+	return exp;
+}
+
+Expr operator+(double a, Var x)
+{
+	return x + a;
+}
+
+Expr operator-(const Expr& x, const Expr& y)
+{
+	Expr exp;
+	exp.coeffs.resize(std::max(x.coeffs.size(), y.coeffs.size()), 0);
+	for (int c = 0; c < exp.coeffs.size(); c++) {
+		if (c < x.coeffs.size() && c < y.coeffs.size()) {
+			exp.coeffs[c] = x.coeffs[c] - y.coeffs[c];
+		}
+		else {
+			exp.coeffs[c] = c < x.coeffs.size() ? x.coeffs[c] : y.coeffs[c];
+		}
+	}
+	exp.constant = x.constant + y.constant;
+	return exp;
+}
+
+Expr operator-(const Expr& x)
+{
+	Expr expr(x);
+	for (int c = 0; c < expr.coeffs.size(); c++) {
+		expr.coeffs[c] = -expr.coeffs[c];
+	}
+	expr.constant = -expr.constant;
+	return expr;
+}
+
+Expr operator-(Var x)
+{
+	return -Expr(x);
+}
+
+Expr operator-(Var x, Var y)
+{
+	Expr exp;
+	exp.coeffs.resize(std::max(x.col + 1, y.col + 1), 0);
+	exp.coeffs[x.col] = x.coeffs;
+	exp.coeffs[y.col] = -y.coeffs;
+	return exp;
+}
+
+Expr operator-(Var x, double a)
+{
+	return x - Var(a);
+}
+
+Expr operator-(double a, Var x)
+{
+	return x - a;
+}
+
+Expr operator*(double a, Var x)
+{
+	Expr exp;
+	exp.coeffs.resize(x.col + 1, 0);
+	exp.coeffs[x.col] = a * x.coeffs;
+	return exp;
+}
+
+Expr operator*(Var x, double a)
+{
+	return a * x;
+}
+
+Expr operator*(const Expr& x, double a)
+{
+	Expr exp = x;
+	for (int c = 0; c < exp.coeffs.size(); c++) {
+		exp.coeffs[c] *= a;
+	}
+	exp.constant *= a;
+	return exp;
+}
+
+Expr operator*(double a, const Expr& x)
+{
+	return x * a;
+}
+
+Expr operator/(Var x, double a)
+{
+	return Expr(x) / a;
+}
+
+Expr operator/(const Expr& x, double a)
+{
+	Expr exp = x;
+	for (int c = 0; c < exp.coeffs.size(); c++) {
+		exp.coeffs[c] /= a;
+	}
+	exp.constant /= a;
+	return exp;
+}
+
+Model::Model():
+	cn(0),
+	bn(0),
+	vars(nullptr)
+{
+
+}
+
+Var* Model::addVars(int num)
+{
+	Var* old = vars;
+	vars = new Var[cn + num];
+	for (int c = cn; c < cn + num; c++) {
+		vars[c].col = c;
+	}
+	memcpy(vars, old, sizeof(Var) * cn);
+	delete[] old;
+	cn += num;
+	return vars;
+}
+
+void Model::addConstr(const Expr& expr, char sense, double rhs)
+{
+	matrix.push_back(vector<double>(1, rhs));			// TODO<44><4F>û<EFBFBD>ȴ<EFBFBD>С<EFBFBD><D0A1>û<EFBFBD><C3BB><EFBFBD>dz<EFBFBD><C7B3><EFBFBD>
+	matrix.back().insert(matrix.back().end(), expr.coeffs.begin(), expr.coeffs.end());
+	for (int c = matrix.back().size(); c <= cn; c++) {
+		matrix.back().push_back(0);
+	}
+}
+
+void Model::setObjective(Expr obje, int sense)
+{
+	matrix.insert(matrix.begin(), obje.coeffs);
+	matrix.front().insert(matrix.front().begin(), obje.constant);		// TODO<44><4F>д<EFBFBD>IJ<EFBFBD><C4B2><EFBFBD><EFBFBD><EFBFBD>
+	for (int c = matrix.front().size(); c <= cn; c++) {
+		matrix.front().push_back(0);
+	}
+}
+
+void Model::print()
+{
+	for (size_t i = 0; i < matrix.size(); i++) {
+		for (size_t j = 0; j < matrix[0].size(); j++) {
+			cout << matrix[i][j] << "\t";
+		}
+		cout << endl;
+	}
+}
+
+Rtn Model::optimize(bool isMax)
+{
+	bn = matrix.size();
+	for (int row = 1; row < bn; row++) {
+		matrix.front().push_back(0);
+		for (int col = 1; col < bn; col++) {
+			matrix[row].push_back(col == row ? 1 : 0);
+		}
+	}
+	cn = matrix.front().size();
+	if (!isMax) {
+		for (int i = 0; i < cn; i++) {
+			matrix.front()[i] = -matrix.front()[i];
+		}
+	}
+
+	for (size_t i = 1; i < bn; i++) {
+		basic.push_back(cn - bn + i);
+	}
+
+	// === <20>жϳ<D0B6>ʼ<EFBFBD><CABC><EFBFBD>Ƿ<EFBFBD>Ϊ<EFBFBD><CEAA><EFBFBD>н<EFBFBD> ===
+	bool initial_feasible = true;
+	for (int row = 1; row < bn; row++) {
+		if (matrix[row].front() < 0) {
+			initial_feasible = false;
+			break;
+		}
+	}
+
+	// === <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ʼ<EFBFBD><CABC><EFBFBD>н<EFBFBD> === 
+	if (!initial_feasible) {
+		vector<double> coeff = matrix.front();
+		matrix.front() = vector<double>(cn, .0);
+		matrix.front().push_back(1);
+		pair<size_t, size_t> t = { -1 ,cn };
+
+		for (int row = 1; row < bn; row++) {
+			matrix[row].push_back(-1);
+			if (t.first == -1 || matrix[row].front() < matrix[t.first].front()) {
+				t.first = row;
+			}
+		}
+
+		_gaussian(t);
+		if (fabs(_simplex()) > 1e-10) {
+			rtn = Rtn::INFEASIBLE;
+			return rtn;
+		}
+
+		// if the x0 in B, we should pivot it.
+		auto iter = find(basic.begin(), basic.end(), cn);
+		if (iter != basic.end()) {
+			for (int col = 1; col < matrix.front().size(); col++) {
+				if (matrix.front()[col] != 1) {
+					t = make_pair(iter - basic.begin(), col);
+					_pivot(t);
+					break;
+				}
+			}
+		}
+
+		for (int row = 0; row < bn; row++) {
+			matrix[row].pop_back();
+		}
+
+		// recover the coefficient line
+		for (int col = 0; col < cn; col++) {
+			matrix.front()[col] = coeff[col];
+		}
+		for (int row = 1; row <= basic.size(); row++) {
+			int norm = matrix.front()[basic[row - 1]];
+			for (int col = 0; col < cn; col++) {
+				matrix.front()[col] -= norm * matrix[row][col];
+			}
+		}
+	}
+
+	res = isMax ? _simplex() : _simplex();
+	return rtn;
+}
+
+double Model::get()
+{
+	return res;
+}
+
+double Model::_simplex()
+{
+	pair<size_t, size_t> t;
+	while (1) {
+		rtn = _pivot(t);
+		if (rtn == Rtn::OPTIMAL || rtn == Rtn::UNBOUNDED) {
+			break;
+		}
+		_gaussian(t);
+	}
+
+	return matrix.front().front();
+}
+
+Rtn Model::_pivot(pair<size_t, size_t>& p)
+{
+	p = make_pair(0, 0);
+	double cmin = INT_MAX;
+	vector<double> coef = matrix.front();
+
+	// === <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD>Ԫ<EFBFBD><D4AA><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>Сֵ ===
+	for (size_t i = 1; i < coef.size(); i++) {
+		if (cmin > coef[i] && find(basic.begin(), basic.end(), i) == basic.end()) {
+			cmin = coef[i];
+			p.second = i;
+		}
+	}
+	if (cmin >= 0) {
+		return Rtn::OPTIMAL;
+	}
+	double bmin = INT_MAX;
+	for (size_t row = 1; row < bn; row++) {
+		double tmp = matrix[row].front() / matrix[row][p.second];
+		if (matrix[row][p.second] > 0 && bmin > tmp) {
+			bmin = tmp;
+			p.first = row;
+		}
+	}
+
+	if (abs(bmin - INT_MAX) < 1e-10) {
+		return Rtn::UNBOUNDED;
+	}
+
+	for (auto iter = basic.begin(); iter != basic.end(); iter++) {
+		if (matrix[p.first][*iter] != 0) {
+			*iter = p.second;
+			break;
+		}
+	}
+	assert(basic[p.first - 1] == p.second);
+	return Rtn::PIVOT;
+}
+
+void Model::_gaussian(pair<size_t, size_t> p)
+{
+	size_t x = p.first, y = p.second;
+
+	// === <20><><EFBFBD>й<EFBFBD>һ<EFBFBD><D2BB> ===
+	double norm = matrix[x][y];
+	for (size_t col = 0; col < matrix[x].size(); col++) {
+		matrix[x][col] /= norm;
+	}
+
+	// === <20><><EFBFBD><EFBFBD><EFBFBD>б任 ===
+	for (size_t row = 0; row < bn; row++) {
+		if (row == x) {
+			continue;
+		}
+		if (matrix[row][y] != 0) {
+			double norm = matrix[row][y];
+			for (size_t col = 0; col < matrix[x].size(); col++) {
+				matrix[row][col] = matrix[row][col] - norm * matrix[x][col];
+			}
+		}
+	}
+	basic[x - 1] = y;	// <20><>Ԫ
+}
+
+
+Expr::Expr(double constant)
+	:constant(constant)
+{
+}
+
+Expr::Expr(Var var, double coeff)
+{
+	this->coeffs.resize(var.col + 1);
+	this->coeffs[var.col] = coeff;
+	this->constant = 0;
+}
+
+Expr Expr::operator=(const Expr& rhs)
+{
+	return *this;
+}
+
+void Expr::operator+=(const Expr& expr)
+{
+	coeffs.resize(std::max(coeffs.size(), expr.coeffs.size()));
+	for (int c = 0; c < expr.coeffs.size(); c++) {
+		coeffs[c] += expr.coeffs[c];
+	}
+}
+
+void Expr::operator-=(const Expr& expr)
+{
+	coeffs.resize(std::max(coeffs.size(), expr.coeffs.size()));
+	for (int c = 0; c < expr.coeffs.size(); c++) {
+		coeffs[c] -= expr.coeffs[c];
+	}
+}
+
+void Expr::operator*=(double mult)
+{
+	for (int c = 0; c < coeffs.size(); c++) {
+		coeffs[c] *= mult;
+	}
+}
+
+void Expr::operator/=(double a)
+{
+	for (int c = 0; c < coeffs.size(); c++) {
+		coeffs[c] /= a;
+	}
+}
+
+Expr Expr::operator+(const Expr& rhs)
+{
+	coeffs.resize(std::max(coeffs.size(), rhs.coeffs.size()));
+	for (int c = 0; c < rhs.coeffs.size(); c++) {
+		coeffs[c] += rhs.coeffs[c];
+	}
+	return *this;
+}
+
+Expr Expr::operator-(const Expr& rhs)
+{
+	coeffs.resize(std::max(coeffs.size(), rhs.coeffs.size()));
+	for (int c = 0; c < rhs.coeffs.size(); c++) {
+		coeffs[c] -= rhs.coeffs[c];
+	}
+	return *this;
+}