Add artificial bound for free variables

Author: filikat

highs/ipm/hipo/ipm/Model.cpp

@@ -277,15 +277,12 @@ void Model::print(const Logger& logger) const {
   // compute max and min for bounds intervals
   double boundmin = kHighsInf;
   double boundmax = 0.0;
-  Int free_vars = 0;
   for (Int i = 0; i < n_; ++i) {
     if (std::isfinite(lower_[i]) && std::isfinite(upper_[i])) {
       const double diff = std::abs(upper_[i] - lower_[i]);
       boundmin = std::min(boundmin, diff);
       boundmax = std::max(boundmax, diff);
     }
-
-    if (!std::isfinite(lower_[i]) && !std::isfinite(upper_[i])) free_vars++;
   }
   if (std::isinf(boundmin)) boundmin = 0.0;
 
@@ -345,7 +342,6 @@ void Model::print(const Logger& logger) const {
 
   if (logger.debug(1)) {
     preprocessor_.print(log_stream);
-    log_stream << textline("Free variables:") << integer(free_vars) << '\n';
 
     log_stream << "Expected nnz: ";
     log_stream << "AS " << sci(AS_nz_, 0, 1) << "; ";
@@ -421,4 +417,32 @@ void Model::printDense() const {
   printf("offset %6.2f\n", offset_);
 }
 
+void Model::adjustFreeVars(std::vector<double>& x, std::vector<double>& xl,
+                           std::vector<double>& xu, const Logger& logger) {
+  for (Int i = 0; i < n_; ++i) {
+    if (is_free_[i]) {
+      if (x[i] < lower_[i] * kFreeVarsCloseRatio) {
+        // getting close to lower bound
+        const double new_lower = std::min(lower_[i] * kFreeVarsIncreaseBound,
+                                          x[i] / kFreeVarsCloseRatio);
+        logger.printDetailed(
+            "Free var %d is at %.1e with lb %.1e, lb changed to %.1e\n", i,
+            x[i], lower_[i], new_lower);
+        lower_[i] = new_lower;
+        xl[i] = x[i] - lower_[i];
+      }
+      if (x[i] > upper_[i] * kFreeVarsCloseRatio) {
+        // getting close to upper bound
+        const double new_upper = std::max(upper_[i] * kFreeVarsIncreaseBound,
+                                          x[i] / kFreeVarsCloseRatio);
+        logger.printDetailed(
+            "Free var %d is at %.1e with ub %.1e, ub changed to %.1e\n", i,
+            x[i], upper_[i], new_upper);
+        upper_[i] = new_upper;
+        xu[i] = upper_[i] - x[i];
+      }
+    }
+  }
+}
+
 }  // namespace hipo

highs/ipm/hipo/ipm/Model.h

@@ -5,9 +5,9 @@
 #include <string>
 #include <vector>
 
-#include "ipm/hipo/auxiliary/Logger.h"
 #include "PreProcess.h"
 #include "ipm/hipo/auxiliary/IntConfig.h"
+#include "ipm/hipo/auxiliary/Logger.h"
 #include "ipm/ipx/lp_solver.h"
 #include "lp_data/HighsLp.h"
 #include "model/HighsHessian.h"
@@ -65,6 +65,8 @@ class Model {
   std::vector<double> one_norm_cols_, one_norm_rows_, inf_norm_cols_,
       inf_norm_rows_;
 
+  std::vector<bool> is_free_;
+
   void preprocess();
   Int checkData() const;
   void computeNorms();
@@ -90,6 +92,9 @@ class Model {
 
   void nzBounds();
 
+  void adjustFreeVars(std::vector<double>& x, std::vector<double>& xl,
+                      std::vector<double>& xu, const Logger& logger);
+
   // Check if variable has finite lower/upper bound
   bool hasLb(Int j) const { return std::isfinite(lower_[j]); }
   bool hasUb(Int j) const { return std::isfinite(upper_[j]); }
@@ -123,6 +128,7 @@ class Model {
   Int64 nzAS() const { return AS_nz_; }
   Int64 nzNElb() const { return NE_nz_lb_; }
   Int64 nzNEub() const { return NE_nz_ub_; }
+  bool free(Int i) const { return is_free_[i]; }
 
   Int loadIntoIpx(ipx::LpSolver& lps) const;
 
@@ -131,6 +137,7 @@ class Model {
   friend struct PreprocessFixedVars;
   friend struct PreprocessScaling;
   friend struct PreprocessFormulation;
+  friend struct PreprocessFreeVars;
 };
 
 }  // namespace hipo

highs/ipm/hipo/ipm/Parameters.h

@@ -51,6 +51,11 @@ const double kSmallScalingCoeff = 1e-4;
 const double kLargeScalingCoeff = 1e4;
 const double kSmallBoundDiff = 1e-3;
 
+// parameters for free variables
+const double kFreeVarsInitialBound = 1e4;
+const double kFreeVarsCloseRatio = 0.5;
+const double kFreeVarsIncreaseBound = 10.0;
+
 // static regularisation
 struct Regularisation {
   double primal = 1e-12;

highs/ipm/hipo/ipm/PreProcess.cpp

@@ -483,11 +483,11 @@ void PreprocessScaling::undo(PreprocessorPoint& point, const Model& model,
 
     // set variables that were ignored
     for (Int i = 0; i < n_post; ++i) {
-      if (!model.hasLb(i)) {
+      if (!model.hasLb(i) || model.is_free_[i]) {
         point.xl[i] = kHighsInf;
         point.zl[i] = 0.0;
       }
-      if (!model.hasUb(i)) {
+      if (!model.hasUb(i) || model.is_free_[i]) {
         point.xu[i] = kHighsInf;
         point.zu[i] = 0.0;
       }
@@ -561,11 +561,11 @@ void PreprocessFormulation::undo(PreprocessorPoint& point, const Model& model,
 
   // force unused entries to have correct value
   for (Int i = 0; i < n_pre; ++i) {
-    if (!model.hasLb(i)) {
+    if (!model.hasLb(i) || model.is_free_[i]) {
       point.xl[i] = kHighsInf;
       point.zl[i] = 0.0;
     }
-    if (!model.hasUb(i)) {
+    if (!model.hasUb(i) || model.is_free_[i]) {
       point.xu[i] = kHighsInf;
       point.zu[i] = 0.0;
     }
@@ -617,6 +617,49 @@ void PreprocessFormulation::print(std::stringstream& stream) const {
   stream << "Added " << n_post - n_pre << " slacks\n";
 }
 
+void PreprocessFreeVars::apply(Model& model) {
+  Int& n = model.n_;
+  Int& m = model.m_;
+  HighsSparseMatrix& A = model.A_;
+  std::vector<double>& b = model.b_;
+  std::vector<double>& c = model.c_;
+  std::vector<double>& lower = model.lower_;
+  std::vector<double>& upper = model.upper_;
+  std::vector<char>& constraints = model.constraints_;
+  HighsHessian& Q = model.Q_;
+  std::vector<bool>& is_free = model.is_free_;
+
+  n_pre = n;
+  m_pre = m;
+
+  is_free.assign(n, false);
+
+  for (Int i = 0; i < n; ++i) {
+    if (!std::isfinite(lower[i]) && !std::isfinite(upper[i]) &&
+        lower[i] != upper[i]) {
+      // free variable
+      is_free[i] = true;
+      ++free_vars;
+      lower[i] = -kFreeVarsInitialBound;
+      upper[i] = kFreeVarsInitialBound;
+    }
+  }
+
+  n_post = n;
+  m_post = m;
+}
+
+void PreprocessFreeVars::undo(PreprocessorPoint& point, const Model& model,
+                              const Iterate& it) const {
+  point.assertConsistency(n_post, m_post);
+  //
+  point.assertConsistency(n_pre, m_pre);
+}
+
+void PreprocessFreeVars::print(std::stringstream& stream) const {
+  if (free_vars > 0) stream << "Found " << free_vars << " free variables\n";
+}
+
 #define APPLY_ACTION(T)                                      \
   stack.push_back(std::unique_ptr<PreprocessAction>(new T)); \
   stack.back()->apply(model);
@@ -629,6 +672,7 @@ void Preprocessor::apply(Model& model) {
   APPLY_ACTION(PreprocessEmptyRows);
   APPLY_ACTION(PreprocessScaling);
   APPLY_ACTION(PreprocessFormulation);
+  APPLY_ACTION(PreprocessFreeVars);
 }
 void Preprocessor::undo(PreprocessorPoint& point, const Model& model,
                         const Iterate& it) const {

highs/ipm/hipo/ipm/PreProcess.h

@@ -80,6 +80,15 @@ struct PreprocessFormulation : public PreprocessAction {
   void print(std::stringstream& stream) const override;
 };
 
+struct PreprocessFreeVars : public PreprocessAction {
+  Int free_vars{};
+
+  void apply(Model& model) override;
+  void undo(PreprocessorPoint& point, const Model& model,
+            const Iterate& it) const override;
+  void print(std::stringstream& stream) const override;
+};
+
 struct Preprocessor {
   std::vector<std::unique_ptr<PreprocessAction>> stack;
 

highs/ipm/hipo/ipm/Solver.cpp

@@ -191,6 +191,8 @@ bool Solver::prepareIter() {
 
   ++iter_;
 
+  model_.adjustFreeVars(it_->x, it_->xl, it_->xu, logger_);
+
   // Clear Newton direction
   it_->clearDir();