Improve bounds for Pow() and quadratics

ampl · Mar 1, 2023 · 2615da3 · 2615da3 · glebbelov · Mar 1, 2023
1 parent 970405f
commit 2615da3
Show file tree

Hide file tree

Showing 3 changed files with 32 additions and 18 deletions.
diff --git a/include/mp/flat/constr_prepro.h b/include/mp/flat/constr_prepro.h
@@ -55,19 +55,24 @@ class ConstraintPreprocessors {
     auto& m = MP_DISPATCH( GetModel() );
     auto lb = std::pow(m.lb(arg), pwr),
         ub = std::pow(m.ub(arg), pwr);
-    if (MPD( is_integer_value(pwr) ) && pwr>=0.0) {
-      // result integer if arg is and integer, >=0 exponent
-      prepro.set_result_type( m.var_type(arg) );
-      if (MPD( is_integer_value(pwr / 2.0) )) {  // exponent is even, >=0
-        bool lb_neg = m.lb(arg)<0.0;
-        bool ub_pos = m.ub(arg)>0.0;
-        if (lb_neg && ub_pos) {
-          ub = std::max(lb, ub); lb = 0.0;
-        } else if (lb_neg) {
-          std::swap(lb, ub);
-        }
-      }
-    }
+		if (MPD( is_integer_value(pwr) )) {
+			if (pwr>=0.0) {
+				// result integer if arg is and integer, >=0 exponent
+				prepro.set_result_type( m.var_type(arg) );
+				if (MPD( is_integer_value(pwr / 2.0) )) {  // exponent is even, >=0
+					bool lb_neg = m.lb(arg)<0.0;
+					bool ub_pos = m.ub(arg)>0.0;
+					if (lb_neg && ub_pos) {
+						ub = std::max(lb, ub); lb = 0.0;
+					} else if (lb_neg) {
+						std::swap(lb, ub);
+					}
+				}
+			}
+		} else {                      // fractional power
+			if (lb<0.0)
+				lb = 0.0;
+		}
     prepro.narrow_result_bounds( std::min(lb, ub),
                           std::max(lb, ub) );
   }

diff --git a/include/mp/flat/converter.h b/include/mp/flat/converter.h
@@ -194,8 +194,10 @@ class FlatConverter :
 	void FixUnusedDefinedVars() {
 		for (auto i=num_vars(); i--; ) {
 			if (HasInitExpression(i) &&
-					! VarUsageRef(i))
-				set_var_ub(i, lb(i));
+					! VarUsageRef(i)) {
+				set_var_lb(i, 0.0);      // fix to 0
+				set_var_ub(i, 0.0);
+			}
 		}
 	}
 

diff --git a/include/mp/flat/expr_bounds.h b/include/mp/flat/expr_bounds.h
@@ -91,9 +91,16 @@ class BoundComputations {
   std::pair<double, double> ProductBounds(Var x, Var y) const {
     const auto& m = MPCD(GetModel());
     auto lx=m.lb(x), ly=m.lb(y), ux=m.ub(x), uy=m.ub(y);
-    std::array<double, 4> pb{lx*ly, lx*uy, ux*ly, ux*uy};
-    return { *std::min_element(pb.begin(), pb.end()),
-          *std::max_element(pb.begin(), pb.end()) };
+		if (x!=y) {                        // different vars
+			std::array<double, 4> pb{lx*ly, lx*uy, ux*ly, ux*uy};
+			return { *std::min_element(pb.begin(), pb.end()),
+						*std::max_element(pb.begin(), pb.end()) };
+		}
+		return {
+			lx<=0.0 && ux>=0.0 ?   // lb_result: 0 included?
+						0.0 : std::min(lx*lx, ux*ux),
+			std::max(lx*lx, ux*ux)
+		};
   }
 
   /// Add / merge bounds and type