From 39f4392aa8e4cc3b76067ce98ec3dd1c25dcc5bd Mon Sep 17 00:00:00 2001 From: Arnaud Baguet Date: Thu, 20 Mar 2025 19:57:56 -0400 Subject: [PATCH 01/21] add qp support for highs --- pyomo/contrib/appsi/solvers/highs.py | 2 +- .../solvers/tests/test_highs_persistent.py | 91 +++++++++++-------- 2 files changed, 55 insertions(+), 38 deletions(-) diff --git a/pyomo/contrib/appsi/solvers/highs.py b/pyomo/contrib/appsi/solvers/highs.py index edaafa5d1d4..83a45d68be3 100644 --- a/pyomo/contrib/appsi/solvers/highs.py +++ b/pyomo/contrib/appsi/solvers/highs.py @@ -616,7 +616,7 @@ def _set_objective(self, obj): ) repn = generate_standard_repn( - obj.expr, quadratic=False, compute_values=False + obj.expr, quadratic=True, compute_values=False ) if repn.nonlinear_expr is not None: raise DegreeError( diff --git a/pyomo/contrib/appsi/solvers/tests/test_highs_persistent.py b/pyomo/contrib/appsi/solvers/tests/test_highs_persistent.py index 2df8f9da390..c5abd21e801 100644 --- a/pyomo/contrib/appsi/solvers/tests/test_highs_persistent.py +++ b/pyomo/contrib/appsi/solvers/tests/test_highs_persistent.py @@ -13,12 +13,11 @@ import sys import pyomo.common.unittest as unittest -import pyomo.environ as pe +import pyomo.environ as pyo from pyomo.common.log import LoggingIntercept from pyomo.common.tee import capture_output from pyomo.contrib.appsi.solvers.highs import Highs -from pyomo.contrib.appsi.base import TerminationCondition opt = Highs() @@ -28,16 +27,16 @@ class TestBugs(unittest.TestCase): def test_mutable_params_with_remove_cons(self): - m = pe.ConcreteModel() - m.x = pe.Var(bounds=(-10, 10)) - m.y = pe.Var() + m = pyo.ConcreteModel() + m.x = pyo.Var(bounds=(-10, 10)) + m.y = pyo.Var() - m.p1 = pe.Param(mutable=True) - m.p2 = pe.Param(mutable=True) + m.p1 = pyo.Param(mutable=True) + m.p2 = pyo.Param(mutable=True) - m.obj = pe.Objective(expr=m.y) - m.c1 = pe.Constraint(expr=m.y >= m.x + m.p1) - m.c2 = pe.Constraint(expr=m.y >= -m.x + m.p2) + m.obj = pyo.Objective(expr=m.y) + m.c1 = pyo.Constraint(expr=m.y >= m.x + m.p1) + m.c2 = pyo.Constraint(expr=m.y >= -m.x + m.p2) m.p1.value = 1 m.p2.value = 1 @@ -52,19 +51,19 @@ def test_mutable_params_with_remove_cons(self): self.assertAlmostEqual(res.best_feasible_objective, -8) def test_mutable_params_with_remove_vars(self): - m = pe.ConcreteModel() - m.x = pe.Var() - m.y = pe.Var() + m = pyo.ConcreteModel() + m.x = pyo.Var() + m.y = pyo.Var() - m.p1 = pe.Param(mutable=True) - m.p2 = pe.Param(mutable=True) + m.p1 = pyo.Param(mutable=True) + m.p2 = pyo.Param(mutable=True) m.y.setlb(m.p1) m.y.setub(m.p2) - m.obj = pe.Objective(expr=m.y) - m.c1 = pe.Constraint(expr=m.y >= m.x + 1) - m.c2 = pe.Constraint(expr=m.y >= -m.x + 1) + m.obj = pyo.Objective(expr=m.y) + m.c1 = pyo.Constraint(expr=m.y >= m.x + 1) + m.c2 = pyo.Constraint(expr=m.y >= -m.x + 1) m.p1.value = -10 m.p2.value = 10 @@ -83,16 +82,16 @@ def test_mutable_params_with_remove_vars(self): def test_fix_and_unfix(self): # Tests issue https://github.com/Pyomo/pyomo/issues/3127 - m = pe.ConcreteModel() - m.x = pe.Var(domain=pe.Binary) - m.y = pe.Var(domain=pe.Binary) - m.fx = pe.Var(domain=pe.NonNegativeReals) - m.fy = pe.Var(domain=pe.NonNegativeReals) - m.c1 = pe.Constraint(expr=m.fx <= m.x) - m.c2 = pe.Constraint(expr=m.fy <= m.y) - m.c3 = pe.Constraint(expr=m.x + m.y <= 1) + m = pyo.ConcreteModel() + m.x = pyo.Var(domain=pyo.Binary) + m.y = pyo.Var(domain=pyo.Binary) + m.fx = pyo.Var(domain=pyo.NonNegativeReals) + m.fy = pyo.Var(domain=pyo.NonNegativeReals) + m.c1 = pyo.Constraint(expr=m.fx <= m.x) + m.c2 = pyo.Constraint(expr=m.fy <= m.y) + m.c3 = pyo.Constraint(expr=m.x + m.y <= 1) - m.obj = pe.Objective(expr=m.fx * 0.5 + m.fy * 0.4, sense=pe.maximize) + m.obj = pyo.Objective(expr=m.fx * 0.5 + m.fy * 0.4, sense=pyo.maximize) opt = Highs() @@ -157,21 +156,21 @@ def test_capture_highs_output(self): self.assertEqual(ref, OUT.getvalue()[-len(ref) :]) def test_warm_start(self): - m = pe.ConcreteModel() + m = pyo.ConcreteModel() # decision variables - m.x1 = pe.Var(domain=pe.Integers, name="x1", bounds=(0, 10)) - m.x2 = pe.Var(domain=pe.Reals, name="x2", bounds=(0, 10)) - m.x3 = pe.Var(domain=pe.Binary, name="x3") + m.x1 = pyo.Var(domain=pyo.Integers, name="x1", bounds=(0, 10)) + m.x2 = pyo.Var(domain=pyo.Reals, name="x2", bounds=(0, 10)) + m.x3 = pyo.Var(domain=pyo.Binary, name="x3") # objective function - m.OBJ = pe.Objective(expr=(3 * m.x1 + 2 * m.x2 + 4 * m.x3), sense=pe.maximize) + m.OBJ = pyo.Objective(expr=(3 * m.x1 + 2 * m.x2 + 4 * m.x3), sense=pyo.maximize) # constraints - m.C1 = pe.Constraint(expr=m.x1 + m.x2 <= 9) - m.C2 = pe.Constraint(expr=3 * m.x1 + m.x2 <= 18) - m.C3 = pe.Constraint(expr=m.x1 <= 7) - m.C4 = pe.Constraint(expr=m.x2 <= 6) + m.C1 = pyo.Constraint(expr=m.x1 + m.x2 <= 9) + m.C2 = pyo.Constraint(expr=3 * m.x1 + m.x2 <= 18) + m.C3 = pyo.Constraint(expr=m.x1 <= 7) + m.C4 = pyo.Constraint(expr=m.x2 <= 6) # MIP start m.x1 = 4 @@ -180,6 +179,24 @@ def test_warm_start(self): # solving process with capture_output() as output: - pe.SolverFactory("appsi_highs").solve(m, tee=True, warmstart=True) + pyo.SolverFactory("appsi_highs").solve(m, tee=True, warmstart=True) log = output.getvalue() self.assertIn("MIP start solution is feasible, objective value is 25", log) + + def test_qp(self): + # test issue #3381 + m = pyo.ConcreteModel() + + m.x1 = pyo.Var(name='x1', domain=pyo.Reals) + m.x2 = pyo.Var(name='x2', domain=pyo.Reals) + + # Quadratic Objective function + m.obj = pyo.Objective(expr=m.x1 * m.x1 + m.x2 * m.x2, sense=pyo.minimize) + + m.con1 = pyo.Constraint(expr=m.x1 >= 1) + m.con2 = pyo.Constraint(expr=m.x2 >= 1) + + opt = Highs() + opt.solve(m) + self.assertAlmostEqual(m.x1.value, 1, places=5) + self.assertAlmostEqual(m.x2.value, 1, places=5) From 823d882e189a1eaa764c61c9a82eea252c31a639 Mon Sep 17 00:00:00 2001 From: Arnaud Baguet Date: Mon, 24 Mar 2025 20:59:53 -0400 Subject: [PATCH 02/21] wrong result in test qp --- .../solvers/tests/test_highs_persistent.py | 90 +++++------ pyomo/contrib/solver/solvers/highs.py | 145 +++++++++++++++++- .../solver/tests/solvers/test_highs.py | 134 ++++++++++++++++ 3 files changed, 313 insertions(+), 56 deletions(-) create mode 100644 pyomo/contrib/solver/tests/solvers/test_highs.py diff --git a/pyomo/contrib/appsi/solvers/tests/test_highs_persistent.py b/pyomo/contrib/appsi/solvers/tests/test_highs_persistent.py index c5abd21e801..0bdd09038c7 100644 --- a/pyomo/contrib/appsi/solvers/tests/test_highs_persistent.py +++ b/pyomo/contrib/appsi/solvers/tests/test_highs_persistent.py @@ -13,7 +13,7 @@ import sys import pyomo.common.unittest as unittest -import pyomo.environ as pyo +import pyomo.environ as pe from pyomo.common.log import LoggingIntercept from pyomo.common.tee import capture_output @@ -27,16 +27,16 @@ class TestBugs(unittest.TestCase): def test_mutable_params_with_remove_cons(self): - m = pyo.ConcreteModel() - m.x = pyo.Var(bounds=(-10, 10)) - m.y = pyo.Var() + m = pe.ConcreteModel() + m.x = pe.Var(bounds=(-10, 10)) + m.y = pe.Var() - m.p1 = pyo.Param(mutable=True) - m.p2 = pyo.Param(mutable=True) + m.p1 = pe.Param(mutable=True) + m.p2 = pe.Param(mutable=True) - m.obj = pyo.Objective(expr=m.y) - m.c1 = pyo.Constraint(expr=m.y >= m.x + m.p1) - m.c2 = pyo.Constraint(expr=m.y >= -m.x + m.p2) + m.obj = pe.Objective(expr=m.y) + m.c1 = pe.Constraint(expr=m.y >= m.x + m.p1) + m.c2 = pe.Constraint(expr=m.y >= -m.x + m.p2) m.p1.value = 1 m.p2.value = 1 @@ -51,19 +51,19 @@ def test_mutable_params_with_remove_cons(self): self.assertAlmostEqual(res.best_feasible_objective, -8) def test_mutable_params_with_remove_vars(self): - m = pyo.ConcreteModel() - m.x = pyo.Var() - m.y = pyo.Var() + m = pe.ConcreteModel() + m.x = pe.Var() + m.y = pe.Var() - m.p1 = pyo.Param(mutable=True) - m.p2 = pyo.Param(mutable=True) + m.p1 = pe.Param(mutable=True) + m.p2 = pe.Param(mutable=True) m.y.setlb(m.p1) m.y.setub(m.p2) - m.obj = pyo.Objective(expr=m.y) - m.c1 = pyo.Constraint(expr=m.y >= m.x + 1) - m.c2 = pyo.Constraint(expr=m.y >= -m.x + 1) + m.obj = pe.Objective(expr=m.y) + m.c1 = pe.Constraint(expr=m.y >= m.x + 1) + m.c2 = pe.Constraint(expr=m.y >= -m.x + 1) m.p1.value = -10 m.p2.value = 10 @@ -82,16 +82,16 @@ def test_mutable_params_with_remove_vars(self): def test_fix_and_unfix(self): # Tests issue https://github.com/Pyomo/pyomo/issues/3127 - m = pyo.ConcreteModel() - m.x = pyo.Var(domain=pyo.Binary) - m.y = pyo.Var(domain=pyo.Binary) - m.fx = pyo.Var(domain=pyo.NonNegativeReals) - m.fy = pyo.Var(domain=pyo.NonNegativeReals) - m.c1 = pyo.Constraint(expr=m.fx <= m.x) - m.c2 = pyo.Constraint(expr=m.fy <= m.y) - m.c3 = pyo.Constraint(expr=m.x + m.y <= 1) + m = pe.ConcreteModel() + m.x = pe.Var(domain=pe.Binary) + m.y = pe.Var(domain=pe.Binary) + m.fx = pe.Var(domain=pe.NonNegativeReals) + m.fy = pe.Var(domain=pe.NonNegativeReals) + m.c1 = pe.Constraint(expr=m.fx <= m.x) + m.c2 = pe.Constraint(expr=m.fy <= m.y) + m.c3 = pe.Constraint(expr=m.x + m.y <= 1) - m.obj = pyo.Objective(expr=m.fx * 0.5 + m.fy * 0.4, sense=pyo.maximize) + m.obj = pe.Objective(expr=m.fx * 0.5 + m.fy * 0.4, sense=pe.maximize) opt = Highs() @@ -156,21 +156,21 @@ def test_capture_highs_output(self): self.assertEqual(ref, OUT.getvalue()[-len(ref) :]) def test_warm_start(self): - m = pyo.ConcreteModel() + m = pe.ConcreteModel() # decision variables - m.x1 = pyo.Var(domain=pyo.Integers, name="x1", bounds=(0, 10)) - m.x2 = pyo.Var(domain=pyo.Reals, name="x2", bounds=(0, 10)) - m.x3 = pyo.Var(domain=pyo.Binary, name="x3") + m.x1 = pe.Var(domain=pe.Integers, name="x1", bounds=(0, 10)) + m.x2 = pe.Var(domain=pe.Reals, name="x2", bounds=(0, 10)) + m.x3 = pe.Var(domain=pe.Binary, name="x3") # objective function - m.OBJ = pyo.Objective(expr=(3 * m.x1 + 2 * m.x2 + 4 * m.x3), sense=pyo.maximize) + m.OBJ = pe.Objective(expr=(3 * m.x1 + 2 * m.x2 + 4 * m.x3), sense=pe.maximize) # constraints - m.C1 = pyo.Constraint(expr=m.x1 + m.x2 <= 9) - m.C2 = pyo.Constraint(expr=3 * m.x1 + m.x2 <= 18) - m.C3 = pyo.Constraint(expr=m.x1 <= 7) - m.C4 = pyo.Constraint(expr=m.x2 <= 6) + m.C1 = pe.Constraint(expr=m.x1 + m.x2 <= 9) + m.C2 = pe.Constraint(expr=3 * m.x1 + m.x2 <= 18) + m.C3 = pe.Constraint(expr=m.x1 <= 7) + m.C4 = pe.Constraint(expr=m.x2 <= 6) # MIP start m.x1 = 4 @@ -179,24 +179,6 @@ def test_warm_start(self): # solving process with capture_output() as output: - pyo.SolverFactory("appsi_highs").solve(m, tee=True, warmstart=True) + pe.SolverFactory("appsi_highs").solve(m, tee=True, warmstart=True) log = output.getvalue() self.assertIn("MIP start solution is feasible, objective value is 25", log) - - def test_qp(self): - # test issue #3381 - m = pyo.ConcreteModel() - - m.x1 = pyo.Var(name='x1', domain=pyo.Reals) - m.x2 = pyo.Var(name='x2', domain=pyo.Reals) - - # Quadratic Objective function - m.obj = pyo.Objective(expr=m.x1 * m.x1 + m.x2 * m.x2, sense=pyo.minimize) - - m.con1 = pyo.Constraint(expr=m.x1 >= 1) - m.con2 = pyo.Constraint(expr=m.x2 >= 1) - - opt = Highs() - opt.solve(m) - self.assertAlmostEqual(m.x1.value, 1, places=5) - self.assertAlmostEqual(m.x2.value, 1, places=5) diff --git a/pyomo/contrib/solver/solvers/highs.py b/pyomo/contrib/solver/solvers/highs.py index 3facfc00d13..b34c4079303 100644 --- a/pyomo/contrib/solver/solvers/highs.py +++ b/pyomo/contrib/solver/solvers/highs.py @@ -96,6 +96,20 @@ def update(self): self.highs.changeColCost(col_ndx, value(self.expr)) +class _MutableQuadraticObjectiveCoefficient: + def __init__(self, expr, highs, value_updater): + self.expr = expr + self.highs = highs + self.value_updater = value_updater # Function to call when value changes + self._last_value = value(expr) + + def update(self): + new_value = value(self.expr) + if new_value != self._last_value: + self.value_updater(self._last_value, new_value) + self._last_value = new_value + + class _MutableObjectiveOffset: def __init__(self, expr, highs): self.expr = expr @@ -472,6 +486,8 @@ def update_parameters(self): self._sol = None if self._last_results_object is not None: self._last_results_object.solution_loader.invalidate() + self._rebuild_hessian = False + for con, helpers in self._mutable_helpers.items(): for helper in helpers: helper.update() @@ -480,6 +496,10 @@ def update_parameters(self): for helper in self._objective_helpers: helper.update() + # Rebuild Hessian if needed + if self._rebuild_hessian: + self._build_and_pass_hessian() + def _set_objective(self, obj): self._sol = None if self._last_results_object is not None: @@ -488,9 +508,19 @@ def _set_objective(self, obj): indices = np.arange(n) costs = np.zeros(n, dtype=np.double) self._objective_helpers = [] + self._quad_coefs = {} + if obj is None: sense = highspy.ObjSense.kMinimize self._solver_model.changeObjectiveOffset(0) + self._solver_model.passHessian( + 0, + 0, + highspy.HessianFormat.kTriangular, + np.array([], dtype=np.int32), + np.array([], dtype=np.int32), + np.array([], dtype=np.double), + ) else: if obj.sense == minimize: sense = highspy.ObjSense.kMinimize @@ -500,9 +530,9 @@ def _set_objective(self, obj): raise ValueError(f'Objective sense is not recognized: {obj.sense}') repn = generate_standard_repn( - obj.expr, quadratic=False, compute_values=False + obj.expr, quadratic=True, compute_values=False ) - if repn.nonlinear_expr is not None: + if repn.nonlinear_expr is not None and repn.polynomial_degree() > 2: raise IncompatibleModelError( f'Highs interface does not support expressions of degree {repn.polynomial_degree()}' ) @@ -530,6 +560,117 @@ def _set_objective(self, obj): self._solver_model.changeObjectiveSense(sense) self._solver_model.changeColsCost(n, indices, costs) + # Process quadratic terms if present and HiGHS has passHessian method + + if repn.quadratic_vars and len(repn.quadratic_vars) > 0: + # Dictionary to collect quadratic coefficients + quad_coefs = {} # (row, col) -> coefficient + + for ndx, (v1, v2) in enumerate(repn.quadratic_vars): + v1_ndx = self._pyomo_var_to_solver_var_map[id(v1)] + v2_ndx = self._pyomo_var_to_solver_var_map[id(v2)] + + # Ensure we're storing the lower triangular part + row = max(v1_ndx, v2_ndx) + col = min(v1_ndx, v2_ndx) + + coef = repn.quadratic_coefs[ndx] + coef_val = value(coef) + + # For different variables, divide by 2 since HiGHS expects the lower triangular part + # of Q where obj = 0.5*x^T*Q*x + c^T*x + if v1_ndx != v2_ndx: + coef_val /= 2.0 + + # Add to the dictionary + if (row, col) in quad_coefs: + quad_coefs[(row, col)] += coef_val + else: + quad_coefs[(row, col)] = coef_val + + # Handle mutable coefficients + if not is_constant(coef): + # Create a value updater function for this specific coefficient + def make_updater(row, col): + def update_quad_coef(old_val, new_val): + # Calculate the delta to add to the existing value + if v1_ndx != v2_ndx: + delta = (new_val - old_val) / 2.0 + else: + delta = new_val - old_val + + # Update the stored coefficient + self._quad_coefs[(row, col)] += delta + + # Signal that we need to rebuild the Hessian + self._rebuild_hessian = True + + return update_quad_coef + + updater = make_updater(row, col) + + # Store the mutable helper + helper = _MutableQuadraticObjectiveCoefficient( + expr=coef, highs=self._solver_model, value_updater=updater + ) + self._objective_helpers.append(helper) + + # Store quadratic coefficients for future updates + self._quad_coefs = quad_coefs + + # Build the CSC format arrays for HiGHS + self._build_and_pass_hessian() + + def _build_and_pass_hessian(self): + # Skip if no quadratic coefficients + if not self._quad_coefs: + return + + dim = len(self._pyomo_var_to_solver_var_map) + quad_coefs = self._quad_coefs + + # Build CSC format for the lower triangular part + q_value = [] + q_index = [] + q_start = [0] * (dim + 1) + + sorted_entries = sorted(quad_coefs.items(), key=lambda x: (x[0][1], x[0][0])) + + last_col = -1 + for (row, col), val in sorted_entries: + # Update column pointers + while col > last_col: + last_col += 1 + q_start[last_col + 1] = len(q_value) + + # Add the entry + q_index.append(row) + q_value.append(val) + + # Fill in remaining column pointers + while last_col < dim - 1: + last_col += 1 + q_start[last_col + 1] = len(q_value) + + # Create NumPy arrays + np_q_start = np.array(q_start, dtype=np.int32) + np_q_index = np.array(q_index, dtype=np.int32) + np_q_value = np.array(q_value, dtype=np.double) + + # Pass the Hessian to HiGHS + nnz = len(q_value) + self._solver_model.passHessian( + dim, + nnz, + highspy.HessianFormat.kTriangular, + np_q_start, + np_q_index, + np_q_value, + ) + + # Reset the rebuild flag + self._rebuild_hessian = False + def _postsolve(self): config = self._active_config timer = config.timer diff --git a/pyomo/contrib/solver/tests/solvers/test_highs.py b/pyomo/contrib/solver/tests/solvers/test_highs.py new file mode 100644 index 00000000000..08668a71ac8 --- /dev/null +++ b/pyomo/contrib/solver/tests/solvers/test_highs.py @@ -0,0 +1,134 @@ +# ___________________________________________________________________________ +# +# Pyomo: Python Optimization Modeling Objects +# Copyright (c) 2008-2025 +# National Technology and Engineering Solutions of Sandia, LLC +# Under the terms of Contract DE-NA0003525 with National Technology and +# Engineering Solutions of Sandia, LLC, the U.S. Government retains certain +# rights in this software. +# This software is distributed under the 3-clause BSD License. +# ___________________________________________________________________________ + +import subprocess +import sys + +import pyomo.common.unittest as unittest +import pyomo.environ as pe + +from pyomo.contrib.solver.solvers.highs import Highs +from pyomo.common.log import LoggingIntercept +from pyomo.common.tee import capture_output + +opt = Highs() +if not opt.available(): + raise unittest.SkipTest + + +class TestBugs(unittest.TestCase): + def test_mutable_params_with_remove_cons(self): + m = pe.ConcreteModel() + m.x = pe.Var(bounds=(-10, 10)) + m.y = pe.Var() + + m.p1 = pe.Param(mutable=True) + m.p2 = pe.Param(mutable=True) + + m.obj = pe.Objective(expr=m.y) + m.c1 = pe.Constraint(expr=m.y >= m.x + m.p1) + m.c2 = pe.Constraint(expr=m.y >= -m.x + m.p2) + + m.p1.value = 1 + m.p2.value = 1 + + opt = Highs() + res = opt.solve(m) + self.assertAlmostEqual(res.objective_bound, 1) + + del m.c1 + m.p2.value = 2 + res = opt.solve(m) + self.assertAlmostEqual(res.objective_bound, -8) + + def test_mutable_params_with_remove_vars(self): + m = pe.ConcreteModel() + m.x = pe.Var() + m.y = pe.Var() + + m.p1 = pe.Param(mutable=True) + m.p2 = pe.Param(mutable=True) + + m.y.setlb(m.p1) + m.y.setub(m.p2) + + m.obj = pe.Objective(expr=m.y) + m.c1 = pe.Constraint(expr=m.y >= m.x + 1) + m.c2 = pe.Constraint(expr=m.y >= -m.x + 1) + + m.p1.value = -10 + m.p2.value = 10 + + opt = Highs() + res = opt.solve(m) + self.assertAlmostEqual(res.objective_bound, 1) + + del m.c1 + del m.c2 + m.p1.value = -9 + m.p2.value = 9 + res = opt.solve(m) + self.assertAlmostEqual(res.objective_bound, -9) + + def test_fix_and_unfix(self): + # Tests issue https://github.com/Pyomo/pyomo/issues/3127 + + m = pe.ConcreteModel() + m.x = pe.Var(domain=pe.Binary) + m.y = pe.Var(domain=pe.Binary) + m.fx = pe.Var(domain=pe.NonNegativeReals) + m.fy = pe.Var(domain=pe.NonNegativeReals) + m.c1 = pe.Constraint(expr=m.fx <= m.x) + m.c2 = pe.Constraint(expr=m.fy <= m.y) + m.c3 = pe.Constraint(expr=m.x + m.y <= 1) + + m.obj = pe.Objective(expr=m.fx * 0.5 + m.fy * 0.4, sense=pe.maximize) + + opt = Highs() + + # solution 1 has m.x == 1 and m.y == 0 + r = opt.solve(m) + self.assertAlmostEqual(m.fx.value, 1, places=5) + self.assertAlmostEqual(m.fy.value, 0, places=5) + self.assertAlmostEqual(r.objective_bound, 0.5, places=5) + + # solution 2 has m.x == 0 and m.y == 1 + m.y.fix(1) + r = opt.solve(m) + self.assertAlmostEqual(m.fx.value, 0, places=5) + self.assertAlmostEqual(m.fy.value, 1, places=5) + self.assertAlmostEqual(r.objective_bound, 0.4, places=5) + + # solution 3 should be equal solution 1 + m.y.unfix() + m.x.fix(1) + r = opt.solve(m) + self.assertAlmostEqual(m.fx.value, 1, places=5) + self.assertAlmostEqual(m.fy.value, 0, places=5) + self.assertAlmostEqual(r.objective_bound, 0.5, places=5) + + def test_qp(self): + # test issue #3381 + m = pe.ConcreteModel() + + m.x1 = pe.Var(name='x1', domain=pe.Reals) + m.x2 = pe.Var(name='x2', domain=pe.Reals) + + # Quadratic Objective function + m.obj = pe.Objective(expr=m.x1 * m.x1 + m.x2 * m.x2, sense=pe.minimize) + + m.con1 = pe.Constraint(expr=m.x1 >= 1) + m.con2 = pe.Constraint(expr=m.x2 >= 1) + + results = opt.solve(m) + self.assertAlmostEqual(m.x1.value, 1, places=5) + self.assertAlmostEqual(m.x2.value, 1, places=5) + self.assertEqual(results.objective_bound, 2) From 2621b3ee0c476f5e241502b263e973c97c023612 Mon Sep 17 00:00:00 2001 From: Arnaud Baguet Date: Mon, 24 Mar 2025 21:36:11 -0400 Subject: [PATCH 03/21] pass qp tests --- pyomo/contrib/solver/solvers/highs.py | 22 +++++++++++-------- .../solver/tests/solvers/test_highs.py | 22 ++++++++++++++++++- 2 files changed, 34 insertions(+), 10 deletions(-) diff --git a/pyomo/contrib/solver/solvers/highs.py b/pyomo/contrib/solver/solvers/highs.py index b34c4079303..c1d283cd70f 100644 --- a/pyomo/contrib/solver/solvers/highs.py +++ b/pyomo/contrib/solver/solvers/highs.py @@ -577,10 +577,9 @@ def _set_objective(self, obj): coef = repn.quadratic_coefs[ndx] coef_val = value(coef) - # For different variables, divide by 2 since HiGHS expects the lower triangular part - # of Q where obj = 0.5*x^T*Q*x + c^T*x - if v1_ndx != v2_ndx: - coef_val /= 2.0 + # Adjust coefficient values for HiGHS's expected format + if v1_ndx == v2_ndx: + coef_val *= 2.0 # Add to the dictionary if (row, col) in quad_coefs: @@ -632,16 +631,17 @@ def _build_and_pass_hessian(self): # Build CSC format for the lower triangular part q_value = [] q_index = [] - q_start = [0] * (dim + 1) + # Create q_start with length dim (not dim+1) as HiGHS expects, contrary to standard CSC + q_start = [0] * dim sorted_entries = sorted(quad_coefs.items(), key=lambda x: (x[0][1], x[0][0])) last_col = -1 for (row, col), val in sorted_entries: - # Update column pointers while col > last_col: last_col += 1 - q_start[last_col + 1] = len(q_value) + if last_col < dim: + q_start[last_col] = len(q_value) # Add the entry q_index.append(row) @@ -650,7 +650,7 @@ def _build_and_pass_hessian(self): # Fill in remaining column pointers while last_col < dim - 1: last_col += 1 - q_start[last_col + 1] = len(q_value) + q_start[last_col] = len(q_value) # Create NumPy arrays np_q_start = np.array(q_start, dtype=np.int32) @@ -659,7 +659,7 @@ def _build_and_pass_hessian(self): # Pass the Hessian to HiGHS nnz = len(q_value) - self._solver_model.passHessian( + status = self._solver_model.passHessian( dim, nnz, highspy.HessianFormat.kTriangular, @@ -667,6 +667,10 @@ def _build_and_pass_hessian(self): np_q_index, np_q_value, ) + if status != highspy.HighsStatus.kOk: + logger.warning( + f"HiGHS returned non-OK status when passing Hessian: {status}" + ) # Reset the rebuild flag self._rebuild_hessian = False diff --git a/pyomo/contrib/solver/tests/solvers/test_highs.py b/pyomo/contrib/solver/tests/solvers/test_highs.py index 08668a71ac8..5daaf06f5f5 100644 --- a/pyomo/contrib/solver/tests/solvers/test_highs.py +++ b/pyomo/contrib/solver/tests/solvers/test_highs.py @@ -115,7 +115,7 @@ def test_fix_and_unfix(self): self.assertAlmostEqual(m.fy.value, 0, places=5) self.assertAlmostEqual(r.objective_bound, 0.5, places=5) - def test_qp(self): + def test_qp1(self): # test issue #3381 m = pe.ConcreteModel() @@ -132,3 +132,23 @@ def test_qp(self): self.assertAlmostEqual(m.x1.value, 1, places=5) self.assertAlmostEqual(m.x2.value, 1, places=5) self.assertEqual(results.objective_bound, 2) + + def test_qp2(self): + # test issue #3381 + m = pe.ConcreteModel() + + m.x1 = pe.Var(name='x1', domain=pe.Reals) + m.x2 = pe.Var(name='x2', domain=pe.Reals) + + # Quadratic Objective function + m.obj = pe.Objective( + expr=m.x1 * m.x1 + m.x2 * m.x2 - m.x1 * m.x2, sense=pe.minimize + ) + + m.con1 = pe.Constraint(expr=m.x1 >= 1) + m.con2 = pe.Constraint(expr=m.x2 >= 1) + + results = opt.solve(m) + self.assertAlmostEqual(m.x1.value, 1, places=5) + self.assertAlmostEqual(m.x2.value, 1, places=5) + self.assertEqual(results.objective_bound, 1) From 0815a03d2b951e85139d9efc5664c2aae167fd6b Mon Sep 17 00:00:00 2001 From: Arnaud Baguet Date: Mon, 24 Mar 2025 21:42:54 -0400 Subject: [PATCH 04/21] rename to quadratic_coefs for consistency with gurobi --- pyomo/contrib/solver/solvers/highs.py | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/pyomo/contrib/solver/solvers/highs.py b/pyomo/contrib/solver/solvers/highs.py index c1d283cd70f..93ccf398f48 100644 --- a/pyomo/contrib/solver/solvers/highs.py +++ b/pyomo/contrib/solver/solvers/highs.py @@ -508,7 +508,7 @@ def _set_objective(self, obj): indices = np.arange(n) costs = np.zeros(n, dtype=np.double) self._objective_helpers = [] - self._quad_coefs = {} + self.quadratic_coefs = {} if obj is None: sense = highspy.ObjSense.kMinimize @@ -564,7 +564,7 @@ def _set_objective(self, obj): if repn.quadratic_vars and len(repn.quadratic_vars) > 0: # Dictionary to collect quadratic coefficients - quad_coefs = {} # (row, col) -> coefficient + quadratic_coefs = {} # (row, col) -> coefficient for ndx, (v1, v2) in enumerate(repn.quadratic_vars): v1_ndx = self._pyomo_var_to_solver_var_map[id(v1)] @@ -582,10 +582,10 @@ def _set_objective(self, obj): coef_val *= 2.0 # Add to the dictionary - if (row, col) in quad_coefs: - quad_coefs[(row, col)] += coef_val + if (row, col) in quadratic_coefs: + quadratic_coefs[(row, col)] += coef_val else: - quad_coefs[(row, col)] = coef_val + quadratic_coefs[(row, col)] = coef_val # Handle mutable coefficients if not is_constant(coef): @@ -599,7 +599,7 @@ def update_quad_coef(old_val, new_val): delta = new_val - old_val # Update the stored coefficient - self._quad_coefs[(row, col)] += delta + self.quadratic_coefs[(row, col)] += delta # Signal that we need to rebuild the Hessian self._rebuild_hessian = True @@ -615,18 +615,18 @@ def update_quad_coef(old_val, new_val): self._objective_helpers.append(helper) # Store quadratic coefficients for future updates - self._quad_coefs = quad_coefs + self.quadratic_coefs = quadratic_coefs # Build the CSC format arrays for HiGHS self._build_and_pass_hessian() def _build_and_pass_hessian(self): # Skip if no quadratic coefficients - if not self._quad_coefs: + if not self.quadratic_coefs: return dim = len(self._pyomo_var_to_solver_var_map) - quad_coefs = self._quad_coefs + quad_coefs = self.quadratic_coefs # Build CSC format for the lower triangular part q_value = [] From 9bba7d04af802807bf006e81f30dcc3e5bf93520 Mon Sep 17 00:00:00 2001 From: Arnaud Baguet Date: Mon, 24 Mar 2025 21:43:25 -0400 Subject: [PATCH 05/21] rename to quadratic_coefs for consistency with gurobi 2 --- pyomo/contrib/solver/solvers/highs.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/pyomo/contrib/solver/solvers/highs.py b/pyomo/contrib/solver/solvers/highs.py index 93ccf398f48..64cd4d5b911 100644 --- a/pyomo/contrib/solver/solvers/highs.py +++ b/pyomo/contrib/solver/solvers/highs.py @@ -626,7 +626,6 @@ def _build_and_pass_hessian(self): return dim = len(self._pyomo_var_to_solver_var_map) - quad_coefs = self.quadratic_coefs # Build CSC format for the lower triangular part q_value = [] @@ -634,7 +633,9 @@ def _build_and_pass_hessian(self): # Create q_start with length dim (not dim+1) as HiGHS expects, contrary to standard CSC q_start = [0] * dim - sorted_entries = sorted(quad_coefs.items(), key=lambda x: (x[0][1], x[0][0])) + sorted_entries = sorted( + self.quadratic_coefs.items(), key=lambda x: (x[0][1], x[0][0]) + ) last_col = -1 for (row, col), val in sorted_entries: From 1b42058fdaf69c6ed8c43e8de83abb5bc653a9b6 Mon Sep 17 00:00:00 2001 From: Arnaud Baguet Date: Mon, 24 Mar 2025 23:15:34 -0400 Subject: [PATCH 06/21] introduce _MutableObjective --- pyomo/contrib/solver/solvers/highs.py | 299 +++++++++++++------------- 1 file changed, 153 insertions(+), 146 deletions(-) diff --git a/pyomo/contrib/solver/solvers/highs.py b/pyomo/contrib/solver/solvers/highs.py index 64cd4d5b911..3bd083fabf2 100644 --- a/pyomo/contrib/solver/solvers/highs.py +++ b/pyomo/contrib/solver/solvers/highs.py @@ -96,18 +96,136 @@ def update(self): self.highs.changeColCost(col_ndx, value(self.expr)) -class _MutableQuadraticObjectiveCoefficient: - def __init__(self, expr, highs, value_updater): +class _MutableQuadraticCoefficient: + def __init__(self, expr, row_idx, col_idx): self.expr = expr + self.row_idx = row_idx + self.col_idx = col_idx + + +class _MutableObjective: + def __init__(self, highs, constant, linear_coefs, quadratic_coefs): self.highs = highs - self.value_updater = value_updater # Function to call when value changes - self._last_value = value(expr) + self.constant = constant + self.linear_coefs = linear_coefs + self.quadratic_coefs = quadratic_coefs + self.last_quadratic_coef_values = [value(i.expr) for i in self.quadratic_coefs] + # Store the quadratic coefficients in dictionary format + self.quad_coef_dict = {} + self._initialize_quad_coef_dict() + # Flag to force first update of quadratic coefficients + self._first_update = True + + def _initialize_quad_coef_dict(self): + for coef in self.quadratic_coefs: + v1_ndx = coef.row_idx + v2_ndx = coef.col_idx + # Ensure we're storing the lower triangular part + row = max(v1_ndx, v2_ndx) + col = min(v1_ndx, v2_ndx) + + coef_val = value(coef.expr) + # Adjust for diagonal elements + if v1_ndx == v2_ndx: + coef_val *= 2.0 + + self.quad_coef_dict[(row, col)] = coef_val def update(self): - new_value = value(self.expr) - if new_value != self._last_value: - self.value_updater(self._last_value, new_value) - self._last_value = new_value + """ + Update the quadratic objective expression. + """ + needs_quadratic_update = self._first_update + + # Update linear coefficients + for coef in self.linear_coefs: + coef.update() + + # Update constant term + self.constant.update() + + # Check if quadratic coefficients changed + for ndx, coef in enumerate(self.quadratic_coefs): + current_val = value(coef.expr) + if current_val != self.last_quadratic_coef_values[ndx]: + needs_quadratic_update = True + + # Update the dictionary entry + v1_ndx = coef.row_idx + v2_ndx = coef.col_idx + row = max(v1_ndx, v2_ndx) + col = min(v1_ndx, v2_ndx) + + # Calculate the delta + delta = current_val - self.last_quadratic_coef_values[ndx] + if v1_ndx == v2_ndx: + delta *= 2.0 + + # Update the stored coefficient + if (row, col) in self.quad_coef_dict: + self.quad_coef_dict[(row, col)] += delta + else: + self.quad_coef_dict[(row, col)] = delta + + self.last_quadratic_coef_values[ndx] = current_val + + # If anything changed, rebuild and pass the Hessian + if needs_quadratic_update: + self._build_and_pass_hessian() + self._first_update = False + + def _build_and_pass_hessian(self): + """Build and pass the Hessian to HiGHS in CSC format""" + if not self.quad_coef_dict: + return + + dim = self.highs.getNumCol() + + # Build CSC format for the lower triangular part + q_value = [] + q_index = [] + q_start = [0] * dim + + sorted_entries = sorted( + self.quad_coef_dict.items(), key=lambda x: (x[0][1], x[0][0]) + ) + + last_col = -1 + for (row, col), val in sorted_entries: + while col > last_col: + last_col += 1 + if last_col < dim: + q_start[last_col] = len(q_value) + + # Add the entry + q_index.append(row) + q_value.append(val) + + # Fill in remaining column pointers + while last_col < dim - 1: + last_col += 1 + q_start[last_col] = len(q_value) + + # Create NumPy arrays + np_q_start = np.array(q_start, dtype=np.int32) + np_q_index = np.array(q_index, dtype=np.int32) + np_q_value = np.array(q_value, dtype=np.double) + + # Pass the Hessian to HiGHS + nnz = len(q_value) + status = self.highs.passHessian( + dim, + nnz, + highspy.HessianFormat.kTriangular, + np_q_start, + np_q_index, + np_q_value, + ) + + if status != highspy.HighsStatus.kOk: + logger.warning( + f"HiGHS returned non-OK status when passing Hessian: {status}" + ) class _MutableObjectiveOffset: @@ -155,7 +273,6 @@ def __init__(self, **kwds): self._solver_con_to_pyomo_con_map = {} self._mutable_helpers = {} self._mutable_bounds = {} - self._objective_helpers = [] self._last_results_object: Optional[Results] = None self._sol = None @@ -486,41 +603,29 @@ def update_parameters(self): self._sol = None if self._last_results_object is not None: self._last_results_object.solution_loader.invalidate() - self._rebuild_hessian = False for con, helpers in self._mutable_helpers.items(): for helper in helpers: helper.update() for k, (v, helper) in self._mutable_bounds.items(): helper.update() - for helper in self._objective_helpers: - helper.update() - # Rebuild Hessian if needed - if self._rebuild_hessian: - self._build_and_pass_hessian() + self._mutable_objective.update() def _set_objective(self, obj): self._sol = None if self._last_results_object is not None: self._last_results_object.solution_loader.invalidate() n = len(self._pyomo_var_to_solver_var_map) - indices = np.arange(n) costs = np.zeros(n, dtype=np.double) - self._objective_helpers = [] - self.quadratic_coefs = {} + + # Initialize empty lists for all coefficient types + mutable_linear_coefficients = [] + mutable_quadratic_coefficients = [] if obj is None: sense = highspy.ObjSense.kMinimize self._solver_model.changeObjectiveOffset(0) - self._solver_model.passHessian( - 0, - 0, - highspy.HessianFormat.kTriangular, - np.array([], dtype=np.int32), - np.array([], dtype=np.int32), - np.array([], dtype=np.double), - ) else: if obj.sense == minimize: sense = highspy.ObjSense.kMinimize @@ -548,133 +653,35 @@ def _set_objective(self, obj): expr=coef, highs=self._solver_model, ) - self._objective_helpers.append(mutable_objective_coef) - - self._solver_model.changeObjectiveOffset(value(repn.constant)) - if not is_constant(repn.constant): - mutable_objective_offset = _MutableObjectiveOffset( - expr=repn.constant, highs=self._solver_model - ) - self._objective_helpers.append(mutable_objective_offset) - - self._solver_model.changeObjectiveSense(sense) - self._solver_model.changeColsCost(n, indices, costs) + mutable_linear_coefficients.append(mutable_objective_coef) - # Process quadratic terms if present and HiGHS has passHessian method - - if repn.quadratic_vars and len(repn.quadratic_vars) > 0: - # Dictionary to collect quadratic coefficients - quadratic_coefs = {} # (row, col) -> coefficient - - for ndx, (v1, v2) in enumerate(repn.quadratic_vars): - v1_ndx = self._pyomo_var_to_solver_var_map[id(v1)] - v2_ndx = self._pyomo_var_to_solver_var_map[id(v2)] - - # Ensure we're storing the lower triangular part - row = max(v1_ndx, v2_ndx) - col = min(v1_ndx, v2_ndx) - - coef = repn.quadratic_coefs[ndx] - coef_val = value(coef) - - # Adjust coefficient values for HiGHS's expected format - if v1_ndx == v2_ndx: - coef_val *= 2.0 - - # Add to the dictionary - if (row, col) in quadratic_coefs: - quadratic_coefs[(row, col)] += coef_val - else: - quadratic_coefs[(row, col)] = coef_val - - # Handle mutable coefficients - if not is_constant(coef): - # Create a value updater function for this specific coefficient - def make_updater(row, col): - def update_quad_coef(old_val, new_val): - # Calculate the delta to add to the existing value - if v1_ndx != v2_ndx: - delta = (new_val - old_val) / 2.0 - else: - delta = new_val - old_val - - # Update the stored coefficient - self.quadratic_coefs[(row, col)] += delta - - # Signal that we need to rebuild the Hessian - self._rebuild_hessian = True + mutable_constant = _MutableObjectiveOffset( + expr=repn.constant, highs=self._solver_model + ) - return update_quad_coef + if repn.quadratic_vars and len(repn.quadratic_vars) > 0: + for ndx, (v1, v2) in enumerate(repn.quadratic_vars): + v1_id = id(v1) + v2_id = id(v2) + v1_ndx = self._pyomo_var_to_solver_var_map[v1_id] + v2_ndx = self._pyomo_var_to_solver_var_map[v2_id] - updater = make_updater(row, col) + coef = repn.quadratic_coefs[ndx] - # Store the mutable helper - helper = _MutableQuadraticObjectiveCoefficient( - expr=coef, highs=self._solver_model, value_updater=updater + mutable_quadratic_coefficients.append( + _MutableQuadraticCoefficient( + expr=coef, row_idx=v1_ndx, col_idx=v2_ndx + ) ) - self._objective_helpers.append(helper) - - # Store quadratic coefficients for future updates - self.quadratic_coefs = quadratic_coefs - - # Build the CSC format arrays for HiGHS - self._build_and_pass_hessian() - - def _build_and_pass_hessian(self): - # Skip if no quadratic coefficients - if not self.quadratic_coefs: - return - - dim = len(self._pyomo_var_to_solver_var_map) - - # Build CSC format for the lower triangular part - q_value = [] - q_index = [] - # Create q_start with length dim (not dim+1) as HiGHS expects, contrary to standard CSC - q_start = [0] * dim - - sorted_entries = sorted( - self.quadratic_coefs.items(), key=lambda x: (x[0][1], x[0][0]) - ) - - last_col = -1 - for (row, col), val in sorted_entries: - while col > last_col: - last_col += 1 - if last_col < dim: - q_start[last_col] = len(q_value) - - # Add the entry - q_index.append(row) - q_value.append(val) - # Fill in remaining column pointers - while last_col < dim - 1: - last_col += 1 - q_start[last_col] = len(q_value) - - # Create NumPy arrays - np_q_start = np.array(q_start, dtype=np.int32) - np_q_index = np.array(q_index, dtype=np.int32) - np_q_value = np.array(q_value, dtype=np.double) - - # Pass the Hessian to HiGHS - nnz = len(q_value) - status = self._solver_model.passHessian( - dim, - nnz, - highspy.HessianFormat.kTriangular, - np_q_start, - np_q_index, - np_q_value, + self._solver_model.changeObjectiveSense(sense) + self._mutable_objective = _MutableObjective( + self._solver_model, + mutable_constant, + mutable_linear_coefficients, + mutable_quadratic_coefficients, ) - if status != highspy.HighsStatus.kOk: - logger.warning( - f"HiGHS returned non-OK status when passing Hessian: {status}" - ) - - # Reset the rebuild flag - self._rebuild_hessian = False + self._mutable_objective.update() def _postsolve(self): config = self._active_config From c659ccc803d56537c361307dc7a76b2edfafa2eb Mon Sep 17 00:00:00 2001 From: Arnaud Baguet Date: Mon, 24 Mar 2025 23:49:25 -0400 Subject: [PATCH 07/21] fix persistency delta bug --- pyomo/contrib/solver/solvers/highs.py | 31 +++++-------------- .../solver/tests/solvers/test_highs.py | 13 ++++++-- 2 files changed, 18 insertions(+), 26 deletions(-) diff --git a/pyomo/contrib/solver/solvers/highs.py b/pyomo/contrib/solver/solvers/highs.py index 3bd083fabf2..fb1c5667afd 100644 --- a/pyomo/contrib/solver/solvers/highs.py +++ b/pyomo/contrib/solver/solvers/highs.py @@ -137,35 +137,25 @@ def update(self): """ needs_quadratic_update = self._first_update - # Update linear coefficients + self.constant.update() for coef in self.linear_coefs: coef.update() - # Update constant term - self.constant.update() - - # Check if quadratic coefficients changed for ndx, coef in enumerate(self.quadratic_coefs): current_val = value(coef.expr) if current_val != self.last_quadratic_coef_values[ndx]: needs_quadratic_update = True - # Update the dictionary entry v1_ndx = coef.row_idx v2_ndx = coef.col_idx row = max(v1_ndx, v2_ndx) col = min(v1_ndx, v2_ndx) - # Calculate the delta - delta = current_val - self.last_quadratic_coef_values[ndx] + # Adjust the diagonal to match Highs' expected format if v1_ndx == v2_ndx: - delta *= 2.0 + current_val *= 2.0 - # Update the stored coefficient - if (row, col) in self.quad_coef_dict: - self.quad_coef_dict[(row, col)] += delta - else: - self.quad_coef_dict[(row, col)] = delta + self.quad_coef_dict[(row, col)] = current_val self.last_quadratic_coef_values[ndx] = current_val @@ -201,25 +191,18 @@ def _build_and_pass_hessian(self): q_index.append(row) q_value.append(val) - # Fill in remaining column pointers while last_col < dim - 1: last_col += 1 q_start[last_col] = len(q_value) - # Create NumPy arrays - np_q_start = np.array(q_start, dtype=np.int32) - np_q_index = np.array(q_index, dtype=np.int32) - np_q_value = np.array(q_value, dtype=np.double) - - # Pass the Hessian to HiGHS nnz = len(q_value) status = self.highs.passHessian( dim, nnz, highspy.HessianFormat.kTriangular, - np_q_start, - np_q_index, - np_q_value, + np.array(q_start, dtype=np.int32), + np.array(q_index, dtype=np.int32), + np.array(q_value, dtype=np.double), ) if status != highspy.HighsStatus.kOk: diff --git a/pyomo/contrib/solver/tests/solvers/test_highs.py b/pyomo/contrib/solver/tests/solvers/test_highs.py index 5daaf06f5f5..3cfe289645f 100644 --- a/pyomo/contrib/solver/tests/solvers/test_highs.py +++ b/pyomo/contrib/solver/tests/solvers/test_highs.py @@ -140,15 +140,24 @@ def test_qp2(self): m.x1 = pe.Var(name='x1', domain=pe.Reals) m.x2 = pe.Var(name='x2', domain=pe.Reals) - # Quadratic Objective function + m.p = pe.Param(initialize=1, mutable=True) + m.obj = pe.Objective( - expr=m.x1 * m.x1 + m.x2 * m.x2 - m.x1 * m.x2, sense=pe.minimize + expr=m.p * m.x1 * m.x1 + m.x2 * m.x2 - m.x1 * m.x2, sense=pe.minimize ) m.con1 = pe.Constraint(expr=m.x1 >= 1) m.con2 = pe.Constraint(expr=m.x2 >= 1) + opt.set_instance(m) results = opt.solve(m) self.assertAlmostEqual(m.x1.value, 1, places=5) self.assertAlmostEqual(m.x2.value, 1, places=5) self.assertEqual(results.objective_bound, 1) + + m.p.value = 2.0 + opt.update_parameters() + results = opt.solve(m) + self.assertAlmostEqual(m.x1.value, 1, places=5) + self.assertAlmostEqual(m.x2.value, 1, places=5) + self.assertEqual(results.objective_bound, 2) From 38a0a85269170b1cd9b8d9e8a3e03eb06ebf5fca Mon Sep 17 00:00:00 2001 From: Arnaud Baguet Date: Tue, 25 Mar 2025 00:00:33 -0400 Subject: [PATCH 08/21] or polynomial degree > 2 --- pyomo/contrib/solver/solvers/highs.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyomo/contrib/solver/solvers/highs.py b/pyomo/contrib/solver/solvers/highs.py index fb1c5667afd..4c720ac0925 100644 --- a/pyomo/contrib/solver/solvers/highs.py +++ b/pyomo/contrib/solver/solvers/highs.py @@ -620,7 +620,7 @@ def _set_objective(self, obj): repn = generate_standard_repn( obj.expr, quadratic=True, compute_values=False ) - if repn.nonlinear_expr is not None and repn.polynomial_degree() > 2: + if repn.nonlinear_expr is not None or repn.polynomial_degree() > 2: raise IncompatibleModelError( f'Highs interface does not support expressions of degree {repn.polynomial_degree()}' ) From 6c4b7b11cb4c7dd87848c8fa50fc886459b87cd7 Mon Sep 17 00:00:00 2001 From: Arnaud Baguet Date: Tue, 25 Mar 2025 00:09:09 -0400 Subject: [PATCH 09/21] remove changes to legacy appsi_highs --- pyomo/contrib/appsi/solvers/highs.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyomo/contrib/appsi/solvers/highs.py b/pyomo/contrib/appsi/solvers/highs.py index 83a45d68be3..edaafa5d1d4 100644 --- a/pyomo/contrib/appsi/solvers/highs.py +++ b/pyomo/contrib/appsi/solvers/highs.py @@ -616,7 +616,7 @@ def _set_objective(self, obj): ) repn = generate_standard_repn( - obj.expr, quadratic=True, compute_values=False + obj.expr, quadratic=False, compute_values=False ) if repn.nonlinear_expr is not None: raise DegreeError( From 891483e7155e93562d45b697d1e1cc07ca53eb61 Mon Sep 17 00:00:00 2001 From: Arnaud Baguet Date: Tue, 25 Mar 2025 21:10:52 -0400 Subject: [PATCH 10/21] add back mutable vs non mutable handling --- pyomo/contrib/solver/solvers/highs.py | 2 ++ pyomo/contrib/solver/tests/solvers/test_highs.py | 5 ----- 2 files changed, 2 insertions(+), 5 deletions(-) diff --git a/pyomo/contrib/solver/solvers/highs.py b/pyomo/contrib/solver/solvers/highs.py index 4c720ac0925..21ca808afb7 100644 --- a/pyomo/contrib/solver/solvers/highs.py +++ b/pyomo/contrib/solver/solvers/highs.py @@ -600,6 +600,7 @@ def _set_objective(self, obj): if self._last_results_object is not None: self._last_results_object.solution_loader.invalidate() n = len(self._pyomo_var_to_solver_var_map) + indices = np.arange(n) costs = np.zeros(n, dtype=np.double) # Initialize empty lists for all coefficient types @@ -658,6 +659,7 @@ def _set_objective(self, obj): ) self._solver_model.changeObjectiveSense(sense) + self._solver_model.changeColsCost(n, indices, costs) self._mutable_objective = _MutableObjective( self._solver_model, mutable_constant, diff --git a/pyomo/contrib/solver/tests/solvers/test_highs.py b/pyomo/contrib/solver/tests/solvers/test_highs.py index 3cfe289645f..be1a2cf056d 100644 --- a/pyomo/contrib/solver/tests/solvers/test_highs.py +++ b/pyomo/contrib/solver/tests/solvers/test_highs.py @@ -9,15 +9,10 @@ # This software is distributed under the 3-clause BSD License. # ___________________________________________________________________________ -import subprocess -import sys - import pyomo.common.unittest as unittest import pyomo.environ as pe from pyomo.contrib.solver.solvers.highs import Highs -from pyomo.common.log import LoggingIntercept -from pyomo.common.tee import capture_output opt = Highs() if not opt.available(): From 6d6c2d2e0524db813da209f0ef44dd9c9d78539a Mon Sep 17 00:00:00 2001 From: Arnaud Baguet Date: Wed, 26 Mar 2025 21:20:28 -0400 Subject: [PATCH 11/21] fix more unit tests --- pyomo/contrib/solver/solvers/highs.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyomo/contrib/solver/solvers/highs.py b/pyomo/contrib/solver/solvers/highs.py index 21ca808afb7..c554682ca9a 100644 --- a/pyomo/contrib/solver/solvers/highs.py +++ b/pyomo/contrib/solver/solvers/highs.py @@ -609,7 +609,7 @@ def _set_objective(self, obj): if obj is None: sense = highspy.ObjSense.kMinimize - self._solver_model.changeObjectiveOffset(0) + mutable_constant = _MutableObjectiveOffset(expr=0, highs=self._solver_model) else: if obj.sense == minimize: sense = highspy.ObjSense.kMinimize From 3d71b39dc7e0c96d1d42d55edcc400cf064d7cb2 Mon Sep 17 00:00:00 2001 From: Arnaud Baguet Date: Thu, 10 Apr 2025 20:00:14 -0400 Subject: [PATCH 12/21] standardize pyomo.environ --- .../solver/tests/solvers/test_highs.py | 80 +++++++++---------- 1 file changed, 40 insertions(+), 40 deletions(-) diff --git a/pyomo/contrib/solver/tests/solvers/test_highs.py b/pyomo/contrib/solver/tests/solvers/test_highs.py index be1a2cf056d..ef557fbdaed 100644 --- a/pyomo/contrib/solver/tests/solvers/test_highs.py +++ b/pyomo/contrib/solver/tests/solvers/test_highs.py @@ -10,7 +10,7 @@ # ___________________________________________________________________________ import pyomo.common.unittest as unittest -import pyomo.environ as pe +import pyomo.environ as pyo from pyomo.contrib.solver.solvers.highs import Highs @@ -21,16 +21,16 @@ class TestBugs(unittest.TestCase): def test_mutable_params_with_remove_cons(self): - m = pe.ConcreteModel() - m.x = pe.Var(bounds=(-10, 10)) - m.y = pe.Var() + m = pyo.ConcreteModel() + m.x = pyo.Var(bounds=(-10, 10)) + m.y = pyo.Var() - m.p1 = pe.Param(mutable=True) - m.p2 = pe.Param(mutable=True) + m.p1 = pyo.Param(mutable=True) + m.p2 = pyo.Param(mutable=True) - m.obj = pe.Objective(expr=m.y) - m.c1 = pe.Constraint(expr=m.y >= m.x + m.p1) - m.c2 = pe.Constraint(expr=m.y >= -m.x + m.p2) + m.obj = pyo.Objective(expr=m.y) + m.c1 = pyo.Constraint(expr=m.y >= m.x + m.p1) + m.c2 = pyo.Constraint(expr=m.y >= -m.x + m.p2) m.p1.value = 1 m.p2.value = 1 @@ -45,19 +45,19 @@ def test_mutable_params_with_remove_cons(self): self.assertAlmostEqual(res.objective_bound, -8) def test_mutable_params_with_remove_vars(self): - m = pe.ConcreteModel() - m.x = pe.Var() - m.y = pe.Var() + m = pyo.ConcreteModel() + m.x = pyo.Var() + m.y = pyo.Var() - m.p1 = pe.Param(mutable=True) - m.p2 = pe.Param(mutable=True) + m.p1 = pyo.Param(mutable=True) + m.p2 = pyo.Param(mutable=True) m.y.setlb(m.p1) m.y.setub(m.p2) - m.obj = pe.Objective(expr=m.y) - m.c1 = pe.Constraint(expr=m.y >= m.x + 1) - m.c2 = pe.Constraint(expr=m.y >= -m.x + 1) + m.obj = pyo.Objective(expr=m.y) + m.c1 = pyo.Constraint(expr=m.y >= m.x + 1) + m.c2 = pyo.Constraint(expr=m.y >= -m.x + 1) m.p1.value = -10 m.p2.value = 10 @@ -76,16 +76,16 @@ def test_mutable_params_with_remove_vars(self): def test_fix_and_unfix(self): # Tests issue https://github.com/Pyomo/pyomo/issues/3127 - m = pe.ConcreteModel() - m.x = pe.Var(domain=pe.Binary) - m.y = pe.Var(domain=pe.Binary) - m.fx = pe.Var(domain=pe.NonNegativeReals) - m.fy = pe.Var(domain=pe.NonNegativeReals) - m.c1 = pe.Constraint(expr=m.fx <= m.x) - m.c2 = pe.Constraint(expr=m.fy <= m.y) - m.c3 = pe.Constraint(expr=m.x + m.y <= 1) + m = pyo.ConcreteModel() + m.x = pyo.Var(domain=pyo.Binary) + m.y = pyo.Var(domain=pyo.Binary) + m.fx = pyo.Var(domain=pyo.NonNegativeReals) + m.fy = pyo.Var(domain=pyo.NonNegativeReals) + m.c1 = pyo.Constraint(expr=m.fx <= m.x) + m.c2 = pyo.Constraint(expr=m.fy <= m.y) + m.c3 = pyo.Constraint(expr=m.x + m.y <= 1) - m.obj = pe.Objective(expr=m.fx * 0.5 + m.fy * 0.4, sense=pe.maximize) + m.obj = pyo.Objective(expr=m.fx * 0.5 + m.fy * 0.4, sense=pyo.maximize) opt = Highs() @@ -112,16 +112,16 @@ def test_fix_and_unfix(self): def test_qp1(self): # test issue #3381 - m = pe.ConcreteModel() + m = pyo.ConcreteModel() - m.x1 = pe.Var(name='x1', domain=pe.Reals) - m.x2 = pe.Var(name='x2', domain=pe.Reals) + m.x1 = pyo.Var(name='x1', domain=pyo.Reals) + m.x2 = pyo.Var(name='x2', domain=pyo.Reals) # Quadratic Objective function - m.obj = pe.Objective(expr=m.x1 * m.x1 + m.x2 * m.x2, sense=pe.minimize) + m.obj = pyo.Objective(expr=m.x1 * m.x1 + m.x2 * m.x2, sense=pyo.minimize) - m.con1 = pe.Constraint(expr=m.x1 >= 1) - m.con2 = pe.Constraint(expr=m.x2 >= 1) + m.con1 = pyo.Constraint(expr=m.x1 >= 1) + m.con2 = pyo.Constraint(expr=m.x2 >= 1) results = opt.solve(m) self.assertAlmostEqual(m.x1.value, 1, places=5) @@ -130,19 +130,19 @@ def test_qp1(self): def test_qp2(self): # test issue #3381 - m = pe.ConcreteModel() + m = pyo.ConcreteModel() - m.x1 = pe.Var(name='x1', domain=pe.Reals) - m.x2 = pe.Var(name='x2', domain=pe.Reals) + m.x1 = pyo.Var(name='x1', domain=pyo.Reals) + m.x2 = pyo.Var(name='x2', domain=pyo.Reals) - m.p = pe.Param(initialize=1, mutable=True) + m.p = pyo.Param(initialize=1, mutable=True) - m.obj = pe.Objective( - expr=m.p * m.x1 * m.x1 + m.x2 * m.x2 - m.x1 * m.x2, sense=pe.minimize + m.obj = pyo.Objective( + expr=m.p * m.x1 * m.x1 + m.x2 * m.x2 - m.x1 * m.x2, sense=pyo.minimize ) - m.con1 = pe.Constraint(expr=m.x1 >= 1) - m.con2 = pe.Constraint(expr=m.x2 >= 1) + m.con1 = pyo.Constraint(expr=m.x1 >= 1) + m.con2 = pyo.Constraint(expr=m.x2 >= 1) opt.set_instance(m) results = opt.solve(m) From 04eb8474543e2f0ed9045277abfe0b74a193a9e2 Mon Sep 17 00:00:00 2001 From: Arnaud Baguet Date: Mon, 14 Apr 2025 18:27:52 -0400 Subject: [PATCH 13/21] rename q_value to hessian_value --- pyomo/contrib/solver/solvers/highs.py | 22 +++++++++++----------- 1 file changed, 11 insertions(+), 11 deletions(-) diff --git a/pyomo/contrib/solver/solvers/highs.py b/pyomo/contrib/solver/solvers/highs.py index c554682ca9a..4c5814752d0 100644 --- a/pyomo/contrib/solver/solvers/highs.py +++ b/pyomo/contrib/solver/solvers/highs.py @@ -172,9 +172,9 @@ def _build_and_pass_hessian(self): dim = self.highs.getNumCol() # Build CSC format for the lower triangular part - q_value = [] - q_index = [] - q_start = [0] * dim + hessian_value = [] + hessian_index = [] + hessian_start = [0] * dim sorted_entries = sorted( self.quad_coef_dict.items(), key=lambda x: (x[0][1], x[0][0]) @@ -185,24 +185,24 @@ def _build_and_pass_hessian(self): while col > last_col: last_col += 1 if last_col < dim: - q_start[last_col] = len(q_value) + hessian_start[last_col] = len(hessian_value) # Add the entry - q_index.append(row) - q_value.append(val) + hessian_index.append(row) + hessian_value.append(val) while last_col < dim - 1: last_col += 1 - q_start[last_col] = len(q_value) + hessian_start[last_col] = len(hessian_value) - nnz = len(q_value) + nnz = len(hessian_value) status = self.highs.passHessian( dim, nnz, highspy.HessianFormat.kTriangular, - np.array(q_start, dtype=np.int32), - np.array(q_index, dtype=np.int32), - np.array(q_value, dtype=np.double), + np.array(hessian_start, dtype=np.int32), + np.array(hessian_index, dtype=np.int32), + np.array(hessian_value, dtype=np.double), ) if status != highspy.HighsStatus.kOk: From c1422b70a0123d8d7cd4c16d0ca20d6ed968e919 Mon Sep 17 00:00:00 2001 From: Arnaud Baguet Date: Tue, 24 Jun 2025 17:22:43 -0400 Subject: [PATCH 14/21] move qp tests to test_solvers --- .../solver/tests/solvers/test_highs.py | 47 ------------------- .../solver/tests/solvers/test_solvers.py | 41 +++++++++++++++- 2 files changed, 40 insertions(+), 48 deletions(-) diff --git a/pyomo/contrib/solver/tests/solvers/test_highs.py b/pyomo/contrib/solver/tests/solvers/test_highs.py index ef557fbdaed..f59a0bfa42d 100644 --- a/pyomo/contrib/solver/tests/solvers/test_highs.py +++ b/pyomo/contrib/solver/tests/solvers/test_highs.py @@ -109,50 +109,3 @@ def test_fix_and_unfix(self): self.assertAlmostEqual(m.fx.value, 1, places=5) self.assertAlmostEqual(m.fy.value, 0, places=5) self.assertAlmostEqual(r.objective_bound, 0.5, places=5) - - def test_qp1(self): - # test issue #3381 - m = pyo.ConcreteModel() - - m.x1 = pyo.Var(name='x1', domain=pyo.Reals) - m.x2 = pyo.Var(name='x2', domain=pyo.Reals) - - # Quadratic Objective function - m.obj = pyo.Objective(expr=m.x1 * m.x1 + m.x2 * m.x2, sense=pyo.minimize) - - m.con1 = pyo.Constraint(expr=m.x1 >= 1) - m.con2 = pyo.Constraint(expr=m.x2 >= 1) - - results = opt.solve(m) - self.assertAlmostEqual(m.x1.value, 1, places=5) - self.assertAlmostEqual(m.x2.value, 1, places=5) - self.assertEqual(results.objective_bound, 2) - - def test_qp2(self): - # test issue #3381 - m = pyo.ConcreteModel() - - m.x1 = pyo.Var(name='x1', domain=pyo.Reals) - m.x2 = pyo.Var(name='x2', domain=pyo.Reals) - - m.p = pyo.Param(initialize=1, mutable=True) - - m.obj = pyo.Objective( - expr=m.p * m.x1 * m.x1 + m.x2 * m.x2 - m.x1 * m.x2, sense=pyo.minimize - ) - - m.con1 = pyo.Constraint(expr=m.x1 >= 1) - m.con2 = pyo.Constraint(expr=m.x2 >= 1) - - opt.set_instance(m) - results = opt.solve(m) - self.assertAlmostEqual(m.x1.value, 1, places=5) - self.assertAlmostEqual(m.x2.value, 1, places=5) - self.assertEqual(results.objective_bound, 1) - - m.p.value = 2.0 - opt.update_parameters() - results = opt.solve(m) - self.assertAlmostEqual(m.x1.value, 1, places=5) - self.assertAlmostEqual(m.x2.value, 1, places=5) - self.assertEqual(results.objective_bound, 2) diff --git a/pyomo/contrib/solver/tests/solvers/test_solvers.py b/pyomo/contrib/solver/tests/solvers/test_solvers.py index 808c6675529..b75f562d21c 100644 --- a/pyomo/contrib/solver/tests/solvers/test_solvers.py +++ b/pyomo/contrib/solver/tests/solvers/test_solvers.py @@ -58,6 +58,7 @@ ] nlp_solvers = [('ipopt', Ipopt)] qcp_solvers = [('gurobi', GurobiPersistent), ('ipopt', Ipopt)] +qp_solvers = qcp_solvers + [("highs", Highs)] miqcqp_solvers = [('gurobi', GurobiPersistent)] nl_solvers = [('ipopt', Ipopt)] nl_solvers_set = {i[0] for i in nl_solvers} @@ -1097,7 +1098,7 @@ def test_mutable_quadratic_coefficient( self.assertAlmostEqual(m.y.value, 0.0869525991355825, 4) @parameterized.expand(input=_load_tests(qcp_solvers)) - def test_mutable_quadratic_objective( + def test_mutable_quadratic_objective_qcp( self, name: str, opt_class: Type[SolverBase], use_presolve: bool ): opt: SolverBase = opt_class() @@ -1128,6 +1129,44 @@ def test_mutable_quadratic_objective( self.assertAlmostEqual(m.x.value, 0.6962249634573562, 4) self.assertAlmostEqual(m.y.value, 0.09227926676152151, 4) + @parameterized.expand(input=_load_tests(qp_solvers)) + def test_mutable_quadratic_objective_qp( + self, name: str, opt_class: Type[SolverBase], use_presolve: bool + ): + opt: SolverBase = opt_class() + if not opt.available(): + raise unittest.SkipTest(f'Solver {opt.name} not available.') + if any(name.startswith(i) for i in nl_solvers_set): + if use_presolve: + opt.config.writer_config.linear_presolve = True + else: + opt.config.writer_config.linear_presolve = False + # test issue #3381 + m = pyo.ConcreteModel() + + m.x1 = pyo.Var(name='x1', domain=pyo.Reals) + m.x2 = pyo.Var(name='x2', domain=pyo.Reals) + + m.p = pyo.Param(initialize=1, mutable=True) + + m.obj = pyo.Objective( + expr=m.p * m.x1 * m.x1 + m.x2 * m.x2 - m.x1 * m.x2, sense=pyo.minimize + ) + + m.con1 = pyo.Constraint(expr=m.x1 >= 1) + m.con2 = pyo.Constraint(expr=m.x2 >= 1) + + results = opt.solve(m) + self.assertAlmostEqual(m.x1.value, 1, places=4) + self.assertAlmostEqual(m.x2.value, 1, places=4) + self.assertAlmostEqual(results.incumbent_objective, 1, 4) + + m.p.value = 2.0 + results = opt.solve(m) + self.assertAlmostEqual(m.x1.value, 1, places=4) + self.assertAlmostEqual(m.x2.value, 1, places=4) + self.assertAlmostEqual(results.incumbent_objective, 2, 4) + @parameterized.expand(input=_load_tests(all_solvers)) def test_fixed_vars( self, name: str, opt_class: Type[SolverBase], use_presolve: bool From be02bba263899528188df2648bc74bd3c86caaec Mon Sep 17 00:00:00 2001 From: Arnaud Baguet Date: Tue, 24 Jun 2025 17:41:30 -0400 Subject: [PATCH 15/21] delete variable x in test --- pyomo/contrib/solver/tests/solvers/test_solvers.py | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/pyomo/contrib/solver/tests/solvers/test_solvers.py b/pyomo/contrib/solver/tests/solvers/test_solvers.py index b75f562d21c..de0dae62aea 100644 --- a/pyomo/contrib/solver/tests/solvers/test_solvers.py +++ b/pyomo/contrib/solver/tests/solvers/test_solvers.py @@ -1167,6 +1167,11 @@ def test_mutable_quadratic_objective_qp( self.assertAlmostEqual(m.x2.value, 1, places=4) self.assertAlmostEqual(results.incumbent_objective, 2, 4) + del m.x1 + results = opt.solve(m) + self.assertAlmostEqual(m.x2.value, 1, places=4) + self.assertAlmostEqual(results.incumbent_objective, 1, 4) + @parameterized.expand(input=_load_tests(all_solvers)) def test_fixed_vars( self, name: str, opt_class: Type[SolverBase], use_presolve: bool From 582358385097bbd26b1e36df909308bf1ffb4a0c Mon Sep 17 00:00:00 2001 From: Arnaud Baguet Date: Tue, 24 Jun 2025 18:45:24 -0400 Subject: [PATCH 16/21] addressing michael's comments MutableQuadraticCoefficient --- pyomo/contrib/solver/solvers/highs.py | 79 +++++++++---------- .../solver/tests/solvers/test_solvers.py | 4 +- 2 files changed, 39 insertions(+), 44 deletions(-) diff --git a/pyomo/contrib/solver/solvers/highs.py b/pyomo/contrib/solver/solvers/highs.py index 1f505112078..b6ff0a8a007 100644 --- a/pyomo/contrib/solver/solvers/highs.py +++ b/pyomo/contrib/solver/solvers/highs.py @@ -97,39 +97,36 @@ def update(self): class _MutableQuadraticCoefficient: - def __init__(self, expr, row_idx, col_idx): + def __init__(self, expr, v1_id, v2_id): self.expr = expr - self.row_idx = row_idx - self.col_idx = col_idx + self.v1_id = v1_id + self.v2_id = v2_id class _MutableObjective: - def __init__(self, highs, constant, linear_coefs, quadratic_coefs): + def __init__( + self, + highs, + constant, + linear_coefs, + quadratic_coefs, + pyomo_var_to_solver_var_map, + ): self.highs = highs self.constant = constant self.linear_coefs = linear_coefs self.quadratic_coefs = quadratic_coefs - self.last_quadratic_coef_values = [value(i.expr) for i in self.quadratic_coefs] + self._pyomo_var_to_solver_var_map = pyomo_var_to_solver_var_map # Store the quadratic coefficients in dictionary format - self.quad_coef_dict = {} - self._initialize_quad_coef_dict() + self._initialize_quad_coef_dicts() # Flag to force first update of quadratic coefficients self._first_update = True - def _initialize_quad_coef_dict(self): + def _initialize_quad_coef_dicts(self): + self.quad_coef_dict = {} for coef in self.quadratic_coefs: - v1_ndx = coef.row_idx - v2_ndx = coef.col_idx - # Ensure we're storing the lower triangular part - row = max(v1_ndx, v2_ndx) - col = min(v1_ndx, v2_ndx) - - coef_val = value(coef.expr) - # Adjust for diagonal elements - if v1_ndx == v2_ndx: - coef_val *= 2.0 - - self.quad_coef_dict[(row, col)] = coef_val + self.quad_coef_dict[(coef.v1_id, coef.v2_id)] = value(coef.expr) + self.previous_quad_coef_dict = self.quad_coef_dict.copy() def update(self): """ @@ -141,23 +138,13 @@ def update(self): for coef in self.linear_coefs: coef.update() - for ndx, coef in enumerate(self.quadratic_coefs): + for coef in self.quadratic_coefs: current_val = value(coef.expr) - if current_val != self.last_quadratic_coef_values[ndx]: + previous_val = self.previous_quad_coef_dict.get((coef.v1_id, coef.v2_id)) + if previous_val is not None and current_val != previous_val: needs_quadratic_update = True - - v1_ndx = coef.row_idx - v2_ndx = coef.col_idx - row = max(v1_ndx, v2_ndx) - col = min(v1_ndx, v2_ndx) - - # Adjust the diagonal to match Highs' expected format - if v1_ndx == v2_ndx: - current_val *= 2.0 - - self.quad_coef_dict[(row, col)] = current_val - - self.last_quadratic_coef_values[ndx] = current_val + self.quad_coef_dict[(coef.v1_id, coef.v2_id)] = current_val + self.previous_quad_coef_dict[(coef.v1_id, coef.v2_id)] = current_val # If anything changed, rebuild and pass the Hessian if needs_quadratic_update: @@ -176,8 +163,20 @@ def _build_and_pass_hessian(self): hessian_index = [] hessian_start = [0] * dim + quad_coef_idx_dict = {} + for (v1_id, v2_id), coef in self.quad_coef_dict.items(): + v1_ndx = self._pyomo_var_to_solver_var_map[v1_id] + v2_ndx = self._pyomo_var_to_solver_var_map[v2_id] + # Ensure we're storing the lower triangular part + row = max(v1_ndx, v2_ndx) + col = min(v1_ndx, v2_ndx) + # Adjust the diagonal to match Highs' expected format + if v1_ndx == v2_ndx: + coef *= 2.0 + quad_coef_idx_dict[(row, col)] = coef + sorted_entries = sorted( - self.quad_coef_dict.items(), key=lambda x: (x[0][1], x[0][0]) + quad_coef_idx_dict.items(), key=lambda x: (x[0][1], x[0][0]) ) last_col = -1 @@ -645,16 +644,11 @@ def _set_objective(self, obj): if repn.quadratic_vars and len(repn.quadratic_vars) > 0: for ndx, (v1, v2) in enumerate(repn.quadratic_vars): - v1_id = id(v1) - v2_id = id(v2) - v1_ndx = self._pyomo_var_to_solver_var_map[v1_id] - v2_ndx = self._pyomo_var_to_solver_var_map[v2_id] - coef = repn.quadratic_coefs[ndx] mutable_quadratic_coefficients.append( _MutableQuadraticCoefficient( - expr=coef, row_idx=v1_ndx, col_idx=v2_ndx + expr=coef, v1_id=id(v1), v2_id=id(v2) ) ) @@ -665,6 +659,7 @@ def _set_objective(self, obj): mutable_constant, mutable_linear_coefficients, mutable_quadratic_coefficients, + self._pyomo_var_to_solver_var_map, ) self._mutable_objective.update() diff --git a/pyomo/contrib/solver/tests/solvers/test_solvers.py b/pyomo/contrib/solver/tests/solvers/test_solvers.py index de0dae62aea..99e5ec3209c 100644 --- a/pyomo/contrib/solver/tests/solvers/test_solvers.py +++ b/pyomo/contrib/solver/tests/solvers/test_solvers.py @@ -1167,10 +1167,10 @@ def test_mutable_quadratic_objective_qp( self.assertAlmostEqual(m.x2.value, 1, places=4) self.assertAlmostEqual(results.incumbent_objective, 2, 4) - del m.x1 + # del m.x1 results = opt.solve(m) self.assertAlmostEqual(m.x2.value, 1, places=4) - self.assertAlmostEqual(results.incumbent_objective, 1, 4) + self.assertAlmostEqual(results.incumbent_objective, 2, 4) @parameterized.expand(input=_load_tests(all_solvers)) def test_fixed_vars( From 3fb1476671a0e89596c581f777f2d27ca6ca88e3 Mon Sep 17 00:00:00 2001 From: Arnaud Baguet Date: Tue, 24 Jun 2025 18:54:55 -0400 Subject: [PATCH 17/21] remove commented bit in test --- pyomo/contrib/solver/tests/solvers/test_solvers.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyomo/contrib/solver/tests/solvers/test_solvers.py b/pyomo/contrib/solver/tests/solvers/test_solvers.py index 99e5ec3209c..bacb6e7204e 100644 --- a/pyomo/contrib/solver/tests/solvers/test_solvers.py +++ b/pyomo/contrib/solver/tests/solvers/test_solvers.py @@ -1167,7 +1167,7 @@ def test_mutable_quadratic_objective_qp( self.assertAlmostEqual(m.x2.value, 1, places=4) self.assertAlmostEqual(results.incumbent_objective, 2, 4) - # del m.x1 + del m.x1 results = opt.solve(m) self.assertAlmostEqual(m.x2.value, 1, places=4) self.assertAlmostEqual(results.incumbent_objective, 2, 4) From a5317fad27ded2a2f58d2e088b40526414bd10ab Mon Sep 17 00:00:00 2001 From: Arnaud Baguet Date: Thu, 3 Jul 2025 19:24:02 -0400 Subject: [PATCH 18/21] add back TerminationCondition --- pyomo/contrib/appsi/solvers/tests/test_highs_persistent.py | 1 + 1 file changed, 1 insertion(+) diff --git a/pyomo/contrib/appsi/solvers/tests/test_highs_persistent.py b/pyomo/contrib/appsi/solvers/tests/test_highs_persistent.py index 5814a8b3eac..9c8a8b7c0f4 100644 --- a/pyomo/contrib/appsi/solvers/tests/test_highs_persistent.py +++ b/pyomo/contrib/appsi/solvers/tests/test_highs_persistent.py @@ -18,6 +18,7 @@ from pyomo.common.log import LoggingIntercept from pyomo.common.tee import capture_output from pyomo.contrib.appsi.solvers.highs import Highs +from pyomo.contrib.appsi.base import TerminationCondition from pyomo.contrib.solver.tests.solvers import instances From c4b08aa15f073726d88d7604b92d14f6c7d78300 Mon Sep 17 00:00:00 2001 From: Michael Bynum Date: Tue, 29 Jul 2025 05:54:27 -0600 Subject: [PATCH 19/21] update mutable qp test --- .../solver/tests/solvers/test_solvers.py | 72 ++++++++++++++----- 1 file changed, 53 insertions(+), 19 deletions(-) diff --git a/pyomo/contrib/solver/tests/solvers/test_solvers.py b/pyomo/contrib/solver/tests/solvers/test_solvers.py index 971fb894f0c..02d6c1fa5fc 100644 --- a/pyomo/contrib/solver/tests/solvers/test_solvers.py +++ b/pyomo/contrib/solver/tests/solvers/test_solvers.py @@ -53,14 +53,14 @@ ('highs', Highs), ] mip_solvers = [ - ('gurobi', GurobiPersistent), + ('gurobi_persistent', GurobiPersistent), ('gurobi_direct', GurobiDirect), ('highs', Highs), ] nlp_solvers = [('ipopt', Ipopt)] -qcp_solvers = [('gurobi', GurobiPersistent), ('ipopt', Ipopt)] +qcp_solvers = [('gurobi_persistent', GurobiPersistent), ('ipopt', Ipopt)] qp_solvers = qcp_solvers + [("highs", Highs)] -miqcqp_solvers = [('gurobi', GurobiPersistent)] +miqcqp_solvers = [('gurobi_persistent', GurobiPersistent)] nl_solvers = [('ipopt', Ipopt)] nl_solvers_set = {i[0] for i in nl_solvers} @@ -1145,33 +1145,67 @@ def test_mutable_quadratic_objective_qp( # test issue #3381 m = pyo.ConcreteModel() - m.x1 = pyo.Var(name='x1', domain=pyo.Reals) - m.x2 = pyo.Var(name='x2', domain=pyo.Reals) + m.x1 = pyo.Var() + m.x2 = pyo.Var() - m.p = pyo.Param(initialize=1, mutable=True) + m.p1 = pyo.Param(initialize=1, mutable=True) + m.p2 = pyo.Param(initialize=1, mutable=True) + m.p3 = pyo.Param(initialize=4, mutable=True) m.obj = pyo.Objective( - expr=m.p * m.x1 * m.x1 + m.x2 * m.x2 - m.x1 * m.x2, sense=pyo.minimize + expr=m.p1*(m.x1 - 1)**2 + m.p2*(m.x2 - 6)**2 - m.p3*m.x2 ) - m.con1 = pyo.Constraint(expr=m.x1 >= 1) - m.con2 = pyo.Constraint(expr=m.x2 >= 1) + m.con = pyo.Constraint(expr=m.x1 >= m.x2) + + results = opt.solve(m) + self.assertAlmostEqual(m.x1.value, 4.5, places=4) + self.assertAlmostEqual(m.x2.value, 4.5, places=4) + self.assertAlmostEqual(results.incumbent_objective, -3.5, 4) + m.p2.value = 2.0 results = opt.solve(m) - self.assertAlmostEqual(m.x1.value, 1, places=4) - self.assertAlmostEqual(m.x2.value, 1, places=4) - self.assertAlmostEqual(results.incumbent_objective, 1, 4) + self.assertAlmostEqual(m.x1.value, 5, places=4) + self.assertAlmostEqual(m.x2.value, 5, places=4) + self.assertAlmostEqual(results.incumbent_objective, -2, 4) + + m.x3 = pyo.Var() + del m.obj + m.obj = pyo.Objective( + expr=m.p2*(m.x2 - 6)**2 - m.p3*m.x2 + m.p1*(m.x3 - 1)**2 + ) + m.con2 = pyo.Constraint(expr=m.x3 >= m.x1) - m.p.value = 2.0 results = opt.solve(m) - self.assertAlmostEqual(m.x1.value, 1, places=4) - self.assertAlmostEqual(m.x2.value, 1, places=4) - self.assertAlmostEqual(results.incumbent_objective, 2, 4) + self.assertAlmostEqual(m.x1.value, 5, places=4) + self.assertAlmostEqual(m.x2.value, 5, places=4) + self.assertAlmostEqual(m.x3.value, 5, places=4) + self.assertAlmostEqual(results.incumbent_objective, -2, 4) + + if opt_class is Highs: + """ + This assertions is not important by iteself. + We just need it to make sure that removing the + variable below is actually testing what we think + (which is that the mutable quadratic coefficients + work correctly even when the column changes) + """ + self.assertIn(opt._pyomo_var_to_solver_var_map[id(m.x1)], {0, 1}) + self.assertIn(opt._pyomo_var_to_solver_var_map[id(m.x2)], {0, 1}) + self.assertEqual(opt._pyomo_var_to_solver_var_map[id(m.x3)], 2) + + del m.con + del m.con2 + m.p1.value = 2 + m.con = pyo.Constraint(expr=m.x3 >= m.x2) - del m.x1 results = opt.solve(m) - self.assertAlmostEqual(m.x2.value, 1, places=4) - self.assertAlmostEqual(results.incumbent_objective, 2, 4) + self.assertAlmostEqual(m.x2.value, 4, places=4) + self.assertAlmostEqual(m.x3.value, 4, places=4) + self.assertAlmostEqual(results.incumbent_objective, 10, 4) + + if opt_class is Highs: + self.assertIn(opt._pyomo_var_to_solver_var_map[id(m.x3)], {0, 1}) @parameterized.expand(input=_load_tests(all_solvers)) def test_fixed_vars( From de53637454c4d8c18e7ae08f0363b3f5e58ad98b Mon Sep 17 00:00:00 2001 From: Michael Bynum Date: Tue, 29 Jul 2025 05:59:08 -0600 Subject: [PATCH 20/21] run black --- .../contrib/solver/tests/solvers/test_solvers.py | 16 +++++++--------- 1 file changed, 7 insertions(+), 9 deletions(-) diff --git a/pyomo/contrib/solver/tests/solvers/test_solvers.py b/pyomo/contrib/solver/tests/solvers/test_solvers.py index 02d6c1fa5fc..a1fa41f2415 100644 --- a/pyomo/contrib/solver/tests/solvers/test_solvers.py +++ b/pyomo/contrib/solver/tests/solvers/test_solvers.py @@ -1153,7 +1153,7 @@ def test_mutable_quadratic_objective_qp( m.p3 = pyo.Param(initialize=4, mutable=True) m.obj = pyo.Objective( - expr=m.p1*(m.x1 - 1)**2 + m.p2*(m.x2 - 6)**2 - m.p3*m.x2 + expr=m.p1 * (m.x1 - 1) ** 2 + m.p2 * (m.x2 - 6) ** 2 - m.p3 * m.x2 ) m.con = pyo.Constraint(expr=m.x1 >= m.x2) @@ -1172,7 +1172,7 @@ def test_mutable_quadratic_objective_qp( m.x3 = pyo.Var() del m.obj m.obj = pyo.Objective( - expr=m.p2*(m.x2 - 6)**2 - m.p3*m.x2 + m.p1*(m.x3 - 1)**2 + expr=m.p2 * (m.x2 - 6) ** 2 - m.p3 * m.x2 + m.p1 * (m.x3 - 1) ** 2 ) m.con2 = pyo.Constraint(expr=m.x3 >= m.x1) @@ -1183,13 +1183,11 @@ def test_mutable_quadratic_objective_qp( self.assertAlmostEqual(results.incumbent_objective, -2, 4) if opt_class is Highs: - """ - This assertions is not important by iteself. - We just need it to make sure that removing the - variable below is actually testing what we think - (which is that the mutable quadratic coefficients - work correctly even when the column changes) - """ + # This assertions is not important by iteself. + # We just need it to make sure that removing the + # variable below is actually testing what we think + # (which is that the mutable quadratic coefficients + # work correctly even when the column changes) self.assertIn(opt._pyomo_var_to_solver_var_map[id(m.x1)], {0, 1}) self.assertIn(opt._pyomo_var_to_solver_var_map[id(m.x2)], {0, 1}) self.assertEqual(opt._pyomo_var_to_solver_var_map[id(m.x3)], 2) From 0d4502500fcc8ce7a51a2bee21a25ae9aaf57ae1 Mon Sep 17 00:00:00 2001 From: Michael Bynum Date: Tue, 29 Jul 2025 06:12:55 -0600 Subject: [PATCH 21/21] typo --- pyomo/contrib/solver/tests/solvers/test_solvers.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyomo/contrib/solver/tests/solvers/test_solvers.py b/pyomo/contrib/solver/tests/solvers/test_solvers.py index a1fa41f2415..5ab36554061 100644 --- a/pyomo/contrib/solver/tests/solvers/test_solvers.py +++ b/pyomo/contrib/solver/tests/solvers/test_solvers.py @@ -1183,7 +1183,7 @@ def test_mutable_quadratic_objective_qp( self.assertAlmostEqual(results.incumbent_objective, -2, 4) if opt_class is Highs: - # This assertions is not important by iteself. + # This assertions is not important by itself. # We just need it to make sure that removing the # variable below is actually testing what we think # (which is that the mutable quadratic coefficients