diff --git a/FEATURES.md b/FEATURES.md
index 034eca4232..76f4c4997d 100644
--- a/FEATURES.md
+++ b/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
 
 
 
diff --git a/check/CMakeLists.txt b/check/CMakeLists.txt
index f8bb62d133..700b5eee36 100644
--- a/check/CMakeLists.txt
+++ b/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
diff --git a/check/TestCAPI.c b/check/TestCAPI.c
index e5848f67a1..d75cc71305 100644
--- a/check/TestCAPI.c
+++ b/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();
diff --git a/check/TestHighsCDouble.cpp b/check/TestHighsCDouble.cpp
new file mode 100644
index 0000000000..2f17a40598
--- /dev/null
+++ b/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);
+  }
+}
diff --git a/check/TestHighsIntegers.cpp b/check/TestHighsIntegers.cpp
index b56ec69908..c08b288610 100644
--- a/check/TestHighsIntegers.cpp
+++ b/check/TestHighsIntegers.cpp
@@ -1,8 +1,8 @@
 #include "HCheckConfig.h"
 #include "catch.hpp"
-#include "util/HighsCDouble.h"
+// #include "util/HighsCDouble.h"
 #include "util/HighsIntegers.h"
-#include "util/HighsRandom.h"
+// #include "util/HighsRandom.h"
 
 const bool dev_run = false;
 
@@ -40,107 +40,3 @@ TEST_CASE("HighsIntegers", "[util]") {
 
   if (dev_run) printf("integral scalar is %g\n", integralscalar);
 }
-
-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);
-  }
-}
diff --git a/check/TestMultiObjective.cpp b/check/TestMultiObjective.cpp
new file mode 100644
index 0000000000..1a15e8d9d2
--- /dev/null
+++ b/check/TestMultiObjective.cpp
@@ -0,0 +1,196 @@
+#include "Highs.h"
+#include "catch.hpp"
+
+const bool dev_run = false;
+
+bool smallDoubleDifference(double v0, double v1) {
+  double difference = std::fabs(v0 - v1);
+  //  printf("smallDoubleDifference = %g\n", difference);
+  return difference < 1e-4;
+}
+
+TEST_CASE("multi-objective", "[util]") {
+  HighsLp lp;
+  lp.num_col_ = 2;
+  lp.num_row_ = 3;
+  lp.col_cost_ = {0, 0};
+  lp.col_lower_ = {0, 0};
+  lp.col_upper_ = {kHighsInf, kHighsInf};
+  lp.row_lower_ = {-kHighsInf, -kHighsInf, -kHighsInf};
+  lp.row_upper_ = {18, 8, 14};
+  lp.a_matrix_.start_ = {0, 3, 6};
+  lp.a_matrix_.index_ = {0, 1, 2, 0, 1, 2};
+  lp.a_matrix_.value_ = {3, 1, 1, 1, 1, 2};
+  Highs h;
+  h.setOptionValue("output_flag", dev_run);
+
+  for (HighsInt k = 0; k < 2; k++) {
+    // Pass 0 is continuous; pass 1 integer
+    if (dev_run)
+      printf(
+          "\n******************\nPass %d: var type is %s\n******************\n",
+          int(k), k == 0 ? "continuous" : "integer");
+    for (HighsInt l = 0; l < 2; l++) {
+      // Pass 0 is with unsigned weights and coefficients
+      double obj_mu = l == 0 ? 1 : -1;
+      if (dev_run)
+        printf(
+            "\n******************\nPass %d: objective multiplier is "
+            "%g\n******************\n",
+            int(l), obj_mu);
+
+      if (k == 0) {
+        lp.integrality_.clear();
+      } else if (k == 1) {
+        lp.integrality_ = {HighsVarType::kInteger, HighsVarType::kInteger};
+      }
+      h.passModel(lp);
+
+      h.setOptionValue("blend_multi_objectives", true);
+
+      HighsLinearObjective linear_objective;
+      std::vector<HighsLinearObjective> linear_objectives;
+      REQUIRE(h.clearLinearObjectives() == HighsStatus::kOk);
+
+      // Begin with an illegal linear objective
+      if (dev_run) printf("\nPass illegal linear objective\n");
+      linear_objective.weight = -obj_mu;
+      linear_objective.offset = -obj_mu;
+      linear_objective.coefficients = {obj_mu * 1, obj_mu * 1, obj_mu * 0};
+      linear_objective.abs_tolerance = 0.0;
+      linear_objective.rel_tolerance = 0.0;
+      REQUIRE(h.addLinearObjective(linear_objective) == HighsStatus::kError);
+      // Now legalise the linear objective so LP has nonunique optimal
+      // solutions on the line joining (2, 6) and (5, 3)
+      if (dev_run) printf("\nPass legal linear objective\n");
+      linear_objective.coefficients = {obj_mu * 1, obj_mu * 1};
+      REQUIRE(h.addLinearObjective(linear_objective) == HighsStatus::kOk);
+
+      REQUIRE(h.run() == HighsStatus::kOk);
+      h.writeSolution("", kSolutionStylePretty);
+      REQUIRE(smallDoubleDifference(h.getInfo().objective_function_value, -7));
+      // Save the linear objective for the next
+      linear_objectives.push_back(linear_objective);
+
+      // Add a second linear objective with a very small minimization
+      // weight that should push the optimal solution to (2, 6)
+      if (dev_run) printf("\nPass second linear objective\n");
+      linear_objective.weight = obj_mu * 1e-4;
+      linear_objective.offset = 0;
+      linear_objective.coefficients = {obj_mu * 1, obj_mu * 0};
+      REQUIRE(h.addLinearObjective(linear_objective) == HighsStatus::kOk);
+
+      REQUIRE(h.run() == HighsStatus::kOk);
+      h.writeSolution("", kSolutionStylePretty);
+      REQUIRE(smallDoubleDifference(h.getSolution().col_value[0], 2));
+      REQUIRE(smallDoubleDifference(h.getSolution().col_value[1], 6));
+      linear_objectives.push_back(linear_objective);
+
+      if (dev_run) printf("\nClear and pass two linear objectives\n");
+      REQUIRE(h.clearLinearObjectives() == HighsStatus::kOk);
+      REQUIRE(h.passLinearObjectives(2, linear_objectives.data()) ==
+              HighsStatus::kOk);
+      REQUIRE(h.run() == HighsStatus::kOk);
+      h.writeSolution("", kSolutionStylePretty);
+      REQUIRE(smallDoubleDifference(h.getSolution().col_value[0], 2));
+      REQUIRE(smallDoubleDifference(h.getSolution().col_value[1], 6));
+
+      // Set illegal priorities - that can be passed OK since
+      // blend_multi_objectives = true
+      if (dev_run)
+        printf(
+            "\nSetting priorities that will be illegal when using "
+            "lexicographic "
+            "optimization\n");
+      linear_objectives[0].priority = 0;
+      linear_objectives[1].priority = 0;
+      REQUIRE(h.passLinearObjectives(2, linear_objectives.data()) ==
+              HighsStatus::kOk);
+
+      // Now test lexicographic optimization
+      h.setOptionValue("blend_multi_objectives", false);
+
+      if (dev_run) printf("\nLexicographic using illegal priorities\n");
+      REQUIRE(h.run() == HighsStatus::kError);
+
+      if (dev_run)
+        printf(
+            "\nSetting priorities that are illegal now blend_multi_objectives "
+            "= "
+            "false\n");
+      REQUIRE(h.passLinearObjectives(2, linear_objectives.data()) ==
+              HighsStatus::kError);
+
+      if (dev_run)
+        printf(
+            "\nSetting legal priorities for blend_multi_objectives = false\n");
+      linear_objectives[0].priority = 10;
+      REQUIRE(h.passLinearObjectives(2, linear_objectives.data()) ==
+              HighsStatus::kOk);
+
+      if (dev_run)
+        printf("\nLexicographic using existing multi objective data\n");
+      REQUIRE(h.run() == HighsStatus::kOk);
+      h.writeSolution("", kSolutionStylePretty);
+      REQUIRE(smallDoubleDifference(h.getSolution().col_value[0], 2));
+      REQUIRE(smallDoubleDifference(h.getSolution().col_value[1], 6));
+
+      // Back to blending
+      h.setOptionValue("blend_multi_objectives", true);
+      //  h.setOptionValue("output_flag", true);
+      REQUIRE(h.clearLinearObjectives() == HighsStatus::kOk);
+      linear_objectives[0].coefficients = {obj_mu * 1.0001, obj_mu * 1};
+      linear_objectives[0].abs_tolerance = 1e-5;
+      linear_objectives[0].rel_tolerance = 0.05;
+      linear_objectives[1].weight = obj_mu * 1e-3;
+      if (dev_run)
+        printf(
+            "\nBlending: first solve objective just giving unique optimal "
+            "solution\n");
+      REQUIRE(h.passLinearObjectives(1, linear_objectives.data()) ==
+              HighsStatus::kOk);
+
+      REQUIRE(h.run() == HighsStatus::kOk);
+      h.writeSolution("", kSolutionStylePretty);
+
+      REQUIRE(h.passLinearObjectives(2, linear_objectives.data()) ==
+              HighsStatus::kOk);
+
+      REQUIRE(h.run() == HighsStatus::kOk);
+      h.writeSolution("", kSolutionStylePretty);
+
+      // Back to lexicographic optimization
+      h.setOptionValue("blend_multi_objectives", false);
+
+      if (dev_run) printf("\nLexicographic using non-trivial tolerances\n");
+      REQUIRE(h.run() == HighsStatus::kOk);
+      h.writeSolution("", kSolutionStylePretty);
+
+      if (k == 0) {
+        REQUIRE(smallDoubleDifference(h.getSolution().col_value[0], 4.9));
+        REQUIRE(smallDoubleDifference(h.getSolution().col_value[1], 3.1));
+      } else {
+        REQUIRE(smallDoubleDifference(h.getSolution().col_value[0], 5));
+        REQUIRE(smallDoubleDifference(h.getSolution().col_value[1], 3));
+      }
+
+      linear_objectives[0].abs_tolerance = kHighsInf;
+
+      REQUIRE(h.passLinearObjectives(2, linear_objectives.data()) ==
+              HighsStatus::kOk);
+
+      REQUIRE(h.run() == HighsStatus::kOk);
+      h.writeSolution("", kSolutionStylePretty);
+
+      //  printf("Solution = [%23.18g, %23.18g]\n",
+      //  h.getSolution().col_value[0], h.getSolution().col_value[1]);
+      if (k == 0) {
+        REQUIRE(smallDoubleDifference(h.getSolution().col_value[0], 1.30069));
+        REQUIRE(smallDoubleDifference(h.getSolution().col_value[1], 6.34966));
+      } else {
+        REQUIRE(smallDoubleDifference(h.getSolution().col_value[0], 2));
+        REQUIRE(smallDoubleDifference(h.getSolution().col_value[1], 6));
+      }
+    }
+  }
+}
diff --git a/docs/make.jl b/docs/make.jl
index c1f6bc70b9..93947c773c 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -63,10 +63,14 @@ Documenter.makedocs(
                 "structures/classes/HighsSparseMatrix.md",
                 "structures/classes/HighsLp.md",
                 "structures/classes/HighsHessian.md",
-                "structures/classes/HighsModel.md",
-                "structures/classes/HighsSolution.md",
-                "structures/classes/HighsBasis.md",
-                "structures/classes/HighsInfo.md",
+                "structures/classes/HighsModel.md"
+            ],
+	    "Structures" => Any[
+                "structures/structs/index.md",
+                "structures/structs/HighsSolution.md",
+                "structures/structs/HighsBasis.md",
+                "structures/structs/HighsInfo.md",
+                "structures/structs/HighsLinearObjective.md"
             ],
 	],
         "Callbacks" => "callbacks.md",
diff --git a/docs/src/guide/further.md b/docs/src/guide/further.md
index 3597928f5b..91a4dc75ca 100644
--- a/docs/src/guide/further.md
+++ b/docs/src/guide/further.md
@@ -104,3 +104,80 @@ HiGHS. Note that this does not affect how the incumbent model is
 solved. There are two corresponding [postsolve](@ref Presolve/postsolve)
 methods, according to whether there are just solution values, or also
 a basis.
+
+## [Multi-objective optimization](@id guide-multi-objective-optimization)
+
+Users can specify multiple linear objectives with respect to which
+HiGHS will optimize by either blending them, or by performing
+lexicographic optimization according to the truth of the
+[blend\_multi\_objectives](@ref blend_multi_objectives) option. Each
+linear objective is represented by the following data, held in the
+[HighsLinearObjective](@ref HighsLinearObjective) structure
+
+- weight: Scalar of type double - The weight of this objective when blending 
+- offset: Scalar of type double - The offset of this objective
+- coefficients: Vector of type double - The coefficients of this objective
+- abs\_tolerance: Scalar of type double - The absolute tolerance on this objective when performing lexicographic optimization 
+- rel\_tolerance: Scalar of type double - The relative tolerance on this objective when performing lexicographic optimization 
+- priority: Scalar of type HighsInt - The priority of this objective when performing lexicographic optimization
+
+### Methods
+
+Multi-objective optimization in HiGHS is defined by the following methods
+
+- [passLinearObjectives](@ref Multi-objective-optimization] - Pass multiple linear objectives as their number `num_linear_objective` and pointer to a vector of `HighsLinearObjective` instances, overwriting any previous linear objectives
+- [addLinearObjective](@ref Multi-objective-optimization] - Add a single `HighsLinearObjective` instance to any already stored in HiGHS
+- [clearLinearObjectives](@ref Multi-objective-optimization] - Clears any linear objectives stored in HiGHS
+
+When there is at least one `HighsLinearObjective` instance in HiGHS,
+the `col_cost_` data in the incumbent model is ignored.
+
+### Blending multiple linear objectives
+
+When [blend\_multi\_objectives](@ref blend_multi_objectives) is `true`,
+as it is by default, any `HighsLinearObjective` instances will be
+combined according to the `weight` values, and the resulting objective
+will be minimized. Hence, any objectives that should be maximized
+within the combination must have a negative `weight` value.
+
+### Lexicographic optimization of multiple linear objectives
+
+When [blend\_multi\_objectives](@ref blend_multi_objectives) is `false`,
+HiGHS will optimize lexicographically with respect to any
+`HighsLinearObjective` instances. This is carried out as follows, according to the
+`priority` values in `HighsLinearObjective` instances. Note that _all
+priority values must be distinct_.
+
+* Minimize/maximize with respect to the linear objective of highest priority value, according to whether its `weight` is positive/negative
+
+* Add a constraint to the model so that the value of the linear objective of highest priority satsifies a bound given by the values of `abs_tolerance` and/or `rel_tolerance`.
+
+    + If the objective was minimized to a value ``f^*>=0``, then the constraint ensures that the this objective value is no greater than
+```math
+
+\min(f^*+`abs_tolerance`,~f^*[1+`rel_tolerance`]).
+```
+
+    + If the objective was minimized to a value ``f^*<0``, then the constraint ensures that the this objective value is no greater than
+```math
+\min(f^*+`abs_tolerance`,~f^*[1-`rel_tolerance`]).
+```
+
+    + If the objective was maximized to a value ``f^*>=0``, then the constraint ensures that the this objective value is no less than
+```math
+\max(f^*-`abs_tolerance`,~f^*[1-`rel_tolerance`]).
+```
+
+    + If the objective was maximized to a value ``f^*<0``, then the constraint ensures that the this objective value is no less than
+```math
+\max(f^*-`abs_tolerance`,~f^*[1+`rel_tolerance`]).
+```
+
+* Minimize/maximize with respect to the linear objective of next highest priority, and then add a corresponding objective constraint to the model, repeating until optimization with respect to the linear objective of lowest priority has taken place.
+
+Note
+
+* Negative values of `abs_tolerance` and `rel_tolerance` will be ignored. This is a convenient way of "switching off" a bounding technique that is not of interest.
+* When the model is continuous, no dual information will be returned if there is more than one linear objective.
+
+
diff --git a/docs/src/interfaces/python/example-py.md b/docs/src/interfaces/python/example-py.md
index 50d40c59e7..fd4eef038b 100644
--- a/docs/src/interfaces/python/example-py.md
+++ b/docs/src/interfaces/python/example-py.md
@@ -237,5 +237,11 @@ print('Basis validity = ', h.basisValidityToString(info.basis_validity))
  * `presolveRuleTypeToString`
  * `postsolve`
  
- 
- 
+## Multi-objective optimization
+
+* `addLinearObjective`
+* `clearLinearObjectives`
+
+
+
+
diff --git a/docs/src/options/definitions.md b/docs/src/options/definitions.md
index 8055067601..f4e932523c 100644
--- a/docs/src/options/definitions.md
+++ b/docs/src/options/definitions.md
@@ -353,3 +353,8 @@
 - Range: {0, 2147483647}
 - Default: 4000
 
+## blend\_multi\_objectives
+- Blend multiple objectives or apply lexicographically
+- Type: boolean
+- Default: "true"
+
diff --git a/docs/src/structures/classes/index.md b/docs/src/structures/classes/index.md
index f2f5ed4326..9277bc5ef8 100644
--- a/docs/src/structures/classes/index.md
+++ b/docs/src/structures/classes/index.md
@@ -6,8 +6,5 @@ The data members of fundamental classes in HiGHS are defined in this section.
  * [HighsLp](@ref)
  * [HighsHessian](@ref)
  * [HighsModel](@ref)
- * [HighsSolution](@ref)
- * [HighsBasis](@ref)
- * [HighsInfo](@ref)
 
 Class data members for internal use only are not documented.
diff --git a/docs/src/structures/classes/HighsBasis.md b/docs/src/structures/structs/HighsBasis.md
similarity index 96%
rename from docs/src/structures/classes/HighsBasis.md
rename to docs/src/structures/structs/HighsBasis.md
index bfc71c7efd..6bac592bd5 100644
--- a/docs/src/structures/classes/HighsBasis.md
+++ b/docs/src/structures/structs/HighsBasis.md
@@ -1,6 +1,6 @@
 # HighsBasis
 
-The basis of a model is communicated via an instance of the HighsBasis class
+The basis of a model is communicated via an instance of the HighsBasis structure
 
 - valid: Scalar of type bool - Indicates whether the basis is valid
 - col\_status: Vector of type [HighsBasisStatus](@ref) - Comparison with [HighsBasisStatus](@ref) gives the basis status of a column
diff --git a/docs/src/structures/classes/HighsInfo.md b/docs/src/structures/structs/HighsInfo.md
similarity index 98%
rename from docs/src/structures/classes/HighsInfo.md
rename to docs/src/structures/structs/HighsInfo.md
index f0ddebee4f..9a3217adc3 100644
--- a/docs/src/structures/classes/HighsInfo.md
+++ b/docs/src/structures/structs/HighsInfo.md
@@ -1,6 +1,6 @@
 # HighsInfo
 
-Scalar information about a solved model is communicated via an instance of the HighsInfo class
+Scalar information about a solved model is communicated via an instance of the HighsInfo structure
 
 ## valid
 - Indicates whether the values in a HighsInfo instance are valid
diff --git a/docs/src/structures/structs/HighsLinearObjective.md b/docs/src/structures/structs/HighsLinearObjective.md
new file mode 100644
index 0000000000..9cc9f09d12
--- /dev/null
+++ b/docs/src/structures/structs/HighsLinearObjective.md
@@ -0,0 +1,11 @@
+# HighsLinearObjective
+
+A linear objective for a model is communicated via an instance of the HighsLinearObjective structure
+
+- weight: Scalar of type double - The weight of this objective when blending 
+- offset: Scalar of type double - The offset of this objective
+- coefficients: Vector of type double - The coefficients of this objective
+- abs\_tolerance: Scalar of type double - The absolute tolerance on this objective when performing lexicographic optimization 
+- rel\_tolerance: Scalar of type double - The relative tolerance on this objective when performing lexicographic optimization 
+- priority: Scalar of type HighsInt - The priority of this objective when performing lexicographic optimization
+
diff --git a/docs/src/structures/classes/HighsSolution.md b/docs/src/structures/structs/HighsSolution.md
similarity index 96%
rename from docs/src/structures/classes/HighsSolution.md
rename to docs/src/structures/structs/HighsSolution.md
index 78fc48a959..5b9dbf6a73 100644
--- a/docs/src/structures/classes/HighsSolution.md
+++ b/docs/src/structures/structs/HighsSolution.md
@@ -1,6 +1,6 @@
 # HighsSolution
 
-The solution of a model is communicated via an instance of the HighsSolution class
+The solution of a model is communicated via an instance of the HighsSolution structure
 
 - value\_valid: Scalar of type bool - Indicates whether the column and row values are valid
 - dual\_valid: Scalar of type bool - Indicates whether the column and row [duals](@ref Dual-values) are valid
diff --git a/docs/src/structures/structs/index.md b/docs/src/structures/structs/index.md
new file mode 100644
index 0000000000..d8c763c2b3
--- /dev/null
+++ b/docs/src/structures/structs/index.md
@@ -0,0 +1,10 @@
+# [Overview](@id structs-overview)
+
+The data members of fundamental structs in HiGHS are defined in this section.
+
+ * [HighsSolution](@ref)
+ * [HighsBasis](@ref)
+ * [HighsInfo](@ref)
+ * [HighsLinearObjective](@ref)
+
+Structure data members for internal use only are not documented.
diff --git a/examples/call_highs_from_python.py b/examples/call_highs_from_python.py
index 2ff88f281b..61e6356f3e 100644
--- a/examples/call_highs_from_python.py
+++ b/examples/call_highs_from_python.py
@@ -493,4 +493,8 @@ def user_interrupt_callback(
 print("row_bound:", iis.row_bound)
 
 print("col_index:", iis.col_index)
-print("col_bound:", iis.col_bound)
\ No newline at end of file
+print("col_bound:", iis.col_bound)
+
+# ~~~
+# Clear so that incumbent model is empty
+h.clear()
diff --git a/src/Highs.h b/src/Highs.h
index 34d1044f5f..04a4a3ce25 100644
--- a/src/Highs.h
+++ b/src/Highs.h
@@ -150,6 +150,24 @@ class Highs {
   HighsStatus passHessian(const HighsInt dim, const HighsInt num_nz,
                           const HighsInt format, const HighsInt* start,
                           const HighsInt* index, const double* value);
+  /**
+   * @brief Pass multiple linear objectives for the incumbent model
+   */
+  HighsStatus passLinearObjectives(
+      const HighsInt num_linear_objective,
+      const HighsLinearObjective* linear_objective);
+
+  /**
+   * @brief Add a linear objective for the incumbent model
+   */
+  HighsStatus addLinearObjective(const HighsLinearObjective& linear_objective,
+                                 const HighsInt iObj = -1);
+
+  /**
+   * @brief Clear the multiple linear objective data
+   */
+  HighsStatus clearLinearObjectives();
+
   /**
    * @brief Pass a column name to the incumbent model
    */
@@ -183,7 +201,7 @@ class Highs {
   HighsStatus presolve();
 
   /**
-   * @brief Solve the incumbent model according to the specified options
+   * @brief Run the solver, accounting for any multiple objective
    */
   HighsStatus run();
 
@@ -1389,6 +1407,8 @@ class Highs {
   ICrashInfo icrash_info_;
 
   HighsModel model_;
+  std::vector<HighsLinearObjective> multi_linear_objective_;
+
   HighsModel presolved_model_;
   HighsTimer timer_;
 
@@ -1397,7 +1417,6 @@ class Highs {
   HighsInfo info_;
   HighsRanging ranging_;
   HighsIis iis_;
-
   std::vector<HighsObjectiveSolution> saved_objective_and_solution_;
 
   HighsPresolveStatus model_presolve_status_ =
@@ -1424,6 +1443,8 @@ class Highs {
 
   bool written_log_header = false;
 
+  HighsStatus solve();
+
   void exactResizeModel() {
     this->model_.lp_.exactResize();
     this->model_.hessian_.exactResize();
@@ -1600,6 +1621,10 @@ class Highs {
                          HighsInt* iis_col_index, HighsInt* iis_row_index,
                          HighsInt* iis_col_bound, HighsInt* iis_row_bound);
 
+  HighsStatus returnFromLexicographicOptimization(
+      const HighsStatus return_status, HighsInt original_lp_num_row);
+  HighsStatus multiobjectiveSolve();
+
   bool aFormatOk(const HighsInt num_nz, const HighsInt format);
   bool qFormatOk(const HighsInt num_nz, const HighsInt format);
   void clearZeroHessian();
@@ -1623,6 +1648,10 @@ class Highs {
                               const bool constraint,
                               const double ill_conditioning_bound);
   bool infeasibleBoundsOk();
+  bool validLinearObjective(const HighsLinearObjective& linear_objective,
+                            const HighsInt iObj) const;
+  bool hasRepeatedLinearObjectivePriorities(
+      const HighsLinearObjective* linear_objective = nullptr) const;
 };
 
 // Start of deprecated methods not in the Highs class
diff --git a/src/highs_bindings.cpp b/src/highs_bindings.cpp
index b8f7e249cf..64db9172d3 100644
--- a/src/highs_bindings.cpp
+++ b/src/highs_bindings.cpp
@@ -159,6 +159,10 @@ HighsStatus highs_passHessianPointers(Highs* h, const HighsInt dim,
                         q_value_ptr);
 }
 
+HighsStatus highs_addLinearObjective(Highs* h, const HighsLinearObjective& linear_objective) {
+  return h->addLinearObjective(linear_objective, -1);
+}
+
 HighsStatus highs_postsolve(Highs* h, const HighsSolution& solution,
                             const HighsBasis& basis) {
   return h->postsolve(solution, basis);
@@ -939,6 +943,8 @@ PYBIND11_MODULE(_core, m) {
       .def("passModel", &highs_passLpPointers)
       .def("passHessian", &highs_passHessian)
       .def("passHessian", &highs_passHessianPointers)
+      .def("addLinearObjective", &highs_addLinearObjective)
+      .def("clearLinearObjectives", &Highs::clearLinearObjectives)
       .def("passColName", &Highs::passColName)
       .def("passRowName", &Highs::passRowName)
       .def("readModel", &Highs::readModel)
@@ -1148,6 +1154,14 @@ PYBIND11_MODULE(_core, m) {
       .def(py::init<>())
       .def_readwrite("simplex_time", &HighsIisInfo::simplex_time)
       .def_readwrite("simplex_iterations", &HighsIisInfo::simplex_iterations);
+  py::class_<HighsLinearObjective>(m, "HighsLinearObjective")
+      .def(py::init<>())
+      .def_readwrite("weight", &HighsLinearObjective::weight)
+      .def_readwrite("offset", &HighsLinearObjective::offset)
+      .def_readwrite("coefficients", &HighsLinearObjective::coefficients)
+      .def_readwrite("abs_tolerance", &HighsLinearObjective::abs_tolerance)
+      .def_readwrite("rel_tolerance", &HighsLinearObjective::rel_tolerance)
+      .def_readwrite("priority", &HighsLinearObjective::priority);
   // constants
   m.attr("kHighsInf") = kHighsInf;
   m.attr("kHighsIInf") = kHighsIInf;
diff --git a/src/highspy/__init__.py b/src/highspy/__init__.py
index 9e3c03a0e1..e3044e8ccb 100644
--- a/src/highspy/__init__.py
+++ b/src/highspy/__init__.py
@@ -29,6 +29,7 @@
     HighsRanging,
     kHighsInf,
     kHighsIInf,
+    HighsLinearObjective,
     HIGHS_VERSION_MAJOR,
     HIGHS_VERSION_MINOR,
     HIGHS_VERSION_PATCH,
diff --git a/src/interfaces/highs_c_api.cpp b/src/interfaces/highs_c_api.cpp
index 8dad0e5216..655c0fce30 100644
--- a/src/interfaces/highs_c_api.cpp
+++ b/src/interfaces/highs_c_api.cpp
@@ -291,6 +291,57 @@ HighsInt Highs_passHessian(void* highs, const HighsInt dim,
       ->passHessian(dim, num_nz, format, start, index, value);
 }
 
+HighsInt Highs_passLinearObjectives(const void* highs,
+                                    const HighsInt num_linear_objective,
+                                    const double* weight, const double* offset,
+                                    const double* coefficients,
+                                    const double* abs_tolerance,
+                                    const double* rel_tolerance,
+                                    const HighsInt* priority) {
+  HighsInt status = Highs_clearLinearObjectives(highs);
+  if (status != kHighsStatusOk) return status;
+  HighsInt num_col = Highs_getNumCol(highs);
+  for (HighsInt iObj = 0; iObj < num_linear_objective; iObj++) {
+    status = Highs_addLinearObjectiveWithIndex(
+        highs, weight[iObj], offset[iObj], &coefficients[iObj * num_col],
+        abs_tolerance[iObj], rel_tolerance[iObj], priority[iObj], iObj);
+    if (status != kHighsStatusOk) return status;
+  }
+  return kHighsStatusOk;
+}
+
+HighsInt Highs_addLinearObjective(const void* highs, const double weight,
+                                  const double offset,
+                                  const double* coefficients,
+                                  const double abs_tolerance,
+                                  const double rel_tolerance,
+                                  const HighsInt priority) {
+  return Highs_addLinearObjectiveWithIndex(highs, weight, offset, coefficients,
+                                           abs_tolerance, rel_tolerance,
+                                           priority, -1);
+}
+
+HighsInt Highs_addLinearObjectiveWithIndex(
+    const void* highs, const double weight, const double offset,
+    const double* coefficients, const double abs_tolerance,
+    const double rel_tolerance, const HighsInt priority, const HighsInt iObj) {
+  HighsLinearObjective linear_objective;
+  HighsInt num_col = Highs_getNumCol(highs);
+  linear_objective.weight = weight;
+  linear_objective.offset = offset;
+  for (HighsInt iCol = 0; iCol < num_col; iCol++)
+    linear_objective.coefficients.push_back(coefficients[iCol]);
+  linear_objective.abs_tolerance = abs_tolerance;
+  linear_objective.rel_tolerance = rel_tolerance;
+  linear_objective.priority = priority;
+  linear_objective.weight = weight;
+  return HighsInt(((Highs*)highs)->addLinearObjective(linear_objective));
+}
+
+HighsInt Highs_clearLinearObjectives(const void* highs) {
+  return HighsInt(((Highs*)highs)->clearLinearObjectives());
+}
+
 HighsInt Highs_passRowName(const void* highs, const HighsInt row,
                            const char* name) {
   return (HighsInt)((Highs*)highs)->passRowName(row, std::string(name));
diff --git a/src/interfaces/highs_c_api.h b/src/interfaces/highs_c_api.h
index a59ec64a17..db8764fae6 100644
--- a/src/interfaces/highs_c_api.h
+++ b/src/interfaces/highs_c_api.h
@@ -557,6 +557,27 @@ HighsInt Highs_passHessian(void* highs, const HighsInt dim,
                            const HighsInt* start, const HighsInt* index,
                            const double* value);
 
+HighsInt Highs_passLinearObjectives(const void* highs,
+                                    const HighsInt num_linear_objective,
+                                    const double* weight, const double* offset,
+                                    const double* coefficients,
+                                    const double* abs_tolerance,
+                                    const double* rel_tolerance,
+                                    const HighsInt* priority);
+
+HighsInt Highs_addLinearObjective(const void* highs, const double weight,
+                                  const double offset,
+                                  const double* coefficients,
+                                  const double abs_tolerance,
+                                  const double rel_tolerance,
+                                  const HighsInt priority);
+
+HighsInt Highs_addLinearObjectiveWithIndex(
+    const void* highs, const double weight, const double offset,
+    const double* coefficients, const double abs_tolerance,
+    const double rel_tolerance, const HighsInt priority, const HighsInt iObj);
+
+HighsInt Highs_clearLinearObjectives(const void* highs);
 /**
  * Pass the name of a row.
  *
@@ -945,7 +966,7 @@ HighsInt Highs_getDualRay(const void* highs, HighsInt* has_dual_ray,
  *                                                filled with the unboundedness
  *                                                direction.
  */
-HighsInt getDualUnboundednessDirection(
+HighsInt Highs_getDualUnboundednessDirection(
     const void* highs, HighsInt* has_dual_unboundedness_direction,
     double* dual_unboundedness_direction_value);
 
diff --git a/src/lp_data/HStruct.h b/src/lp_data/HStruct.h
index e54657d406..3963d9df1c 100644
--- a/src/lp_data/HStruct.h
+++ b/src/lp_data/HStruct.h
@@ -144,4 +144,14 @@ struct HighsIllConditioning {
   void clear();
 };
 
+struct HighsLinearObjective {
+  double weight;
+  double offset;
+  std::vector<double> coefficients;
+  double abs_tolerance;
+  double rel_tolerance;
+  HighsInt priority;
+  void clear();
+};
+
 #endif /* LP_DATA_HSTRUCT_H_ */
diff --git a/src/lp_data/Highs.cpp b/src/lp_data/Highs.cpp
index 67c6d44e33..6352a68a09 100644
--- a/src/lp_data/Highs.cpp
+++ b/src/lp_data/Highs.cpp
@@ -59,6 +59,7 @@ HighsStatus Highs::clear() {
 
 HighsStatus Highs::clearModel() {
   model_.clear();
+  multi_linear_objective_.clear();
   return clearSolver();
 }
 
@@ -577,6 +578,36 @@ HighsStatus Highs::passHessian(const HighsInt dim, const HighsInt num_nz,
   return passHessian(hessian);
 }
 
+HighsStatus Highs::passLinearObjectives(
+    const HighsInt num_linear_objective,
+    const HighsLinearObjective* linear_objective) {
+  if (num_linear_objective < 0) return HighsStatus::kOk;
+  this->multi_linear_objective_.clear();
+  for (HighsInt iObj = 0; iObj < num_linear_objective; iObj++)
+    if (this->addLinearObjective(linear_objective[iObj], iObj) !=
+        HighsStatus::kOk)
+      return HighsStatus::kError;
+  return HighsStatus::kOk;
+}
+
+HighsStatus Highs::addLinearObjective(
+    const HighsLinearObjective& linear_objective, const HighsInt iObj) {
+  if (model_.isQp()) {
+    highsLogUser(options_.log_options, HighsLogType::kError,
+                 "Cannot define additional linear objective for QP\n");
+    return HighsStatus::kError;
+  }
+  if (!this->validLinearObjective(linear_objective, iObj))
+    return HighsStatus::kError;
+  this->multi_linear_objective_.push_back(linear_objective);
+  return HighsStatus::kOk;
+}
+
+HighsStatus Highs::clearLinearObjectives() {
+  this->multi_linear_objective_.clear();
+  return HighsStatus::kOk;
+}
+
 HighsStatus Highs::passColName(const HighsInt col, const std::string& name) {
   const HighsInt num_col = this->model_.lp_.num_col_;
   if (col < 0 || col >= num_col) {
@@ -876,9 +907,15 @@ HighsStatus Highs::presolve() {
   return returnFromHighs(return_status);
 }
 
+HighsStatus Highs::run() {
+  HighsInt num_linear_objective = this->multi_linear_objective_.size();
+  if (num_linear_objective == 0) return this->solve();
+  return this->multiobjectiveSolve();
+}
+
 // Checks the options calls presolve and postsolve if needed. Solvers are called
 // with callSolveLp(..)
-HighsStatus Highs::run() {
+HighsStatus Highs::solve() {
   HighsInt min_highs_debug_level = kHighsDebugLevelMin;
   // kHighsDebugLevelCostly;
   // kHighsDebugLevelMax;
diff --git a/src/lp_data/HighsInfo.cpp b/src/lp_data/HighsInfo.cpp
index ce01a585b8..2792f64550 100644
--- a/src/lp_data/HighsInfo.cpp
+++ b/src/lp_data/HighsInfo.cpp
@@ -292,7 +292,7 @@ void reportInfo(FILE* file, const InfoRecordInt64& info,
     fprintf(file, "\n# %s\n# [type: int64_t]\n%s = %" PRId64 "\n",
             info.description.c_str(), info.name.c_str(), *info.value);
   } else {
-    fprintf(file, "%s = %" PRId64 "\n", info.name.c_str(), *info.value);
+    fprintf(file, "%-30s = %" PRId64 "\n", info.name.c_str(), *info.value);
   }
 }
 
@@ -306,7 +306,7 @@ void reportInfo(FILE* file, const InfoRecordInt& info,
     fprintf(file, "\n# %s\n# [type: HighsInt]\n%s = %" HIGHSINT_FORMAT "\n",
             info.description.c_str(), info.name.c_str(), *info.value);
   } else {
-    fprintf(file, "%s = %" HIGHSINT_FORMAT "\n", info.name.c_str(),
+    fprintf(file, "%-30s = %" HIGHSINT_FORMAT "\n", info.name.c_str(),
             *info.value);
   }
 }
@@ -321,6 +321,6 @@ void reportInfo(FILE* file, const InfoRecordDouble& info,
     fprintf(file, "\n# %s\n# [type: double]\n%s = %g\n",
             info.description.c_str(), info.name.c_str(), *info.value);
   } else {
-    fprintf(file, "%s = %g\n", info.name.c_str(), *info.value);
+    fprintf(file, "%-30s = %g\n", info.name.c_str(), *info.value);
   }
 }
diff --git a/src/lp_data/HighsInterface.cpp b/src/lp_data/HighsInterface.cpp
index a40813d1fe..1c7ad16b9f 100644
--- a/src/lp_data/HighsInterface.cpp
+++ b/src/lp_data/HighsInterface.cpp
@@ -1755,7 +1755,6 @@ HighsStatus Highs::getDualRayInterface(bool& has_dual_ray,
       this->setOptionValue("presolve", kHighsOffString);
       this->setOptionValue("solve_relaxation", true);
       HighsStatus call_status = this->run();
-      this->writeSolution("", kSolutionStylePretty);
       if (call_status != HighsStatus::kOk) return_status = call_status;
       has_dual_ray = ekk_instance_.status_.has_dual_ray;
       has_invert = ekk_instance_.status_.has_invert;
@@ -3584,3 +3583,356 @@ bool Highs::infeasibleBoundsOk() {
                  int(num_true_infeasible_bound));
   return num_true_infeasible_bound == 0;
 }
+
+bool Highs::validLinearObjective(const HighsLinearObjective& linear_objective,
+                                 const HighsInt iObj) const {
+  HighsInt linear_objective_coefficients_size =
+      linear_objective.coefficients.size();
+  if (linear_objective_coefficients_size != this->model_.lp_.num_col_) {
+    highsLogUser(
+        options_.log_options, HighsLogType::kError,
+        "Coefficient vector for linear objective %s has size %d != %d = "
+        "lp.num_col_\n",
+        iObj >= 0 ? std::to_string(iObj).c_str() : "",
+        int(linear_objective_coefficients_size),
+        int(this->model_.lp_.num_col_));
+    return false;
+  }
+  if (!options_.blend_multi_objectives &&
+      hasRepeatedLinearObjectivePriorities(&linear_objective)) {
+    highsLogUser(
+        options_.log_options, HighsLogType::kError,
+        "Repeated priorities for lexicographic optimization is illegal\n");
+    return false;
+  }
+  return true;
+}
+
+bool Highs::hasRepeatedLinearObjectivePriorities(
+    const HighsLinearObjective* linear_objective) const {
+  // Look for repeated values in the linear objective priorities, also
+  // comparing linear_objective if it's not a null pointer. Cost is
+  // O(n^2), but who will have more than O(1) linear objectives!
+  HighsInt num_linear_objective = this->multi_linear_objective_.size();
+  if (num_linear_objective <= 0 ||
+      num_linear_objective <= 1 && !linear_objective)
+    return false;
+  for (HighsInt iObj0 = 0; iObj0 < num_linear_objective; iObj0++) {
+    HighsInt priority0 = this->multi_linear_objective_[iObj0].priority;
+    for (HighsInt iObj1 = iObj0 + 1; iObj1 < num_linear_objective; iObj1++) {
+      HighsInt priority1 = this->multi_linear_objective_[iObj1].priority;
+      if (priority1 == priority0) return true;
+    }
+    if (linear_objective) {
+      if (linear_objective->priority == priority0) return true;
+    }
+  }
+  return false;
+}
+
+bool comparison(std::pair<HighsInt, HighsInt> x1,
+                std::pair<HighsInt, HighsInt> x2) {
+  return x1.first >= x2.first;
+}
+
+HighsStatus Highs::returnFromLexicographicOptimization(
+    HighsStatus return_status, HighsInt original_lp_num_row) {
+  // Save model_status_ and info_ since they are cleared by calling
+  // deleteRows
+  HighsModelStatus model_status = this->model_status_;
+  HighsInfo info = this->info_;
+  HighsInt num_linear_objective = this->multi_linear_objective_.size();
+  if (num_linear_objective > 1) {
+    this->deleteRows(original_lp_num_row, this->model_.lp_.num_row_ - 1);
+    // Recover model_status_ and info_, and then account for lack of basis or
+    // dual solution
+    this->model_status_ = model_status;
+    this->info_ = info;
+    info_.objective_function_value = 0;
+    info_.basis_validity = kBasisValidityInvalid;
+    info_.dual_solution_status = kSolutionStatusNone;
+    info_.num_dual_infeasibilities = kHighsIllegalInfeasibilityCount;
+    info_.max_dual_infeasibility = kHighsIllegalInfeasibilityMeasure;
+    info_.sum_dual_infeasibilities = kHighsIllegalInfeasibilityMeasure;
+    info_.max_complementarity_violation = kHighsIllegalComplementarityViolation;
+    info_.sum_complementarity_violations =
+        kHighsIllegalComplementarityViolation;
+    this->solution_.value_valid = true;
+    this->model_.lp_.col_cost_.assign(this->model_.lp_.num_col_, 0);
+  }
+  return return_status;
+}
+
+HighsStatus Highs::multiobjectiveSolve() {
+  const HighsInt coeff_logging_size_limit = 10;
+  HighsInt num_linear_objective = this->multi_linear_objective_.size();
+
+  assert(num_linear_objective > 0);
+  HighsLp& lp = this->model_.lp_;
+  for (HighsInt iObj = 0; iObj < num_linear_objective; iObj++) {
+    HighsLinearObjective& multi_linear_objective =
+        this->multi_linear_objective_[iObj];
+    if (multi_linear_objective.coefficients.size() !=
+        static_cast<size_t>(lp.num_col_)) {
+      highsLogUser(options_.log_options, HighsLogType::kError,
+                   "Multiple linear objective coefficient vector %d has size "
+                   "incompatible with model\n",
+                   int(iObj));
+      return HighsStatus::kError;
+    }
+  }
+
+  std::unique_ptr<std::stringstream> multi_objective_log;
+  highsLogUser(options_.log_options, HighsLogType::kInfo,
+               "Solving with %d multiple linear objectives, %s\n",
+               int(num_linear_objective),
+               this->options_.blend_multi_objectives
+                   ? "blending objectives by weight"
+                   : "using lexicographic optimization by priority");
+  highsLogUser(
+      options_.log_options, HighsLogType::kInfo,
+      "Ix      weight      offset     abs_tol     rel_tol    priority%s\n",
+      lp.num_col_ < coeff_logging_size_limit ? "   coefficients" : "");
+  for (HighsInt iObj = 0; iObj < num_linear_objective; iObj++) {
+    HighsLinearObjective& linear_objective =
+        this->multi_linear_objective_[iObj];
+    multi_objective_log =
+        std::unique_ptr<std::stringstream>(new std::stringstream());
+    *multi_objective_log << highsFormatToString(
+        "%2d %11.6g %11.6g %11.6g %11.6g %11d  ", int(iObj),
+        linear_objective.weight, linear_objective.offset,
+        linear_objective.abs_tolerance, linear_objective.rel_tolerance,
+        linear_objective.priority);
+    if (lp.num_col_ < coeff_logging_size_limit) {
+      for (HighsInt iCol = 0; iCol < lp.num_col_; iCol++)
+        *multi_objective_log << highsFormatToString(
+            "%s c_{%1d} = %g", iCol == 0 ? "" : ",", int(iCol),
+            linear_objective.coefficients[iCol]);
+    }
+    *multi_objective_log << "\n";
+    highsLogUser(options_.log_options, HighsLogType::kInfo, "%s",
+                 multi_objective_log->str().c_str());
+  }
+  this->clearSolver();
+  if (this->options_.blend_multi_objectives) {
+    // Objectives are blended by weight and minimized
+    lp.offset_ = 0;
+    lp.col_cost_.assign(lp.num_col_, 0);
+    for (HighsInt iObj = 0; iObj < num_linear_objective; iObj++) {
+      HighsLinearObjective& multi_linear_objective =
+          this->multi_linear_objective_[iObj];
+      lp.offset_ +=
+          multi_linear_objective.weight * multi_linear_objective.offset;
+      for (HighsInt iCol = 0; iCol < lp.num_col_; iCol++)
+        lp.col_cost_[iCol] += multi_linear_objective.weight *
+                              multi_linear_objective.coefficients[iCol];
+    }
+    lp.sense_ = ObjSense::kMinimize;
+
+    multi_objective_log =
+        std::unique_ptr<std::stringstream>(new std::stringstream());
+    *multi_objective_log << highsFormatToString(
+        "Solving with blended objective");
+    if (lp.num_col_ < coeff_logging_size_limit) {
+      *multi_objective_log << highsFormatToString(
+          ": %s %g", lp.sense_ == ObjSense::kMinimize ? "min" : "max",
+          lp.offset_);
+      for (HighsInt iCol = 0; iCol < lp.num_col_; iCol++)
+        *multi_objective_log << highsFormatToString(
+            " + (%g) x[%d]", lp.col_cost_[iCol], int(iCol));
+    }
+    *multi_objective_log << "\n";
+    highsLogUser(options_.log_options, HighsLogType::kInfo, "%s",
+                 multi_objective_log->str().c_str());
+    return this->solve();
+  }
+
+  // Objectives are applied lexicographically
+  if (model_.isQp() && num_linear_objective > 1) {
+    // Lexicographic optimization with a single linear objective is
+    // trivially standard optimization, so is OK
+    highsLogUser(
+        options_.log_options, HighsLogType::kError,
+        "Cannot perform non-trivial lexicographic optimization for QP\n");
+    return HighsStatus::kError;
+  }
+  // Check whether there are repeated linear objective priorities
+  if (hasRepeatedLinearObjectivePriorities()) {
+    highsLogUser(
+        options_.log_options, HighsLogType::kError,
+        "Repeated priorities for lexicographic optimization is illegal\n");
+    return HighsStatus::kError;
+  }
+  std::vector<std::pair<HighsInt, HighsInt>> priority_objective;
+
+  for (HighsInt iObj = 0; iObj < num_linear_objective; iObj++)
+    priority_objective.push_back(
+        std::make_pair(this->multi_linear_objective_[iObj].priority, iObj));
+  std::sort(priority_objective.begin(), priority_objective.end(), comparison);
+  // Clear LP objective
+  lp.offset_ = 0;
+  lp.col_cost_.assign(lp.num_col_, 0);
+  const HighsInt original_lp_num_row = lp.num_row_;
+  std::vector<HighsInt> index(lp.num_col_);
+  std::vector<double> value(lp.num_col_);
+  // Use the solution of one MIP to provide an integer feasible
+  // solution of the next
+  HighsSolution solution;
+  for (HighsInt iIx = 0; iIx < num_linear_objective; iIx++) {
+    HighsInt priority = priority_objective[iIx].first;
+    HighsInt iObj = priority_objective[iIx].second;
+    // Use this objective
+    HighsLinearObjective& linear_objective =
+        this->multi_linear_objective_[iObj];
+    lp.offset_ = linear_objective.offset;
+    lp.col_cost_ = linear_objective.coefficients;
+    lp.sense_ =
+        linear_objective.weight > 0 ? ObjSense::kMinimize : ObjSense::kMaximize;
+    if (lp.isMip() && solution.value_valid) {
+      HighsStatus set_solution_status = this->setSolution(solution);
+      if (set_solution_status == HighsStatus::kError) {
+        highsLogUser(options_.log_options, HighsLogType::kError,
+                     "Failure to use one MIP to provide an integer feasible "
+                     "solution of the next\n");
+        return returnFromLexicographicOptimization(HighsStatus::kError,
+                                                   original_lp_num_row);
+      }
+      bool valid, integral, feasible;
+      HighsStatus assess_primal_solution =
+          assessPrimalSolution(valid, integral, feasible);
+      if (!valid || !integral || !feasible) {
+        highsLogUser(options_.log_options, HighsLogType::kWarning,
+                     "Failure to use one MIP to provide an integer feasible "
+                     "solution of the next: "
+                     "status is valid = %s, integral = %s, feasible = %s\n",
+                     highsBoolToString(valid).c_str(),
+                     highsBoolToString(integral).c_str(),
+                     highsBoolToString(feasible).c_str());
+      }
+    }
+    multi_objective_log =
+        std::unique_ptr<std::stringstream>(new std::stringstream());
+    *multi_objective_log << highsFormatToString("Solving with objective %d",
+                                                int(iObj));
+    if (lp.num_col_ < coeff_logging_size_limit) {
+      *multi_objective_log << highsFormatToString(
+          ": %s %g", lp.sense_ == ObjSense::kMinimize ? "min" : "max",
+          lp.offset_);
+      for (HighsInt iCol = 0; iCol < lp.num_col_; iCol++)
+        *multi_objective_log << highsFormatToString(
+            " + (%g) x[%d]", lp.col_cost_[iCol], int(iCol));
+    }
+    *multi_objective_log << "\n";
+    highsLogUser(options_.log_options, HighsLogType::kInfo, "%s",
+                 multi_objective_log->str().c_str());
+    HighsStatus solve_status = this->solve();
+    if (solve_status == HighsStatus::kError)
+      return returnFromLexicographicOptimization(HighsStatus::kError,
+                                                 original_lp_num_row);
+    if (model_status_ != HighsModelStatus::kOptimal) {
+      highsLogUser(options_.log_options, HighsLogType::kWarning,
+                   "After priority %d solve, model status is %s\n",
+                   int(priority), modelStatusToString(model_status_).c_str());
+      return returnFromLexicographicOptimization(HighsStatus::kWarning,
+                                                 original_lp_num_row);
+    }
+    if (iIx == num_linear_objective - 1) break;
+    if (lp.isMip()) {
+      // Save the solution to provide an integer feasible solution of
+      // the next MIP
+      solution.col_value = this->solution_.col_value;
+      solution.value_valid = true;
+    }
+    // Add the constraint
+    HighsInt nnz = 0;
+    for (HighsInt iCol = 0; iCol < lp.num_col_; iCol++) {
+      if (lp.col_cost_[iCol]) {
+        index[nnz] = iCol;
+        value[nnz] = lp.col_cost_[iCol];
+        nnz++;
+      }
+    }
+    double objective = info_.objective_function_value;
+    HighsStatus add_row_status;
+    double lower_bound = -kHighsInf;
+    double upper_bound = kHighsInf;
+    if (lp.sense_ == ObjSense::kMinimize) {
+      // Minimizing, so set a greater upper bound than the objective
+      if (linear_objective.abs_tolerance >= 0)
+        upper_bound = objective + linear_objective.abs_tolerance;
+      if (linear_objective.rel_tolerance >= 0) {
+        if (objective >= 0) {
+          // Guarantees objective of at least (1+t).f^*
+          //
+          // so ((1+t).f^*-f^*)/f^* = t
+          upper_bound = std::min(
+              objective * (1.0 + linear_objective.rel_tolerance), upper_bound);
+        } else if (objective < 0) {
+          // Guarantees objective of at least (1-t).f^*
+          //
+          // so ((1-t).f^*-f^*)/f^* = -t
+          upper_bound = std::min(
+              objective * (1.0 - linear_objective.rel_tolerance), upper_bound);
+        }
+      }
+      upper_bound -= lp.offset_;
+    } else {
+      // Maximizing, so set a lesser lower bound than the objective
+      if (linear_objective.abs_tolerance >= 0)
+        lower_bound = objective - linear_objective.abs_tolerance;
+      if (linear_objective.rel_tolerance >= 0) {
+        if (objective >= 0) {
+          // Guarantees objective of at most (1-t).f^*
+          //
+          // so ((1-t).f^*-f^*)/f^* = -t
+          lower_bound = std::max(
+              objective * (1.0 - linear_objective.rel_tolerance), lower_bound);
+        } else if (objective < 0) {
+          // Guarantees objective of at least (1+t).f^*
+          //
+          // so ((1+t).f^*-f^*)/f^* = t
+          lower_bound = std::max(
+              objective * (1.0 + linear_objective.rel_tolerance), lower_bound);
+        }
+      }
+      lower_bound -= lp.offset_;
+    }
+    if (lower_bound == -kHighsInf && upper_bound == kHighsInf)
+      highsLogUser(options_.log_options, HighsLogType::kWarning,
+                   "After priority %d solve, no objective constraint due to "
+                   "absolute tolerance being %g < 0,"
+                   " and relative tolerance being %g < 0\n",
+                   int(priority), linear_objective.abs_tolerance,
+                   linear_objective.rel_tolerance);
+    multi_objective_log =
+        std::unique_ptr<std::stringstream>(new std::stringstream());
+    *multi_objective_log << highsFormatToString(
+        "Add constraint for objective %d: ", int(iObj));
+    if (nnz < coeff_logging_size_limit) {
+      *multi_objective_log << highsFormatToString("%g <= ", lower_bound);
+      for (HighsInt iEl = 0; iEl < nnz; iEl++)
+        *multi_objective_log << highsFormatToString(
+            "%s(%g) x[%d]", iEl > 0 ? " + " : "", value[iEl], int(index[iEl]));
+      *multi_objective_log << highsFormatToString(" <= %g\n", upper_bound);
+    } else {
+      *multi_objective_log << highsFormatToString("Bounds [%g, %g]\n",
+                                                  lower_bound, upper_bound);
+    }
+    highsLogUser(options_.log_options, HighsLogType::kInfo, "%s",
+                 multi_objective_log->str().c_str());
+    add_row_status =
+        this->addRow(lower_bound, upper_bound, nnz, index.data(), value.data());
+    assert(add_row_status == HighsStatus::kOk);
+  }
+  return returnFromLexicographicOptimization(HighsStatus::kOk,
+                                             original_lp_num_row);
+}
+
+void HighsLinearObjective::clear() {
+  this->weight = 0.0;
+  this->offset = 0.0;
+  this->coefficients.clear();
+  this->abs_tolerance = 0.0;
+  this->rel_tolerance = 0.0;
+  this->priority = 0;
+}
diff --git a/src/lp_data/HighsOptions.h b/src/lp_data/HighsOptions.h
index cfb395a3fc..82424a0adc 100644
--- a/src/lp_data/HighsOptions.h
+++ b/src/lp_data/HighsOptions.h
@@ -348,6 +348,9 @@ struct HighsOptionsStruct {
   // Options for IIS calculation
   HighsInt iis_strategy;
 
+  // Option for multi-objective optimization
+  bool blend_multi_objectives;
+
   // Advanced options
   HighsInt log_dev_level;
   bool log_githash;
@@ -485,6 +488,8 @@ struct HighsOptionsStruct {
         pdlp_d_gap_tol(0.0),
         qp_iteration_limit(0),
         qp_nullspace_limit(0),
+        iis_strategy(0),
+        blend_multi_objectives(false),
         log_dev_level(0),
         log_githash(false),
         solve_relaxation(false),
@@ -1120,6 +1125,12 @@ class HighsOptions : public HighsOptionsStruct {
         kIisStrategyMax);
     records.push_back(record_int);
 
+    record_bool = new OptionRecordBool(
+        "blend_multi_objectives",
+        "Blend multiple objectives or apply lexicographically: Default = true",
+        advanced, &blend_multi_objectives, true);
+    records.push_back(record_bool);
+
     // Fix the number of user settable options
     num_user_settable_options_ = static_cast<HighsInt>(records.size());