Merge pull request #2053 from ERGO-Code/fix-1873

Added blending and lexicographic handling of multiple linear objectives

FEATURES.md

@@ -2,39 +2,13 @@
 
 ## Code changes
 
-When primal infeasiblity is detected in presolve, no dual ray is available so, previously, the `has_dual_ray` parameter of `Highs::getDualRay` returned false and that was it. Now, if a null pointer is not passed for `dual_ray_value`, `Highs::getDualRay` will compute a dual ray - at the cost of solving the feasiblility LP without presolve. The same is now true for `Highs::getPrimalRay`. `Highs::getDualUnboundednessDirection` has been introduced to determine the product between the constraint matrix and the dual ray, forcing the calculation of the latter if necessary. Once a dual ray is known for the incumbent model in HiGHS, subsequent calls to `Highs::getDualRay` and `Highs::getDualUnboundednessDirection` will be vastly cheaper
-
-The method `Highs::getDualObjectiveValue` now exitsts to compute the dual objective value, returning `HighsStatus::kError` if it is not possible.
-
-The method `Highs::getStandardFormLp` now exists to return the incumbent LP in standard form - overlooking any integrality or Hessian. To determine the sizes of the vectors of data, the method is called without specifying pointers to the data arrays.
-
-Added documentation on the use of presolve when solving an incumbent model, and clarifying the use of the method `Highs::presolve`.
-
-HiGHS will now read a `MIPLIB` solution file
-
-Added time limit check to `HPresolve::strengthenInequalities`
-
-Added `getColIntegrality` to `highspy`
-
-Now computing the primal-dual integral, reporting it, and making it available as `HighsInfo::primal_dual_integral`
-
-Trivial primal heuristics "all zero", "all lower bound", "all upper bound", and "lock point" added to the MIP solver
-
-# Build changes
-
-Added wheels for Python 3.13 
-
-Updated command line options and moved them out of the library and into the executable 
-
-
-
-
-
-
-
+HiGHS now handles multiple linear objectives by either blending using weights, or performing lexicographic optimization: see https://ergo-code.github.io/HiGHS/stable/guide/further/#guide-multi-objective-optimization
 
+Fixed minor bug in bound checking in presolve
 
+Fixed bug in `floor(HighsCDouble x)` and `ceil(HighsCDouble x)` when argument is small
 
+Added some sanity checks to Highs::writeLocalModel to prevent segfaults if called directly by a user
 
 
 

check/CMakeLists.txt

@@ -48,6 +48,7 @@ if ((NOT FAST_BUILD OR ALL_TESTS) AND NOT (BUILD_EXTRA_UNIT_ONLY))
       TestBasis.cpp
       TestBasisSolves.cpp
       TestCrossover.cpp
+      TestHighsCDouble.cpp
       TestHighsHash.cpp
       TestHighsIntegers.cpp
       TestHighsParallel.cpp
@@ -66,6 +67,7 @@ if ((NOT FAST_BUILD OR ALL_TESTS) AND NOT (BUILD_EXTRA_UNIT_ONLY))
       TestLpModification.cpp
       TestLpOrientation.cpp
       TestModelProperties.cpp
+      TestMultiObjective.cpp
       TestPdlp.cpp
       TestPresolve.cpp
       TestQpSolver.cpp

check/TestCAPI.c

@@ -1855,6 +1855,118 @@ void test_getModel() {
   Highs_destroy(highs);
 }
 
+void test_multiObjective() {
+  void* highs;
+  highs = Highs_create();
+  const double inf = Highs_getInfinity(highs);
+
+  HighsInt num_col = 2;
+  HighsInt num_row = 3;
+  HighsInt num_nz = num_col * num_row;
+  HighsInt a_format = kHighsMatrixFormatColwise;
+  HighsInt sense = kHighsObjSenseMaximize;
+  double offset = -1;
+  double col_cost[2] = {1, 1};
+  double col_lower[2] = {0, 0};
+  double col_upper[2] = {inf, inf};
+  double row_lower[3] = {-inf, -inf, -inf};
+  double row_upper[3] = {18, 8, 14};
+  HighsInt a_start[3] = {0, 3, 6};
+  HighsInt a_index[6] = {0, 1, 2, 0, 1, 2};
+  double a_value[6] = {3, 1, 1, 1, 1, 2};
+  HighsInt integrality[2] = {kHighsVarTypeInteger, kHighsVarTypeInteger};
+
+  Highs_setBoolOptionValue(highs, "output_flag", dev_run);
+  HighsInt return_status = Highs_passLp(highs, num_col, num_row, num_nz, a_format, sense,
+					offset, col_cost, col_lower, col_upper,
+					row_lower, row_upper, a_start, a_index, a_value);
+  assert(return_status == kHighsStatusOk);
+
+  return_status = Highs_clearLinearObjectives(highs);
+  assert(return_status == kHighsStatusOk);
+
+  double weight = -1;
+  double linear_objective_offset = -1;
+  double coefficients[2] = {1, 1};
+  double abs_tolerance = 0;
+  double rel_tolerance = 0;
+  HighsInt priority = 10;
+  return_status = Highs_addLinearObjective(highs, weight, linear_objective_offset, coefficients, abs_tolerance, rel_tolerance, priority);
+  assert(return_status == kHighsStatusOk);
+
+  weight = 1e-4;
+  linear_objective_offset = 0;
+  coefficients[0] = 1;
+  coefficients[1] = 0;
+  priority = 0;
+  return_status = Highs_addLinearObjective(highs, weight, linear_objective_offset, coefficients, abs_tolerance, rel_tolerance, priority);
+  assert(return_status == kHighsStatusOk);
+
+  return_status = Highs_run(highs);
+  assert(return_status == kHighsStatusOk);
+  HighsInt model_status = Highs_getModelStatus(highs);
+  assert(model_status == kHighsModelStatusOptimal);
+
+  Highs_writeSolutionPretty(highs, "");
+  double* col_value = (double*)malloc(sizeof(double) * num_col);
+  return_status =
+    Highs_getSolution(highs, col_value, NULL, NULL, NULL);
+  assertDoubleValuesEqual("col_value[0]", col_value[0], 2);
+  assertDoubleValuesEqual("col_value[1]", col_value[1], 6);
+
+  Highs_setBoolOptionValue(highs, "blend_multi_objectives", 0);
+
+  if (dev_run) printf("\n***************\nLexicographic 1\n***************\n");
+  double weight2[2] = {-1, 1e-4};
+  double linear_objective_offset2[2] = {-1, 0};
+  double coefficients2[4] = {1, 1, 1, 0};
+  double abs_tolerance2[2] = {0, -1};
+  double rel_tolerance2[2] = {0, -1};
+  HighsInt priority2[2] = {10, 0};
+  return_status = Highs_passLinearObjectives(highs, 2, weight2, linear_objective_offset2, coefficients2, abs_tolerance2, rel_tolerance2, priority2);
+  return_status = Highs_run(highs);
+  assert(return_status == kHighsStatusOk);
+  model_status = Highs_getModelStatus(highs);
+  assert(model_status == kHighsModelStatusOptimal);
+  Highs_writeSolutionPretty(highs, "");
+  return_status =
+    Highs_getSolution(highs, col_value, NULL, NULL, NULL);
+  assertDoubleValuesEqual("col_value[0]", col_value[0], 2);
+  assertDoubleValuesEqual("col_value[1]", col_value[1], 6);
+
+  //  weight2[1] = 1e-5;
+  coefficients2[0] = 1.0001;
+  abs_tolerance2[0] = 1e-5;
+  rel_tolerance2[0] = 0.05;
+  return_status = Highs_passLinearObjectives(highs, 2, weight2, linear_objective_offset2, coefficients2, abs_tolerance2, rel_tolerance2, priority2);
+  return_status = Highs_run(highs);
+  assert(return_status == kHighsStatusOk);
+  model_status = Highs_getModelStatus(highs);
+  assert(model_status == kHighsModelStatusOptimal);
+  Highs_writeSolutionPretty(highs, "");
+  return_status =
+    Highs_getSolution(highs, col_value, NULL, NULL, NULL);
+  assertDoubleValuesEqual("col_value[0]", col_value[0], 4.9);
+  assertDoubleValuesEqual("col_value[1]", col_value[1], 3.1);
+
+  if (dev_run) printf("\n***************\nLexicographic 2\n***************\n");
+  abs_tolerance2[0] = -1;
+
+  return_status = Highs_passLinearObjectives(highs, 2, weight2, linear_objective_offset2, coefficients2, abs_tolerance2, rel_tolerance2, priority2);
+  return_status = Highs_run(highs);
+  assert(return_status == kHighsStatusOk);
+  model_status = Highs_getModelStatus(highs);
+  assert(model_status == kHighsModelStatusOptimal);
+  Highs_writeSolutionPretty(highs, "");
+  return_status =
+    Highs_getSolution(highs, col_value, NULL, NULL, NULL);
+  assertDoubleValuesEqual("col_value[0]", col_value[0], 1.30069);
+  assertDoubleValuesEqual("col_value[1]", col_value[1], 6.34966);
+
+  Highs_destroy(highs);
+  free(col_value);
+}
+
 /*
 The horrible C in this causes problems in some of the CI tests,
 so suppress thius test until the C has been improved
@@ -1919,6 +2031,7 @@ int main() {
   test_ranging();
   test_feasibilityRelaxation();
   test_getModel();
+  test_multiObjective();
   return 0;
 }
 //  test_setSolution();

check/TestHighsCDouble.cpp

@@ -0,0 +1,109 @@
+#include "HCheckConfig.h"
+#include "catch.hpp"
+#include "util/HighsCDouble.h"
+#include "util/HighsRandom.h"
+
+void testCeil(HighsCDouble x) {
+  double ceil_x;
+  double double_x;
+  ceil_x = double(ceil(x));
+  double_x = double(x);
+  REQUIRE(ceil_x >= double_x);
+  REQUIRE(ceil(x) >= x);
+}
+
+void testFloor(HighsCDouble x) {
+  double floor_x;
+  double double_x;
+  floor_x = double(floor(x));
+  double_x = double(x);
+  REQUIRE(floor_x <= double_x);
+  REQUIRE(floor(x) <= x);
+}
+
+TEST_CASE("HighsCDouble-ceil", "[util]") {
+  HighsCDouble x;
+  x = -1e-34;
+  testCeil(x);
+  x = -1e-32;
+  testCeil(x);
+  x = -1e-30;
+  testCeil(x);
+  x = -1e-23;
+  testCeil(x);
+  x = -1e-12;
+  testCeil(x);
+  x = -1e-1;
+  testCeil(x);
+  x = -0.99;
+  testCeil(x);
+
+  x = 0.99;
+  testCeil(x);
+  x = 1e-1;
+  testCeil(x);
+  x = 1e-12;
+  testCeil(x);
+  // This and rest failed in #2041
+  x = 1e-23;
+  testCeil(x);
+  x = 1e-30;
+  testCeil(x);
+  x = 1e-32;
+  testCeil(x);
+  x = 1e-34;
+  testCeil(x);
+
+  HighsRandom rand;
+  for (HighsInt k = 0; k < 1000; k++) {
+    double man = rand.fraction();
+    HighsInt power = 2 - rand.integer(5);
+    double exp = std::pow(10, power);
+    x = man * exp;
+    testCeil(x);
+  }
+}
+
+TEST_CASE("HighsCDouble-floor", "[util]") {
+  HighsCDouble x;
+
+  x = 1e-34;
+  testFloor(x);
+  x = 1e-32;
+  testFloor(x);
+  x = 1e-30;
+  testFloor(x);
+  x = 1e-23;
+  testFloor(x);
+  x = 1e-12;
+  testFloor(x);
+  x = 1e-1;
+  testFloor(x);
+  x = 0.99;
+  testFloor(x);
+
+  x = -0.99;
+  testFloor(x);
+  x = -1e-1;
+  testFloor(x);
+  x = -1e-12;
+  testFloor(x);
+  // This and rest failed in #2041
+  x = -1e-23;
+  testFloor(x);
+  x = -1e-30;
+  testFloor(x);
+  x = -1e-32;
+  testFloor(x);
+  x = -1e-34;
+  testFloor(x);
+
+  HighsRandom rand;
+  for (HighsInt k = 0; k < 1000; k++) {
+    double man = rand.fraction();
+    HighsInt power = 2 - rand.integer(5);
+    double exp = std::pow(10, power);
+    x = man * exp;
+    testFloor(x);
+  }
+}