diff --git a/linopy/io.py b/linopy/io.py index 9ab0d678..069eca44 100644 --- a/linopy/io.py +++ b/linopy/io.py @@ -745,6 +745,61 @@ def to_mosek( task.putobjsense(mosek.objsense.maximize) else: task.putobjsense(mosek.objsense.minimize) + + ## Quadratic Constraints + if len(m.quadratic_constraints): + # Get the number of quadratic constraints + n_qcons = len(M.qclabels) + + # Append quadratic constraints to the task + # In MOSEK, quadratic constraints are added to regular constraints + # with quadratic terms via putqconk + qc_start_idx = m.ncons # Start after linear constraints + task.appendcons(n_qcons) + + # Get matrices for QC + Qc_list = M.Qc + qc_linear = M.qc_linear + qc_sense = M.qc_sense + qc_rhs = M.qc_rhs + + for i, label in enumerate(M.qclabels): + con_idx = qc_start_idx + i + + # Set constraint name + task.putconname(con_idx, f"qc{label}") + + # Set constraint bound based on sense + sense = qc_sense[i] + rhs = qc_rhs[i] + if sense == "<=": + bk = mosek.boundkey.up + bl, bu = 0.0, rhs + elif sense == ">=": + bk = mosek.boundkey.lo + bl, bu = rhs, 0.0 + else: # "=" + bk = mosek.boundkey.fx + bl, bu = rhs, rhs + task.putconbound(con_idx, bk, bl, bu) + + # Add linear terms if any + if qc_linear is not None: + row = qc_linear.getrow(i).tocoo() + if row.nnz > 0: + task.putarow(con_idx, list(row.col), list(row.data)) + + # Add quadratic terms + # MOSEK expects lower triangular part only + Q = Qc_list[i] + if Q.nnz > 0: + # Get lower triangular part (MOSEK requirement) + Q_lower = tril(Q).tocoo() + # MOSEK uses 0.5 * x'Qx convention, but our Q is already doubled + # So we need to divide by 2 + task.putqconk(con_idx, list(Q_lower.row), list(Q_lower.col), + list(Q_lower.data / 2)) + return task diff --git a/linopy/matrices.py b/linopy/matrices.py index 27f9e196..838971d1 100644 --- a/linopy/matrices.py +++ b/linopy/matrices.py @@ -12,7 +12,6 @@ import numpy as np import pandas as pd -import scipy.sparse from numpy import ndarray from pandas.core.indexes.base import Index from pandas.core.series import Series @@ -330,6 +329,10 @@ def qc_linear(self) -> csc_matrix | None: con_map = pd.Series(index=labels, data=np.arange(n_cons)) # Filter to linear terms only + # Note: 'vars' column may not exist if there are no linear terms + if "vars" not in df.columns: + return csc_matrix((n_cons, n_vars)) + linear_df = df[~df["is_quadratic"] & (df["vars"] >= 0)] if linear_df.empty: diff --git a/test/test_quadratic_constraint.py b/test/test_quadratic_constraint.py index 7073a1a3..35e6c296 100644 --- a/test/test_quadratic_constraint.py +++ b/test/test_quadratic_constraint.py @@ -301,31 +301,82 @@ def test_lp_file_with_multidimensional_constraint(self) -> None: class TestSolverValidation: """Tests for solver validation with quadratic constraints.""" - def test_solver_validation_warning( + def test_highs_rejects_quadratic_constraints( self, m: Model, x: linopy.Variable, y: linopy.Variable ) -> None: - """Test that unsupported solvers raise an error for QC problems.""" - m.add_objective(x + y) + """Test that HiGHS raises an error for quadratic constraints.""" + if "highs" not in linopy.available_solvers: + pytest.skip("HiGHS not available") + + m.add_objective(x + y) # Linear objective m.add_quadratic_constraints(x * x + y * y, "<=", 100, name="qc1") - # HiGHS doesn't support quadratic constraints + # HiGHS supports QP (quadratic objective) but not QCP (quadratic constraints) from linopy.solvers import quadratic_constraint_solvers if "highs" not in quadratic_constraint_solvers: - with pytest.raises(ValueError, match="does not support quadratic"): + with pytest.raises(ValueError, match="does not support quadratic constraints"): m.solve(solver_name="highs") - def test_to_highspy_raises_on_quadratic_constraints( + def test_highs_accepts_quadratic_objective( self, m: Model, x: linopy.Variable, y: linopy.Variable ) -> None: - """Test that to_highspy raises ValueError for quadratic constraints.""" - m.add_objective(x + y) - m.add_quadratic_constraints(x * x + y * y, "<=", 100, name="qc1") + """Test that HiGHS accepts quadratic objectives (but not QC).""" + if "highs" not in linopy.available_solvers: + pytest.skip("HiGHS not available") + + # Quadratic objective, no quadratic constraints + m.add_objective(x * x + y * y) + m.add_constraints(x + y >= 1, name="c1") + + # This should work - HiGHS supports QP + status, _ = m.solve(solver_name="highs") + assert status == "ok" + + def test_supported_solver_accepts_quadratic_constraints( + self, m: Model, x: linopy.Variable, y: linopy.Variable + ) -> None: + """Test that supported solvers accept quadratic constraints.""" + from linopy.solvers import quadratic_constraint_solvers - from linopy.io import to_highspy + # Find a solver that supports QC + available_qc_solvers = [ + s for s in quadratic_constraint_solvers if s in linopy.available_solvers + ] + if not available_qc_solvers: + pytest.skip("No QC-supporting solver available") + + solver = available_qc_solvers[0] + + m.add_objective(x + y, sense="max") + m.add_constraints(x + y <= 10, name="budget") + m.add_quadratic_constraints(x * x + y * y, "<=", 25, name="circle") + + # Should succeed + status, _ = m.solve(solver_name=solver) + assert status == "ok" + + def test_is_quadratic_with_qc_only( + self, m: Model, x: linopy.Variable, y: linopy.Variable + ) -> None: + """Test that is_quadratic is True when only QC are present.""" + m.add_objective(x + y) # Linear objective + m.add_quadratic_constraints(x * x, "<=", 10, name="qc") + + assert m.has_quadratic_constraints is True + assert m.objective.is_quadratic is False + assert m.is_quadratic is True # True because of QC + + def test_is_quadratic_with_quadratic_objective_only( + self, m: Model, x: linopy.Variable, y: linopy.Variable + ) -> None: + """Test that is_quadratic is True when only quadratic objective.""" + m.add_objective(x * x + y * y) # Quadratic objective + m.add_constraints(x + y <= 10, name="c") # Linear constraint - with pytest.raises(ValueError, match="HiGHS does not support quadratic"): - to_highspy(m) + assert m.has_quadratic_constraints is False + assert m.objective.is_quadratic is True + assert m.is_quadratic is True # True because of objective class TestQuadraticConstraintRepr: @@ -521,6 +572,45 @@ def test_netcdf_roundtrip_multidimensional(self) -> None: fn.unlink() +class TestMOSEKExport: + """Tests for MOSEK direct API export with quadratic constraints.""" + + def test_to_mosek_with_quadratic_constraints( + self, m: Model, x: linopy.Variable, y: linopy.Variable + ) -> None: + """Test that to_mosek works with quadratic constraints.""" + if "mosek" not in linopy.available_solvers: + pytest.skip("MOSEK 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") + + from linopy.io import to_mosek + + task = to_mosek(m) + # If we got here without error, the export worked + assert task is not None + + def test_to_mosek_multidimensional(self) -> None: + """Test MOSEK export with multi-dimensional quadratic constraints.""" + if "mosek" not in linopy.available_solvers: + pytest.skip("MOSEK 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") + + from linopy.io import to_mosek + + task = to_mosek(m) + assert task is not None + + class TestDualValues: """Tests for dual value retrieval for quadratic constraints."""