Feature/dual-value-retrieval

linopy/constants.py

@@ -209,6 +209,7 @@ class Solution:
     primal: pd.Series = field(default_factory=_pd_series_float)
     dual: pd.Series = field(default_factory=_pd_series_float)
     objective: float = field(default=np.nan)
+    qc_dual: pd.Series = field(default_factory=_pd_series_float)
 
 
 @dataclass

linopy/model.py

@@ -1461,6 +1461,20 @@ def solve(
                     vals = dual.reindex(idx).values.reshape(con.labels.shape)
                 con.dual = xr.DataArray(vals, con.labels.coords)
 
+        # Assign quadratic constraint dual values
+        if not result.solution.qc_dual.empty:
+            qc_dual = result.solution.qc_dual.copy()
+            qc_dual = set_int_index(qc_dual)
+            qc_dual.loc[-1] = nan
+
+            for name, qcon in self.quadratic_constraints.items():
+                idx = np.ravel(qcon.labels)
+                try:
+                    vals = qc_dual[idx].values.reshape(qcon.labels.shape)
+                except KeyError:
+                    vals = qc_dual.reindex(idx).values.reshape(qcon.labels.shape)
+                qcon.dual = xr.DataArray(vals, qcon.labels.coords)
+
         return result.status.status.value, result.status.termination_condition.value
 
     def compute_infeasibilities(self) -> list[int]:

linopy/solvers.py

@@ -1195,7 +1195,19 @@ def get_solver_solution() -> Solution:
                 logger.warning("Dual values of MILP couldn't be parsed")
                 dual = pd.Series(dtype=float)
 
-            return Solution(sol, dual, objective)
+            # Retrieve quadratic constraint dual values
+            qc_dual = pd.Series(dtype=float)
+            try:
+                qcs = m.getQConstrs()
+                if qcs:
+                    qc_dual = pd.Series(
+                        {qc.QCName: qc.QCPi for qc in qcs}, dtype=float
+                    )
+            except (AttributeError, gurobipy.GurobiError):
+                # QCPi not available (non-convex or QCPDual=0)
+                pass
+
+            return Solution(sol, dual, objective, qc_dual)
 
         solution = self.safe_get_solution(status=status, func=get_solver_solution)
         solution = solution = maybe_adjust_objective_sign(solution, io_api, sense)

test/test_quadratic_constraint.py

@@ -469,7 +469,6 @@ def test_empty_quadratic_constraints(self, m: Model) -> None:
         assert m.matrices.qc_linear is None
 
 
-
 class TestNetCDFSerialization:
     """Tests for netCDF serialization of quadratic constraints."""
 
@@ -520,3 +519,47 @@ def test_netcdf_roundtrip_multidimensional(self) -> None:
         assert len(m2.quadratic_constraints["circles"].labels.values.ravel()) == 3
 
         fn.unlink()
+
+
+class TestDualValues:
+    """Tests for dual value retrieval for quadratic constraints."""
+
+    def test_qc_dual_with_gurobi(
+        self, m: Model, x: linopy.Variable, y: linopy.Variable
+    ) -> None:
+        """Test that dual values can be retrieved for convex QC with Gurobi."""
+        if "gurobi" not in linopy.available_solvers:
+            pytest.skip("Gurobi not available")
+
+        m.add_constraints(x + y <= 8, name="budget")
+        m.add_quadratic_constraints(x * x + y * y, "<=", 25, name="circle")
+        m.add_objective(x + 2 * y, sense="max")
+
+        # Solve with QCPDual enabled
+        m.solve(solver_name="gurobi", QCPDual=1)
+
+        # Check dual values exist
+        dual = m.quadratic_constraints["circle"].dual
+        assert dual is not None
+        assert not dual.isnull().all()
+        # Dual should be positive for binding <= constraint
+        assert float(dual.values) > 0
+
+    def test_qc_dual_multidimensional(self) -> None:
+        """Test dual values for multi-dimensional quadratic constraints."""
+        if "gurobi" not in linopy.available_solvers:
+            pytest.skip("Gurobi not available")
+
+        m = Model()
+        x = m.add_variables(lower=0, coords=[range(3)], name="x")
+        y = m.add_variables(lower=0, coords=[range(3)], name="y")
+
+        m.add_constraints(x + y <= 8, name="budget")
+        m.add_quadratic_constraints(x * x + y * y, "<=", 25, name="circles")
+        m.add_objective((x + 2 * y).sum(), sense="max")
+
+        m.solve(solver_name="gurobi", QCPDual=1)
+
+        dual = m.quadratic_constraints["circles"].dual
+        assert dual.shape == (3,)
+        assert not dual.isnull().all()