appsi_highs: allow solving quadratic problems

pyomo/contrib/appsi/solvers/highs.py

@@ -147,6 +147,20 @@ def update(self):
         ub = value(self.upper_expr)
         self.highs.changeRowBounds(row_ndx, lb, ub)
 
+class _MutableQuadraticCoefficient(object):
+    def __init__(self, pyomo_con, pyomo_var1_id, pyomo_var2_id, con_map, var_map, expr, highs):
+        self.expr = expr
+        self.highs = highs
+        self.pyomo_var1_id = pyomo_var1_id
+        self.pyomo_var2_id = pyomo_var2_id
+        self.pyomo_con = pyomo_con
+        self.con_map = con_map
+        self.var_map = var_map
+
+    def update(self):
+        row_ndx = self.con_map[self.pyomo_con]
+        col_ndx = self.var_map[self.pyomo_var_id]
+        self.highs.changeCoeff(row_ndx, col_ndx, value(self.expr))
 
 class Highs(PersistentBase, PersistentSolver):
     """
@@ -397,12 +411,8 @@ def _add_constraints(self, cons: List[ConstraintData]):
 
         for con in cons:
             repn = generate_standard_repn(
-                con.body, quadratic=False, compute_values=False
+                con.body, quadratic=True, compute_values=False
             )
-            if repn.nonlinear_expr is not None:
-                raise DegreeError(
-                    f'Highs interface does not support expressions of degree {repn.polynomial_degree()}'
-                )
 
             starts.append(len(coef_values))
             for ndx, coef in enumerate(repn.linear_coefs):
@@ -426,6 +436,29 @@ def _add_constraints(self, cons: List[ConstraintData]):
                 var_indices.append(self._pyomo_var_to_solver_var_map[v_id])
                 coef_values.append(coef_val)
 
+            # TODO fix quadratic expression
+            for coef, v in zip(repn.quadratic_coefs, repn.quadratic_vars):
+                x, y = v
+                if not is_constant(coef):
+                    mutable_quadratic_coefficient = _MutableQuadraticCoefficient(
+                        pyomo_con=con,
+                        pyomo_var1_id=id(x),
+                        pyomo_var2_id=id(y),
+                        con_map=self._pyomo_con_to_solver_con_map,
+                        var_map=self._pyomo_var_to_solver_var_map,
+                        expr=coef,
+                        highs=self._solver_model,
+                    )
+                    if con not in self._mutable_helpers:
+                        self._mutable_helpers[con] = list()
+                    self._mutable_helpers[con].append(mutable_quadratic_coefficient)
+                    if coef_val == 0:
+                        continue
+                coef_val = value(coef)
+                var_indices.append(self._pyomo_var_to_solver_var_map[id(x)])
+                var_indices.append(self._pyomo_var_to_solver_var_map[id(y)])
+                coef_values.append(coef_val)
+
             if con.has_lb():
                 lb = con.lower - repn.constant
             else:
@@ -589,12 +622,12 @@ 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(
-                    f'Highs interface does not support expressions of degree {repn.polynomial_degree()}'
-                )
+
+            degree = repn.polynomial_degree()
+            if (degree is None) or (degree > 2):
+                raise DegreeError('Highs does not support expressions of degree {0}.'.format(degree))
 
             for coef, v in zip(repn.linear_coefs, repn.linear_vars):
                 v_id = id(v)
@@ -609,6 +642,18 @@ def _set_objective(self, obj):
                     )
                     self._objective_helpers.append(mutable_objective_coef)
 
+            # TODO fix representation of quadratic vars
+            for coef, v in zip(repn.quadratic_coefs, repn.quadratic_vars):
+                x, y = v
+                if not is_constant(coef):
+                    mutable_objective_coef = _MutableObjectiveCoefficient(
+                        pyomo_var_id=id(x),
+                        var_map=self._pyomo_var_to_solver_var_map,
+                        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(

pyomo/contrib/appsi/solvers/tests/test_persistent_solvers.py

@@ -50,6 +50,7 @@
     ('gurobi', Gurobi),
     ('ipopt', Ipopt),
     ('cplex', Cplex),
+    ('highs', Highs),
     ('maingo', MAiNGO),
 ]
 miqcqp_solvers = [('gurobi', Gurobi), ('cplex', Cplex), ('maingo', MAiNGO)]