Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
55 changes: 55 additions & 0 deletions linopy/io.py
Original file line number Diff line number Diff line change
Expand Up @@ -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))
Comment on lines +798 to +801
Copy link
Copy Markdown

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

⚠️ Potential issue | 🔴 Critical

🧩 Analysis chain

Script executed:

#!/bin/bash
# Find and examine matrices.py for Qc property
fd -t f "matrices.py" --exec wc -l {} \;

Repository: FBumann/linopy


<shell_result>
Repository: FBumann/linopy
Exit code: 0

stdout:

357 ./linopy/matrices.py
79 ./test/test_matrices.py

</shell_result>


Script executed:

#!/bin/bash
# Find matrices.py and read the Qc property
fd -t f "matrices.py" | head -5

Repository: FBumann/linopy


<shell_result>
Repository: FBumann/linopy
Exit code: 0

stdout:

linopy/matrices.py
test/test_matrices.py

</shell_result>


Script executed:

#!/bin/bash
# Read the relevant lines from io.py with context
sed -n '790,810p' linopy/io.py

Repository: FBumann/linopy


<shell_result>
Repository: FBumann/linopy
Exit code: 0

stdout:

                    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


def to_gurobipy(
    m: Model, env: Any | None = None, explicit_coordinate_names: bool = False
) -> Any:
    """
    Export the model to gurobipy.

</shell_result>


Script executed:

#!/bin/bash
# Search for Qc property definition in matrices.py
rg -n "def Qc" --type py

Repository: FBumann/linopy


<shell_result>
Repository: FBumann/linopy
Exit code: 0

stdout:

linopy/matrices.py:241:    def Qc(self) -> list[csc_matrix]:

</shell_result>


Script executed:

#!/bin/bash
# Look for Q matrix handling and MOSEK references in io.py around the target lines
sed -n '780,820p' linopy/io.py

Repository: FBumann/linopy


<shell_result>
Repository: FBumann/linopy
Exit code: 0

stdout:

                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


def to_gurobipy(
    m: Model, env: Any | None = None, explicit_coordinate_names: bool = False
) -> Any:
    """
    Export the model to gurobipy.

    This function does not write the model to intermediate files but directly
    passes it to gurobipy. Note that for large models this is not
    computationally efficient.

    Parameters
    ----------
    m : linopy.Model
    env : gurobipy.Env

</shell_result>


Script executed:

#!/bin/bash
# Search for Q_lower construction and MOSEK references
rg -n "Q_lower" linopy/io.py -B5 -A5 | head -50

Repository: FBumann/linopy


<shell_result>
Repository: FBumann/linopy
Exit code: 0

stdout:

792-            # Add quadratic terms
793-            # MOSEK expects lower triangular part only
794-            Q = Qc_list[i]
795-            if Q.nnz > 0:
796-                # Get lower triangular part (MOSEK requirement)
797:                Q_lower = tril(Q).tocoo()
798-                # MOSEK uses 0.5 * x'Qx convention, but our Q is already doubled
799-                # So we need to divide by 2
800:                task.putqconk(con_idx, list(Q_lower.row), list(Q_lower.col),
801:                              list(Q_lower.data / 2))
802-
803-    return task
804-
805-
806-def to_gurobipy(

</shell_result>


Script executed:

#!/bin/bash
# Read the Qc property implementation
sed -n '241,280p' linopy/matrices.py

Repository: FBumann/linopy


<shell_result>
Repository: FBumann/linopy
Exit code: 0

stdout:

    def Qc(self) -> list[csc_matrix]:
        """
        List of Q matrices for quadratic constraints.

        Returns a list where each element is a sparse matrix representing the
        quadratic terms of one constraint. The matrix follows the convention
        x'Qx, where Q is symmetric with doubled diagonal terms.
        """
        m = self._parent
        if not len(m.quadratic_constraints):
            return []

        df = self.flat_qcons
        labels = self.qclabels
        n_vars = len(self.vlabels)

        # Build variable label to index mapping
        var_map = pd.Series(index=self.vlabels, data=np.arange(n_vars))

        matrices = []
        for label in labels:
            label_df = df[(df["labels"] == label) & df["is_quadratic"]]

            if label_df.empty:
                # No quadratic terms - empty matrix
                matrices.append(csc_matrix((n_vars, n_vars)))
                continue

            rows = []
            cols = []
            data = []

            for _, row in label_df.iterrows():
                var1 = int(row["vars1"])
                var2 = int(row["vars2"])
                coeff = row["coeffs"]

                if var1 < 0 or var2 < 0:
                    continue

</shell_result>


Script executed:

#!/bin/bash
# Search for tests related to MOSEK quadratic constraints
rg -n "putqconk|to_mosek" --type py -B2 -A2 | head -60

Repository: FBumann/linopy


<shell_result>
Repository: FBumann/linopy
Exit code: 0

stdout:

test/test_quadratic_constraint.py-576-    """Tests for MOSEK direct API export with quadratic constraints."""
test/test_quadratic_constraint.py-577-
test/test_quadratic_constraint.py:578:    def test_to_mosek_with_quadratic_constraints(
test/test_quadratic_constraint.py-579-        self, m: Model, x: linopy.Variable, y: linopy.Variable
test/test_quadratic_constraint.py-580-    ) -> None:
test/test_quadratic_constraint.py:581:        """Test that to_mosek works with quadratic constraints."""
test/test_quadratic_constraint.py-582-        if "mosek" not in linopy.available_solvers:
test/test_quadratic_constraint.py-583-            pytest.skip("MOSEK not available")
--
test/test_quadratic_constraint.py-587-        m.add_objective(x + 2 * y, sense="max")
test/test_quadratic_constraint.py-588-
test/test_quadratic_constraint.py:589:        from linopy.io import to_mosek
test/test_quadratic_constraint.py-590-
test/test_quadratic_constraint.py:591:        task = to_mosek(m)
test/test_quadratic_constraint.py-592-        # If we got here without error, the export worked
test/test_quadratic_constraint.py-593-        assert task is not None
test/test_quadratic_constraint.py-594-
test/test_quadratic_constraint.py:595:    def test_to_mosek_multidimensional(self) -> None:
test/test_quadratic_constraint.py-596-        """Test MOSEK export with multi-dimensional quadratic constraints."""
test/test_quadratic_constraint.py-597-        if "mosek" not in linopy.available_solvers:
--
test/test_quadratic_constraint.py-606-        m.add_objective((x + 2 * y).sum(), sense="max")
test/test_quadratic_constraint.py-607-
test/test_quadratic_constraint.py:608:        from linopy.io import to_mosek
test/test_quadratic_constraint.py-609-
test/test_quadratic_constraint.py:610:        task = to_mosek(m)
test/test_quadratic_constraint.py-611-        assert task is not None
test/test_quadratic_constraint.py-612-
--
linopy/solvers.py-1738-            )
linopy/solvers.py-1739-        with mosek.Task() as m:
linopy/solvers.py:1740:            m = model.to_mosek(m, explicit_coordinate_names=explicit_coordinate_names)
linopy/solvers.py-1741-
linopy/solvers.py-1742-            return self._solve(
--
linopy/model.py-60-    to_gurobipy,
linopy/model.py-61-    to_highspy,
linopy/model.py:62:    to_mosek,
linopy/model.py-63-    to_netcdf,
linopy/model.py-64-)
--
linopy/model.py-1688-    to_gurobipy = to_gurobipy
linopy/model.py-1689-
linopy/model.py:1690:    to_mosek = to_mosek
linopy/model.py-1691-
linopy/model.py-1692-    to_highspy = to_highspy
--
linopy/io.py-637-
linopy/io.py-638-
linopy/io.py:639:def to_mosek(
linopy/io.py-640-    m: Model, task: Any | None = None, explicit_coordinate_names: bool = False
linopy/io.py-641-) -> Any:
--
linopy/io.py-754-        # Append quadratic constraints to the task
linopy/io.py-755-        # In MOSEK, quadratic constraints are added to regular constraints
linopy/io.py:756:        # with quadratic terms via putqconk
linopy/io.py-757-        qc_start_idx = m.ncons  # Start after linear constraints
linopy/io.py-758-        task.appendcons(n_qcons)
--
linopy/io.py-798-                # MOSEK uses 0.5 * x'Qx convention, but our Q is already doubled

</shell_result>


Script executed:

#!/bin/bash
# Look for how Qc matrices are constructed - are diagonals doubled?
sed -n '241,320p' linopy/matrices.py

Repository: FBumann/linopy


<shell_result>
Repository: FBumann/linopy
Exit code: 0

stdout:

    def Qc(self) -> list[csc_matrix]:
        """
        List of Q matrices for quadratic constraints.

        Returns a list where each element is a sparse matrix representing the
        quadratic terms of one constraint. The matrix follows the convention
        x'Qx, where Q is symmetric with doubled diagonal terms.
        """
        m = self._parent
        if not len(m.quadratic_constraints):
            return []

        df = self.flat_qcons
        labels = self.qclabels
        n_vars = len(self.vlabels)

        # Build variable label to index mapping
        var_map = pd.Series(index=self.vlabels, data=np.arange(n_vars))

        matrices = []
        for label in labels:
            label_df = df[(df["labels"] == label) & df["is_quadratic"]]

            if label_df.empty:
                # No quadratic terms - empty matrix
                matrices.append(csc_matrix((n_vars, n_vars)))
                continue

            rows = []
            cols = []
            data = []

            for _, row in label_df.iterrows():
                var1 = int(row["vars1"])
                var2 = int(row["vars2"])
                coeff = row["coeffs"]

                if var1 < 0 or var2 < 0:
                    continue

                # Map to matrix indices
                i = var_map.get(var1, -1)
                j = var_map.get(var2, -1)

                if i < 0 or j < 0:
                    continue

                if i == j:
                    # Diagonal term - double it for x'Qx convention
                    rows.append(i)
                    cols.append(j)
                    data.append(coeff * 2)
                else:
                    # Off-diagonal - add symmetric entries
                    rows.extend([i, j])
                    cols.extend([j, i])
                    data.extend([coeff, coeff])

            Q = csc_matrix(
                (data, (rows, cols)), shape=(n_vars, n_vars)
            )
            matrices.append(Q)

        return matrices

    @property
    def qc_linear(self) -> csc_matrix | None:
        """
        Matrix of linear coefficients for quadratic constraints.

        Returns a sparse matrix of shape (n_qconstraints, n_variables) where
        each row contains the linear coefficients for one quadratic constraint.
        """
        m = self._parent
        if not len(m.quadratic_constraints):
            return None

        df = self.flat_qcons
        labels = self.qclabels
        n_cons = len(labels)

</shell_result>


Script executed:

#!/bin/bash
# Look for how quadratic terms are handled in to_gurobipy
sed -n '900,1000p' linopy/io.py | grep -A10 -B10 "quadratic\|Q_"

Repository: FBumann/linopy


<shell_result>
Repository: FBumann/linopy
Exit code: 0

stdout:

                elif sign == ">=":
                    sense = gurobipy.GRB.GREATER_EQUAL
                else:
                    sense = gurobipy.GRB.EQUAL

                # Determine constraint name
                if explicit_coordinate_names:
                    label_pos = m.quadratic_constraints.get_label_position(label)
                    if label_pos:
                        con_name, coord = label_pos
                        qc_name = f"{clean_name(con_name)}{print_coord(coord)}#{label}"
                    else:
                        qc_name = f"{clean_name(name)}#{label}"
                else:
                    qc_name = f"qc{label}"

                model.addQConstr(qexpr, sense, rhs, name=qc_name)

</shell_result>


Script executed:

#!/bin/bash
# Find to_gurobipy and how it handles quadratic constraints
rg -n "to_gurobipy" linopy/io.py -A100 | grep -A5 -B5 "Qc\|quadratic" | head -50

Repository: FBumann/linopy


<shell_result>
Repository: FBumann/linopy
Exit code: 0

stdout:

850-    if len(m.constraints):
851-        names = np.vectorize(print_constraint)(M.clabels).astype(object)
852-        c = model.addMConstr(M.A, x, M.sense, M.b)  # type: ignore
853-        c.setAttr("ConstrName", list(names))  # type: ignore
854-
855-    # Add quadratic constraints
856-    if len(m.quadratic_constraints):
857-        for name in m.quadratic_constraints:
858-            qcon = m.quadratic_constraints[name]
859-            df = qcon.to_polars()
860-
861-            # Build label -> (sign, rhs) lookup for multi-dimensional constraints
862-            labels_flat = qcon.labels.values.ravel()
863-            signs_flat = np.broadcast_to(qcon.sign.values, qcon.labels.shape).ravel()
--
870-
871-            # Build QuadExpr for each constraint label
872-            for label in df["labels"].unique().to_list():
873-                label_df = df.filter(pl.col("labels") == label)
874-
875-                # Build the quadratic expression
876-                qexpr = gurobipy.QuadExpr()
877-
878-                # Add linear terms
879-                linear_terms = label_df.filter(~pl.col("is_quadratic"))
880-                for row in linear_terms.iter_rows(named=True):
881-                    coeff = row["coeffs"]
882-                    var_idx = int(row["vars"])
883-                    if var_idx >= 0:
884-                        qexpr.addTerms(coeff, x[var_idx].item())
885-
886-                # Add quadratic terms
887-                quad_terms = label_df.filter(pl.col("is_quadratic"))
888-                for row in quad_terms.iter_rows(named=True):
889-                    coeff = row["coeffs"]
890-                    var1 = int(row["vars1"])
891-                    var2 = int(row["vars2"])
892-                    qexpr.addTerms(coeff, x[var1].item(), x[var2].item())

</shell_result>


Script executed:

#!/bin/bash
# Get more context around the quadratic constraint handling in to_mosek
sed -n '700,810p' linopy/io.py

Repository: FBumann/linopy


<shell_result>
Repository: FBumann/linopy
Exit code: 0

stdout:

    task.putvarboundslice(0, m.nvars, bkx, blx, bux)

    if len(m.binaries.labels) + len(m.integers.labels) > 0:
        idx = [i for (i, v) in enumerate(M.vtypes) if v in ["B", "I"]]
        task.putvartypelist(idx, [mosek.variabletype.type_int] * len(idx))
        if len(m.binaries.labels) > 0:
            bidx = [i for (i, v) in enumerate(M.vtypes) if v == "B"]
            task.putvarboundlistconst(bidx, mosek.boundkey.ra, 0.0, 1.0)

    ## Constraints

    if len(m.constraints) > 0:
        names = np.vectorize(print_constraint)(M.clabels).astype(object)
        for i, n in enumerate(names):
            task.putconname(i, n)
        bkc = [
            (
                (mosek.boundkey.up if b < np.inf else mosek.boundkey.fr)
                if s == "<"
                else (
                    (mosek.boundkey.lo if b > -np.inf else mosek.boundkey.up)
                    if s == ">"
                    else mosek.boundkey.fx
                )
            )
            for s, b in zip(M.sense, M.b)
        ]
        blc = [b if b > -np.inf else 0.0 for b in M.b]
        buc = [b if b < np.inf else 0.0 for b in M.b]
        # blc = M.b
        # buc = M.b
        if M.A is not None:
            A = M.A.tocsr()
            task.putarowslice(
                0, m.ncons, A.indptr[:-1], A.indptr[1:], A.indices, A.data
            )
            task.putconboundslice(0, m.ncons, bkc, blc, buc)

    ## Objective
    if M.Q is not None:
        Q = (0.5 * tril(M.Q + M.Q.transpose())).tocoo()
        task.putqobj(Q.row, Q.col, Q.data)
    task.putclist(list(np.arange(m.nvars)), M.c)

    if m.objective.sense == "max":
        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


def to_gurobipy(
    m: Model, env: Any | None = None, explicit_coordinate_names: bool = False
) -> Any:
    """
    Export the model to gurobipy.

</shell_result>


Script executed:

#!/bin/bash
# Search for any MOSEK documentation or comments about the 0.5 convention
rg -n "0\.5.*x.*Qx\|putqobj\|putqconk" linopy/io.py -B3 -A3

Repository: FBumann/linopy


<shell_result>
Repository: FBumann/linopy
Exit code: 1

</shell_result>


Script executed:

#!/bin/bash
# Search for related issues or TODOs about quadratic constraints
rg -n "TODO|FIXME|BUG|quadratic" linopy/io.py | head -20

Repository: FBumann/linopy


<shell_result>
Repository: FBumann/linopy
Exit code: 0

stdout:

142:def objective_write_quadratic_terms(
195:        objective_write_quadratic_terms(f, quads, print_variable)
425:def quadratic_constraints_to_file(
432:    Write out quadratic constraints of a model to an LP file.
434:    LP format for quadratic constraints (Gurobi/CPLEX style):
437:    if not len(m.quadratic_constraints):
445:        logger.info("Writing quadratic constraints.")
447:    names = list(m.quadratic_constraints)
451:            desc="Writing quadratic constraints.",
456:        qcon = m.quadratic_constraints[name]
482:                label_pos = m.quadratic_constraints.get_label_position(label)
496:            linear_terms = label_df.filter(~pl.col("is_quadratic"))
505:            # Write quadratic terms in brackets
506:            quad_terms = label_df.filter(pl.col("is_quadratic"))
553:        quadratic_constraints_to_file(
750:    if len(m.quadratic_constraints):
751:        # Get the number of quadratic constraints
754:        # Append quadratic constraints to the task
755:        # In MOSEK, quadratic constraints are added to regular constraints
756:        # with quadratic terms via putqconk

</shell_result>


Script executed:

#!/bin/bash
# Check the imports and understand tril function
head -50 linopy/io.py | grep -E "^import|^from"

Repository: FBumann/linopy


<shell_result>
Repository: FBumann/linopy
Exit code: 0

stdout:

from __future__ import annotations
import logging
import shutil
import time
from collections.abc import Callable
from io import BufferedWriter
from pathlib import Path
from tempfile import TemporaryDirectory
from typing import TYPE_CHECKING, Any
import numpy as np
import pandas as pd
import polars as pl
import xarray as xr
from numpy import ones_like, zeros_like
from scipy.sparse import tril, triu
from tqdm import tqdm
from linopy import solvers
from linopy.constants import CONCAT_DIM
from linopy.objective import Objective

</shell_result>


Script executed:

#!/bin/bash
# Look at the quadratic constraint test to see if it verifies correctness
cat linopy/test/test_quadratic_constraint.py | sed -n '576,620p'

Repository: FBumann/linopy


<shell_result>
Repository: FBumann/linopy
Exit code: 0

stderr:

cat: linopy/test/test_quadratic_constraint.py: No such file or directory

</shell_result>


Script executed:

#!/bin/bash
# Look at the to_mosek function docstring and full implementation
sed -n '639,750p' linopy/io.py

Repository: FBumann/linopy


<shell_result>
Repository: FBumann/linopy
Exit code: 0

stdout:

def to_mosek(
    m: Model, task: Any | None = None, explicit_coordinate_names: bool = False
) -> Any:
    """
    Export model to MOSEK.

    Export the model directly to MOSEK without writing files.

    Parameters
    ----------
    m : linopy.Model
    task : empty MOSEK task

    Returns
    -------
    task : MOSEK Task object
    """

    import mosek

    print_variable, print_constraint = get_printers_scalar(
        m, explicit_coordinate_names=explicit_coordinate_names
    )

    if task is None:
        task = mosek.Task()

    task.appendvars(m.nvars)
    task.appendcons(m.ncons)

    M = m.matrices
    # for j, n in enumerate(("x" + M.vlabels.astype(str).astype(object))):
    #    task.putvarname(j, n)

    labels = np.vectorize(print_variable)(M.vlabels).astype(object)
    task.generatevarnames(
        np.arange(0, len(labels)), "%0", [len(labels)], None, [0], list(labels)
    )

    ## Variables

    # MOSEK uses bound keys (free, bounded below or above, ranged and fixed)
    # plus bound values (lower and upper), and it is considered an error to
    # input an infinite value for a finite bound.
    # bkx and bkc define the boundkeys based on upper and lower bound, and blx,
    # bux, blc and buc define the finite bounds. The numerical value of a bound
    # indicated to be infinite by the bound key is ignored by MOSEK.
    bkx = [
        (
            (
                (mosek.boundkey.ra if lb < ub else mosek.boundkey.fx)
                if ub < np.inf
                else mosek.boundkey.lo
            )
            if (lb > -np.inf)
            else (mosek.boundkey.up if (ub < np.inf) else mosek.boundkey.fr)
        )
        for (lb, ub) in zip(M.lb, M.ub)
    ]
    blx = [b if b > -np.inf else 0.0 for b in M.lb]
    bux = [b if b < np.inf else 0.0 for b in M.ub]
    task.putvarboundslice(0, m.nvars, bkx, blx, bux)

    if len(m.binaries.labels) + len(m.integers.labels) > 0:
        idx = [i for (i, v) in enumerate(M.vtypes) if v in ["B", "I"]]
        task.putvartypelist(idx, [mosek.variabletype.type_int] * len(idx))
        if len(m.binaries.labels) > 0:
            bidx = [i for (i, v) in enumerate(M.vtypes) if v == "B"]
            task.putvarboundlistconst(bidx, mosek.boundkey.ra, 0.0, 1.0)

    ## Constraints

    if len(m.constraints) > 0:
        names = np.vectorize(print_constraint)(M.clabels).astype(object)
        for i, n in enumerate(names):
            task.putconname(i, n)
        bkc = [
            (
                (mosek.boundkey.up if b < np.inf else mosek.boundkey.fr)
                if s == "<"
                else (
                    (mosek.boundkey.lo if b > -np.inf else mosek.boundkey.up)
                    if s == ">"
                    else mosek.boundkey.fx
                )
            )
            for s, b in zip(M.sense, M.b)
        ]
        blc = [b if b > -np.inf else 0.0 for b in M.b]
        buc = [b if b < np.inf else 0.0 for b in M.b]
        # blc = M.b
        # buc = M.b
        if M.A is not None:
            A = M.A.tocsr()
            task.putarowslice(
                0, m.ncons, A.indptr[:-1], A.indptr[1:], A.indices, A.data
            )
            task.putconboundslice(0, m.ncons, bkc, blc, buc)

    ## Objective
    if M.Q is not None:
        Q = (0.5 * tril(M.Q + M.Q.transpose())).tocoo()
        task.putqobj(Q.row, Q.col, Q.data)
    task.putclist(list(np.arange(m.nvars)), M.c)

    if m.objective.sense == "max":
        task.putobjsense(mosek.objsense.maximize)
    else:
        task.putobjsense(mosek.objsense.minimize)

    ## Quadratic Constraints
    if len(m.quadratic_constraints):

</shell_result>


Script executed:

#!/bin/bash
# Search for any issues or comments about quadratic constraints in git history or issues
find . -name "*.md" -o -name "*.txt" | xargs grep -l "quadratic\|putqconk" 2>/dev/null | head -5

Repository: FBumann/linopy


<shell_result>
Repository: FBumann/linopy
Exit code: 0

stdout:

./dev-scripts/quadratic_constraints_plan.md
./dev-scripts/quadratic_constraints_remaining_tasks.md
./dev-scripts/quadratic_constraints_status.md

</shell_result>


Script executed:

#!/bin/bash
python3 << 'EOF'
import numpy as np
from scipy.sparse import csc_matrix, tril

# Simulate how Qc creates a matrix for a constraint like: x1^2 + 2*x1*x2 + x2^2
# User coefficients: diag (x1^2): 1, off-diag (x1*x2): 2, diag (x2^2): 1

# As per Qc code:
# Diagonal: coeff * 2 => [2, 2]
# Off-diagonal: symmetric pair with original coeff => [(1,2):2, (2,1):2]

rows = [0, 1, 0, 1]  # diagonal at (0,0) and (1,1), off-diagonal at (0,1) and (1,0)
cols = [0, 1, 1, 0]
data = [2, 2, 2, 2]  # [diag*2, diag*2, coeff, coeff]

Q_full = csc_matrix((data, (rows, cols)), shape=(2, 2))

print("Full symmetric Q matrix (as created by Qc):")
print(Q_full.toarray())
print()

# Now apply what the code does for quadratic constraints:
Q_lower = tril(Q_full).tocoo()
print("Lower triangular part (tril):")
print("rows:", Q_lower.row, "cols:", Q_lower.col, "data:", Q_lower.data)
print()

# After division by 2:
Q_lower_div2 = Q_lower.data / 2
print("After division by 2:", Q_lower_div2)
print()

# What the code sends to MOSEK:
print("What putqconk receives:")
print(f"  At (0,0): {Q_lower_div2[0]} (diagonal)")
print(f"  At (1,0): {Q_lower_div2[1]} (off-diagonal below diagonal)")
print()

# What should be received (reconstructing full then taking lower):
Q_reconstructed = 0.5 * tril(Q_full + Q_full.T).tocoo()
print("What OBJECTIVE sends to MOSEK (0.5 * tril(Q + Q.T)):")
print("rows:", Q_reconstructed.row, "cols:", Q_reconstructed.col, "data:", Q_reconstructed.data)

EOF

Repository: FBumann/linopy


<shell_result>
Repository: FBumann/linopy
Script execution failed with error: Request timed out after 120000ms (requestId=cc11bc26-5f94-45c7-b7b4-40c7a21ecd8b)
</shell_result>


Script executed:

#!/bin/bash
# Read the quadratic constraints documentation
cat ./dev-scripts/quadratic_constraints_plan.md

Repository: FBumann/linopy


<shell_result>
Repository: FBumann/linopy
Exit code: 0

stdout:


<!-- blackdoc:off -->

# Quadratic Constraints and Expressions in Linopy

## Design Document and Implementation Plan

### Executive Summary

This document outlines a plan to extend linopy with support for **quadratic constraints** (QCP/QCQP). Currently, linopy supports:
- Linear constraints (`Constraint` class)
- Linear expressions (`LinearExpression` class)
- Quadratic expressions (`QuadraticExpression` class) — **only for objectives**

The goal is to enable quadratic constraints of the form:

```
x'Qx + a'x ≤ b   (or ≥, =)
```

This feature would make linopy suitable for a broader class of optimization problems including convex QCPs, non-convex QCQPs (where supported by solvers), and Second-Order Cone Programs (SOCPs).

---

## 1. Current Architecture Analysis

### 1.1 Linear Expression (`expressions.py`)

The `LinearExpression` class stores:
```python
Dataset {
    'coeffs': DataArray[float]  # shape: (..., _term)
    'vars': DataArray[int]      # shape: (..., _term), variable labels
    'const': DataArray[float]   # shape: (...)
}
```

Key method for constraint creation (`expressions.py:843-866`):
```python
def to_constraint(self, sign: SignLike, rhs: ConstantLike) -> Constraint:
    all_to_lhs = (self - rhs).data
    data = assign_multiindex_safe(
        all_to_lhs[["coeffs", "vars"]],
        sign=sign,
        rhs=-all_to_lhs.const
    )
    return constraints.Constraint(data, model=self.model)
```

### 1.2 Quadratic Expression (`expressions.py`)

The `QuadraticExpression` class adds a `_factor` dimension (size=2) for storing two variables per quadratic term:
```python
Dataset {
    'coeffs': DataArray[float]  # shape: (..., _factor, _term)
    'vars': DataArray[int]      # shape: (..., _factor, _term)
    'const': DataArray[float]   # shape: (...)
}
```

Currently, `QuadraticExpression.to_constraint()` raises `NotImplementedError` (`expressions.py:1805-1808`).

### 1.3 Linear Constraint (`constraints.py`)

The `Constraint` class stores:
```python
Dataset {
    'coeffs': DataArray[float]  # LHS coefficients
    'vars': DataArray[int]      # Variable labels
    'sign': DataArray[str]      # '=', '<=', '>='
    'rhs': DataArray[float]     # Right-hand side
    'labels': DataArray[int]    # Constraint labels (-1 if masked)
    'dual': DataArray[float]    # [OPTIONAL] Dual values
}
```

### 1.4 Solver Support for Quadratic Constraints

| Solver | QCP Support | QCQP Support | Non-Convex Support |
|--------|-------------|--------------|-------------------|
| Gurobi | ✅ Yes | ✅ Yes | ✅ Yes (v9.0+) |
| CPLEX | ✅ Yes | ✅ Yes | ⚠️ Limited |
| MOSEK | ✅ Yes | ✅ Yes | ❌ Convex only |
| Xpress | ✅ Yes | ✅ Yes | ⚠️ Limited |
| COPT | ✅ Yes | ✅ Yes | ⚠️ Limited |
| SCIP | ✅ Yes | ✅ Yes | ✅ Yes |
| HiGHS | ❌ No | ❌ No | ❌ No |
| GLPK | ❌ No | ❌ No | ❌ No |
| CBC | ❌ No | ❌ No | ❌ No |

**Key Insight**: HiGHS (a common default solver) does NOT support quadratic constraints. This has implications for default behavior and error handling.

---

## 2. Proposed Design

### 2.1 New Class: `QuadraticConstraint`

Create a new `QuadraticConstraint` class parallel to `Constraint`:

```python
class QuadraticConstraint:
    """
    A quadratic constraint of the form: x'Qx + a'x ≤ b (or ≥, =)

    Dataset structure:
    {
        'quad_coeffs': DataArray[float]   # shape: (..., _factor, _qterm)
        'quad_vars': DataArray[int]       # shape: (..., _factor, _qterm)
        'lin_coeffs': DataArray[float]    # shape: (..., _term)
        'lin_vars': DataArray[int]        # shape: (..., _term)
        'sign': DataArray[str]            # '=', '<=', '>='
        'rhs': DataArray[float]           # Right-hand side constant
        'labels': DataArray[int]          # Constraint labels
        'dual': DataArray[float]          # [OPTIONAL] Dual values (only for convex)
    }
    """
```

**Design Rationale**:
- Separate `quad_*` and `lin_*` arrays to allow efficient handling of purely linear terms
- Use `_qterm` dimension (distinct from `_term`) for quadratic terms
- Maintain API consistency with `Constraint` class

### 2.2 Container Class: `QuadraticConstraints`

Add a container class analogous to `Constraints`:

```python
class QuadraticConstraints:
    """
    Container for multiple QuadraticConstraint objects.
    Provides dict-like access and aggregation properties.
    """
```

### 2.3 Model Integration

Extend the `Model` class:

```python
class Model:
    def __init__(self, ...):
        self.constraints = Constraints()      # Linear constraints
        self.quadratic_constraints = QuadraticConstraints()  # NEW

    def add_quadratic_constraints(
        self,
        lhs: QuadraticExpression | Callable,
        sign: SignLike,
        rhs: ConstantLike,
        name: str | None = None,
        coords: CoordsLike | None = None,
        mask: MaskLike | None = None,
    ) -> QuadraticConstraint:
        """Add quadratic constraint(s) to the model."""

    @property
    def has_quadratic_constraints(self) -> bool:
        """Return True if model has any quadratic constraints."""

    @property
    def type(self) -> str:
        """Return problem type: 'LP', 'QP', 'MILP', 'MIQP', 'QCP', 'QCQP', etc."""
```

### 2.4 Expression API Changes

Implement `QuadraticExpression.to_constraint()`:

```python
def to_constraint(self, sign: SignLike, rhs: ConstantLike) -> QuadraticConstraint:
    """
    Convert quadratic expression to a quadratic constraint.

    Parameters
    ----------
    sign : str
        Constraint sense: '<=', '>=', or '='
    rhs : float or array-like
        Right-hand side constant

    Returns
    -------
    QuadraticConstraint
    """
```

Enable comparison operators on `QuadraticExpression`:
```python
# These would create QuadraticConstraint objects
quad_expr <= 10  # Works (returns QuadraticConstraint)
quad_expr >= 5   # Works
quad_expr == 0   # Works
```

---

## 3. Implementation Details

### 3.1 Data Storage for Quadratic Constraints

**Option A: Unified Storage** (simpler, less efficient)
```python
# Store everything with _factor dimension, linear terms have vars[_factor=1] = -1
Dataset {
    'coeffs': DataArray[float]  # shape: (..., _factor, _term)
    'vars': DataArray[int]      # shape: (..., _factor, _term)
    'sign': DataArray[str]
    'rhs': DataArray[float]
    'labels': DataArray[int]
}
```

**Option B: Split Storage** (recommended, more efficient)
```python
# Separate linear and quadratic terms
Dataset {
    'quad_coeffs': DataArray[float]  # shape: (..., _factor, _qterm)
    'quad_vars': DataArray[int]      # shape: (..., _factor, _qterm)
    'lin_coeffs': DataArray[float]   # shape: (..., _term)
    'lin_vars': DataArray[int]       # shape: (..., _term)
    'sign': DataArray[str]
    'rhs': DataArray[float]
    'labels': DataArray[int]
}
```

**Recommendation**: Option B provides clearer separation, easier debugging, and more efficient matrix construction for solvers that handle linear and quadratic parts separately.

### 3.2 Matrix Representation

Add to `MatrixAccessor`:

```python
@property
def Qc(self) -> list[tuple[csc_matrix, ndarray, float, str]]:
    """
    List of quadratic constraint matrices.

    Returns list of tuples: (Q_i, a_i, b_i, sense_i)
    where constraint i is: x'Q_i x + a_i'x {sense_i} b_i
    """

@property
def qc_labels(self) -> ndarray:
    """Labels of quadratic constraints."""
```

### 3.3 Solver Export Functions

#### LP File Format

The LP file format supports quadratic constraints in the `QCROWS` section:

```
Subject To
 c1: x + y <= 10

QCROWS
 qc1: [ x^2 + 2 x * y + y^2 ] + x + y <= 5
End
```

Add function:
```python
def quadratic_constraints_to_file(
    m: Model,
    f: BufferedWriter,
    progress: bool = False,
    explicit_coordinate_names: bool = False,
) -> None:
    """Write quadratic constraints to LP file."""
```

#### Direct API Export

**Gurobi** (`addQConstr` or matrix interface):
```python
def to_gurobipy(m: Model, env=None, ...):
    # ... existing code ...

    # Add quadratic constraints
    for qc in m.quadratic_constraints:
        model.addQConstr(Q, sense, rhs, name)
```

**MOSEK** (`putqconk`):
```python
def to_mosek(m: Model, task=None, ...):
    # ... existing code ...

    # Add quadratic constraints
    for k, (Q, a, b, sense) in enumerate(M.Qc):
        task.putqconk(k, Q.row, Q.col, Q.data)
        task.putarow(k, a.nonzero()[0], a[a.nonzero()])
```

### 3.4 Solution Handling

Quadratic constraints may have dual values (for convex problems):

```python
class QuadraticConstraint:
    @property
    def dual(self) -> DataArray | None:
        """
        Dual values for the quadratic constraint.

        Note: Only available for convex quadratic constraints
        and when the solver provides them.
        """
```

---

## 4. API Design Considerations

### 4.1 Consistency with Existing API

The API should feel natural to existing linopy users:

```python
import linopy as lp

m = lp.Model()
x = m.add_variables(coords=[range(3)], name='x')
y = m.add_variables(name='y')

# Linear constraint (existing)
m.add_constraints(x.sum() <= 10, name='linear_budget')

# Quadratic constraint (new - Option A: via add_constraints)
m.add_constraints(x @ x + y <= 5, name='quad_con')

# Quadratic constraint (new - Option B: via add_quadratic_constraints)
m.add_quadratic_constraints(x @ x + y <= 5, name='quad_con')
```

**Question for discussion**: Should quadratic constraints be added via:
- **Option A**: Same `add_constraints()` method (auto-detect based on expression type)
- **Option B**: Separate `add_quadratic_constraints()` method

**Recommendation**: Start with **Option B** for clarity, with Option A as a future enhancement. This makes the API explicit about what type of constraint is being added.

### 4.2 Operator Overloading

Enable natural syntax on `QuadraticExpression`:

```python
# All should return QuadraticConstraint
x * x <= 10
(x @ x) + y >= 5
2 * x * y == 0
```

### 4.3 Error Handling

```python
# Clear error for unsupported solvers
m.solve(solver='highs')
# Raises: "Solver 'highs' does not support quadratic constraints.
#          Use one of: ['gurobi', 'cplex', 'mosek', 'xpress', 'copt', 'scip']"

# Warning for non-convex constraints with convex-only solvers
m.solve(solver='mosek')
# Warning: "MOSEK requires convex quadratic constraints.
#           Non-convex constraints may cause solver failure."
```

---

## 5. File Structure Changes

### 5.1 New Files

None required - extend existing modules.

### 5.2 Modified Files

| File | Changes |
|------|---------|
| `expressions.py` | Implement `QuadraticExpression.to_constraint()` |
| `constraints.py` | Add `QuadraticConstraint` and `QuadraticConstraints` classes |
| `model.py` | Add `add_quadratic_constraints()`, `quadratic_constraints` property |
| `io.py` | Add LP file export for quadratic constraints |
| `solvers.py` | Add `QUADRATIC_CONSTRAINT_SOLVERS` list |
| `matrices.py` | Add `Qc` property for quadratic constraint matrices |
| `constants.py` | Add any new constants (e.g., `QTERM_DIM = "_qterm"`) |

### 5.3 New Test Files

- `test/test_quadratic_constraint.py` — Unit tests for QuadraticConstraint class
- `test/test_quadratic_optimization.py` — Integration tests with solvers

---

## 6. Implementation Phases

### Phase 1: Core Data Structures (Week 1-2)

1. Add `QTERM_DIM` constant
2. Implement `QuadraticConstraint` class with basic functionality:
   - `__init__`, `__repr__`, `__getitem__`
   - Properties: `labels`, `sign`, `rhs`, `lhs`, `mask`, etc.
   - Methods: `to_polars()`, `flat`
3. Implement `QuadraticConstraints` container
4. Add `QuadraticExpression.to_constraint()` method

### Phase 2: Model Integration (Week 2-3)

1. Add `Model.quadratic_constraints` property
2. Implement `Model.add_quadratic_constraints()` method
3. Update `Model.type` property for QCP/QCQP detection
4. Add `Model.has_quadratic_constraints` property
5. Update constraint label management

### Phase 3: Solver Export (Week 3-4)

1. Extend LP file writer with `QCROWS` section
2. Update `to_gurobipy()` for quadratic constraints
3. Update `to_mosek()` for quadratic constraints
4. Update other direct-API solvers (CPLEX, Xpress, COPT)
5. Add solver compatibility checks

### Phase 4: Solution Handling (Week 4-5)

1. Parse quadratic constraint duals from solver results
2. Map duals back to constraint coordinates
3. Add `QuadraticConstraint.dual` property

### Phase 5: Testing & Documentation (Week 5-6)

1. Comprehensive unit tests
2. Integration tests with each supported solver
3. Update documentation and examples
4. Add tutorial notebook

---

## 7. Code Examples

### 7.1 Basic Usage

```python
import linopy as lp

# Create model
m = lp.Model()

# Add variables
x = m.add_variables(lower=0, upper=10, name='x')
y = m.add_variables(lower=0, upper=10, name='y')

# Quadratic objective (already supported)
m.add_objective(x**2 + y**2 + x + y)

# Linear constraint (already supported)
m.add_constraints(x + y >= 1, name='sum_bound')

# NEW: Quadratic constraint
m.add_quadratic_constraints(x**2 + y**2 <= 25, name='circle')

# Solve with quadratic-capable solver
m.solve(solver='gurobi')

print(f"x = {x.solution.values}")
print(f"y = {y.solution.values}")
```

### 7.2 Multi-dimensional Quadratic Constraints

```python
import linopy as lp
import pandas as pd

m = lp.Model()

# Index set
times = pd.Index(range(24), name='time')

# Variables
power = m.add_variables(coords=[times], name='power')
reserve = m.add_variables(coords=[times], name='reserve')

# Quadratic constraint at each time step
# power²[t] + reserve²[t] <= capacity[t]
capacity = [100] * 24
qc = m.add_quadratic_constraints(
    power**2 + reserve**2 <= capacity,
    name='capacity_limit'
)

print(qc)
# QuadraticConstraint `capacity_limit` [time: 24]:
# -------------------------------------------
# 0: power[0]² + reserve[0]² <= 100
# 1: power[1]² + reserve[1]² <= 100
# ...
```

### 7.3 Rule-based Quadratic Constraints

```python
def capacity_rule(m, t):
    """Quadratic capacity constraint at time t."""
    return m['power'][t]**2 + m['reserve'][t]**2 <= capacity[t]

m.add_quadratic_constraints(capacity_rule, coords=[times], name='capacity')
```

---

## 8. Design Decisions (Resolved)

### Q1: Unified vs Separate `add_constraints` method?

**Decision**: Use separate `add_quadratic_constraints()` method.

```python
m.add_quadratic_constraints(x**2 <= 10)  # Quadratic
m.add_constraints(x <= 10)               # Linear
```

**Rationale**: Explicit API is clearer. Auto-detection can be added later by routing `add_constraints()` to `add_quadratic_constraints()` when a `QuadraticExpression` is detected.

### Q2: Storage location for quadratic constraints?

**Decision**: Separate containers.

```python
m.constraints           # Linear only
m.quadratic_constraints # Quadratic only
```

**Rationale**: Simpler implementation, matches current pattern for variables/constraints, avoids complexity of mixed container.

### Q3: How to handle mixed linear+quadratic in same named constraint group?

**Decision**: Each named constraint should be uniformly linear or quadratic. Mixed cases require two separate constraints.

### Q4: Convexity checking?

**Decision**: Defer to solver.

**Rationale**: Avoids computational overhead and complex eigenvalue analysis. Solvers like Gurobi and MOSEK provide clear error messages for non-convex constraints. We can add clear documentation about convexity requirements.

---

## 9. References

### Solver Documentation

- [Gurobi Quadratic Constraints](https://docs.gurobi.com/projects/optimizer/en/current/concepts/modeling/constraints.html)
- [MOSEK Conic Constraints](https://docs.mosek.com/latest/pythonapi/tutorial-cqo-shared.html)
- [CPLEX Quadratic Constraints](https://www.ibm.com/docs/en/icos/latest?topic=programming-adding-quadratic-constraints)

### Related Issues/PRs

- Current `NotImplementedError` in `expressions.py:1805-1808`
- Test showing limitation: `test/test_quadratic_expression.py:331`

### Mathematical Background

A quadratically constrained quadratic program (QCQP) has the form:

```
minimize    (1/2)x'Q₀x + c'x
subject to  (1/2)x'Qᵢx + aᵢ'x ≤ bᵢ    for i = 1,...,m
            Ax = b
            l ≤ x ≤ u
```

Where:
- Q₀ is the objective quadratic matrix (may be zero for QCP)
- Qᵢ are constraint quadratic matrices
- aᵢ are constraint linear coefficient vectors
- bᵢ are constraint right-hand sides

For **convex** QCPs, all Qᵢ must be positive semi-definite (PSD).

---

## 10. Summary

Adding quadratic constraint support to linopy is a significant but feasible enhancement. The key design decisions are:

1. **New class**: `QuadraticConstraint` parallel to `Constraint`
2. **Split storage**: Separate `quad_*` and `lin_*` arrays for efficiency
3. **Explicit API**: `add_quadratic_constraints()` method
4. **Solver filtering**: Clear error messages for unsupported solvers
5. **Phased implementation**: Core → Model → Export → Tests

This enhancement would expand linopy's capabilities to cover:
- Convex QCPs (portfolio optimization, geometric programming)
- QCQPs (facility location, engineering design)
- SOCPs via quadratic constraint reformulation

---

*Document Version: 1.0*
*Date: 2025-01-25*
*Status: Draft for Discussion*

</shell_result>


Script executed:

#!/bin/bash
# Also check the status file
cat ./dev-scripts/quadratic_constraints_status.md

Repository: FBumann/linopy


<shell_result>
Repository: FBumann/linopy
Exit code: 0

stdout:

# Quadratic Constraints Implementation Status

**Date:** 2024-11-25
**Status:** Phase 1 Complete ✅

## Related Documents

- **Original Plan:** `dev-scripts/quadratic_constraints_plan.md` - Detailed design document
- **Remaining Tasks:** `dev-scripts/quadratic_constraints_remaining_tasks.md` - What's left to do

---

## Summary

Quadratic constraints (QCP/QCQP) have been successfully added to linopy. Users can now create and solve optimization problems with constraints of the form:

```
x'Qx + a'x ≤ b   (or ≥, =)
```

## Usage Example

```python
import linopy

m = linopy.Model()

# Variables
x = m.add_variables(lower=0, name='x')
y = m.add_variables(lower=0, name='y')

# Linear constraint (existing)
m.add_constraints(x + y <= 10, name='budget')

# Quadratic constraint (NEW!)
m.add_quadratic_constraints(
    x*x + 2*x*y + y*y + 3*x + 4*y,
    "<=",
    100,
    name='quadratic_budget'
)

# Objective
m.add_objective(x + 2*y)

# Solve (with Gurobi, MOSEK, CPLEX, etc.)
m.solve(solver_name='gurobi')
```

## What Was Implemented

### Core Components

| Component | File | Status |
|-----------|------|--------|
| `QTERM_DIM` constant | `constants.py` | ✅ |
| `QuadraticConstraint` class | `constraints.py` | ✅ |
| `QuadraticConstraints` container | `constraints.py` | ✅ |
| `QuadraticExpression.to_constraint()` | `expressions.py` | ✅ |
| `Model.add_quadratic_constraints()` | `model.py` | ✅ |
| `Model.quadratic_constraints` property | `model.py` | ✅ |
| `Model.has_quadratic_constraints` | `model.py` | ✅ |
| Model type detection (QCLP/QCQP) | `model.py` | ✅ |
| Solver validation | `model.py` | ✅ |
| `QUADRATIC_CONSTRAINT_SOLVERS` list | `solvers.py` | ✅ |
| LP file export | `io.py` | ✅ |
| Gurobi direct export | `io.py` | ✅ |
| Unit tests | `test/test_quadratic_constraint.py` | ✅ |

### Supported Solvers

| Solver | Support | Notes |
|--------|---------|-------|
| Gurobi | ✅ | Full support via `addQConstr()` |
| CPLEX | ✅ | Via LP file |
| MOSEK | ✅ | Via LP file (direct API pending) |
| Xpress | ✅ | Via LP file |
| COPT | ✅ | Via LP file |
| SCIP | ✅ | Via LP file |
| HiGHS | ❌ | Does not support QC - validation error raised |

### Model Type Strings

The `Model.type` property now returns:

| Type | Meaning |
|------|---------|
| `LP` | Linear constraints, linear objective |
| `QP` | Linear constraints, quadratic objective |
| `QCLP` | Quadratic constraints, linear objective |
| `QCQP` | Quadratic constraints, quadratic objective |
| `MILP` | Mixed-integer linear |
| `MIQP` | Mixed-integer quadratic objective |
| `MIQCLP` | Mixed-integer with quadratic constraints |
| `MIQCQP` | Mixed-integer QC with quadratic objective |

---

## What's NOT Yet Implemented

See `dev-scripts/quadratic_constraints_remaining_tasks.md` for full details.

### High Priority
1. Matrix accessor (`Qc` property in `matrices.py`)
2. MOSEK direct API support
3. HiGHS explicit validation in `to_highspy()`
4. Solution/dual retrieval for quadratic constraints

### Medium Priority
5. Multi-dimensional quadratic constraints (with coordinates)
6. netCDF serialization
7. Constraint modification methods

### Low Priority
8. Convexity checking (optional)
9. Documentation

---

## Testing

All 23 unit tests pass:

```bash
pytest test/test_quadratic_constraint.py -v
# 23 passed
```

Existing quadratic expression tests updated:
```bash
pytest test/test_quadratic_expression.py -v
# 32 passed (test_quadratic_to_constraint updated)
```

---

## Design Decisions Made

1. **Separate method**: Using `add_quadratic_constraints()` instead of overloading `add_constraints()` - clearer API, can add auto-detection later

2. **Separate container**: `Model.quadratic_constraints` is separate from `Model.constraints` - cleaner code, explicit handling

3. **Defer convexity**: Convexity checking deferred to solver - avoids false positives, solvers handle it better

4. **LP format**: Using Gurobi-style LP format with `[ ]` brackets for quadratic terms

---

## Files Modified

```
linopy/
├── constants.py      # +QTERM_DIM
├── constraints.py    # +QuadraticConstraint, +QuadraticConstraints (~500 lines)
├── expressions.py    # +to_constraint() for QuadraticExpression
├── model.py          # +add_quadratic_constraints(), properties, validation
├── solvers.py        # +QUADRATIC_CONSTRAINT_SOLVERS
└── io.py             # +quadratic_constraints_to_file(), updated to_gurobipy()

test/
├── test_quadratic_constraint.py  # NEW - 23 tests
└── test_quadratic_expression.py  # Updated 1 test
```

---

## Next Steps for Code Agent

1. Pick a task from `quadratic_constraints_remaining_tasks.md`
2. Tasks are independent - any order works
3. Recommended first: Matrix Accessor (enables other features)
4. Run `pytest test/test_quadratic_constraint.py` to verify nothing breaks

</shell_result>


Script executed:

#!/bin/bash
# Read the remaining tasks file
cat ./dev-scripts/quadratic_constraints_remaining_tasks.md

Repository: FBumann/linopy


<shell_result>
Repository: FBumann/linopy
Exit code: 0

stdout:

# Quadratic Constraints - Remaining Tasks

## Implementation Status

### ✅ Completed (Phase 1 - Core Implementation)

1. **Core Data Structures**
   - `QTERM_DIM` constant in `constants.py`
   - `QuadraticConstraint` class in `constraints.py`
   - `QuadraticConstraints` container class
   - `QuadraticExpression.to_constraint()` method

2. **Model Integration**
   - `Model.quadratic_constraints` property
   - `Model.add_quadratic_constraints()` method
   - `Model.has_quadratic_constraints` property
   - `Model.type` property updated for QCP/QCQP detection

3. **Solver Support**
   - `QUADRATIC_CONSTRAINT_SOLVERS` list in `solvers.py`
   - Solver validation in `Model.solve()` - rejects unsupported solvers

4. **Export Functionality**
   - LP file export via `quadratic_constraints_to_file()` in `io.py`
   - Gurobi direct API export via updated `to_gurobipy()`

5. **Tests**
   - `test/test_quadratic_constraint.py` - 23 unit tests
   - Updated `test_quadratic_to_constraint` in `test_quadratic_expression.py`

---

## 🔲 Remaining Tasks (Phase 2)

### High Priority

#### 1. Matrix Accessor for Quadratic Constraints
**File:** `linopy/matrices.py`

Add `Qc` property to `MatrixAccessor` class for quadratic constraint matrices.

```python
@property
def Qc(self) -> list[scipy.sparse.csc_matrix]:
    """Return list of Q matrices for quadratic constraints."""
    # Each quadratic constraint has its own Q matrix
    pass

@property
def qc_linear(self) -> scipy.sparse.csc_matrix:
    """Return linear coefficients for quadratic constraints."""
    pass

@property
def qc_sense(self) -> np.ndarray:
    """Return sense array for quadratic constraints."""
    pass

@property
def qc_rhs(self) -> np.ndarray:
    """Return RHS values for quadratic constraints."""
    pass
```

#### 2. MOSEK Direct API Support
**File:** `linopy/io.py` - `to_mosek()` function

Add quadratic constraint support to MOSEK export:
```python
# After linear constraints section
if len(m.quadratic_constraints):
    for name in m.quadratic_constraints:
        qcon = m.quadratic_constraints[name]
        # Use task.appendcone() or task.putqconk() for quadratic constraints
```

Reference: [MOSEK Python API - Quadratic Constraints](https://docs.mosek.com/latest/pythonapi/tutorial-qcqo.html)

#### 3. HiGHSpy Validation
**File:** `linopy/io.py` - `to_highspy()` function

Add explicit error if model has quadratic constraints:
```python
if len(m.quadratic_constraints):
    raise ValueError(
        "HiGHS does not support quadratic constraints. "
        "Use a solver that supports QCP: gurobi, cplex, mosek, xpress, copt, scip"
    )
```

#### 4. Solution Retrieval for Quadratic Constraints
**Files:** `linopy/solvers.py`, `linopy/constraints.py`

- Add dual value retrieval for quadratic constraints (where supported)
- Store duals in `QuadraticConstraint.dual` property
- Update solver result parsing

### Medium Priority

#### 5. Multi-dimensional Quadratic Constraints
**File:** `linopy/constraints.py`, `linopy/model.py`

Currently, quadratic constraints are primarily scalar. Add support for:
- Broadcasting over coordinates (like linear constraints)
- `iterate_slices()` support for memory-efficient processing
- Coordinate-based indexing

Example API:
```python
# Should work with coordinates
m.add_quadratic_constraints(
    x * x + y * y,  # where x, y have dims=['time', 'node']
    "<=",
    100,
    name="qc"
)
```

#### 6. Constraint Modification Methods
**File:** `linopy/constraints.py`

Add methods to `QuadraticConstraint`:
```python
def modify_rhs(self, new_rhs: ConstantLike) -> None:
    """Modify the right-hand side of the constraint."""

def modify_coeffs(self, new_coeffs: xr.DataArray) -> None:
    """Modify coefficients of the constraint."""
```

#### 7. netCDF Serialization
**File:** `linopy/io.py` - `to_netcdf()` and `read_netcdf()`

Add quadratic constraints to model serialization:
```python
# In to_netcdf()
qcons = [
    with_prefix(qcon.data, f"quadratic_constraints-{name}")
    for name, qcon in m.quadratic_constraints.items()
]

# In read_netcdf()
# Parse quadratic_constraints-* prefixed datasets
```

### Low Priority

#### 8. Convexity Checking (Optional)
**File:** `linopy/constraints.py` or new `linopy/analysis.py`

Add optional convexity verification:
```python
def check_convexity(self) -> bool:
    """
    Check if quadratic constraint is convex.

    A quadratic constraint x'Qx + a'x <= b is convex if Q is
    positive semidefinite.
    """
    # Extract Q matrix
    # Check eigenvalues or use Cholesky decomposition
    pass
```

#### 9. Constraint Printing Improvements
**File:** `linopy/constraints.py`

Enhance `_format_single_constraint()` for better display:
- Handle large constraints with truncation
- Add option for matrix form display
- Support LaTeX output

#### 10. Documentation
**Files:** `doc/` directory

- Add quadratic constraints section to user guide
- Document supported solvers and their limitations
- Add examples for common QCP formulations (portfolio optimization, etc.)

---

## Testing Tasks

### Unit Tests to Add

1. **Multi-dimensional constraints** - Test with coordinates
2. **Edge cases** - Empty constraints, single term, all linear terms
3. **Numerical precision** - Very small/large coefficients
4. **Memory efficiency** - Large constraint sets with `iterate_slices`

### Integration Tests

1. **Solver round-trip** - Create model, solve, verify solution
2. **File format round-trip** - Write LP, read back, compare
3. **Cross-solver consistency** - Same problem, multiple solvers

### Solver-Specific Tests

```python
@pytest.mark.parametrize("solver", ["gurobi", "mosek", "cplex"])
def test_qcp_solve(solver):
    """Test solving QCP with different solvers."""
    if solver not in available_solvers:
        pytest.skip(f"{solver} not available")
    # ... test code
```

---

## Code Quality Tasks

1. **Type hints** - Ensure all new functions have complete type annotations
2. **Docstrings** - Add NumPy-style docstrings to all public methods
3. **Linting** - Run `ruff check --fix` on all modified files
4. **MyPy** - Fix any type errors in new code

---

## Architecture Notes

### Data Structure

```
QuadraticConstraint.data (xarray.Dataset):
├── quad_coeffs: (_qterm, _factor) float64  # Quadratic term coefficients
├── quad_vars: (_qterm, _factor) int64      # Variable indices for quad terms
├── lin_coeffs: (_term) float64             # Linear term coefficients
├── lin_vars: (_term) int64                 # Variable indices for linear terms
├── sign: str                               # "<=", ">=", or "="
├── rhs: float64                            # Right-hand side value
└── labels: int64                           # Constraint label/index
```

### LP File Format

```
qc0:
+3.0 x0
+4.0 x1
+ [
+1.0 x0 ^ 2
+2.0 x0 * x1
+1.0 x1 ^ 2
]
<= +100.0
```

### Solver Compatibility Matrix

| Solver | QCP Support | API Method |
|--------|-------------|------------|
| Gurobi | ✅ | `addQConstr()` |
| CPLEX  | ✅ | `add_quadratic_constraint()` |
| MOSEK  | ✅ | `putqconk()` |
| Xpress | ✅ | `addConstraint()` with quadratic |
| COPT   | ✅ | `addQConstr()` |
| SCIP   | ✅ | LP file import |
| HiGHS  | ❌ | Not supported |

---

## Suggested Task Order for Next Agent

1. **Matrix Accessor** (`matrices.py`) - Enables programmatic access to constraint data
2. **MOSEK support** (`io.py`) - Important solver for QCP
3. **Multi-dimensional constraints** - Core functionality improvement
4. **netCDF serialization** - Model persistence
5. **Documentation** - User-facing docs

Each task is relatively independent and can be completed in a single session.

</shell_result>


Script executed:

#!/bin/bash
# Let's trace through what the Qc property creates and what tril does
# by examining the code flow more carefully

# Show the complete Qc property one more time to verify our understanding
sed -n '241,320p' linopy/matrices.py | grep -A30 "if i == j"

Repository: FBumann/linopy


<shell_result>
Repository: FBumann/linopy
Exit code: 0

stdout:

                if i == j:
                    # Diagonal term - double it for x'Qx convention
                    rows.append(i)
                    cols.append(j)
                    data.append(coeff * 2)
                else:
                    # Off-diagonal - add symmetric entries
                    rows.extend([i, j])
                    cols.extend([j, i])
                    data.extend([coeff, coeff])

            Q = csc_matrix(
                (data, (rows, cols)), shape=(n_vars, n_vars)
            )
            matrices.append(Q)

        return matrices

    @property
    def qc_linear(self) -> csc_matrix | None:
        """
        Matrix of linear coefficients for quadratic constraints.

        Returns a sparse matrix of shape (n_qconstraints, n_variables) where
        each row contains the linear coefficients for one quadratic constraint.
        """
        m = self._parent
        if not len(m.quadratic_constraints):
            return None

        df = self.flat_qcons

</shell_result>


Script executed:

#!/bin/bash
# Let me trace the mathematical logic:
# User input constraint: x1^2 + 2*x1*x2 + x2^2 + 3*x1 + 4*x2 <= 5
# This represents: 1*x1^2 + 2*x1*x2 + 1*x2^2 (quadratic) + 3*x1 + 4*x2 (linear) <= 5

# In standard form for MOSEK: x'Qx + a'x <= b
# where Q = [[1, 1], [1, 1]] (symmetric, with off-diag=1)
#       a = [3, 4]
#       b = 5

# How Qc creates the matrix:
# - Diagonal entry (0,0): 1 * 2 = 2
# - Diagonal entry (1,1): 1 * 2 = 2  
# - Off-diagonal entries (0,1) and (1,0): 2 (as symmetric pair)
# So Q_full = [[2, 2], [2, 2]]

# What tril extracts:
# Q_lower = [[2, 0], [2, 2]] in dense form
# As COO: rows=[0,1,1], cols=[0,0,1], data=[2,2,2]

# After division by 2:
# data = [1, 1, 1]
# So what MOSEK receives:
# At (0,0): 1, at (1,0): 1, at (1,1): 1
# MOSEK reconstructs full: Q = [[1, 1], [1, 1]]

# This is WRONG! The off-diagonal should be 2, not 1
# Because we want x'Qx = 2*x1*x2 = [x1 x2] [[?, ?], [?, ?]] [x1; x2]
# For off-diagonal off-diag element, contribution is 2*off_diag*x1*x2
# So we need off_diag = 1 in matrix to get 2*1*x1*x2 = 2*x1*x2

# Let me verify MOSEK convention...
echo "Analysis complete - waiting for MOSEK docs"

Repository: FBumann/linopy


<shell_result>
Repository: FBumann/linopy
Exit code: 0

stdout:

Analysis complete - waiting for MOSEK docs

</shell_result>


Script executed:

#!/bin/bash
# Look for test files related to MOSEK
find test -name "*mosek*" -o -name "*quadratic*" | xargs ls -la 2>/dev/null

Repository: FBumann/linopy


<shell_result>
Repository: FBumann/linopy
Exit code: 0

stdout:

-rw-r--r-- 1 jailuser jailuser 23403 Nov 25 22:27 test/test_quadratic_constraint.py
-rw-r--r-- 1 jailuser jailuser 10983 Nov 25 22:27 test/test_quadratic_expression.py

</shell_result>


Script executed:

#!/bin/bash
# Look at the MOSEK test
sed -n '576,620p' test/test_quadratic_constraint.py

Repository: FBumann/linopy


<shell_result>
Repository: FBumann/linopy
Exit code: 0

stdout:

    """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."""

    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."""

</shell_result>


Script executed:

#!/bin/bash
# Search for all MOSEK related tests
rg -n "to_mosek\|def test.*mosek" test/test_quadratic_constraint.py -A10

Repository: FBumann/linopy


<shell_result>
Repository: FBumann/linopy
Exit code: 1

</shell_result>


Script executed:

#!/bin/bash
# Check the objective Q handling in detail
sed -n '730,750p' linopy/io.py | cat -n

Repository: FBumann/linopy


<shell_result>
Repository: FBumann/linopy
Exit code: 0

stdout:

     1	        # buc = M.b
     2	        if M.A is not None:
     3	            A = M.A.tocsr()
     4	            task.putarowslice(
     5	                0, m.ncons, A.indptr[:-1], A.indptr[1:], A.indices, A.data
     6	            )
     7	            task.putconboundslice(0, m.ncons, bkc, blc, buc)
     8	
     9	    ## Objective
    10	    if M.Q is not None:
    11	        Q = (0.5 * tril(M.Q + M.Q.transpose())).tocoo()
    12	        task.putqobj(Q.row, Q.col, Q.data)
    13	    task.putclist(list(np.arange(m.nvars)), M.c)
    14	
    15	    if m.objective.sense == "max":
    16	        task.putobjsense(mosek.objsense.maximize)
    17	    else:
    18	        task.putobjsense(mosek.objsense.minimize)
    19	
    20	    ## Quadratic Constraints
    21	    if len(m.quadratic_constraints):

</shell_result>


Script executed:

#!/bin/bash
# Search for M.Q property or Q property in matrices.py related to objective
rg -n "def Q\|@property" linopy/matrices.py | head -30

Repository: FBumann/linopy


<shell_result>
Repository: FBumann/linopy
Exit code: 0

</shell_result>


Script executed:

#!/bin/bash
# Find the Q property for objective
sed -n '1,200p' linopy/matrices.py | grep -n "def Q\|@property" | head -10

Repository: FBumann/linopy


<shell_result>
Repository: FBumann/linopy
Exit code: 0

stdout:

70:    @property
76:    @property
95:    @property
125:    @property
131:    @property
139:    @property
148:    @property
154:    @property
160:    @property
173:    @property

</shell_result>


Script executed:

#!/bin/bash
# Find Q property more specifically
rg -n "^\s*def Q\(" linopy/matrices.py -A2

Repository: FBumann/linopy


<shell_result>
Repository: FBumann/linopy
Exit code: 0

stdout:

174:    def Q(self) -> csc_matrix | None:
175-        """Matrix objective coefficients of quadratic terms of all non-missing variables."""
176-        m = self._parent

</shell_result>


Script executed:

#!/bin/bash
# Look for Q in property lines with context
rg -n "@property" linopy/matrices.py -A1 | grep -B1 -A1 "\bQ\b"

Repository: FBumann/linopy


<shell_result>
Repository: FBumann/linopy
Exit code: 0

stdout:

173:    @property
174-    def Q(self) -> csc_matrix | None:
--

</shell_result>


Script executed:

#!/bin/bash
# Read the Q property implementation
sed -n '174,220p' linopy/matrices.py

Repository: FBumann/linopy


<shell_result>
Repository: FBumann/linopy
Exit code: 0

stdout:

    def Q(self) -> csc_matrix | None:
        """Matrix objective coefficients of quadratic terms of all non-missing variables."""
        m = self._parent
        expr = m.objective.expression
        if not isinstance(expr, expressions.QuadraticExpression):
            return None
        return expr.to_matrix()[self.vlabels][:, self.vlabels]

    # Quadratic constraint accessors

    @cached_property
    def flat_qcons(self) -> pd.DataFrame:
        """Flat DataFrame of all quadratic constraints."""
        m = self._parent
        return m.quadratic_constraints.flat

    @property
    def qclabels(self) -> ndarray:
        """Vector of labels of all non-missing quadratic constraints."""
        df: pd.DataFrame = self.flat_qcons
        if df.empty:
            return np.array([], dtype=int)
        return np.sort(df["labels"].unique())

    @property
    def qc_sense(self) -> ndarray:
        """Vector of senses of all non-missing quadratic constraints."""
        m = self._parent
        if not len(m.quadratic_constraints):
            return np.array([], dtype="<U2")

        labels = self.qclabels
        senses = np.empty(len(labels), dtype="<U2")

        for name in m.quadratic_constraints:
            qcon = m.quadratic_constraints[name]
            qc_labels = qcon.labels.values.ravel()
            qc_signs = np.broadcast_to(qcon.sign.values, qcon.labels.shape).ravel()
            for i, lab in enumerate(qc_labels):
                if lab != -1:
                    idx = np.searchsorted(labels, lab)
                    senses[idx] = str(qc_signs[i])

        return senses

    @property
    def qc_rhs(self) -> ndarray:

</shell_result>


Perfect! Now I have the complete picture. The Q property delegates to expr.to_matrix() from the QuadraticExpression. This confirms the pattern. Let me finalize my analysis and provide the corrected review comment.

Verify the coefficient scaling for MOSEK quadratic constraints.

The Qc matrices from matrices.py (line 241) have only diagonal terms doubled (documented: "Q is symmetric with doubled diagonal terms"). The code divides all entries uniformly by 2, which incorrectly scales off-diagonal terms.

For example, a constraint like x*y with coefficient 1:

  • Qc stores symmetrically as (0,1)=1, (1,0)=1 (not doubled)
  • tril(Q) extracts (1,0)=1
  • Division by 2 gives 0.5
  • MOSEK receives coefficient 0.5 for the off-diagonal term, which is incorrect (should be 1)

The objective export (line 745) handles this correctly by reconstructing the full symmetric matrix: 0.5 * tril(M.Q + M.Q.transpose()). Apply the same pattern to quadratic constraints rather than dividing all entries uniformly.

                # 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))

The uniform division by 2 is incorrect for off-diagonal entries. Consider reconstructing the symmetric matrix like the objective does to ensure proper coefficient scaling.

🤖 Prompt for AI Agents
In linopy/io.py around lines 798-801, the code divides all entries of the
lower-triangular Q by 2 which incorrectly scales off-diagonal terms; instead
reconstruct the symmetric Q the same way the objective export does (i.e., form
0.5 * (tril_Q + tril_Q.T) so off-diagonals keep full value while diagonals are
halved), then pass the rows, cols and data of that reconstructed symmetric
lower-triangular matrix to task.putqconk (ensure conversion to lists as done
elsewhere).


return task


Expand Down
5 changes: 4 additions & 1 deletion linopy/matrices.py
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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:
Expand Down
114 changes: 102 additions & 12 deletions test/test_quadratic_constraint.py
Original file line number Diff line number Diff line change
Expand Up @@ -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:
Expand Down Expand Up @@ -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."""

Expand Down