分支定界法基本框架

CMakeLists.txt

@@ -2,7 +2,7 @@ cmake_minimum_required(VERSION 3.16.0)
 
 project(Solver)
 
-set(CMAKE_CXX_STANDARD 17)
+set(CMAKE_CXX_STANDARD 20)
 set(CMAKE_CXX_STANDARD_REQUIRED ON)
 
 # include(./FindGUROBI.cmake)
@@ -15,9 +15,9 @@ set(INCLUDE_DIR "${PROJECT_SOURCE_DIR}/inc")
 file(GLOB SOURCES "${SOURCE_DIR}/*.cpp")
 file(GLOB HEADERS "${INCLUDE_DIR}/*.hpp")
 
-#include_directories(
+# include_directories(
-# ${GUROBI_INCLUDE_DIRS}
+#  ${GUROBI_INCLUDE_DIRS}
-#)
+# )
 
 aux_source_directory(${SOURCE_DIR} src)
 aux_source_directory(${INCLUDE_DIR} inc)

inc/common.hpp

@@ -87,17 +87,28 @@ namespace sv {
 		EQUAL
 	};
 
+
+	enum class VarType {
+		CONTINUOUS,
+		BINARY,
+		INTEGER,
+	};
+
 	class Var {
 	public:
-		Var(double coef = 1) :col(0), val(0), coeffs(coef) {};
+		friend class Model;
+		friend class LinSolver;
+		friend class Expr;
+
+		Var(double coef = 1, VarType type_ = VarType::CONTINUOUS);
 		double get(DoubleAttr attr);
 		int get(IntAttr attr);
-		friend class Model;
+
-		friend class Expr;
 	private:
 		double coeffs;
 		double val;
 		int col;
+		VarType type;
 	};
 
 	Expr operator+(const Expr& x, const Expr& y);
@@ -131,7 +142,7 @@ namespace sv {
 		Expr(double constant = 0.0);
 		Expr(Var var, double coeff = 1.0);
 
-		friend class Model;
+		friend class LinSolver;
 
 		friend Expr operator+(const Expr& x, const Expr& y);
 		friend Expr operator+(const Expr& x);

inc/node.hpp

@@ -1,5 +1,17 @@
 #pragma once
+#include <vector>
 
-class Node {
+namespace sv {
+	class Var;
+	class Model;
+	class Node
+	{
+	public:
+		Model* mdl;
+		double lower_bound;
+		double upper_bound;
 
-};
+		bool is_integer;
+		std::vector<Var> branch_list;
+	};
+}

inc/solver.hpp

@@ -3,31 +3,59 @@
 
 namespace sv {
 
-	class Model
+	class LinSolver {
-	{
 	public:
-		Model();
+		friend class Model;
-		Var* addVars(int col);
+		LinSolver();
+		~LinSolver();
+
+		LinSolver(const LinSolver& solver);
+		LinSolver& operator=(const LinSolver& solver);
+
+		Var* addVars(int col, VarType type);
+		const Var& getVar(int idx);
 		void addConstr(const Expr& expr, ConstrOper sense, double rhs);
 		void setObjective(Expr obje, int sense = MDL_MAXIMIZE);
 		void print();
-		rtn optimize();
+
 		double get(DoubleAttr attr);
 		int get(IntAttr attr);
-	private:
+
+		rtn optimize();
+	protected:
+
 		double _simplex();
 		rtn _pivot(std::pair<size_t, size_t>& p);
+		rtn feasible_solution();
 		void _gaussian(std::pair<size_t, size_t> p);
 
 		Var* vars;
 
 		matrix mt;
+		matrix mt_cvt;
 		size_t cn, bn;
 		std::vector<int> basic;
 		rtn rtn_;
-		double res_;
+		double obj_;
 		int sense;
 	};
+
+	class Model
+	{
+	public:
+		Model();
+		rtn optimize();
+		
+		Var* addVars(int col, VarType type);
+		void addConstr(const Expr& expr, ConstrOper sense, double rhs);
+		void setObjective(Expr obje, int sense = MDL_MAXIMIZE);
+
+		double get(DoubleAttr attr);
+		int get(IntAttr attr);
+
+	private:
+		LinSolver solver;
+	};
 }
 
 

src/common.cpp

@@ -1,6 +1,16 @@
 #include "common.hpp"
+#include <iostream>
 using namespace sv;
 
+
+Var::Var(double coef, VarType type_) :
+	col(0), 
+	val(0), 
+	coeffs(coef), 
+	type(type_) 
+{
+};
+
 double Var::get(DoubleAttr attr)
 {
 	switch (attr) {

src/main.cpp

@@ -6,18 +6,24 @@
 using namespace std;
 using namespace sv;
 
+
 int main(int argc, char* argv[])
 {
 	Model mdl;
 
-	Var* vars = mdl.addVars(2);
+	Var* vars = mdl.addVars(2, VarType::INTEGER);
 
-	mdl.addConstr(vars[0] + vars[1], ConstrOper::LESS_EQUAL, 2);
-	mdl.addConstr(vars[0] + vars[1], ConstrOper::GREATER_EQUAL, 1);
-	mdl.addConstr(5 * vars[0] - 2 * vars[1], ConstrOper::GREATER_EQUAL, -2);
 
-	mdl.setObjective(vars[0] + 2 * vars[1], MDL_MAXIMIZE);
+	//mdl.addConstr(vars[0] + vars[1], ConstrOper::LESS_EQUAL, 2);
+	//mdl.addConstr(vars[0] + vars[1], ConstrOper::GREATER_EQUAL, 1);
+	//mdl.addConstr(5 * vars[0] - 2 * vars[1], ConstrOper::GREATER_EQUAL, -2);
 
+	//mdl.setObjective(vars[0] + 2 * vars[1], MDL_MAXIMIZE);
+
+	mdl.addConstr(2 * vars[0] + vars[1], ConstrOper::LESS_EQUAL, 10);
+	mdl.addConstr(3 * vars[0] + 6 * vars[1], ConstrOper::LESS_EQUAL, 40);
+
+	mdl.setObjective(100 * vars[0] + 150 * vars[1], MDL_MAXIMIZE);
 	switch(mdl.optimize()) {
 	case OPTIMAL:			
 		cout << "OPTIMAL SOLUTION: " << mdl.get(DoubleAttr::Obj) << endl; 

src/solver.cpp

@@ -4,7 +4,8 @@
 #include <iostream>
 #include <cassert>
 #include <algorithm>
-
+#include <stack>
+#include <cassert>
 using namespace sv;
 
 using std::make_pair;
@@ -13,35 +14,107 @@ using std::cout;
 using std::endl;
 using std::vector;
 
-Model::Model():
+
+struct Node
+{
+	LinSolver solver;
+	double lower_bound;
+	double upper_bound;
+
+};
+
+LinSolver::LinSolver() :
 	cn(0),
-	bn(0),
+	bn(1),
 	sense(0),
 	vars(nullptr)
 {
-
 }
 
-Var* Model::addVars(int num)
+sv::LinSolver::~LinSolver()
+{
+	delete[] vars;
+}
+
+sv::LinSolver::LinSolver(const LinSolver& solver)
+	: vars(nullptr),
+	cn(solver.cn),
+	bn(solver.bn),
+	mt(solver.mt),
+	basic(solver.basic),
+	rtn_(solver.rtn_),
+	obj_(solver.obj_),
+	sense(solver.sense) 
+{
+	delete[] vars;
+	vars = nullptr;
+	if (solver.vars != nullptr && cn > 0) {
+		vars = new Var[cn];
+		for (int i = 0; i < cn; i++) {
+			vars[i] = solver.vars[i];
+		}
+	}
+}
+
+LinSolver& sv::LinSolver::operator=(const LinSolver& solver)
+{
+	if (this == &solver) {
+		return *this;
+	}
+	delete[] vars; // <20>ͷ<EFBFBD>ԭ<EFBFBD><D4AD><EFBFBD>ڴ<EFBFBD>
+	vars = nullptr;
+
+	cn = solver.cn;
+	if (cn > 0) {
+		vars = new Var[cn];
+		for (int i = 0; i < cn; i++) {
+			vars[i] = solver.vars[i]; // <20><EFBFBD><EEBFBD>
+		}
+	}
+
+	mt = solver.mt;
+	cn = solver.cn, bn = solver.bn;
+	basic = solver.basic;
+	rtn_ = solver.rtn_;
+	obj_ = solver.obj_;
+	sense = solver.sense;
+
+	return *this;
+}
+
+Var* LinSolver::addVars(int num, VarType type)
 {
 	Var* old = vars;
 	vars = new Var[cn + num];
 	for (int c = cn; c < cn + num; c++) {
 		vars[c].col = c;
+		vars[c].type = type;
 	}
+
+	//for (int i = 0; i < cn; i++) {
+	//	*(vars + i) = *(old + i);
+	//}
 	memcpy(vars, old, sizeof(Var) * cn);
 	delete[] old;
 	cn += num;
 	return vars;
 }
 
-void Model::addConstr(const Expr& expr, ConstrOper sense, double rhs)
+const Var& sv::LinSolver::getVar(int idx)
+{
+	assert(idx >= 0 && idx < cn); // <20><><EFBFBD><EFBFBD>Խ<EFBFBD><D4BD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
+	return vars[idx];
+}
+
+void LinSolver::addConstr(const Expr& expr, ConstrOper sense, double rhs)
 {
 	if (sense == ConstrOper::LESS_EQUAL) {
+		bn++;
 		mt.push_back(vector<double>(1, rhs - expr.constant));
 		mt.back().insert(mt.back().end(), expr.coeffs.begin(), expr.coeffs.end());
 	}
 	else if (sense == ConstrOper::GREATER_EQUAL) {
+		bn++;
 		mt.push_back(vector<double>(1, expr.constant - rhs));
 		for (int coeff : expr.coeffs) {
 			mt.back().push_back(-coeff);
@@ -57,7 +130,7 @@ void Model::addConstr(const Expr& expr, ConstrOper sense, double rhs)
 	}
 }
 
-void Model::setObjective(Expr obje, int _sense)
+void LinSolver::setObjective(Expr obje, int _sense)
 {
 	assert(_sense == 1 || _sense == -1);
 	if (sense == 0) {
@@ -78,33 +151,178 @@ void Model::setObjective(Expr obje, int _sense)
 			}
 		}
 	}
+	for (int i = 0; i < mt.front().size(); i++) {
+		mt.front()[i] = _sense * mt.front()[i];
+	}
 	sense = _sense;
 }
 
-void Model::print()
+rtn LinSolver::optimize()
 {
-	for (size_t i = 0; i < mt.size(); i++) {
+	assert(sense);
-		for (size_t j = 0; j < mt[0].size(); j++) {
+	mt_cvt = mt;
-			cout << mt[i][j] << "\t";
+	rtn_ = feasible_solution();
+	if (rtn_ == LOADED) {
+		obj_ = _simplex();
+	}
+
+	cn = mt_cvt.front().size() - bn;
+	for (int row = 1; row < bn; row++) {
+		if (basic[row - 1] - 1 <= cn) {
+			vars[basic[row - 1] - 1].val = mt_cvt[row].front();
+		}
+	}
+	
+	return rtn_;
+}
+
+void LinSolver::print()
+{
+	for (size_t i = 0; i < mt_cvt.size(); i++) {
+		for (size_t j = 0; j < mt_cvt[0].size(); j++) {
+			cout << mt_cvt[i][j] << "\t";
 		}
 		cout << endl;
 	}
 }
 
+Model::Model()
+{
+
+}
+
 rtn Model::optimize()
 {
-	bn = mt.size();
+	solver.optimize();
-	for (int row = 1; row < bn; row++) {
+	if (solver.rtn_ != OPTIMAL) {
-		mt.front().push_back(0);
+		return solver.rtn_;
-		for (int col = 1; col < bn; col++) {
-			mt[row].push_back(col == row ? 1 : 0);
-		}
-	}
-	cn = mt.front().size();
-	for (int i = 0; i < cn; i++) {
-		mt.front()[i] = sense * mt.front()[i];
 	}
 
+	double global_upper_bound = solver.obj_, global_lower_bound = 0;
+
+	std::stack<Node> list_;
+	
+	Node incumbent_node, root_node;
+	root_node.solver = solver;
+	root_node.upper_bound = global_upper_bound, root_node.lower_bound = global_lower_bound;
+
+	list_.push(root_node);
+	int cnt = 0;
+	while (list_.size() && global_upper_bound - global_lower_bound > 1e-10) {
+		cout << ++cnt << endl;
+		Node current_node = list_.top();
+		list_.pop();
+		current_node.solver.optimize();
+		
+		if (current_node.solver.get(IntAttr::Status) == OPTIMAL) {
+			int branch_var_index = -1;
+
+			for (int i = 0; i < current_node.solver.get(IntAttr::NumVars); i++) {
+				if (current_node.solver.vars[i].type == VarType::INTEGER) {
+					if (fabs(int(current_node.solver.vars[i].val) - current_node.solver.vars[i].val) > 1e-10) {
+						branch_var_index = i;
+						break;
+					}
+				}
+			}
+
+			if (branch_var_index == -1) {
+				current_node.lower_bound = current_node.solver.obj_;
+				current_node.upper_bound = current_node.solver.obj_;
+				if (current_node.lower_bound > global_lower_bound) {
+					global_lower_bound = current_node.lower_bound;
+					incumbent_node = current_node;
+				}
+			}
+			else {
+				if (current_node.upper_bound >= global_lower_bound) {
+					const Var& branch_var = current_node.solver.getVar(branch_var_index);
+					int left_var_bound = branch_var.val;
+					int right_var_bound = branch_var.val + 1;
+
+					Node left_node = current_node;
+					left_node.solver.addConstr(branch_var, ConstrOper::LESS_EQUAL, left_var_bound);
+					list_.push(left_node);
+
+					Node right_node = current_node;
+					right_node.solver.addConstr(branch_var, ConstrOper::GREATER_EQUAL, right_var_bound);
+					list_.push(right_node);
+				}
+			}
+		}
+	}
+	solver.rtn_ = incumbent_node.solver.rtn_;
+	solver.obj_ = incumbent_node.solver.obj_;
+	for (int i = 0; i < solver.cn; i++) {
+		solver.vars[i].val = incumbent_node.solver.vars[i].val;
+	}
+	return solver.rtn_;
+}
+
+Var* sv::Model::addVars(int col, VarType type)
+{
+	return solver.addVars(col, type);
+}
+
+void sv::Model::addConstr(const Expr& expr, ConstrOper sense, double rhs)
+{
+	return solver.addConstr(expr, sense, rhs);
+}
+
+void sv::Model::setObjective(Expr obje, int sense)
+{
+	return solver.setObjective(obje, sense);
+}
+
+double sv::Model::get(DoubleAttr attr)
+{
+	return solver.get(attr);
+}
+
+int sv::Model::get(IntAttr attr)
+{
+	return solver.get(attr);
+}
+
+double LinSolver::get(DoubleAttr attr)
+{
+	return -sense * obj_;
+}
+
+int LinSolver::get(IntAttr attr)
+{
+	switch (attr) {
+	case IntAttr::NumVars:
+		return cn;
+	case IntAttr::Status:
+		return rtn_;
+	}
+	return -1;
+}
+
+double LinSolver::_simplex()
+{
+	pair<size_t, size_t> t;
+	while (1) {
+		rtn_ = _pivot(t);
+		if (rtn_ == OPTIMAL || rtn_ == UNBOUNDED) {
+			break;
+		}
+		_gaussian(t);
+	}	
+	return obj_ = mt_cvt.front().front();
+}
+
+rtn LinSolver::feasible_solution()
+{
+	for (int row = 1; row < bn; row++) {
+		mt_cvt.front().push_back(0);
+		for (int col = 1; col < bn; col++) {
+			mt_cvt[row].push_back(col == row ? 1 : 0);
+		}
+	}
+	cn = mt_cvt.front().size();
+	basic.clear();
 	for (size_t i = 1; i < bn; i++) {
 		basic.push_back(cn - bn + i);
 	}
@@ -112,7 +330,7 @@ rtn Model::optimize()
 	// === <20>жϳ<D0B6>ʼ<EFBFBD><CABC><EFBFBD>Ƿ<EFBFBD>Ϊ<EFBFBD><CEAA><EFBFBD>н<EFBFBD> ===
 	bool initial_feasible = true;
 	for (int row = 1; row < bn; row++) {
-		if (mt[row].front() < 0) {
+		if (mt_cvt[row].front() < 0) {
 			initial_feasible = false;
 			break;
 		}
@@ -120,14 +338,14 @@ rtn Model::optimize()
 
 	// === <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ʼ<EFBFBD><CABC><EFBFBD>н<EFBFBD> === 
 	if (!initial_feasible) {
-		vector<double> coeff = mt.front();
+		vector<double> coeff = mt_cvt.front();
-		mt.front() = vector<double>(cn, .0);
+		mt_cvt.front() = vector<double>(cn, .0);
-		mt.front().push_back(1);
+		mt_cvt.front().push_back(1);
 		pair<size_t, size_t> t = { -1 ,cn };
 
 		for (int row = 1; row < bn; row++) {
-			mt[row].push_back(-1);
+			mt_cvt[row].push_back(-1);
-			if (t.first == -1 || mt[row].front() < mt[t.first].front()) {
+			if (t.first == -1 || mt_cvt[row].front() < mt_cvt[t.first].front()) {
 				t.first = row;
 			}
 		}
@@ -141,8 +359,8 @@ rtn Model::optimize()
 		// 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 < mt.front().size(); col++) {
+			for (int col = 1; col < mt_cvt.front().size(); col++) {
-				if (fabs(mt.front()[col]) > 1e-10) {
+				if (fabs(mt_cvt.front()[col]) > 1e-10) {
 					t = make_pair(iter - basic.begin() + 1, col);
 					_gaussian(t);
 					break;
@@ -151,58 +369,28 @@ rtn Model::optimize()
 		}
 
 		for (int row = 0; row < bn; row++) {
-			mt[row].pop_back();
+			mt_cvt[row].pop_back();
 		}
 
 		// recover the coefficient line
 		for (int col = 0; col < cn; col++) {
-			mt.front()[col] = coeff[col];
+			mt_cvt.front()[col] = coeff[col];
 		}
 		for (int row = 1; row <= basic.size(); row++) {
-			int norm = mt.front()[basic[row - 1]];
+			int norm = mt_cvt.front()[basic[row - 1]];
 			for (int col = 0; col < cn; col++) {
-				mt.front()[col] -= norm * mt[row][col];
+				mt_cvt.front()[col] -= norm * mt_cvt[row][col];
 			}
 		}
 	}
-
+	return LOADED;
-	res_ = -sense * _simplex();
-	for (int row = 1; row < bn; row++) {
-		if (basic[row - 1] - 1 < cn - bn) {
-			vars[basic[row - 1] - 1].val = mt[row].front();
-		}
-	}
-	return rtn_;
 }
 
-double Model::get(DoubleAttr attr)
+rtn LinSolver::_pivot(pair<size_t, size_t>& p)
-{
-	return res_;
-}
-
-int Model::get(IntAttr attr)
-{
-	return cn - bn;
-}
-
-double Model::_simplex()
-{
-	pair<size_t, size_t> t;
-	while (1) {
-		rtn_ = _pivot(t);
-		if (rtn_ == OPTIMAL || rtn_ == UNBOUNDED) {
-			break;
-		}
-		_gaussian(t);
-	}
-	return mt.front().front();
-}
-
-rtn Model::_pivot(pair<size_t, size_t>& p)
 {
 	p = make_pair(0, 0);
 	double cmin = INT_MAX;
-	vector<double> coef = mt.front();
+	vector<double> coef = mt_cvt.front();
 
 	// === <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD>Ԫ<EFBFBD><D4AA><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>Сֵ ===
 	for (size_t i = 1; i < coef.size(); i++) {
@@ -216,8 +404,8 @@ rtn Model::_pivot(pair<size_t, size_t>& p)
 	}
 	double bmin = INT_MAX;
 	for (size_t row = 1; row < bn; row++) {
-		double tmp = mt[row].front() / mt[row][p.second];
+		double tmp = mt_cvt[row].front() / mt_cvt[row][p.second];
-		if (mt[row][p.second] > 0 && bmin > tmp) {
+		if (mt_cvt[row][p.second] > 0 && bmin > tmp) {
 			bmin = tmp;
 			p.first = row;
 		}
@@ -228,7 +416,7 @@ rtn Model::_pivot(pair<size_t, size_t>& p)
 	}
 
 	for (auto iter = basic.begin(); iter != basic.end(); iter++) {
-		if (mt[p.first][*iter] != 0) {
+		if (mt_cvt[p.first][*iter] != 0) {
 			*iter = p.second;
 			break;
 		}
@@ -237,14 +425,14 @@ rtn Model::_pivot(pair<size_t, size_t>& p)
 	return PIVOT;
 }
 
-void Model::_gaussian(pair<size_t, size_t> p)
+void LinSolver::_gaussian(pair<size_t, size_t> p)
 {
 	size_t x = p.first, y = p.second;
 
 	// === <20><><EFBFBD>й<EFBFBD>һ<EFBFBD><D2BB> ===
-	double norm = mt[x][y];
+	double norm = mt_cvt[x][y];
-	for (size_t col = 0; col < mt[x].size(); col++) {
+	for (size_t col = 0; col < mt_cvt[x].size(); col++) {
-		mt[x][col] /= norm;
+		mt_cvt[x][col] /= norm;
 	}
 
 	// === <20><><EFBFBD><EFBFBD><EFBFBD>б任 ===
@@ -252,10 +440,10 @@ void Model::_gaussian(pair<size_t, size_t> p)
 		if (row == x) {
 			continue;
 		}
-		if (mt[row][y] != 0) {
+		if (mt_cvt[row][y] != 0) {
-			double norm = mt[row][y];
+			double norm = mt_cvt[row][y];
-			for (size_t col = 0; col < mt[x].size(); col++) {
+			for (size_t col = 0; col < mt_cvt[x].size(); col++) {
-				mt[row][col] = mt[row][col] - norm * mt[x][col];
+				mt_cvt[row][col] = mt_cvt[row][col] - norm * mt_cvt[x][col];
 			}
 		}
 	}