diff --git a/include/mp/flat/constr_algebraic.h b/include/mp/flat/constr_algebraic.h index c45882f75..337ba6d6b 100644 --- a/include/mp/flat/constr_algebraic.h +++ b/include/mp/flat/constr_algebraic.h @@ -70,11 +70,12 @@ class AlgebraicConstraint : return Body::ComputeValue(x) - RhsOrRange::lb(); } - /// Compute violation + /// Compute violation. + /// Negative if if holds with slack. double ComputeViolation(const ArrayRef& x) const { auto bd = Body::ComputeValue(x); - return std::max( + return std::max( // same for strict cmp? RhsOrRange::lb() - bd, bd - RhsOrRange::ub()); } diff --git a/include/mp/flat/constr_base.h b/include/mp/flat/constr_base.h index 0d7216ef7..3aa8229ce 100644 --- a/include/mp/flat/constr_base.h +++ b/include/mp/flat/constr_base.h @@ -190,14 +190,6 @@ double ComputeViolation( return std::fabs(x[resvar] - ComputeValue(c, x)); } -template -template -double -CustomFunctionalConstraint -::ComputeViolation(const VarVec& x) const -{ return mp::ComputeViolation(*this, x); } - /// A base class for numerical functional constraint. /// It provides default properties of such a constraint @@ -286,6 +278,27 @@ class ConditionalConstraint : return GetConstraint()==cc.GetConstraint(); } + /// Compute violation of a conditional constraint. + /// If the subconstr is violated but should hold, + /// return the exact gap, despite this is a logical constraint. + /// If the subconstr holds but should not, + /// return the opposite gap. Similar to Indicator. + template + double ComputeViolation(const VarVec& x) { + auto viol = GetConstraint().ComputeViolation(x); + bool ccon_valid = viol<=0.0; + bool has_arg = x[GetResultVar()] >= 0.5; + switch (this->GetContext().GetValue()) { + case Context::CTX_MIX: // Viol is non-positive if holds + return has_arg == ccon_valid ? 0.0 : std::fabs(viol); + case Context::CTX_POS: + return has_arg <= ccon_valid ? 0.0 : viol; + case Context::CTX_NEG: + return has_arg >= ccon_valid ? 0.0 : -viol; + default: + return INFINITY; + } + } }; //////////////////////////////////////////////////////////////////////// diff --git a/include/mp/flat/constr_eval.h b/include/mp/flat/constr_eval.h index b1d9c73a3..cf8639e02 100644 --- a/include/mp/flat/constr_eval.h +++ b/include/mp/flat/constr_eval.h @@ -235,29 +235,117 @@ double ComputeValue(const TanhConstraint& con, const VarVec& x) { return std::tanh(x[con.GetArguments()[0]]); } +/// Compute result of the asinh constraint. +template +double ComputeValue(const AsinhConstraint& con, const VarVec& x) { + return std::asinh(x[con.GetArguments()[0]]); +} + +/// Compute result of the acosh constraint. +template +double ComputeValue(const AcoshConstraint& con, const VarVec& x) { + return std::acosh(x[con.GetArguments()[0]]); +} + +/// Compute result of the atanh constraint. +template +double ComputeValue(const AtanhConstraint& con, const VarVec& x) { + return std::atanh(x[con.GetArguments()[0]]); +} + /// Compute result of a conditional constraint. -/// Actually, for violation itself, we could do this: -/// If the subconstr is violated but should hold, -/// return the exact gap, despite this is a logical constraint. -/// If the subconstr holds but should not, return 0.0. +/// This is not used to compute violation. template double ComputeValue( const ConditionalConstraint& con, const VarVec& x) { auto viol = con.GetConstraint().ComputeViolation(x); bool ccon_valid = viol<=0.0; bool has_arg = x[con.GetResultVar()] >= 0.5; - switch (con.GetContext()) { + switch (con.GetContext().GetValue()) { case Context::CTX_MIX: return has_arg == ccon_valid; case Context::CTX_POS: - return has_arg < ccon_valid; + return has_arg <= ccon_valid; case Context::CTX_NEG: - return has_arg > ccon_valid; + return has_arg >= ccon_valid; default: return INFINITY; } } +/// Compute violation of the QuadraticCone constraint. +template +double ComputeViolation( + const QuadraticConeConstraint& con, const VarVec& x) { + const auto& args = con.GetArguments(); + const auto& params = con.GetParameters(); + assert(args.size()==params.size()); + double sum = 0.0; + for (auto i=args.size(); --i; ) + sum += std::pow( params[i]*x[args[i]], 2.0 ); + return std::sqrt(sum) - params[0]*x[args[0]]; +} + +/// Compute violation of the RotatedQuadraticCone constraint. +template +double ComputeViolation( + const RotatedQuadraticConeConstraint& con, const VarVec& x) { + const auto& args = con.GetArguments(); + const auto& params = con.GetParameters(); + assert(args.size()==params.size()); + double sum = 0.0; + for (auto i=args.size(); --i>1; ) + sum += std::pow( params[i]*x[args[i]], 2.0 ); + return sum + - 2.0 * params[0]*x[args[0]] * params[1]*x[args[1]]; +} + +/// Compute violation of the ExponentialCone constraint. +template // ax >= by exp(cz / (by)) +double ComputeViolation( // where ax, by >= 0 + const ExponentialConeConstraint& con, const VarVec& x) { + const auto& args = con.GetArguments(); + const auto& params = con.GetParameters(); + auto ax = params[0]*x[args[0]]; + auto by = params[1]*x[args[1]]; + if (0.0==std::fabs(by)) + return -ax; + auto cz = params[2]*x[args[2]]; + return by * std::exp(cz / by) - ax; +} + +/// Compute result of the PL constraint. +template +double ComputeValue(const PLConstraint& con, const VarVec& x) { + const auto& plp = con.GetParameters().GetPLPoints(); + assert(!plp.empty()); + auto x0 = x[con.GetArguments()[0]]; // position + if (x0plp.x_.back()) + return plp.y_.back() + + plp.PostSlope()*(x0 - plp.x_.back()); + int i0=0; + for ( ; x0 > plp.x_[i0]; ++i0); + return plp.x_[i0]==x0 + ? plp.y_[i0] + : plp.y_[i0] + + (plp.y_[i0+1]-plp.y_[i0]) + * (x0-plp.x_[i0]) / (plp.x_[i0+1]-plp.x_[i0]); +} + +/// Should be here, +/// after ComputeViolation() is specialized +/// for some constraints. +template +template +double +CustomFunctionalConstraint +::ComputeViolation(const VarVec& x) const +{ return mp::ComputeViolation(*this, x); } + } // namespace mp diff --git a/include/mp/flat/constr_functional.h b/include/mp/flat/constr_functional.h index b7f7d5f2d..83622c8de 100644 --- a/include/mp/flat/constr_functional.h +++ b/include/mp/flat/constr_functional.h @@ -359,11 +359,15 @@ class PLConParams { /// Construct from PLPoints PLConParams(PLPoints plp) : plp_(std::move(plp)) { } /// Produce PLPoints, either stored or from PLSlopes - operator PLPoints() const { return plp_.empty() ? pls_ : plp_; } + const PLPoints& GetPLPoints() const { + if (plp_.empty()) + plp_ = pls_; + return plp_; + } private: PLSlopes pls_; - PLPoints plp_; + mutable PLPoints plp_; }; diff --git a/include/mp/flat/constr_general.h b/include/mp/flat/constr_general.h index 49bab9786..41b0d6505 100644 --- a/include/mp/flat/constr_general.h +++ b/include/mp/flat/constr_general.h @@ -223,7 +223,8 @@ DEF_STATIC_CONSTR_WITH_PRM( QuadraticConeConstraint, VarArray, DblParamArray, " with factors p1..pn"); /// Rotated quadratic cone DEF_STATIC_CONSTR_WITH_PRM( RotatedQuadraticConeConstraint, VarArray, DblParamArray, - "Rotated quadratic cone p1*x1*p2*x2 >= sqrt((p3*x3)^2 + ...))," + "Rotated quadratic cone " + "2 * p1*x1*p2*x2 >= (p3*x3)^2 + ...)," " x1, x2 >= 0, with factors p1..pn"); /// Exponential cone DEF_STATIC_CONSTR_WITH_PRM( ExponentialConeConstraint, VarArray3, DblParamArray3, diff --git a/include/mp/flat/redef/MIP/piecewise_linear.h b/include/mp/flat/redef/MIP/piecewise_linear.h index cf90cccda..2cf51cb71 100644 --- a/include/mp/flat/redef/MIP/piecewise_linear.h +++ b/include/mp/flat/redef/MIP/piecewise_linear.h @@ -29,7 +29,7 @@ class PLConverter_MIP : /// Convert in any context void Convert(const ItemType& cc, int ) { - points_ = cc.GetParameters(); // extract PLPoints + points_ = cc.GetParameters().GetPLPoints(); i0=0; // first breakpoint i1=points_.x_.size()-1; // last breakpoint y = cc.GetResultVar(); diff --git a/solvers/cplexmp/cplexmpmodelapi.cc b/solvers/cplexmp/cplexmpmodelapi.cc index cfd55c5e2..1929f5877 100644 --- a/solvers/cplexmp/cplexmpmodelapi.cc +++ b/solvers/cplexmp/cplexmpmodelapi.cc @@ -141,7 +141,7 @@ void CplexModelAPI::AddConstraint(const IndicatorConstraintLinGE &ic) { } void CplexModelAPI::AddConstraint(const PLConstraint& plc) { - PLPoints plp(plc.GetParameters()); + const auto& plp = plc.GetParameters().GetPLPoints(); CPLEX_CALL( CPXaddpwl(env(), lp(), plc.GetResultVar(), plc.GetArguments()[0], plp.PreSlope(), plp.PostSlope(), diff --git a/solvers/gurobi/gurobimodelapi.cc b/solvers/gurobi/gurobimodelapi.cc index b8fafc286..88ae936cc 100644 --- a/solvers/gurobi/gurobimodelapi.cc +++ b/solvers/gurobi/gurobimodelapi.cc @@ -225,7 +225,7 @@ void GurobiModelAPI::AddConstraint(const TanConstraint &cc) { } void GurobiModelAPI::AddConstraint(const PLConstraint& plc) { - PLPoints plp(plc.GetParameters()); + const auto& plp = plc.GetParameters().GetPLPoints(); GRB_CALL( GRBaddgenconstrPWL(model(), plc.name(), plc.GetArguments()[0], plc.GetResultVar(), plp.x_.size(), plp.x_.data(), plp.y_.data()) ); diff --git a/test/end2end/cases/categorized/fast/nonlinear/expA_1.mod b/test/end2end/cases/categorized/fast/nonlinear/expA_1.mod index 529648da0..fc5cf613b 100644 --- a/test/end2end/cases/categorized/fast/nonlinear/expA_1.mod +++ b/test/end2end/cases/categorized/fast/nonlinear/expA_1.mod @@ -4,14 +4,14 @@ # expA_1.mod # ------------------------------------------------------------- -param ubx integer := 10; -param uby integer := 20; +param ubx integer := 100; +param uby integer := 200; -var x >= 0, <= ubx; -var y >= -41, <= uby; +var x >= -100, <= ubx; +var y >= -205, <= uby; maximize TotalSum: - x-y; + -x-y; subj to ExpA: y >= 0.1 ^ x; diff --git a/test/end2end/cases/categorized/fast/nonlinear/modellist.json b/test/end2end/cases/categorized/fast/nonlinear/modellist.json index 2cd4e1abc..0cea157b0 100644 --- a/test/end2end/cases/categorized/fast/nonlinear/modellist.json +++ b/test/end2end/cases/categorized/fast/nonlinear/modellist.json @@ -3,7 +3,7 @@ "name" : "expA_1", "tags" : ["nonlinear"], "values": { - "if abs(_obj[1]-10)<=0.01 then 1 else 0": 1 + "if abs(_obj[1]+0.7977905982)<=0.01 then 1 else 0": 1 } }, { @@ -11,7 +11,7 @@ "tags" : ["nonlinear", "check_pl_approx_expA"], "options": { "ANYSOLVER_options": "acc:expa=1" }, "values": { - "if abs(_obj[1]-10)<=0.01 then 1 else 0": 1 + "if abs(_obj[1]+0.7977905982)<=0.01 then 1 else 0": 1 } }, {