From d4456e9632ca6f9342dc94769251284509fd6355 Mon Sep 17 00:00:00 2001 From: hit_lu Date: Fri, 23 May 2025 17:59:24 +0800 Subject: [PATCH] =?UTF-8?q?=E5=88=86=E6=94=AF=E5=AE=9A=E7=95=8C=E6=B3=95?= =?UTF-8?q?=E5=9F=BA=E6=9C=AC=E6=A1=86=E6=9E=B6?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- CMakeLists.txt | 8 +- inc/common.hpp | 19 ++- inc/node.hpp | 16 ++- inc/solver.hpp | 42 +++++- src/common.cpp | 10 ++ src/main.cpp | 16 ++- src/solver.cpp | 340 ++++++++++++++++++++++++++++++++++++++----------- 7 files changed, 353 insertions(+), 98 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 1ec41fe..152d1ba 100644 --- a/CMakeLists.txt +++ b/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( -# ${GUROBI_INCLUDE_DIRS} -#) +# include_directories( +# ${GUROBI_INCLUDE_DIRS} +# ) aux_source_directory(${SOURCE_DIR} src) aux_source_directory(${INCLUDE_DIR} inc) diff --git a/inc/common.hpp b/inc/common.hpp index 5dc0dbd..d871f09 100644 --- a/inc/common.hpp +++ b/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); diff --git a/inc/node.hpp b/inc/node.hpp index 7c06824..e58aa9d 100644 --- a/inc/node.hpp +++ b/inc/node.hpp @@ -1,5 +1,17 @@ #pragma once +#include -class Node { +namespace sv { + class Var; + class Model; + class Node + { + public: + Model* mdl; + double lower_bound; + double upper_bound; -}; \ No newline at end of file + bool is_integer; + std::vector branch_list; + }; +} diff --git a/inc/solver.hpp b/inc/solver.hpp index ca974c7..5fe27ed 100644 --- a/inc/solver.hpp +++ b/inc/solver.hpp @@ -3,31 +3,59 @@ namespace sv { - class Model - { + class LinSolver { public: - Model(); - Var* addVars(int col); + friend class Model; + 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& p); + rtn feasible_solution(); void _gaussian(std::pair p); Var* vars; matrix mt; + matrix mt_cvt; size_t cn, bn; std::vector 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; + }; } diff --git a/src/common.cpp b/src/common.cpp index 2869025..87bace8 100644 --- a/src/common.cpp +++ b/src/common.cpp @@ -1,6 +1,16 @@ #include "common.hpp" +#include using namespace sv; + +Var::Var(double coef, VarType type_) : + col(0), + val(0), + coeffs(coef), + type(type_) +{ +}; + double Var::get(DoubleAttr attr) { switch (attr) { diff --git a/src/main.cpp b/src/main.cpp index 37e8956..72cb57e 100644 --- a/src/main.cpp +++ b/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; diff --git a/src/solver.cpp b/src/solver.cpp index 4953777..b44ca18 100644 --- a/src/solver.cpp +++ b/src/solver.cpp @@ -4,7 +4,8 @@ #include #include #include - +#include +#include 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; // 释放原有内存 + vars = nullptr; + + cn = solver.cn; + if (cn > 0) { + vars = new Var[cn]; + for (int i = 0; i < cn; i++) { + vars[i] = solver.vars[i]; // 深拷贝 + } + } + + 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); // 添加越界检查 + return vars[idx]; +} + +void LinSolver::addConstr(const Expr& expr, ConstrOper sense, double rhs) { if (sense == ConstrOper::LESS_EQUAL) { + bn++; mt.push_back(vector(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(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++) { - for (size_t j = 0; j < mt[0].size(); j++) { - cout << mt[i][j] << "\t"; + assert(sense); + mt_cvt = mt; + 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(); - for (int row = 1; row < bn; row++) { - mt.front().push_back(0); - for (int col = 1; col < bn; col++) { - mt[row].push_back(col == row ? 1 : 0); + solver.optimize(); + if (solver.rtn_ != OPTIMAL) { + return solver.rtn_; + } + + double global_upper_bound = solver.obj_, global_lower_bound = 0; + + std::stack 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); + } + } } } - cn = mt.front().size(); - for (int i = 0; i < cn; i++) { - mt.front()[i] = sense * mt.front()[i]; + 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 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() // === 判断初始解是否为可行解 === 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() // === 构造初始可行解 === if (!initial_feasible) { - vector coeff = mt.front(); - mt.front() = vector(cn, .0); - mt.front().push_back(1); + vector coeff = mt_cvt.front(); + mt_cvt.front() = vector(cn, .0); + mt_cvt.front().push_back(1); pair t = { -1 ,cn }; for (int row = 1; row < bn; row++) { - mt[row].push_back(-1); - if (t.first == -1 || mt[row].front() < mt[t.first].front()) { + mt_cvt[row].push_back(-1); + 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++) { - if (fabs(mt.front()[col]) > 1e-10) { + for (int col = 1; col < mt_cvt.front().size(); col++) { + 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]; } } } - - 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_; + return LOADED; } -double Model::get(DoubleAttr attr) -{ - return res_; -} - -int Model::get(IntAttr attr) -{ - return cn - bn; -} - -double Model::_simplex() -{ - pair t; - while (1) { - rtn_ = _pivot(t); - if (rtn_ == OPTIMAL || rtn_ == UNBOUNDED) { - break; - } - _gaussian(t); - } - return mt.front().front(); -} - -rtn Model::_pivot(pair& p) +rtn LinSolver::_pivot(pair& p) { p = make_pair(0, 0); double cmin = INT_MAX; - vector coef = mt.front(); + vector coef = mt_cvt.front(); // === 非主轴元素中找最小值 === for (size_t i = 1; i < coef.size(); i++) { @@ -216,8 +404,8 @@ rtn Model::_pivot(pair& p) } double bmin = INT_MAX; for (size_t row = 1; row < bn; row++) { - double tmp = mt[row].front() / mt[row][p.second]; - if (mt[row][p.second] > 0 && bmin > tmp) { + double tmp = mt_cvt[row].front() / mt_cvt[row][p.second]; + if (mt_cvt[row][p.second] > 0 && bmin > tmp) { bmin = tmp; p.first = row; } @@ -228,7 +416,7 @@ rtn Model::_pivot(pair& 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& p) return PIVOT; } -void Model::_gaussian(pair p) +void LinSolver::_gaussian(pair p) { size_t x = p.first, y = p.second; // === 主行归一化 === - double norm = mt[x][y]; - for (size_t col = 0; col < mt[x].size(); col++) { - mt[x][col] /= norm; + double norm = mt_cvt[x][y]; + for (size_t col = 0; col < mt_cvt[x].size(); col++) { + mt_cvt[x][col] /= norm; } // === 其余行变换 === @@ -252,10 +440,10 @@ void Model::_gaussian(pair p) if (row == x) { continue; } - if (mt[row][y] != 0) { - double norm = mt[row][y]; - for (size_t col = 0; col < mt[x].size(); col++) { - mt[row][col] = mt[row][col] - norm * mt[x][col]; + if (mt_cvt[row][y] != 0) { + double norm = mt_cvt[row][y]; + for (size_t col = 0; col < mt_cvt[x].size(); col++) { + mt_cvt[row][col] = mt_cvt[row][col] - norm * mt_cvt[x][col]; } } }