diff --git a/pyomo/contrib/solver/solvers/highs.py b/pyomo/contrib/solver/solvers/highs.py
index 5a85fa8814d..2fdac4942c8 100644
--- a/pyomo/contrib/solver/solvers/highs.py
+++ b/pyomo/contrib/solver/solvers/highs.py
@@ -97,6 +97,120 @@ def update(self):
         self.highs.changeColCost(col_ndx, value(self.expr))
 
 
+class _MutableQuadraticCoefficient:
+    def __init__(self, expr, v1_id, v2_id):
+        self.expr = expr
+        self.v1_id = v1_id
+        self.v2_id = v2_id
+
+
+class _MutableObjective:
+    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._pyomo_var_to_solver_var_map = pyomo_var_to_solver_var_map
+        # Store the quadratic coefficients in dictionary format
+        self._initialize_quad_coef_dicts()
+        # Flag to force first update of quadratic coefficients
+        self._first_update = True
+
+    def _initialize_quad_coef_dicts(self):
+        self.quad_coef_dict = {}
+        for coef in self.quadratic_coefs:
+            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):
+        """
+        Update the quadratic objective expression.
+        """
+        needs_quadratic_update = self._first_update
+
+        self.constant.update()
+        for coef in self.linear_coefs:
+            coef.update()
+
+        for coef in self.quadratic_coefs:
+            current_val = value(coef.expr)
+            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
+                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:
+            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
+        hessian_value = []
+        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(
+            quad_coef_idx_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:
+                    hessian_start[last_col] = len(hessian_value)
+
+            # Add the entry
+            hessian_index.append(row)
+            hessian_value.append(val)
+
+        while last_col < dim - 1:
+            last_col += 1
+            hessian_start[last_col] = len(hessian_value)
+
+        nnz = len(hessian_value)
+        status = self.highs.passHessian(
+            dim,
+            nnz,
+            highspy.HessianFormat.kTriangular,
+            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:
+            logger.warning(
+                f"HiGHS returned non-OK status when passing Hessian: {status}"
+            )
+
+
 class _MutableObjectiveOffset:
     def __init__(self, expr, highs):
         self.expr = expr
@@ -142,7 +256,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
 
@@ -473,13 +586,14 @@ def update_parameters(self):
         self._sol = None
         if self._last_results_object is not None:
             self._last_results_object.solution_loader.invalidate()
+
         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()
+
+        self._mutable_objective.update()
 
     def _set_objective(self, obj):
         self._sol = None
@@ -488,10 +602,14 @@ def _set_objective(self, obj):
         n = len(self._pyomo_var_to_solver_var_map)
         indices = np.arange(n)
         costs = np.zeros(n, dtype=np.double)
-        self._objective_helpers = []
+
+        # 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)
+            mutable_constant = _MutableObjectiveOffset(expr=0, highs=self._solver_model)
         else:
             if obj.sense == minimize:
                 sense = highspy.ObjSense.kMinimize
@@ -501,9 +619,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 or repn.polynomial_degree() > 2:
                 raise IncompatibleModelError(
                     f'Highs interface does not support expressions of degree {repn.polynomial_degree()}'
                 )
@@ -519,17 +637,32 @@ def _set_objective(self, obj):
                         expr=coef,
                         highs=self._solver_model,
                     )
-                    self._objective_helpers.append(mutable_objective_coef)
+                    mutable_linear_coefficients.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)
+            mutable_constant = _MutableObjectiveOffset(
+                expr=repn.constant, highs=self._solver_model
+            )
+
+            if repn.quadratic_vars and len(repn.quadratic_vars) > 0:
+                for ndx, (v1, v2) in enumerate(repn.quadratic_vars):
+                    coef = repn.quadratic_coefs[ndx]
+
+                    mutable_quadratic_coefficients.append(
+                        _MutableQuadraticCoefficient(
+                            expr=coef, v1_id=id(v1), v2_id=id(v2)
+                        )
+                    )
 
         self._solver_model.changeObjectiveSense(sense)
         self._solver_model.changeColsCost(n, indices, costs)
+        self._mutable_objective = _MutableObjective(
+            self._solver_model,
+            mutable_constant,
+            mutable_linear_coefficients,
+            mutable_quadratic_coefficients,
+            self._pyomo_var_to_solver_var_map,
+        )
+        self._mutable_objective.update()
 
     def _postsolve(self):
         config = self._active_config
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..f59a0bfa42d
--- /dev/null
+++ b/pyomo/contrib/solver/tests/solvers/test_highs.py
@@ -0,0 +1,111 @@
+#  ___________________________________________________________________________
+#
+#  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 pyomo.common.unittest as unittest
+import pyomo.environ as pyo
+
+from pyomo.contrib.solver.solvers.highs import Highs
+
+opt = Highs()
+if not opt.available():
+    raise unittest.SkipTest
+
+
+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.p1 = pyo.Param(mutable=True)
+        m.p2 = pyo.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.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 = pyo.ConcreteModel()
+        m.x = pyo.Var()
+        m.y = pyo.Var()
+
+        m.p1 = pyo.Param(mutable=True)
+        m.p2 = pyo.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.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 = 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 = pyo.Objective(expr=m.fx * 0.5 + m.fy * 0.4, sense=pyo.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)
diff --git a/pyomo/contrib/solver/tests/solvers/test_solvers.py b/pyomo/contrib/solver/tests/solvers/test_solvers.py
index 9615f1cc598..5ab36554061 100644
--- a/pyomo/contrib/solver/tests/solvers/test_solvers.py
+++ b/pyomo/contrib/solver/tests/solvers/test_solvers.py
@@ -53,13 +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)]
-miqcqp_solvers = [('gurobi', GurobiPersistent)]
+qcp_solvers = [('gurobi_persistent', GurobiPersistent), ('ipopt', Ipopt)]
+qp_solvers = qcp_solvers + [("highs", Highs)]
+miqcqp_solvers = [('gurobi_persistent', GurobiPersistent)]
 nl_solvers = [('ipopt', Ipopt)]
 nl_solvers_set = {i[0] for i in nl_solvers}
 
@@ -1098,7 +1099,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()
@@ -1129,6 +1130,81 @@ 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()
+        m.x2 = pyo.Var()
+
+        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.p1 * (m.x1 - 1) ** 2 + m.p2 * (m.x2 - 6) ** 2 - m.p3 * m.x2
+        )
+
+        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, 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)
+
+        results = opt.solve(m)
+        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 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
+            # 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)
+
+        results = opt.solve(m)
+        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(
         self, name: str, opt_class: Type[SolverBase], use_presolve: bool