diff --git a/CMakeLists.txt b/CMakeLists.txt index 152d1ba..d106132 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,6 +1,6 @@ cmake_minimum_required(VERSION 3.16.0) -project(Solver) +project(solver) set(CMAKE_CXX_STANDARD 20) set(CMAKE_CXX_STANDARD_REQUIRED ON) @@ -22,10 +22,10 @@ file(GLOB HEADERS "${INCLUDE_DIR}/*.hpp") aux_source_directory(${SOURCE_DIR} src) aux_source_directory(${INCLUDE_DIR} inc) -add_executable(Solver ${SOURCES} ${HEADERS}) +add_executable(solver ${SOURCES} ${HEADERS}) -#target_link_libraries(Solver +#target_link_libraries(solver # ${GUROBI_LIBRARIES} #) -target_include_directories(Solver PRIVATE ${SOURCE_DIR} ${INCLUDE_DIR}) +target_include_directories(solver PRIVATE ${SOURCE_DIR} ${INCLUDE_DIR}) diff --git a/TodoList.md b/TodoList.md new file mode 100644 index 0000000..6327dc3 --- /dev/null +++ b/TodoList.md @@ -0,0 +1,19 @@ +## P0 基础修复 + +1. 补齐变量约束能力：完善变量 LB/UB 上下界定义与校验；对 BINARY 类型严格强制 0-1 取值 +2. 增加测试与回归用例：搭建基础测试基线，覆盖整数、二进制、简单 MIP 算例，建立回归机制 + +## P1 核心能力 + +3. 完善终止条件与工程可控性：实现 TimeLimit 时间限制、MIPGap 收敛间隙终止条件；完善运行日志输出 +4. 迭代 MIP 分支与节点选择策略：按顺序落地：best-bound 选节点 → pseudocost 分支策略 + +## P2 进阶算法 + +5. 增加启发式、割平面与预处理能力：实现一种rounding舍入启发式、root层Gomory割平面、简单预处理逻辑 + +## P3 底层求解架构 + +6. 学习稀疏单纯形，对接HiGHS等成熟求解器架构：调研学习稀疏单纯形、稀疏基分解、对偶单纯形热启动机制 + +## P4 工程收口与标准化（待定） diff --git a/inc/common.hpp b/inc/common.hpp index d871f09..cd045bd 100644 --- a/inc/common.hpp +++ b/inc/common.hpp @@ -4,171 +4,163 @@ #define MDL_MINIMIZE 1 #define MDL_MAXIMIZE -1 -#define LOADED 0 -#define OPTIMAL 1 -#define INFEASIBLE 2 -#define INF_OR_UNBD 3 -#define UNBOUNDED 4 -#define CUTOFF 5 -#define ITERATION_LIMIT 6 -#define NODE_LIMIT 7 -#define TIME_LIMIT 8 -#define SOLUTION_LIMIT 9 -#define INTERRUPTED 10 -#define NUMERIC 11 -#define SUBOPTIMAL 12 -#define INPROGRESS 13 -#define USER_OBJ_LIMIT 14 -#define WORK_LIMIT 15 -#define MEM_LIMIT 16 -#define PIVOT 17 +#define LOADED 0 +#define OPTIMAL 1 +#define INFEASIBLE 2 +#define INF_OR_UNBD 3 +#define UNBOUNDED 4 +#define CUTOFF 5 +#define ITERATION_LIMIT 6 +#define NODE_LIMIT 7 +#define TIME_LIMIT 8 +#define SOLUTION_LIMIT 9 +#define INTERRUPTED 10 +#define NUMERIC 11 +#define SUBOPTIMAL 12 +#define INPROGRESS 13 +#define USER_OBJ_LIMIT 14 +#define WORK_LIMIT 15 +#define MEM_LIMIT 16 +#define PIVOT 17 namespace sv { - using matrix = std::vector>; - using rtn = int; - class Expr; - class Var; +using matrix = std::vector>; +using rtn = int; +class Expr; +class Var; - enum class IntAttr { - NumConstrs, - NumVars, - NumIntVars, - NumBinVars, - ModelSense, - IsMIP, - IsMultiObj, - Status, - SolCount, - Lazy, - NumObj, - NumCol - }; +enum class IntAttr { + NumConstrs, + NumVars, + NumIntVars, + NumBinVars, + ModelSense, + IsMIP, + IsMultiObj, + Status, + SolCount, + Lazy, + NumObj, + NumCol +}; - enum class DoubleAttr { - Runtime, - Work, - ObjCon, - LB, - UB, - Obj, - Start, - RHS, - Coeff, - MaxCoeff, - MinCoeff, - MaxBound, - MinBound, - ObjVal, - MIPGap, - IterCount, - NodeCount, - X, - Slack, - }; +enum class DoubleAttr { + Runtime, + Work, + ObjCon, + LB, + UB, + Obj, + Start, + RHS, + Coeff, + MaxCoeff, + MinCoeff, + MaxBound, + MinBound, + ObjVal, + MIPGap, + IterCount, + NodeCount, + X, + Slack, +}; - enum class StringAttr { - ModelName, - VarName, - ConstrName, - QCName, - GenConstrName, - ObjNName, - ScenNName, - BatchID, - VTag, - CTag, - QCTag, - BatchErrorMessage - }; +enum class StringAttr { + ModelName, + VarName, + ConstrName, + QCName, + GenConstrName, + ObjNName, + ScenNName, + BatchID, + VTag, + CTag, + QCTag, + BatchErrorMessage +}; - enum class ConstrOper { - LESS_EQUAL, - GREATER_EQUAL, - EQUAL - }; +enum class ConstrOper { LESS_EQUAL, GREATER_EQUAL, EQUAL }; +enum class VarType { + CONTINUOUS, + BINARY, + INTEGER, +}; - enum class VarType { - CONTINUOUS, - BINARY, - INTEGER, - }; +class Var { +public: + friend class Model; + friend class LinSolver; + friend class Expr; - class Var { - public: - 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); - Var(double coef = 1, VarType type_ = VarType::CONTINUOUS); - double get(DoubleAttr attr); - int get(IntAttr attr); +private: + double coeffs; + double val; + int col; + VarType type; +}; - private: - double coeffs; - double val; - int col; - VarType type; - }; +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); - 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); +class Expr { +private: + double constant; + std::vector coeffs; + std::vector vars; +public: + Expr(const Expr &expr) = default; + Expr(double constant = 0.0); + Expr(Var var, double coeff = 1.0); - class Expr - { - private: - double constant; - std::vector coeffs; - std::vector vars; + friend class LinSolver; - public: - Expr(const Expr& expr) = default; - Expr(double constant = 0.0); - Expr(Var var, double coeff = 1.0); + 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); - friend class LinSolver; - - 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); - }; -} + 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); +}; +} // namespace sv diff --git a/inc/solver.hpp b/inc/solver.hpp index b6315cb..b7eef4f 100644 --- a/inc/solver.hpp +++ b/inc/solver.hpp @@ -3,60 +3,56 @@ namespace sv { - class LinSolver { - public: - friend class Model; - LinSolver(); - ~LinSolver(); +class LinSolver { +public: + friend class Model; + LinSolver(); + ~LinSolver(); - LinSolver(const LinSolver& solver); - LinSolver& operator=(const LinSolver& solver); + LinSolver(const LinSolver &solver); + LinSolver &operator=(const LinSolver &solver); - Var* addVars(int col, VarType type); - Var& getVar(int idx); - void addConstr(const Expr& expr, ConstrOper sense, double rhs); - void setObjective(Expr obje, int sense = MDL_MAXIMIZE); - void print(); + Var *addVars(int col, VarType type); + Var &getVar(int idx); + void addConstr(const Expr &expr, ConstrOper sense, double rhs); + void setObjective(Expr obje, int sense = MDL_MAXIMIZE); + void print(); - double get(DoubleAttr attr); - int get(IntAttr attr); + double get(DoubleAttr attr); + int get(IntAttr attr); - rtn optimize(); - protected: + rtn optimize(); - double _simplex(); - rtn _pivot(std::pair& p); - rtn feasible_solution(); - void _gaussian(std::pair p); +protected: + double _simplex(); + rtn _pivot(std::pair &p); + rtn feasible_solution(); + void _gaussian(std::pair p); - std::vector> vars; + std::vector> vars; - matrix table; - matrix ope_table; - size_t cn, bn; - std::vector basic; - rtn rtn_; - double obj_; - int sense; - }; + matrix table; + matrix ope_table; + size_t cn, bn; + std::vector basic; + rtn rtn_; + 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; - }; -} +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; +}; +} // namespace sv diff --git a/solver b/solver new file mode 100755 index 0000000..0c72b07 Binary files /dev/null and b/solver differ diff --git a/src/common.cpp b/src/common.cpp index b8b02fe..63ed9c8 100644 --- a/src/common.cpp +++ b/src/common.cpp @@ -2,17 +2,9 @@ #include using namespace sv; +Var::Var(double coef, VarType type_) : col(0), val(0), coeffs(coef), type(type_){}; -Var::Var(double coef, VarType type_) : - col(0), - val(0), - coeffs(coef), - type(type_) -{ -}; - -double Var::get(DoubleAttr attr) -{ +double Var::get(DoubleAttr attr) { switch (attr) { case DoubleAttr::Coeff: return coeffs; @@ -22,13 +14,11 @@ double Var::get(DoubleAttr attr) return -1; } -int Var::get(IntAttr attr) -{ +int Var::get(IntAttr attr) { return col; } -Expr sv::operator+(const Expr& x, const Expr& y) -{ +Expr sv::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++) { @@ -43,13 +33,11 @@ Expr sv::operator+(const Expr& x, const Expr& y) return exp; } -Expr sv::operator+(const Expr& x) -{ +Expr sv::operator+(const Expr &x) { return x; } -Expr sv::operator+(Var x, Var y) -{ +Expr sv::operator+(Var x, Var y) { Expr exp; exp.coeffs.resize(std::max(x.get(IntAttr::NumCol) + 1, y.get(IntAttr::NumCol) + 1), 0); exp.coeffs.at(x.get(IntAttr::NumCol)) = x.get(DoubleAttr::Coeff); @@ -57,21 +45,19 @@ Expr sv::operator+(Var x, Var y) return exp; } -Expr sv::operator+(Var x, double a) -{ +Expr sv::operator+(Var x, double a) { Expr exp; - exp.coeffs.resize(x.get(IntAttr::NumCol) + 1); + exp.coeffs.resize(x.get(IntAttr::NumCol) + 1, 0); + exp.coeffs.at(x.get(IntAttr::NumCol)) = x.get(DoubleAttr::Coeff); exp.constant = a; return exp; } -Expr sv::operator+(double a, Var x) -{ +Expr sv::operator+(double a, Var x) { return x + a; } -Expr sv::operator-(const Expr& x, const Expr& y) -{ +Expr sv::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++) { @@ -79,15 +65,14 @@ Expr sv::operator-(const Expr& x, const Expr& y) exp.coeffs.at(c) = x.coeffs.at(c) - y.coeffs.at(c); } else { - exp.coeffs.at(c) = c < x.coeffs.size() ? x.coeffs.at(c) : y.coeffs.at(c); + exp.coeffs.at(c) = c < x.coeffs.size() ? x.coeffs.at(c) : -y.coeffs.at(c); } } - exp.constant = x.constant + y.constant; + exp.constant = x.constant - y.constant; return exp; } -Expr sv::operator-(const Expr& x) -{ +Expr sv::operator-(const Expr &x) { Expr expr(x); for (int c = 0; c < expr.coeffs.size(); c++) { expr.coeffs.at(c) = -expr.coeffs.at(c); @@ -96,13 +81,11 @@ Expr sv::operator-(const Expr& x) return expr; } -Expr sv::operator-(Var x) -{ +Expr sv::operator-(Var x) { return -Expr(x); } -Expr sv::operator-(Var x, Var y) -{ +Expr sv::operator-(Var x, Var y) { Expr exp; exp.coeffs.resize(std::max(x.get(IntAttr::NumCol) + 1, y.get(IntAttr::NumCol) + 1), 0); exp.coeffs.at(x.get(IntAttr::NumCol)) = x.get(DoubleAttr::Coeff); @@ -110,31 +93,30 @@ Expr sv::operator-(Var x, Var y) return exp; } -Expr sv::operator-(Var x, double a) -{ - return x - Var(a); +Expr sv::operator-(Var x, double a) { + Expr exp; + exp.coeffs.resize(x.get(IntAttr::NumCol) + 1, 0); + exp.coeffs.at(x.get(IntAttr::NumCol)) = x.get(DoubleAttr::Coeff); + exp.constant = -a; + return exp; } -Expr sv::operator-(double a, Var x) -{ - return x - a; +Expr sv::operator-(double a, Var x) { + return a + (-x); } -Expr sv::operator*(double a, Var x) -{ +Expr sv::operator*(double a, Var x) { Expr exp; exp.coeffs.resize(x.get(IntAttr::NumCol) + 1, 0); exp.coeffs.at(x.get(IntAttr::NumCol)) = a * x.get(DoubleAttr::Coeff); return exp; } -Expr sv::operator*(Var x, double a) -{ +Expr sv::operator*(Var x, double a) { return a * x; } -Expr sv::operator*(const Expr& x, double a) -{ +Expr sv::operator*(const Expr &x, double a) { Expr exp = x; for (int c = 0; c < exp.coeffs.size(); c++) { exp.coeffs.at(c) *= a; @@ -143,18 +125,15 @@ Expr sv::operator*(const Expr& x, double a) return exp; } -Expr sv::operator*(double a, const Expr& x) -{ +Expr sv::operator*(double a, const Expr &x) { return x * a; } -Expr sv::operator/(Var x, double a) -{ +Expr sv::operator/(Var x, double a) { return Expr(x) / a; } -Expr sv::operator/(const Expr& x, double a) -{ +Expr sv::operator/(const Expr &x, double a) { Expr exp = x; for (int c = 0; c < exp.coeffs.size(); c++) { exp.coeffs.at(c) /= a; @@ -163,25 +142,24 @@ Expr sv::operator/(const Expr& x, double a) return exp; } -Expr::Expr(double constant) - :constant(constant) -{ -} +Expr::Expr(double constant) : constant(constant) {} -Expr::Expr(Var var, double coeff) -{ +Expr::Expr(Var var, double coeff) { this->coeffs.resize(var.col + 1); this->coeffs.at(var.col) = coeff; this->constant = 0; } -Expr Expr::operator=(const Expr& rhs) -{ +Expr &Expr::operator=(const Expr &rhs) { + if (this != &rhs) { + constant = rhs.constant; + coeffs = rhs.coeffs; + vars = rhs.vars; + } return *this; } -void Expr::operator+=(const Expr& expr) -{ +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.at(c) += expr.coeffs.at(c); @@ -189,8 +167,7 @@ void Expr::operator+=(const Expr& expr) constant += expr.constant; } -void Expr::operator-=(const Expr& expr) -{ +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.at(c) -= expr.coeffs.at(c); @@ -198,24 +175,21 @@ void Expr::operator-=(const Expr& expr) constant -= expr.constant; } -void Expr::operator*=(double mult) -{ +void Expr::operator*=(double mult) { for (int c = 0; c < coeffs.size(); c++) { coeffs.at(c) *= mult; } constant *= mult; } -void Expr::operator/=(double a) -{ +void Expr::operator/=(double a) { for (int c = 0; c < coeffs.size(); c++) { coeffs.at(c) /= a; } constant /= a; } -Expr Expr::operator+(const Expr& rhs) -{ +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.at(c) += rhs.coeffs.at(c); @@ -224,8 +198,7 @@ Expr Expr::operator+(const Expr& rhs) return *this; } -Expr Expr::operator-(const Expr& rhs) -{ +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.at(c) -= rhs.coeffs.at(c); diff --git a/src/main.cpp b/src/main.cpp index b1e9b2f..ce86d4e 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -6,36 +6,35 @@ using namespace std; using namespace sv; - -int main(int argc, char* argv[]) -{ +int main(int argc, char *argv[]) { Model mdl; - - Var* int_var = mdl.addVars(3, VarType::INTEGER); - Var* con_var = mdl.addVars(2, VarType::CONTINUOUS); - mdl.addConstr(2 * int_var[0] + int_var[1], ConstrOper::LESS_EQUAL, 10); - mdl.addConstr(3 * int_var[0] + 6 * int_var[1], ConstrOper::LESS_EQUAL, 40); - mdl.addConstr(3 * int_var[0] + 6 * int_var[1] + 4 * int_var[2], ConstrOper::LESS_EQUAL, 50); - mdl.addConstr(2.3 * con_var[0] + 2.6 * con_var[1], ConstrOper::LESS_EQUAL, 80); - mdl.addConstr(con_var[0] + 2 * con_var[1], ConstrOper::LESS_EQUAL, 70); - mdl.setObjective(100 * int_var[0] + 150 * int_var[1] + 120 * int_var[2] + 82.6 * con_var[0] + 90.4 * con_var[1], MDL_MAXIMIZE); + Var *x = mdl.addVars(3, VarType::INTEGER); + Var *y = mdl.addVars(2, VarType::CONTINUOUS); - switch(mdl.optimize()) { - case OPTIMAL: - cout << "OPTIMAL SOLUTION: " << mdl.get(DoubleAttr::Obj) << endl; + mdl.addConstr(3 * x[0] + 6 * x[1], ConstrOper::LESS_EQUAL, 28); + mdl.addConstr(3 * x[0] + 6 * x[1] + 4 * x[2], ConstrOper::LESS_EQUAL, 30); + mdl.addConstr(3 * y[0] + 2 * y[1], ConstrOper::LESS_EQUAL, 80); + mdl.addConstr(y[0] + 2 * y[1], ConstrOper::LESS_EQUAL, 70); + mdl.addConstr(2 * x[1] + x[2], ConstrOper::EQUAL, 7); + mdl.setObjective(100 * x[0] + 150 * x[1] + 120 * x[2] + 80 * y[0] + 90 * y[1] + 15, + MDL_MAXIMIZE); + + switch (mdl.optimize()) { + case OPTIMAL: + cout << "OPTIMAL SOLUTION: " << mdl.get(DoubleAttr::Obj) << endl; for (int i = 0; i < 3; i++) { - cout << "integer var [" << i << "] : " << int_var[i].get(DoubleAttr::X) << endl; + cout << "integer var [" << i << "] : " << x[i].get(DoubleAttr::X) << endl; } for (int i = 0; i < 2; i++) { - cout << "continuous var [" << i << "] : " << con_var[i].get(DoubleAttr::X) << endl; + cout << "continuous var [" << i << "] : " << y[i].get(DoubleAttr::X) << endl; } break; - case INFEASIBLE: - cout << "INFEASIBLE MODEL" << endl; + case INFEASIBLE: + cout << "INFEASIBLE MODEL" << endl; break; - case UNBOUNDED: - cout << "UNBOUNDED SOLUTION" << endl; + case UNBOUNDED: + cout << "UNBOUNDED SOLUTION" << endl; break; default: assert(false); diff --git a/src/solver.cpp b/src/solver.cpp index e169d77..79eb803 100644 --- a/src/solver.cpp +++ b/src/solver.cpp @@ -5,68 +5,53 @@ #include #include #include -#include -#include +#include #include -#include +#include + using namespace sv; -using std::make_pair; -using std::pair; using std::cout; using std::endl; +using std::make_pair; +using std::pair; using std::vector; - -struct Node -{ +struct Node { LinSolver solver; - double lower_bound; - double upper_bound; - + double bound; // LP bound in true-objective space }; -LinSolver::LinSolver() : - obj_(0), - rtn_(LOADED), - cn(0), - bn(1), - sense(0) -{ +LinSolver::LinSolver() : obj_(0), rtn_(LOADED), cn(0), bn(1), sense(0) {} -} - -sv::LinSolver::~LinSolver() -{ - for (auto& var : vars) { +sv::LinSolver::~LinSolver() { + for (auto &var : vars) { delete[] var.second; } } -sv::LinSolver::LinSolver(const LinSolver& solver) -{ +sv::LinSolver::LinSolver(const LinSolver &solver) { *this = solver; } -LinSolver& sv::LinSolver::operator=(const LinSolver& solver) -{ +LinSolver &sv::LinSolver::operator=(const LinSolver &solver) { if (this == &solver) { return *this; } - for (auto& var : vars) { + for (auto &var : vars) { delete[] var.second; } vars.clear(); vars.reserve(solver.vars.size()); - for (auto& var : solver.vars) { + for (auto &var : solver.vars) { vars.push_back(std::make_pair(var.first, new Var[var.first])); for (int i = 0; i < var.first; i++) { vars.back().second[i] = var.second[i]; } } cn = solver.cn; + bn = solver.bn; table = solver.table; - cn = solver.cn, bn = solver.bn; basic = solver.basic; rtn_ = solver.rtn_; obj_ = solver.obj_; @@ -75,85 +60,73 @@ LinSolver& sv::LinSolver::operator=(const LinSolver& solver) return *this; } -Var* LinSolver::addVars(int num, VarType type) -{ - Var* var = new Var[num]; +Var *LinSolver::addVars(int num, VarType type) { + Var *var = new Var[num]; for (int c = 0; c < num; c++) { - var[c].col = c + cn, var[c].type = type; + var[c].col = c + cn; + var[c].type = type; } vars.push_back(std::make_pair(num, var)); cn += num; return vars.back().second; } -Var& sv::LinSolver::getVar(int idx) -{ - assert(idx >= 0 && idx < cn); +Var &sv::LinSolver::getVar(int idx) { + assert(idx >= 0 && idx < static_cast(cn)); static Var err_var; - for (auto& var : vars) { - if (var.first <= idx) { - idx -= var.first; - } - else { - return var.second[idx]; + int offset = 0; + for (auto &var : vars) { + if (idx < offset + var.first) { + return var.second[idx - offset]; } + offset += var.first; } return err_var; } -void LinSolver::addConstr(const Expr& expr, ConstrOper sense, double rhs) -{ +void LinSolver::addConstr(const Expr &expr, ConstrOper sense, double rhs) { + if (sense == ConstrOper::EQUAL) { + addConstr(expr, ConstrOper::LESS_EQUAL, rhs); + addConstr(expr, ConstrOper::GREATER_EQUAL, rhs); + return; + } + + bn++; if (sense == ConstrOper::LESS_EQUAL) { - bn++; table.push_back(vector(1, rhs - expr.constant)); table.back().insert(table.back().end(), expr.coeffs.begin(), expr.coeffs.end()); } - else if (sense == ConstrOper::GREATER_EQUAL) { - bn++; + else { table.push_back(vector(1, expr.constant - rhs)); - for (int coeff : expr.coeffs) { + for (double coeff : expr.coeffs) { table.back().push_back(-coeff); } } - else { - addConstr(expr, ConstrOper::LESS_EQUAL, rhs); - addConstr(expr, ConstrOper::GREATER_EQUAL, rhs); - } - - for (int c = table.back().size(); c <= cn; c++) { + + for (size_t c = table.back().size(); c <= cn; c++) { table.back().push_back(0); } } -void LinSolver::setObjective(Expr obje, int _sense) -{ +void LinSolver::setObjective(Expr obje, int _sense) { assert(_sense == 1 || _sense == -1); if (sense == 0) { table.insert(table.begin(), obje.coeffs); - table.front().insert(table.front().begin(), -obje.constant); - for (int c = obje.coeffs.size() + 1; c <= cn; c++) { - table.front().push_back(0); - } } else { - table.front().front() = -obje.constant; - for (int col = 0; col < cn; col++) { - if (col < obje.coeffs.size()) { - table.front().at(col + 1) = obje.coeffs.at(col); - } - else { - table.front().at(col) = 0; - } - } + table.front() = obje.coeffs; } - for (int row = 0; row < table.front().size(); row++) { - table.front().at(row) = _sense * table.front().at(row); + table.front().insert(table.front().begin(), -obje.constant); + for (size_t c = obje.coeffs.size() + 1; c <= cn; c++) { + table.front().push_back(0); + } + for (size_t col = 0; col < table.front().size(); col++) { + table.front().at(col) = _sense * table.front().at(col); } sense = _sense; } -rtn LinSolver::optimize() -{ +rtn LinSolver::optimize() { assert(sense); ope_table = table; rtn_ = LOADED; @@ -163,18 +136,22 @@ rtn LinSolver::optimize() } if (rtn_ == OPTIMAL) { - cn = ope_table.front().size() - bn; - for (int row = 1; row < bn; row++) { - if (basic.at(row - 1) - 1 < cn) { - getVar(basic.at(row - 1) - 1).val = ope_table.at(row).front(); + const size_t num_vars = ope_table.front().size() - bn; + cn = num_vars; + for (size_t i = 0; i < num_vars; i++) { + getVar(static_cast(i)).val = 0; + } + for (size_t row = 1; row < bn; row++) { + int var_idx = basic.at(row - 1) - 1; + if (var_idx >= 0 && static_cast(var_idx) < num_vars) { + getVar(var_idx).val = ope_table.at(row).front(); } } } return rtn_; } -void LinSolver::print() -{ +void LinSolver::print() { for (size_t row = 0; row < ope_table.size(); row++) { for (size_t col = 0; col < ope_table.front().size(); col++) { cout << ope_table.at(row).at(col) << "\t"; @@ -183,163 +160,206 @@ void LinSolver::print() } } -Model::Model() -{ +Model::Model() {} +namespace { + +bool is_integer_type(VarType type) { + return type == VarType::INTEGER || type == VarType::BINARY; } -rtn Model::optimize() -{ +bool is_fractional(double val, double eps = 1e-6) { + return fabs(val - std::round(val)) > eps; +} + +bool is_better(double candidate, double incumbent, int sense, double eps = 1e-10) { + if (sense == MDL_MAXIMIZE) { + return candidate > incumbent + eps; + } + return candidate < incumbent - eps; +} + +} // namespace + +rtn Model::optimize() { solver.optimize(); if (solver.rtn_ != OPTIMAL) { return solver.rtn_; } - double global_upper_bound = solver.obj_, global_lower_bound = 0; + const int sense = solver.sense; + const int num_vars = solver.get(IntAttr::NumVars); - std::stack list_; - - Node root_node = { solver, 0, solver.obj_ }; - Node incumbent_node = root_node; + bool has_integer = false; + for (int i = 0; i < num_vars; i++) { + if (is_integer_type(solver.getVar(i).type)) { + has_integer = true; + break; + } + } + if (!has_integer) { + return solver.rtn_; + } + + double best_obj = (sense == MDL_MAXIMIZE) ? -std::numeric_limits::infinity() + : std::numeric_limits::infinity(); + bool found_integer = false; + Node incumbent_node{solver, solver.get(DoubleAttr::Obj)}; + + std::stack open_nodes; + open_nodes.push(Node{solver, solver.get(DoubleAttr::Obj)}); + while (!open_nodes.empty()) { + Node current_node = std::move(open_nodes.top()); + open_nodes.pop(); + + if (found_integer && !is_better(current_node.bound, best_obj, sense)) { + continue; + } - list_.push(root_node); - while (list_.size() && global_upper_bound - global_lower_bound > 1e-10) { - Node current_node = list_.top(); - list_.pop(); current_node.solver.optimize(); - - if (current_node.solver.get(IntAttr::Status) == OPTIMAL) { - int branch_var_index = -1; + if (current_node.solver.get(IntAttr::Status) != OPTIMAL) { + continue; + } - for (int i = 0; i < current_node.solver.get(IntAttr::NumVars); i++) { - if (current_node.solver.getVar(i).type == VarType::INTEGER) { - if (fabs(int(current_node.solver.getVar(i).val) - current_node.solver.getVar(i).val) > 1e-10) { - branch_var_index = i; - break; - } - } + const double lp_obj = current_node.solver.get(DoubleAttr::Obj); + current_node.bound = lp_obj; + + if (found_integer && !is_better(lp_obj, best_obj, sense)) { + continue; + } + + int branch_var_index = -1; + for (int i = 0; i < current_node.solver.get(IntAttr::NumVars); i++) { + Var &var = current_node.solver.getVar(i); + if (!is_integer_type(var.type)) { + continue; } - - 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; + if (is_fractional(var.val)) { + if (branch_var_index < 0) { + branch_var_index = i; } } 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); - } + // Snap near-integer values so later checks / output stay clean + var.val = std::round(var.val); } } + + if (branch_var_index == -1) { + if (!found_integer || is_better(lp_obj, best_obj, sense)) { + best_obj = lp_obj; + found_integer = true; + incumbent_node = current_node; + } + continue; + } + + const Var &branch_var = current_node.solver.getVar(branch_var_index); + const int left_bound = static_cast(std::floor(branch_var.val)); + const int right_bound = left_bound + 1; + + Node left_node = current_node; + left_node.solver.addConstr(branch_var, ConstrOper::LESS_EQUAL, left_bound); + open_nodes.push(std::move(left_node)); + + Node right_node = current_node; + right_node.solver.addConstr(branch_var, ConstrOper::GREATER_EQUAL, right_bound); + open_nodes.push(std::move(right_node)); } + + if (!found_integer) { + return solver.rtn_ = INFEASIBLE; + } + solver.rtn_ = incumbent_node.solver.rtn_; solver.obj_ = incumbent_node.solver.obj_; - for (int i = 0; i < solver.cn; i++) { + for (int i = 0; i < solver.get(IntAttr::NumVars); i++) { solver.getVar(i).val = incumbent_node.solver.getVar(i).val; } return solver.rtn_; } -Var* sv::Model::addVars(int col, VarType type) -{ +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::addConstr(const Expr &expr, ConstrOper sense, double rhs) { + solver.addConstr(expr, sense, rhs); } -void sv::Model::setObjective(Expr obje, int sense) -{ - return solver.setObjective(obje, sense); +void sv::Model::setObjective(Expr obje, int sense) { + solver.setObjective(obje, sense); } -double sv::Model::get(DoubleAttr attr) -{ +double sv::Model::get(DoubleAttr attr) { return solver.get(attr); } -int sv::Model::get(IntAttr attr) -{ +int sv::Model::get(IntAttr attr) { return solver.get(attr); } -double LinSolver::get(DoubleAttr attr) -{ +double LinSolver::get(DoubleAttr attr) { return -sense * obj_; } -int LinSolver::get(IntAttr attr) -{ +int LinSolver::get(IntAttr attr) { switch (attr) { case IntAttr::NumVars: - return cn; + return static_cast(cn); case IntAttr::Status: return rtn_; + default: + return -1; } - return -1; } -double LinSolver::_simplex() -{ +double LinSolver::_simplex() { pair t; - while (1) { + while (true) { rtn_ = _pivot(t); if (rtn_ == OPTIMAL || rtn_ == UNBOUNDED) { break; } _gaussian(t); - } + } return obj_ = ope_table.front().front(); } -rtn LinSolver::feasible_solution() -{ - for (int row = 1; row < bn; row++) { +rtn LinSolver::feasible_solution() { + for (size_t row = 1; row < bn; row++) { ope_table.front().push_back(0); - for (int col = 1; col < bn; col++) { - ope_table.at(row).push_back(col == row ? 1 : 0); + for (size_t col = 1; col < bn; col++) { + ope_table.at(row).push_back(col == row ? 1.0 : 0.0); } } cn = ope_table.front().size(); basic.clear(); + basic.reserve(bn - 1); for (size_t i = 1; i < bn; i++) { - basic.push_back(cn - bn + i); + basic.push_back(static_cast(cn - bn + i)); } - // === 判断初始解是否为可行解 === + // Check whether the initial basic solution is feasible bool initial_feasible = true; - for (int row = 1; row < bn; row++) { + for (size_t row = 1; row < bn; row++) { if (ope_table.at(row).front() < 0) { initial_feasible = false; break; } } - // === 构造初始可行解 === + // Two-phase method when the initial basis is infeasible if (!initial_feasible) { vector coeff = ope_table.front(); - ope_table.front() = vector(cn, .0); + ope_table.front() = vector(cn, 0.0); ope_table.front().push_back(1); - pair t = { -1 ,cn }; + pair t = {static_cast(-1), cn}; - for (int row = 1; row < bn; row++) { + for (size_t row = 1; row < bn; row++) { ope_table.at(row).push_back(-1); - if (t.first == -1 || ope_table.at(row).front() < ope_table.at(t.first).front()) { + if (t.first == static_cast(-1) || + ope_table.at(row).front() < ope_table.at(t.first).front()) { t.first = row; } } @@ -349,29 +369,30 @@ rtn LinSolver::feasible_solution() return rtn_ = INFEASIBLE; } rtn_ = LOADED; - // if the x0 in B, we should pivot it. - auto iter = find(basic.begin(), basic.end(), cn); + + // If artificial variable remains basic, pivot it out + auto iter = find(basic.begin(), basic.end(), static_cast(cn)); if (iter != basic.end()) { - for (int col = 1; col < ope_table.front().size(); col++) { + for (size_t col = 1; col < ope_table.front().size(); col++) { if (fabs(ope_table.front().at(col)) > 1e-10) { - t = make_pair(iter - basic.begin() + 1, col); + t = make_pair(static_cast(iter - basic.begin() + 1), col); _gaussian(t); break; } } } - - for (int row = 0; row < bn; row++) { + + for (size_t row = 0; row < bn; row++) { ope_table.at(row).pop_back(); } - // recover the coefficient line - for (int col = 0; col < cn; col++) { + // Restore original objective row and re-express in current basis + for (size_t col = 0; col < cn; col++) { ope_table.front().at(col) = coeff.at(col); } - for (int row = 1; row <= basic.size(); row++) { - int norm = ope_table.front().at(basic.at(row - 1)); - for (int col = 0; col < cn; col++) { + for (size_t row = 1; row <= basic.size(); row++) { + double norm = ope_table.front().at(basic.at(row - 1)); + for (size_t col = 0; col < cn; col++) { ope_table.front().at(col) -= norm * ope_table.at(row).at(col); } } @@ -379,15 +400,15 @@ rtn LinSolver::feasible_solution() return rtn_; } -rtn LinSolver::_pivot(pair& p) -{ +rtn LinSolver::_pivot(pair &p) { p = make_pair(0, 0); double cmin = DBL_MAX; - vector coef = ope_table.front(); + const vector &coef = ope_table.front(); - // === 非主轴元素中找最小值 === + // Entering variable: most negative reduced cost for (size_t col = 1; col < coef.size(); col++) { - if (cmin > coef.at(col) && find(basic.begin(), basic.end(), col) == basic.end()) { + if (cmin > coef.at(col) && + find(basic.begin(), basic.end(), static_cast(col)) == basic.end()) { cmin = coef.at(col); p.second = col; } @@ -395,51 +416,47 @@ rtn LinSolver::_pivot(pair& p) if (cmin >= 0) { return OPTIMAL; } + double bmin = DBL_MAX; for (size_t row = 1; row < bn; row++) { - double tmp = ope_table.at(row).front() / ope_table.at(row).at(p.second); - if (ope_table.at(row).at(p.second) > 0 && bmin > tmp) { - bmin = tmp; - p.first = row; + const double pivot_col = ope_table.at(row).at(p.second); + if (pivot_col > 0) { + const double tmp = ope_table.at(row).front() / pivot_col; + if (bmin > tmp) { + bmin = tmp; + p.first = row; + } } } - if (abs(bmin - DBL_MAX) < 1e-10) { + if (bmin >= DBL_MAX / 2) { return UNBOUNDED; } - for (auto iter = basic.begin(); iter != basic.end(); iter++) { - if (ope_table.at(p.first).at(*iter) != 0) { - *iter = p.second; - break; - } - } - assert(basic.at(p.first - 1) == p.second); + basic.at(p.first - 1) = static_cast(p.second); return PIVOT; } -void LinSolver::_gaussian(pair p) -{ +void LinSolver::_gaussian(pair p) { size_t x = p.first, y = p.second; - // === 主行归一化 === + // Normalize pivot row double norm = ope_table.at(x).at(y); for (size_t col = 0; col < ope_table.at(x).size(); col++) { ope_table.at(x).at(col) /= norm; } - // === 其余行变换 === + // Eliminate pivot column in other rows for (size_t row = 0; row < bn; row++) { if (row == x) { continue; } if (ope_table.at(row).at(y) != 0) { - double norm = ope_table.at(row).at(y); + double row_norm = ope_table.at(row).at(y); for (size_t col = 0; col < ope_table.at(x).size(); col++) { - ope_table.at(row).at(col) = ope_table.at(row).at(col) - norm * ope_table.at(x).at(col); + ope_table.at(row).at(col) -= row_norm * ope_table.at(x).at(col); } } } - basic.at(x - 1) = y; // 换元 + basic.at(x - 1) = static_cast(y); } -