Repository: PyPSA/linopy

Issue:
Support Quadratic terms 
A new class QuadraticExpression could be added which follows the same scheme as LinearExpression but has the data fields
* coeffs 
* vars1
* vars2

Arithmetic operations with lhs have to be checked rigorously. 

Discussion:

Patch:
diff --git a/doc/release_notes.rst b/doc/release_notes.rst
index d1fc0e9..9c98835 100644
--- a/doc/release_notes.rst
+++ b/doc/release_notes.rst
@@ -6,6 +6,8 @@ Upcoming Release
 
 
 **New Features**
+
+* Linopy now supports quadratic programming. Therefore a new class `QuadraticExpression` was created, which can be assigned to the objective function. The `QuadraticExpression` class supports the same arithmetic operations as the `LinearExpression` and can be created by multiplying two `Variable` or `LinearExpression` objects. Note for the latter, the number of stacked terms must be equal to one (`expr.nterm == 1`).
 * `LinearExpression`'s now support constant values. This allows defining linear expressions with numeric constant values, like `x + 5`.
 * When defining constraints, it is not needed to separate variables from constants anymore. Thus, expressions  like `x <= y` or `5 * x + 10 >= y` are supported.
 * The classes `Variable`, `LinearExpression` and `Constraint` now have a `loc` method.
@@ -15,12 +17,15 @@ Upcoming Release
 * The `groupby` function now supports passing `pandas.Dataframe`s as groupby keys. These allows to group by multiple variables at once.
 * The performance of the `groupby` function was strongly increased. In large operations a speedup of 10x was observed.
 
+
 **Deprecations**
+
 * The class `AnonymousConstraint` is now deprecated in the favor of `Constraint`. The latter can now be assigned to a model or not.
 * The `ravel` and `iter_ravel` method of the `Variables` and `Constraints` class is now deprecated in favor of the `flat` method.
 
 
 **Breaking Changes**
+
 * The `data` attribute of Variables and Constraints now returns a `xarray.Dataset` object instead of a `xarray.DataArray` object with the labels only.
 * The deprecated `groupby_sum` function was removed in favor of the `groupby` method.
 * The deprecated `rolling_sum` function was removed in favor of the `rolling` method.
@@ -34,10 +39,12 @@ Upcoming Release
 
 
 **Internal Changes**
+
 * The internal data fields in `Variable` and `Constraint` are now always broadcasted to have aligned indexes. This allows for a more consistent handling of the objects.
 * The inner structure of the `Variable`, `Variables`, `Constraint` and `Constraints` class has changed to a more stable design. All information of the `Variable` and the `Constraint` class is now stored in the `data` field. The `data` field is a `xarray.Dataset` object. The `Variables` and `Constraints` class "simple" containers for the `Variable` and `Constraint` objects, stored in dictionary under the `data` field. This design allows for a more flexible handling of individual variables and constraints.
 
 **Other**
+
 * License changed to MIT license.
 
 
diff --git a/linopy/common.py b/linopy/common.py
index 932bbc4..d614f23 100644
--- a/linopy/common.py
+++ b/linopy/common.py
@@ -18,7 +18,7 @@ from xarray import DataArray, Dataset, align, apply_ufunc, merge
 from xarray.core import indexing, utils
 
 from linopy.config import options
-from linopy.constants import SIGNS, SIGNS_alternative, sign_replace_dict
+from linopy.constants import SIGNS, TERM_DIM, SIGNS_alternative, sign_replace_dict
 
 
 def maybe_replace_sign(sign):
@@ -111,7 +111,7 @@ def fill_missing_coords(ds):
 
     # Fill in missing integer coordinates
     for dim in ds.dims:
-        if dim not in ds.coords and not dim.endswith("_term"):
+        if dim not in ds.coords and dim != TERM_DIM:
             ds.coords[dim] = arange(ds.sizes[dim])
 
     return ds
@@ -225,15 +225,25 @@ def print_single_expression(c, v, const, model):
     # catch case that to many terms would be printed
     def print_line(expr, const):
         res = []
-        for i, (coeff, (name, coord)) in enumerate(expr):
-            coord_string = print_coord(coord)
+        for i, (coeff, var) in enumerate(expr):
+            coeff_string = f"{coeff:+.4g}"
             if i:
-                coeff_string = f"{coeff:+.4g}"
                 # split sign and coefficient
-                res.append(f"{coeff_string[0]} {coeff_string[1:]} {name}{coord_string}")
+                coeff_string = f"{coeff_string[0]} {coeff_string[1:]}"
+
+            if isinstance(var, list):
+                var_string = ""
+                for name, coords in var:
+                    if name is not None:
+                        coord_string = print_coord(coords)
+                        var_string += f" {name}{coord_string}"
             else:
-                coeff_string = f"{coeff:+.4g} " if coeff != 1 else ""
-                res.append(f"{coeff_string}{name}{coord_string}")
+                name, coords = var
+                coord_string = print_coord(coords)
+                var_string = f" {name}{coord_string}"
+
+            res.append(f"{coeff_string}{var_string}")
+
         if const != 0:
             const_string = f"{const:+.4g}"
             if len(res):
@@ -242,14 +252,19 @@ def print_single_expression(c, v, const, model):
                 res.append(const_string)
         return " ".join(res) if len(res) else "None"
 
-    if isinstance(c, np.ndarray):
+    if v.ndim == 1:
         mask = v != -1
         c, v = c[mask], v[mask]
+    else:
+        mask = (v != -1).any(0)
+        c = c[mask]
+        v = v[:, mask]
 
     max_terms = options.get_value("display_max_terms")
     if len(c) > max_terms:
         truncate = max_terms // 2
-        expr = list(zip(c[:truncate], model.variables.get_label_position(v[:truncate])))
+        positions = model.variables.get_label_position(v[..., :truncate])
+        expr = list(zip(c[:truncate], positions))
         res = print_line(expr, const)
         res += " ... "
         expr = list(
diff --git a/linopy/constants.py b/linopy/constants.py
index 7406884..2438448 100644
--- a/linopy/constants.py
+++ b/linopy/constants.py
@@ -33,6 +33,13 @@ sign_replace_dict = {
     short_LESS_EQUAL: LESS_EQUAL,
 }
 
+TERM_DIM = "_term"
+STACKED_TERM_DIM = "_stacked_term"
+GROUPED_TERM_DIM = "_grouped_term"
+FACTOR_DIM = "_factor"
+CONCAT_DIM = "_concat"
+HELPER_DIMS = [TERM_DIM, STACKED_TERM_DIM, GROUPED_TERM_DIM, FACTOR_DIM, CONCAT_DIM]
+
 
 class ModelStatus(Enum):
     """
@@ -75,7 +82,7 @@ class SolverStatus(Enum):
         for status in STATUS_TO_TERMINATION_CONDITION_MAP:
             if termination_condition in STATUS_TO_TERMINATION_CONDITION_MAP[status]:
                 return status
-        return cls.unknown
+        return cls("unknown")
 
 
 class TerminationCondition(Enum):
diff --git a/linopy/constraints.py b/linopy/constraints.py
index cf4363b..759ab07 100644
--- a/linopy/constraints.py
+++ b/linopy/constraints.py
@@ -36,7 +36,14 @@ from linopy.common import (
     save_join,
 )
 from linopy.config import options
-from linopy.constants import EQUAL, GREATER_EQUAL, LESS_EQUAL, SIGNS_pretty
+from linopy.constants import (
+    EQUAL,
+    GREATER_EQUAL,
+    HELPER_DIMS,
+    LESS_EQUAL,
+    TERM_DIM,
+    SIGNS_pretty,
+)
 
 
 def conwrap(method, *default_args, **new_default_kwargs):
@@ -110,9 +117,7 @@ class Constraint:
 
         data = data.assign_attrs(name=name)
 
-        if "_term" in data.dims:
-            data = data.rename({"_term": f"{name}_term"})
-        (data,) = xr.broadcast(data, exclude=[f"{name}_term"])
+        (data,) = xr.broadcast(data, exclude=[TERM_DIM])
 
         self._data = data
         self._model = model
@@ -147,7 +152,7 @@ class Constraint:
 
     @property
     def coord_dims(self):
-        return {k: self.data.dims[k] for k in self.dims if not k.endswith("_term")}
+        return {k: self.data.dims[k] for k in self.dims if k not in HELPER_DIMS}
 
     @property
     def is_assigned(self):
@@ -225,7 +230,7 @@ class Constraint:
         """
         Return the term dimension of the constraint.
         """
-        return self.name + "_term"
+        return TERM_DIM
 
     @property
     def mask(self):
@@ -284,7 +289,7 @@ class Constraint:
         The function raises an error in case no model is set as a
         reference.
         """
-        data = self.data[["coeffs", "vars"]].rename({self.term_dim: "_term"})
+        data = self.data[["coeffs", "vars"]].rename({self.term_dim: TERM_DIM})
         return expressions.LinearExpression(data, self.model)
 
     @lhs.setter
@@ -292,10 +297,8 @@ class Constraint:
         value = expressions.as_expression(
             value, self.model, coords=self.coords, dims=self.coord_dims
         )
-        self._data = (
-            self.data.drop_vars(["coeffs", "vars"])
-            .assign(coeffs=value.coeffs, vars=value.vars, rhs=self.rhs - value.const)
-            .rename(_term=self.name + "_term")
+        self._data = self.data.drop_vars(["coeffs", "vars"]).assign(
+            coeffs=value.coeffs, vars=value.vars, rhs=self.rhs - value.const
         )
 
     @property
@@ -501,7 +504,7 @@ class Constraints:
             coords = " (" + ", ".join(ds.coords) + ")" if ds.coords else ""
             r += f" * {name}{coords}\n"
         if not len(list(self)):
-            r += "<empty>"
+            r += "<empty>\n"
         return r
 
     def __getitem__(
@@ -806,7 +809,10 @@ class Constraints:
         -------
         pd.DataFrame
         """
-        df = pd.concat([self[k].flat for k in self], ignore_index=True)
+        dfs = [self[k].flat for k in self]
+        if not len(dfs):
+            return pd.DataFrame(columns=["coeffs", "vars", "labels", "key"])
+        df = pd.concat(dfs, ignore_index=True)
         unique_labels = df.labels.unique()
         map_labels = pd.Series(np.arange(len(unique_labels)), index=unique_labels)
         df["key"] = df.labels.map(map_labels)
@@ -822,13 +828,16 @@ class Constraints:
         # TODO: rename "filter_missings" to "~labels_as_coordinates"
         cons = self.flat
 
+        if not len(self):
+            return None
+
         if filter_missings:
             vars = self.model.variables.flat
             shape = (cons.key.max() + 1, vars.key.max() + 1)
             cons["vars"] = cons.vars.map(vars.set_index("labels").key)
             return csc_matrix((cons.coeffs, (cons.key, cons.vars)), shape=shape)
         else:
-            shape = (self.model._cCounter, self.model._xCounter)
+            shape = self.model.shape
             return csc_matrix((cons.coeffs, (cons.labels, cons.vars)), shape=shape)
 
     def reset_dual(self):
diff --git a/linopy/expressions.py b/linopy/expressions.py
index ba298c4..8b82367 100644
--- a/linopy/expressions.py
+++ b/linopy/expressions.py
@@ -19,6 +19,7 @@ import xarray.core.groupby
 import xarray.core.rolling
 from deprecation import deprecated
 from numpy import arange, array, nan
+from scipy.sparse import csc_matrix
 from xarray import DataArray, Dataset
 from xarray.core.dataarray import DataArrayCoordinates
 
@@ -34,7 +35,16 @@ from linopy.common import (
     print_single_expression,
 )
 from linopy.config import options
-from linopy.constants import EQUAL, GREATER_EQUAL, LESS_EQUAL
+from linopy.constants import (
+    EQUAL,
+    FACTOR_DIM,
+    GREATER_EQUAL,
+    GROUPED_TERM_DIM,
+    HELPER_DIMS,
+    LESS_EQUAL,
+    STACKED_TERM_DIM,
+    TERM_DIM,
+)
 
 
 def exprwrap(method, *default_args, **new_default_kwargs):
@@ -154,10 +164,12 @@ class LinearExpressionGroupby:
 
             group_dim = group.index.name
             arrays = [group, group.groupby(group).cumcount()]
-            idx = pd.MultiIndex.from_arrays(arrays, names=[group_name, "_grouped_term"])
+            idx = pd.MultiIndex.from_arrays(
+                arrays, names=[group_name, GROUPED_TERM_DIM]
+            )
             ds = self.data.assign_coords({group_dim: idx})
             ds = ds.unstack(group_dim, fill_value=LinearExpression.fill_value)
-            ds = LinearExpression._sum(ds, dims="_grouped_term")
+            ds = LinearExpression._sum(ds, dims=GROUPED_TERM_DIM)
 
             if int_map is not None:
                 index = ds.indexes["group"].map({v: k for k, v in int_map.items()})
@@ -171,7 +183,7 @@ class LinearExpressionGroupby:
 
             def func(ds):
                 ds = LinearExpression._sum(ds, self.groupby._group_dim)
-                ds = ds.assign_coords(_term=np.arange(len(ds._term)))
+                ds = ds.assign_coords({TERM_DIM: np.arange(len(ds._term))})
                 return ds
 
             return self.map(func, **kwargs)
@@ -205,15 +217,14 @@ class LinearExpressionRolling:
         data = self.rolling.construct("_rolling_term", keep_attrs=True)
         ds = (
             data[["coeffs", "vars"]]
-            .rename(_term="_stacked_term")
-            .stack(_term=["_stacked_term", "_rolling_term"])
-            .reset_index("_term", drop=True)
+            .rename({TERM_DIM: STACKED_TERM_DIM})
+            .stack({TERM_DIM: [STACKED_TERM_DIM, "_rolling_term"]}, create_index=False)
         )
         ds["const"] = data.const.sum("_rolling_term")
         return LinearExpression(ds, self.model)
 
 
-@forward_as_properties(data=["attrs", "coords", "indexes", "sizes"])
+@forward_as_properties(data=["attrs", "coords", "indexes", "sizes"], const=["ndim"])
 class LinearExpression:
     """
     A linear expression consisting of terms of coefficients and variables.
@@ -262,12 +273,12 @@ class LinearExpression:
         from linopy.model import Model
 
         if data is None:
-            da = xr.DataArray([], dims=["_term"])
+            da = xr.DataArray([], dims=[TERM_DIM])
             data = Dataset({"coeffs": da, "vars": da, "const": 0})
         elif isinstance(data, DataArray):
             # assume only constant are passed
             const = fill_missing_coords(data)
-            da = xr.DataArray([], dims=["_term"])
+            da = xr.DataArray([], dims=[TERM_DIM])
             data = Dataset({"coeffs": da, "vars": da, "const": const})
         elif not isinstance(data, Dataset):
             raise ValueError(
@@ -286,19 +297,19 @@ class LinearExpression:
 
         data = fill_missing_coords(data)
 
-        term_dims = [d for d in data.dims if d.endswith("_term")]
-        if not term_dims:
+        if TERM_DIM not in data.dims:
             raise ValueError("data must contain one dimension ending with '_term'")
-        term_dim = term_dims[0]
 
         if "const" not in data:
             data = data.assign(const=0)
 
-        data = xr.broadcast(data, exclude=["_term"])[0]
-        data[["coeffs", "vars"]] = xr.broadcast(data[["coeffs", "vars"]])[0]
+        data = xr.broadcast(data, exclude=HELPER_DIMS)[0]
+        data[["coeffs", "vars"]] = xr.broadcast(
+            data[["coeffs", "vars"]], exclude=[FACTOR_DIM]
+        )[0]
 
         # transpose with new Dataset to really ensure correct order
-        data = Dataset(data.transpose(..., term_dim))
+        data = Dataset(data.transpose(..., TERM_DIM))
 
         if not isinstance(model, Model):
             raise ValueError("model must be an instance of linopy.Model")
@@ -311,13 +322,13 @@ class LinearExpression:
         Print the expression arrays.
         """
         max_lines = options["display_max_rows"]
-        dims = list(self.dims)
-        dim_sizes = list(self.sizes.values())[:-1]
+        dims = list(self.coord_dims)
+        dim_sizes = list(self.coord_dims.values())
         size = np.prod(dim_sizes)  # that the number of theoretical printouts
         masked_entries = self.mask.sum().values if self.mask is not None else 0
         lines = []
 
-        header_string = "LinearExpression"
+        header_string = self.type
 
         if size > 1:
             for indices in generate_indices_for_printout(dim_sizes, max_lines):
@@ -361,7 +372,7 @@ class LinearExpression:
         other = as_expression(
             other, model=self.model, coords=self.coords, dims=self.coord_dims
         )
-        return merge(self, other)
+        return merge(self, other, cls=self.__class__)
 
     def __radd__(self, other):
         # This is needed for using python's sum function
@@ -374,7 +385,7 @@ class LinearExpression:
         other = as_expression(
             other, model=self.model, coords=self.coords, dims=self.coord_dims
         )
-        return merge(self, -other)
+        return merge(self, -other, cls=self.__class__)
 
     def __neg__(self):
         """
@@ -386,19 +397,27 @@ class LinearExpression:
         """
         Multiply the expr by a factor.
         """
-        if isinstance(
-            other, (LinearExpression, variables.Variable, variables.ScalarVariable)
-        ):
+        if isinstance(other, QuadraticExpression):
             raise TypeError(
                 "unsupported operand type(s) for *: "
                 f"{type(self)} and {type(other)}. "
-                "Non-linear expressions are not yet supported."
+                "Higher order non-linear expressions are not yet supported."
             )
-        coeffs = self.coeffs * as_dataarray(
-            other, coords=self.coords, dims=self.coord_dims
-        )
-        assert set(coeffs.shape) == set(self.coeffs.shape)
-        return self.assign(coeffs=coeffs)
+        elif isinstance(other, (variables.Variable, variables.ScalarVariable)):
+            other = other.to_linexpr()
+
+        if type(other) is LinearExpression:
+            if other.nterm > 1:
+                raise TypeError("Multiplication of multiple terms is not supported.")
+            ds = other.data[["coeffs", "vars"]].sel(_term=0).broadcast_like(self.data)
+            ds = ds.assign(const=other.const)
+            return merge(self, ds, dim=FACTOR_DIM, cls=QuadraticExpression)
+        else:
+            coeffs = self.coeffs * as_dataarray(
+                other, coords=self.coords, dims=self.coord_dims
+            )
+            assert set(coeffs.shape) == set(self.coeffs.shape)
+            return self.assign(coeffs=coeffs)
 
     def __rmul__(self, other):
         """
@@ -448,6 +467,10 @@ class LinearExpression:
     def fill_value(cls):
         return cls._fill_value
 
+    @property
+    def type(self):
+        return "LinearExpression"
+
     @property
     def data(self):
         return self._data
@@ -468,7 +491,7 @@ class LinearExpression:
 
     @property
     def coord_dims(self):
-        return {k: self.data.dims[k] for k in self.dims if not k.endswith("_term")}
+        return {k: self.data.dims[k] for k in self.dims if k not in HELPER_DIMS}
 
     @property
     def vars(self):
@@ -504,18 +527,17 @@ class LinearExpression:
         data = _expr_unwrap(expr)
 
         if dims is None:
-            vars = DataArray(data.vars.data.ravel(), dims="_term")
-            coeffs = DataArray(data.coeffs.data.ravel(), dims="_term")
+            vars = DataArray(data.vars.data.ravel(), dims=TERM_DIM)
+            coeffs = DataArray(data.coeffs.data.ravel(), dims=TERM_DIM)
             const = data.const.sum()
             ds = xr.Dataset({"vars": vars, "coeffs": coeffs, "const": const})
         else:
-            dims = [d for d in np.atleast_1d(dims) if d != "_term"]
+            dims = [d for d in np.atleast_1d(dims) if d != TERM_DIM]
             ds = (
                 data[["coeffs", "vars"]]
                 .reset_index(dims, drop=True)
-                .rename(_term="_stacked_term")
-                .stack(_term=["_stacked_term"] + dims)
-                .reset_index("_term", drop=True)
+                .rename({TERM_DIM: STACKED_TERM_DIM})
+                .stack({TERM_DIM: [STACKED_TERM_DIM] + dims}, create_index=False)
             )
             ds["const"] = data.const.sum(dims)
 
@@ -682,12 +704,20 @@ class LinearExpression:
         coeffs = coeffs.reshape((nterm, *shape))
         vars = vars.reshape((nterm, *shape))
 
-        coeffs = DataArray(coeffs, coords, dims=("_term", *coords))
-        vars = DataArray(vars, coords, dims=("_term", *coords))
-        ds = Dataset({"coeffs": coeffs, "vars": vars}).transpose(..., "_term")
+        coeffs = DataArray(coeffs, coords, dims=(TERM_DIM, *coords))
+        vars = DataArray(vars, coords, dims=(TERM_DIM, *coords))
+        ds = Dataset({"coeffs": coeffs, "vars": vars}).transpose(..., TERM_DIM)
 
         return cls(ds, model)
 
+    def to_quadexpr(self):
+        """Convert LinearExpression to QuadraticExpression."""
+        vars = self.data.vars.expand_dims(FACTOR_DIM)
+        fill_value = self.fill_value["vars"]
+        vars = xr.concat([vars, xr.full_like(vars, fill_value)], dim=FACTOR_DIM)
+        data = self.data.assign(vars=vars)
+        return QuadraticExpression(data, self.model)
+
     def to_constraint(self, sign, rhs):
         """
         Convert a linear expression to a constraint.
@@ -852,7 +882,6 @@ class LinearExpression:
         """
         Get the total size of the linear expression.
         """
-        assert self.vars.size == self.coeffs.size
         return self.vars.size
 
     def empty(self):
@@ -866,7 +895,7 @@ class LinearExpression:
         Move all non-zero term entries to the front and cut off all-zero
         entries in the term-axis.
         """
-        data = self.data.transpose(..., "_term")
+        data = self.data.transpose(..., TERM_DIM)
 
         cdata = data.coeffs.data
         axis = cdata.ndim - 1
@@ -890,7 +919,7 @@ class LinearExpression:
         cdata[tuple(mod_nnz)] = data.coeffs.data[nnz]
         data.coeffs.data = cdata
 
-        return self.__class__(data.sel(_term=slice(0, nterm)), self.model)
+        return self.__class__(data.sel({TERM_DIM: slice(0, nterm)}), self.model)
 
     def sanitize(self):
         """
@@ -924,9 +953,6 @@ class LinearExpression:
         -------
         df : pandas.DataFrame
         """
-        # if exclude_lhs:
-        #     ds = self.data.drop_vars(['coeffs', 'vars'])
-        # else:
         ds = self.data
         if not ds.sizes:
             # fallback for weird error raised due to missing index
@@ -991,9 +1017,173 @@ class LinearExpression:
     roll = exprwrap(Dataset.roll)
 
 
+@forward_as_properties(data=["attrs", "coords", "indexes", "sizes"])
+class QuadraticExpression(LinearExpression):
+    """
+    A quadratic expression consisting of terms of coefficients and variables.
+
+    The QuadraticExpression class is a subclass of LinearExpression which allows to
+    apply most xarray functions on it.
+    """
+
+    __slots__ = ("_data", "_model")
+    __array_ufunc__ = None
+    __array_priority__ = 10000
+
+    _fill_value = {"vars": -1, "coeffs": np.nan, "const": 0}
+
+    def __init__(self, data, model):
+        super().__init__(data, model)
+
+        if FACTOR_DIM not in data.vars.dims:
+            raise ValueError(f"Data does not include dimension {FACTOR_DIM}")
+        elif data.sizes[FACTOR_DIM] != 2:
+            raise ValueError(f"Size of dimension {FACTOR_DIM} must be 2.")
+
+        # transpose data to have _term as last dimension and _factor as second last
+        data = xr.Dataset(data.transpose(..., FACTOR_DIM, TERM_DIM))
+        self._data = data
+
+    def __mul__(self, other):
+        """
+        Multiply the expr by a factor.
+        """
+        if isinstance(
+            other,
+            (
+                LinearExpression,
+                QuadraticExpression,
+                variables.Variable,
+                variables.ScalarVariable,
+            ),
+        ):
+            raise TypeError(
+                "unsupported operand type(s) for *: "
+                f"{type(self)} and {type(other)}. "
+                "Higher order non-linear expressions are not yet supported."
+            )
+        return super().__mul__(other)
+
+    @property
+    def type(self):
+        return "QuadraticExpression"
+
+    def __add__(self, other):
+        """
+        Add an expression to others.
+        """
+        other = as_expression(
+            other, model=self.model, coords=self.coords, dims=self.coord_dims
+        )
+        if type(other) is LinearExpression:
+            other = other.to_quadexpr()
+        return merge(self, other, cls=self.__class__)
+
+    def __radd__(self, other):
+        """
+        Add others to expression.
+        """
+        if type(other) is LinearExpression:
+            other = other.to_quadexpr()
+            return other.__add__(self)
+        elif other == 0:
+            return self
+        else:
+            return NotImplemented
+
+    def __sub__(self, other):
+        """
+        Subtract others from expression.
+        """
+        other = as_expression(
+            other, model=self.model, coords=self.coords, dims=self.coord_dims
+        )
+        if type(other) is LinearExpression:
+            other = other.to_quadexpr()
+        return merge(self, -other, cls=self.__class__)
+
+    def __rsub__(self, other):
+        """
+        Subtract expression from others.
+        """
+        if type(other) is LinearExpression:
+            other = other.to_quadexpr()
+            return other.__sub__(self)
+        else:
+            NotImplemented
+
+    @classmethod
+    def _sum(cls, expr: "QuadraticExpression", dims=None) -> Dataset:
+        data = _expr_unwrap(expr)
+        dims = dims or list(set(data.dims) - set(HELPER_DIMS))
+        return LinearExpression._sum(expr, dims)
+
+    def to_constraint(self, sign, rhs):
+        raise NotImplementedError(
+            "Quadratic expressions cannot be used in constraints."
+        )
+
+    @property
+    def flat(self):
+        """
+        Return a flattened expression.
+        """
+        vars = self.data.vars.assign_coords(
+            {FACTOR_DIM: ["vars1", "vars2"]}
+        ).to_dataset(FACTOR_DIM)
+        ds = self.data.drop_vars("vars").assign(vars)
+        if not ds.sizes:
+            # fallback for weird error raised due to missing index
+            df = pd.DataFrame({k: ds[k].item() for k in ds}, index=[0])
+        else:
+            df = ds.to_dataframe()
+        if "mask" in df:
+            mask = df.pop("mask")
+            df = df[mask]
+        df = df[((df.vars1 != -1) | (df.vars2 != -1)) & (df.coeffs != 0)]
+        # Group repeated variables in the same constraint
+        df = df.groupby(["vars1", "vars2"], as_index=False).sum()
+
+        any_nan = df.isna().any()
+        if any_nan.any():
+            fields = ", ".join("`" + df.columns[any_nan] + "`")
+            raise ValueError(
+                f"Expression `{self.name}` contains nan's in field(s) {fields}"
+            )
+
+        return df
+
+    def to_matrix(self):
+        """
+        Return a sparse matrix representation of the expression only including
+        quadratic terms.
+
+        Note that the matrix is formulated following the convention of the
+        optimization problem, i.e. the quadratic term is 0.5 x^T Q x.
+        The matrix Q is therefore symmetric and the diagonal terms are doubled.
+
+        """
+        df = self.flat
+        # drop linear terms
+        df = df[(df.vars1 != -1) & (df.vars2 != -1)]
+
+        # symmetrize cross terms and double diagonal terms
+        cross_terms = df.vars1 != df.vars2
+        df.loc[~cross_terms, "coeffs"] *= 2
+        vals = dict(vars1=df.vars2[cross_terms], vars2=df.vars1[cross_terms])
+        df = pd.concat([df, df[cross_terms].assign(**vals)])
+
+        # assign matrix
+        data = df.coeffs
+        row = df.vars1
+        col = df.vars2
+        nvars = self.model.shape[1]
+        return csc_matrix((data, (row, col)), shape=(nvars, nvars))
+
+
 def as_expression(obj, model=None, **kwargs):
     """
-    Convert an object to a LinearExpression.
+    Convert an object to a LinearExpression or QuadraticExpression.
 
     Parameters
     ----------
@@ -1025,13 +1215,13 @@ def as_expression(obj, model=None, **kwargs):
         return LinearExpression(obj, model)
 
 
-def merge(*exprs, dim="_term", cls=LinearExpression, **kwargs):
+def merge(*exprs, dim=TERM_DIM, cls=LinearExpression, **kwargs):
     """
-    Merge multiple linear expression together.
+    Merge multiple expression together.
 
     This function is a bit faster than summing over multiple linear expressions.
     In case a list of LinearExpression with exactly the same shape is passed
-    and the dimension to concatenate on is "_term", the concatenation uses
+    and the dimension to concatenate on is TERM_DIM, the concatenation uses
     the coordinates of the first object as a basis which overrides the
     coordinates of the consecutive objects.
 
@@ -1053,17 +1243,23 @@ def merge(*exprs, dim="_term", cls=LinearExpression, **kwargs):
     -------
     res : linopy.LinearExpression
     """
+    linopy_types = (variables.Variable, LinearExpression, QuadraticExpression)
+
     exprs = exprs[0] if len(exprs) == 1 else list(exprs)  # allow passing a list
     model = exprs[0].model
 
-    if cls == LinearExpression and dim == "_term":
-        override = check_common_keys_values([e.coord_dims for e in exprs])
+    if cls in linopy_types and dim in HELPER_DIMS:
+        coord_dims = [
+            {k: v for k, v in e.dims.items() if k not in HELPER_DIMS} for e in exprs
+        ]
+        override = check_common_keys_values(coord_dims)
     else:
         override = False
 
-    data = [e.data if isinstance(e, cls) else e for e in exprs]
-    if cls == LinearExpression and len({ds.dims["_term"] for ds in data}) > 1:
-        data = [ds.assign_coords(_term=arange(ds.dims["_term"])) for ds in data]
+    data = [e.data if isinstance(e, linopy_types) else e for e in exprs]
+    for d in set(HELPER_DIMS) & set(_ for ds in data for _ in ds.dims):
+        if len({ds.dims[d] for ds in data}) > 1:
+            data = [ds.assign_coords({d: arange(ds.dims[d])}) for ds in data]
 
     if not kwargs:
         kwargs = {
@@ -1074,15 +1270,21 @@ def merge(*exprs, dim="_term", cls=LinearExpression, **kwargs):
         if override:
             kwargs["join"] = "override"
 
-    if dim == "_term":
+    if dim == TERM_DIM:
         ds = xr.concat([d[["coeffs", "vars"]] for d in data], dim, **kwargs)
-        const = xr.concat([d["const"] for d in data], dim, **kwargs).sum("_term")
+        const = xr.concat([d["const"] for d in data], dim, **kwargs).sum(TERM_DIM)
+        ds["const"] = const
+    elif dim == FACTOR_DIM:
+        ds = xr.concat([d[["vars"]] for d in data], dim, **kwargs)
+        coeffs = xr.concat([d["coeffs"] for d in data], dim, **kwargs).prod(FACTOR_DIM)
+        ds["coeffs"] = coeffs
+        const = xr.concat([d["const"] for d in data], dim, **kwargs).prod(FACTOR_DIM)
         ds["const"] = const
     else:
         ds = xr.concat(data, dim, **kwargs)
 
-    if "_term" in ds.coords:
-        ds = ds.reset_index("_term", drop=True)
+    for d in set(HELPER_DIMS) & set(ds.coords):
+        ds = ds.reset_index(d, drop=True)
 
     return cls(ds, model)
 
@@ -1219,7 +1421,7 @@ class ScalarLinearExpression:
         )
 
     def to_linexpr(self):
-        coeffs = xr.DataArray(list(self.coeffs), dims="_term")
-        vars = xr.DataArray(list(self.vars), dims="_term")
+        coeffs = xr.DataArray(list(self.coeffs), dims=TERM_DIM)
+        vars = xr.DataArray(list(self.vars), dims=TERM_DIM)
         ds = xr.Dataset({"coeffs": coeffs, "vars": vars})
         return LinearExpression(ds, self.model)
diff --git a/linopy/io.py b/linopy/io.py
index 8c58ea2..936fa68 100644
--- a/linopy/io.py
+++ b/linopy/io.py
@@ -13,21 +13,79 @@ import numpy as np
 import pandas as pd
 import xarray as xr
 from numpy import asarray, concatenate, ones_like, zeros_like
+from scipy.sparse import csc_matrix, triu
 from tqdm import tqdm
 
 from linopy import solvers
-from linopy.constants import EQUAL, GREATER_EQUAL
+from linopy.constants import CONCAT_DIM, EQUAL, GREATER_EQUAL
 
 logger = logging.getLogger(__name__)
 
 
 ufunc_kwargs = dict(vectorize=True)
-concat_dim = "_concat_dim"
-concat_kwargs = dict(dim=concat_dim, coords="minimal")
+concat_kwargs = dict(dim=CONCAT_DIM, coords="minimal")
 
 TQDM_COLOR = "#80bfff"
 
 
+def handle_batch(batch, f, batch_size):
+    """
+    Write out a batch to a file and reset the batch.
+    """
+    if len(batch) >= batch_size:
+        f.writelines(batch)  # write out a batch
+        batch = []  # reset batch
+    return batch
+
+
+def objective_write_linear_terms(df, f, batch, batch_size):
+    """
+    Write the linear terms of the objective to a file.
+    """
+    coeffs = df.coeffs.values
+    vars = df.vars.values
+    for idx in range(len(coeffs)):
+        coeff = coeffs[idx]
+        var = vars[idx]
+        batch.append(f"{coeff:+.12g} x{var}\n")
+        batch = handle_batch(batch, f, batch_size)
+    return batch
+
+
+def objective_write_cross_terms(quadratic, f, batch, batch_size):
+    """
+    Write the cross terms of the objective to a file.
+    """
+    is_cross = quadratic.vars1 != quadratic.vars2
+    cross = quadratic[is_cross]
+    coeffs = cross.coeffs.values
+    vars1 = cross.vars1.values
+    vars2 = cross.vars2.values
+    for idx in range(len(coeffs)):
+        coeff = coeffs[idx]
+        var1 = vars1[idx]
+        var2 = vars2[idx]
+        batch.append(f"{coeff:+.12g} x{var1} * x{var2}\n")
+        batch = handle_batch(batch, f, batch_size)
+    return batch
+
+
+def objective_write_quad_terms(quadratic, f, batch, batch_size):
+    """
+    Write the quadratic terms of the objective to a file.
+    """
+    is_cross = quadratic.vars1 != quadratic.vars2
+    quad = quadratic[~is_cross]
+    coeffs = quad.coeffs.values
+    vars = quad.vars1.values
+    for idx in range(len(coeffs)):
+        coeff = coeffs[idx]
+        var = vars[idx]
+        batch.append(f"{coeff:+.12g} x{var} ^ 2\n")
+        batch = handle_batch(batch, f, batch_size)
+    return batch
+
+
 def objective_to_file(m, f, log=False, batch_size=10000):
     """
     Write out the objective of a model to a lp file.
@@ -39,29 +97,37 @@ def objective_to_file(m, f, log=False, batch_size=10000):
     df = m.objective.flat
 
     if np.isnan(df.coeffs).any():
-        raise ValueError(
+        logger.warning(
             "Objective coefficients are missing (nan) where variables are not (-1)."
         )
 
-    coeffs = df.coeffs.values
-    vars = df.vars.values
+    if m.is_linear:
+        batch = objective_write_linear_terms(df, f, [], batch_size)
 
-    batch = []  # to store batch of lines
-    for idx in range(len(df)):
-        coeff = coeffs[idx]
-        var = vars[idx]
-
-        batch.append(f"{coeff:+.12g} x{var}\n")
+    elif m.is_quadratic:
+        is_linear = (df.vars1 == -1) | (df.vars2 == -1)
+        linear = df[is_linear]
+        linear = linear.assign(
+            vars=linear.vars1.where(linear.vars1 != -1, linear.vars2)
+        )
+        batch = objective_write_linear_terms(linear, f, [], batch_size)
 
-        if len(batch) >= batch_size:
-            f.writelines(batch)  # write out a batch
-            batch = []  # reset batch
+        if not is_linear.all():
+            batch.append("+ [\n")
+            quadratic = df[~is_linear]
+            quadratic = quadratic.assign(coeffs=2 * quadratic.coeffs)
+            batch = objective_write_cross_terms(quadratic, f, batch, batch_size)
+            batch = objective_write_quad_terms(quadratic, f, batch, batch_size)
+            batch.append("] / 2\n")
 
     if batch:  # write the remaining lines
         f.writelines(batch)
 
 
 def constraints_to_file(m, f, log=False, batch_size=50000):
+    if not len(m.constraints):
+        return
+
     f.write("\n\ns.t.\n\n")
     names = m.constraints
     if log:
@@ -109,9 +175,7 @@ def constraints_to_file(m, f, log=False, batch_size=50000):
             else:
                 batch.append(f"{coeff:+.12g} x{var}\n")
 
-            if len(batch) >= batch_size:
-                f.writelines(batch)
-                batch = []  # reset batch
+            batch = handle_batch(batch, f, batch_size)
 
             prev_label = label
 
@@ -150,10 +214,7 @@ def bounds_to_file(m, f, log=False, batch_size=10000):
             lower = lowers[idx]
             upper = uppers[idx]
             batch.append(f"{lower:+.12g} <= x{label} <= {upper:+.12g}\n")
-
-            if len(batch) >= batch_size:
-                f.writelines(batch)  # write out a batch
-                batch = []  # reset batch
+            batch = handle_batch(batch, f, batch_size)
 
     if batch:  # write the remaining lines
         f.writelines(batch)
@@ -181,10 +242,7 @@ def binaries_to_file(m, f, log=False, batch_size=1000):
 
         for label in df.labels.values:
             batch.append(f"x{label}\n")
-
-            if len(batch) >= batch_size:
-                f.writelines(batch)  # write out a batch
-                batch = []  # reset batch
+            batch = handle_batch(batch, f, batch_size)
 
     if batch:  # write the remaining lines
         f.writelines(batch)
@@ -212,10 +270,7 @@ def integers_to_file(m, f, log=False, batch_size=1000):
 
         for label in df.labels.values:
             batch.append(f"x{label}\n")
-
-            if len(batch) >= batch_size:
-                f.writelines(batch)  # write out a batch
-                batch = []  # reset batch
+            batch = handle_batch(batch, f, batch_size)
 
     if batch:  # write the remaining lines
         f.writelines(batch)
@@ -293,11 +348,15 @@ def to_gurobipy(m):
         kwargs["vtype"] = M.vtypes
     x = model.addMVar(M.vlabels.shape, M.lb, M.ub, name=list(names), **kwargs)
 
-    model.setObjective(M.c @ x)
+    if m.is_quadratic:
+        model.setObjective(0.5 * x.T @ M.Q @ x + M.c @ x)
+    else:
+        model.setObjective(M.c @ x)
 
-    names = "c" + M.clabels.astype(str).astype(object)
-    c = model.addMConstr(M.A, x, M.sense, M.b)
-    c.setAttr("ConstrName", list(names))
+    if len(m.constraints):
+        names = "c" + M.clabels.astype(str).astype(object)
+        c = model.addMConstr(M.A, x, M.sense, M.b)
+        c.setAttr("ConstrName", list(names))
 
     model.update()
     return model
@@ -334,16 +393,33 @@ def to_highspy(m):
             labels = np.arange(len(vtypes))[vtypes == "B"]
             n = len(labels)
             h.changeColsBounds(n, labels, zeros_like(labels), ones_like(labels))
+
+    # linear objective
     h.changeColsCost(len(M.c), np.arange(len(M.c), dtype=np.int32), M.c)
-    A = M.A.tocsr()
-    num_cons = A.shape[0]
-    lower = np.where(M.sense != "<", M.b, -np.inf)
-    upper = np.where(M.sense != ">", M.b, np.inf)
-    h.addRows(num_cons, lower, upper, A.nnz, A.indptr, A.indices, A.data)
+
+    # linear constraints
+    A = M.A
+    if A is not None:
+        A = A.tocsr()
+        num_cons = A.shape[0]
+        lower = np.where(M.sense != "<", M.b, -np.inf)
+        upper = np.where(M.sense != ">", M.b, np.inf)
+        h.addRows(num_cons, lower, upper, A.nnz, A.indptr, A.indices, A.data)
+
     lp = h.getLp()
-    lp.row_names_ = "c" + M.clabels.astype(str).astype(object)
     lp.col_names_ = "x" + M.vlabels.astype(str).astype(object)
+    if len(M.clabels):
+        lp.row_names_ = "c" + M.clabels.astype(str).astype(object)
     h.passModel(lp)
+
+    # quadrative objective
+    Q = M.Q
+    if Q is not None:
+        Q = triu(Q)
+        Q = Q.tocsr()
+        num_vars = Q.shape[0]
+        h.passHessian(num_vars, Q.nnz, 1, Q.indptr, Q.indices, Q.data)
+
     return h
 
 
diff --git a/linopy/matrices.py b/linopy/matrices.py
index e7062e0..f206043 100644
--- a/linopy/matrices.py
+++ b/linopy/matrices.py
@@ -11,10 +11,14 @@ from functools import cached_property
 import numpy as np
 import pandas as pd
 
+from linopy import expressions
 
-def create_vector(indices, values, fill_value=np.nan):
+
+def create_vector(indices, values, fill_value=np.nan, shape=None):
     """Create a vector of a size equal to the maximum index plus one."""
-    vector = np.full(max(indices) + 1, fill_value)
+    if shape is None:
+        shape = max(indices) + 1
+    vector = np.full(shape, fill_value)
     vector[indices] = values
     return vector
 
@@ -60,7 +64,7 @@ class MatrixAccessor:
 
         ds = pd.concat(specs)
         ds = df.set_index("key").labels.map(ds)
-        return create_vector(ds.index, ds.values, "")
+        return create_vector(ds.index, ds.values, fill_value="")
 
     @property
     def lb(self):
@@ -78,21 +82,22 @@ class MatrixAccessor:
     def clabels(self):
         "Vector of labels of all non-missing constraints."
         df = self.flat_cons
-        return create_vector(df.key, df.labels, -1)
+        if df.empty:
+            return np.array([], dtype=int)
+        return create_vector(df.key, df.labels, fill_value=-1)
 
     @property
     def A(self):
         "Constraint matrix of all non-missing constraints and variables."
         m = self._parent
-        return m.constraints.to_matrix(filter_missings=False)[self.clabels][
-            :, self.vlabels
-        ]
+        A = m.constraints.to_matrix(filter_missings=False)
+        return A[self.clabels][:, self.vlabels] if A is not None else None
 
     @property
     def sense(self):
         "Vector of senses of all non-missing constraints."
         df = self.flat_cons
-        return create_vector(df.key, df.sign.astype(np.dtype("<U1")), "")
+        return create_vector(df.key, df.sign.astype(np.dtype("<U1")), fill_value="")
 
     @property
     def b(self):
@@ -105,5 +110,19 @@ class MatrixAccessor:
         "Vector of objective coefficients of all non-missing variables."
         m = self._parent
         ds = m.objective.flat
+        if isinstance(m.objective, expressions.QuadraticExpression):
+            ds = ds[(ds.vars1 == -1) | (ds.vars2 == -1)]
+            ds["vars"] = ds.vars1.where(ds.vars1 != -1, ds.vars2)
+
         vars = ds.vars.map(self.flat_vars.set_index("labels").key)
-        return create_vector(vars, ds.coeffs, 0)
+        shape = self.flat_vars.key.max() + 1
+        return create_vector(vars, ds.coeffs, fill_value=0, shape=shape)
+
+    @property
+    def Q(self):
+        "Matrix objective coefficients of quadratic terms of all non-missing variables."
+        m = self._parent
+        if m.is_linear:
+            return None
+
+        return m.objective.to_matrix()[self.vlabels][:, self.vlabels]
diff --git a/linopy/model.py b/linopy/model.py
index a6bd2b8..813c304 100644
--- a/linopy/model.py
+++ b/linopy/model.py
@@ -26,12 +26,16 @@ from linopy.common import (
     replace_by_map,
     save_join,
 )
-from linopy.constants import ModelStatus, TerminationCondition
+from linopy.constants import TERM_DIM, ModelStatus, TerminationCondition
 from linopy.constraints import AnonymousScalarConstraint, Constraint, Constraints
-from linopy.expressions import LinearExpression, ScalarLinearExpression
+from linopy.expressions import (
+    LinearExpression,
+    QuadraticExpression,
+    ScalarLinearExpression,
+)
 from linopy.io import to_block_files, to_file, to_gurobipy, to_highspy, to_netcdf
 from linopy.matrices import MatrixAccessor
-from linopy.solvers import available_solvers
+from linopy.solvers import available_solvers, quadratic_solvers
 from linopy.variables import ScalarVariable, Variable, Variables
 
 logger = logging.getLogger(__name__)
@@ -286,8 +290,10 @@ class Model:
         """
         var_string = self.variables.__repr__().split("\n", 2)[2]
         con_string = self.constraints.__repr__().split("\n", 2)[2]
+        model_string = f"Linopy {self.type} model"
+
         return (
-            f"Linopy model\n============\n\n"
+            f"{model_string}\n{'=' * len(model_string)}\n\n"
             f"Variables:\n----------\n{var_string}\n"
             f"Constraints:\n------------\n{con_string}\n"
             f"Status:\n-------\n{self.status}"
@@ -603,13 +609,22 @@ class Model:
 
         if isinstance(expr, (list, tuple)):
             expr = self.linexpr(*expr)
-        assert isinstance(expr, LinearExpression)
+
+        if not isinstance(expr, (LinearExpression, QuadraticExpression)):
+            raise ValueError(
+                f"Invalid type of `expr` ({type(expr)})."
+                " Must be a LinearExpression or QuadraticExpression."
+            )
 
         if self.chunk is not None:
             expr = expr.chunk(self.chunk)
 
-        if expr.vars.ndim > 1:
+        if len(expr.coord_dims):
             expr = expr.sum()
+
+        if expr.const != 0:
+            raise ValueError("Constant values in objective function not supported.")
+
         self._objective = expr
         return self._objective
 
@@ -637,7 +652,9 @@ class Model:
             vars = vars.where(~vars.isin(labels), -1)
             self.constraints[k].data["vars"] = vars
 
-        self.objective = self.objective.sel(_term=~self.objective.vars.isin(labels))
+        self.objective = self.objective.sel(
+            {TERM_DIM: ~self.objective.vars.isin(labels)}
+        )
 
     def remove_constraints(self, name):
         """
@@ -656,6 +673,13 @@ class Model:
         """
         self.constraints.remove(name)
 
+    @property
+    def continuous(self):
+        """
+        Get all continuous variables.
+        """
+        return self.variables.continuous
+
     @property
     def binaries(self):
         """
@@ -670,10 +694,33 @@ class Model:
         """
         return self.variables.integers
 
+    @property
+    def is_linear(self):
+        return type(self.objective) is LinearExpression
+
+    @property
+    def is_quadratic(self):
+        return type(self.objective) is QuadraticExpression
+
+    @property
+    def type(self):
+        if (len(self.binaries) or len(self.integers)) and len(self.continuous):
+            variable_type = "MI"
+        elif len(self.binaries) or len(self.integers):
+            variable_type = "I"
+        else:
+            variable_type = ""
+
+        objective_type = "Q" if self.is_quadratic else "L"
+
+        return f"{variable_type}{objective_type}P"
+
     @property
     def nvars(self):
         """
         Get the total number of variables.
+
+        This excludes all variables which are not active.
         """
         return self.variables.nvars
 
@@ -681,9 +728,20 @@ class Model:
     def ncons(self):
         """
         Get the total number of constraints.
+
+        This excludes all constraints which are not active.
         """
         return self.constraints.ncons
 
+    @property
+    def shape(self):
+        """
+        Get the shape of the non-filtered constraint matrix.
+
+        This includes all constraints and variables which are not active.
+        """
+        return (self._cCounter, self._xCounter)
+
     @property
     def blocks(self):
         """
@@ -854,7 +912,7 @@ class Model:
 
     def solve(
         self,
-        solver_name="gurobi",
+        solver_name=None,
         io_api=None,
         problem_fn=None,
         solution_fn=None,
@@ -876,7 +934,7 @@ class Model:
         ----------
         solver_name : str, optional
             Name of the solver to use, this must be in `linopy.available_solvers`.
-            The default is 'gurobi'.
+            Default to the first entry in `linopy.available_solvers`.
         io_api : str, optional
             Api to use for communicating with the solver, must be one of
             {'lp', 'mps', 'direct'}. If set to 'lp'/'mps' the problem is written to an
@@ -946,6 +1004,12 @@ class Model:
                     c.dual = solved.constraints[k].dual
             return self.status, self.termination_condition
 
+        if len(available_solvers) == 0:
+            raise RuntimeError("No solver installed.")
+
+        if solver_name is None:
+            solver_name = available_solvers[0]
+
         logger.info(f" Solve linear problem using {solver_name.title()} solver")
         assert solver_name in available_solvers, f"Solver {solver_name} not installed"
 
@@ -967,6 +1031,11 @@ class Model:
         if sanitize_zeros:
             self.constraints.sanitize_zeros()
 
+        if self.is_quadratic and solver_name not in quadratic_solvers:
+            raise ValueError(
+                f"Solver {solver_name} does not support quadratic problems."
+            )
+
         try:
             func = getattr(solvers, f"run_{solver_name}")
             result = func(
@@ -1001,7 +1070,10 @@ class Model:
 
         for name, var in self.variables.items():
             idx = np.ravel(var.labels)
-            vals = sol[idx].values.reshape(var.labels.shape)
+            try:
+                vals = sol[idx].values.reshape(var.labels.shape)
+            except KeyError:
+                vals = sol.reindex(idx).values.reshape(var.labels.shape)
             var.solution = xr.DataArray(vals, var.coords)
 
         if not result.solution.dual.empty:
@@ -1029,6 +1101,9 @@ class Model:
         labels : xr.DataArray
             Labels of the infeasible constraints. Labels with value -1 are not in the set.
         """
+        if "gurobi" not in available_solvers:
+            raise ImportError("Gurobi is required for this method.")
+
         import gurobipy
 
         solver_model = getattr(self, "solver_model")
diff --git a/linopy/solvers.py b/linopy/solvers.py
index 1206095..066dc3e 100644
--- a/linopy/solvers.py
+++ b/linopy/solvers.py
@@ -11,6 +11,7 @@ import re
 import subprocess as sub
 from pathlib import Path
 
+import numpy as np
 import pandas as pd
 
 from linopy.constants import (
@@ -21,15 +22,13 @@ from linopy.constants import (
     TerminationCondition,
 )
 
+quadratic_solvers = ["gurobi", "xpress", "cplex", "highs"]
+
 available_solvers = []
 
 which = "where" if os.name == "nt" else "which"
-if sub.run([which, "glpsol"], stdout=sub.DEVNULL, stderr=sub.STDOUT).returncode == 0:
-    available_solvers.append("glpk")
-
-if sub.run([which, "cbc"], stdout=sub.DEVNULL, stderr=sub.STDOUT).returncode == 0:
-    available_solvers.append("cbc")
 
+# the first available solver will be the default solver
 try:
     import gurobipy
 
@@ -44,6 +43,12 @@ try:
 except (ModuleNotFoundError, ImportError):
     pass
 
+if sub.run([which, "glpsol"], stdout=sub.DEVNULL, stderr=sub.STDOUT).returncode == 0:
+    available_solvers.append("glpk")
+
+
+if sub.run([which, "cbc"], stdout=sub.DEVNULL, stderr=sub.STDOUT).returncode == 0:
+    available_solvers.append("cbc")
 
 try:
     import cplex
@@ -254,13 +259,17 @@ def run_glpk(
     else:
         p.wait()
 
+    if not os.path.exists(solution_fn):
+        status = Status(SolverStatus.warning, TerminationCondition.unknown)
+        return Result(status, Solution())
+
     f = open(solution_fn)
 
     def read_until_break(f):
-        linebreak = False
-        while not linebreak:
+        while True:
             line = f.readline()
-            linebreak = line == "\n"
+            if line == "\n" or line == "":
+                break
             yield line
 
     info = io.StringIO("".join(read_until_break(f))[:-2])
@@ -445,17 +454,22 @@ def run_cplex(
 
     if warmstart_fn:
         m.start.read_basis(warmstart_fn)
-    m.solve()
+
     is_lp = m.problem_type[m.get_problem_type()] == "LP"
 
-    if log_fn is not None:
-        log_f.close()
+    try:
+        m.solve()
+    except cplex.exceptions.errors.CplexSolverError as e:
+        pass
 
     condition = m.solution.get_status_string()
     termination_condition = CONDITION_MAP.get(condition, condition)
     status = Status.from_termination_condition(termination_condition)
     status.legacy_status = condition
 
+    if log_fn is not None:
+        log_f.close()
+
     def get_solver_solution() -> Solution:
         if basis_fn and is_lp:
             try:
diff --git a/linopy/variables.py b/linopy/variables.py
index 1dbd333..4fd4b09 100644
--- a/linopy/variables.py
+++ b/linopy/variables.py
@@ -33,6 +33,7 @@ from linopy.common import (
     save_join,
 )
 from linopy.config import options
+from linopy.constants import TERM_DIM
 
 logger = logging.getLogger(__name__)
 
@@ -207,7 +208,7 @@ class Variable:
         """
         coefficient = as_dataarray(coefficient, coords=self.coords, dims=self.dims)
         ds = Dataset({"coeffs": coefficient, "vars": self.labels}).expand_dims(
-            "_term", -1
+            TERM_DIM, -1
         )
         return expressions.LinearExpression(ds, self.model)
 
@@ -261,13 +262,10 @@ class Variable:
         """
         Multiply variables with a coefficient.
         """
-        if isinstance(other, (expressions.LinearExpression, Variable)):
-            raise TypeError(
-                "unsupported operand type(s) for *: "
-                f"{type(self)} and {type(other)}. "
-                "Non-linear expressions are not yet supported."
-            )
-        return self.to_linexpr(other)
+        if isinstance(other, (expressions.LinearExpression, Variable, ScalarVariable)):
+            return self.to_linexpr() * other
+        else:
+            return self.to_linexpr(other)
 
     def __rmul__(self, other):
         """
@@ -803,7 +801,7 @@ class Variables:
             coords = " (" + ", ".join(ds.coords) + ")" if ds.coords else ""
             r += f" * {name}{coords}\n"
         if not len(list(self)):
-            r += "<empty>"
+            r += "<empty>\n"
         return r
 
     def __len__(self):
@@ -939,13 +937,10 @@ class Variables:
         """
         Get tuple of name and coordinate for variable labels.
         """
-        coords = []
-        return_list = True
-        if not isinstance(values, Iterable):
-            values = [values]
-            return_list = False
 
-        for value in values:
+        def find_single(value):
+            if value == -1:
+                return None, None
             for name, var in self.items():
                 labels = var.labels
                 start, stop = var.range
@@ -960,9 +955,17 @@ class Variables:
                     }
 
                     # Add the name of the DataArray and the coordinates to the result list
-                    coords.append((name, coord))
-
-        return coords if return_list else coords[0]
+                    return name, coord
+
+        ndim = np.array(values).ndim
+        if ndim == 0:
+            return find_single(values)
+        elif ndim == 1:
+            return [find_single(v) for v in values]
+        elif ndim == 2:
+            return [[find_single(v) for v in _] for _ in values.T]
+        else:
+            raise ValueError("Array's with more than two dimensions is not supported")
 
     @deprecated("0.2", details="Use `to_dataframe` or `flat` instead.")
     def iter_ravel(self, key, filter_missings=False):


Test patch:
diff --git a/test/test_constraint.py b/test/test_constraint.py
index 2d73837..15ef8fb 100644
--- a/test/test_constraint.py
+++ b/test/test_constraint.py
@@ -299,6 +299,10 @@ def test_constraint_sel(c):
     assert isinstance(c.isel(first=[1, 2]), linopy.constraints.Constraint)
 
 
+def test_constraint_flat(c):
+    assert isinstance(c.flat, pd.DataFrame)
+
+
 def test_constraint_assignment_with_anonymous_constraints(m, x, y):
     m.add_constraints(x + y == 0, name="c2")
     assert m.constraints["c2"].vars.notnull().all()
@@ -372,7 +376,7 @@ def test_constraint_matrix_masked_variables():
     assert A.shape == (m.ncons, m.nvars)
 
     A = m.constraints.to_matrix(filter_missings=False)
-    assert A.shape == (m._cCounter, m._xCounter)
+    assert A.shape == m.shape
 
 
 def test_constraint_matrix_masked_constraints():
@@ -390,7 +394,7 @@ def test_constraint_matrix_masked_constraints():
     assert A.shape == (m.ncons, m.nvars)
 
     A = m.constraints.to_matrix(filter_missings=False)
-    assert A.shape == (m._cCounter, m._xCounter)
+    assert A.shape == m.shape
 
 
 def test_constraint_matrix_masked_constraints_and_variables():
@@ -408,7 +412,7 @@ def test_constraint_matrix_masked_constraints_and_variables():
     assert A.shape == (m.ncons, m.nvars)
 
     A = m.constraints.to_matrix(filter_missings=False)
-    assert A.shape == (m._cCounter, m._xCounter)
+    assert A.shape == m.shape
 
 
 def test_get_name_by_label():
diff --git a/test/test_constraints.py b/test/test_constraints.py
index 18a3b93..d0f9ea5 100644
--- a/test/test_constraints.py
+++ b/test/test_constraints.py
@@ -123,3 +123,23 @@ def test_masked_constraints():
     m.add_constraints(1 * x + 10 * y, EQUAL, 0, mask=mask)
     assert (m.constraints.labels.con0[0:5, :] != -1).all()
     assert (m.constraints.labels.con0[5:10, :] == -1).all()
+
+
+def test_constraints_flat():
+    m = Model()
+
+    lower = xr.DataArray(np.zeros((10, 10)), coords=[range(10), range(10)])
+    upper = xr.DataArray(np.ones((10, 10)), coords=[range(10), range(10)])
+    x = m.add_variables(lower, upper)
+    y = m.add_variables()
+
+    assert isinstance(m.constraints.flat, pd.DataFrame)
+    assert m.constraints.flat.empty
+    assert m.constraints.to_matrix() is None
+
+    m.add_constraints(1 * x + 10 * y, EQUAL, 0)
+    m.add_constraints(1 * x + 10 * y, LESS_EQUAL, 0)
+    m.add_constraints(1 * x + 10 * y, GREATER_EQUAL, 0)
+
+    assert isinstance(m.constraints.flat, pd.DataFrame)
+    assert not m.constraints.flat.empty
diff --git a/test/test_linear_expression.py b/test/test_linear_expression.py
index 31f626a..c548b78 100644
--- a/test/test_linear_expression.py
+++ b/test/test_linear_expression.py
@@ -13,6 +13,7 @@ import xarray as xr
 from xarray.testing import assert_equal
 
 from linopy import LinearExpression, Model, merge
+from linopy.constants import TERM_DIM
 from linopy.testing import assert_linequal
 
 
@@ -73,7 +74,7 @@ def test_linexpr_with_wrong_data(m):
     with pytest.raises(ValueError):
         LinearExpression(data, m)
 
-    coords = [pd.Index([0], name="a"), pd.Index([1, 2], name="_term")]
+    coords = [pd.Index([0], name="a"), pd.Index([1, 2], name=TERM_DIM)]
     coeffs = xr.DataArray(np.array([[1, 2]]), coords=coords)
     vars = xr.DataArray(np.array([[1, 2]]), coords=coords)
     data = xr.Dataset({"coeffs": coeffs, "vars": vars})
@@ -83,8 +84,8 @@ def test_linexpr_with_wrong_data(m):
 
 def test_linexpr_with_data_without_coords(m):
     lhs = 1 * m["x"]
-    vars = xr.DataArray(lhs.vars.values, dims=["dim_0", "_term"])
-    coeffs = xr.DataArray(lhs.coeffs.values, dims=["dim_0", "_term"])
+    vars = xr.DataArray(lhs.vars.values, dims=["dim_0", TERM_DIM])
+    coeffs = xr.DataArray(lhs.coeffs.values, dims=["dim_0", TERM_DIM])
     data = xr.Dataset({"vars": vars, "coeffs": coeffs})
     expr = LinearExpression(data, m)
     assert_linequal(expr, lhs)
@@ -102,7 +103,7 @@ def test_fill_value():
 def test_linexpr_with_scalars(m):
     expr = m.linexpr((10, "x"), (1, "y"))
     target = xr.DataArray(
-        [[10, 1], [10, 1]], coords={"dim_0": [0, 1]}, dims=["dim_0", "_term"]
+        [[10, 1], [10, 1]], coords={"dim_0": [0, 1]}, dims=["dim_0", TERM_DIM]
     )
     assert_equal(expr.coeffs, target)
 
@@ -187,18 +188,9 @@ def test_linear_expression_multi_indexed(u):
 
 
 def test_linear_expression_with_errors(m, x):
-    with pytest.raises(TypeError):
-        x.to_linexpr(x)
-
-    with pytest.raises(TypeError):
-        x * x
-
     with pytest.raises(TypeError):
         x / x
 
-    with pytest.raises(TypeError):
-        x * (1 * x)
-
     with pytest.raises(TypeError):
         x / (1 * x)
 
@@ -293,7 +285,7 @@ def test_linear_expression_sum(x, y, z, v):
     assert res.nterm == expr.size
     assert res.data.notnull().all().to_array().all()
 
-    assert_linequal(expr.sum(["dim_0", "_term"]), expr.sum("dim_0"))
+    assert_linequal(expr.sum(["dim_0", TERM_DIM]), expr.sum("dim_0"))
 
     # test special case otherride coords
     expr = v.loc[:9] + v.loc[10:]
@@ -315,7 +307,7 @@ def test_linear_expression_sum_with_const(x, y, z, v):
     assert res.data.notnull().all().to_array().all()
     assert (res.const == 60).item()
 
-    assert_linequal(expr.sum(["dim_0", "_term"]), expr.sum("dim_0"))
+    assert_linequal(expr.sum(["dim_0", TERM_DIM]), expr.sum("dim_0"))
 
     # test special case otherride coords
     expr = v.loc[:9] + v.loc[10:]
@@ -366,13 +358,18 @@ def test_linear_expression_multiplication_invalid(x, y, z):
 
     with pytest.raises(TypeError):
         expr = 10 * x + y + z
-        expr * x
+        expr * expr
 
     with pytest.raises(TypeError):
         expr = 10 * x + y + z
         expr / x
 
 
+def test_linear_expression_loc(x, y):
+    expr = x + y
+    assert expr.loc[0].size < expr.loc[:5].size
+
+
 def test_linear_expression_where(v):
     expr = np.arange(20) * v
     expr = expr.where(expr.coeffs >= 10)
@@ -592,9 +589,9 @@ def test_linear_expression_outer_sum(x, y):
 def test_rename(x, y, z):
     expr = 10 * x + y + z
     renamed = expr.rename({"dim_0": "dim_5"})
-    assert set(renamed.dims) == {"dim_1", "dim_5", "_term"}
+    assert set(renamed.dims) == {"dim_1", "dim_5", TERM_DIM}
     assert renamed.nterm == 3
 
     renamed = expr.rename({"dim_0": "dim_1", "dim_1": "dim_2"})
-    assert set(renamed.dims) == {"dim_1", "dim_2", "_term"}
+    assert set(renamed.dims) == {"dim_1", "dim_2", TERM_DIM}
     assert renamed.nterm == 3
diff --git a/test/test_model.py b/test/test_model.py
index 53e3e0d..59755af 100644
--- a/test/test_model.py
+++ b/test/test_model.py
@@ -76,6 +76,10 @@ def test_objective():
     assert m.objectiverange.min() == 2
     assert m.objectiverange.max() == 2
 
+    # test objective with constant which is not supported
+    with pytest.raises(ValueError):
+        m.objective = m.objective + 3
+
 
 def test_remove_variable():
     m = Model()
diff --git a/test/test_optimization.py b/test/test_optimization.py
index a31cd49..bbf4a50 100644
--- a/test/test_optimization.py
+++ b/test/test_optimization.py
@@ -12,7 +12,8 @@ import pytest
 from xarray.testing import assert_equal
 
 from linopy import GREATER_EQUAL, Model
-from linopy.solvers import available_solvers
+from linopy.constants import SolverStatus, Status, TerminationCondition
+from linopy.solvers import available_solvers, quadratic_solvers
 
 params = [(name, "lp") for name in available_solvers]
 # mps io is only supported via highspy
@@ -152,6 +153,34 @@ def milp_model_r():
     return m
 
 
+@pytest.fixture
+def quadratic_model():
+    m = Model()
+
+    lower = pd.Series(0, range(10))
+    x = m.add_variables(lower, name="x")
+    y = m.add_variables(lower, name="y")
+
+    m.add_constraints(x + y, GREATER_EQUAL, 10)
+
+    m.add_objective(x * x)
+    return m
+
+
+@pytest.fixture
+def quadratic_model_cross_terms():
+    m = Model()
+
+    lower = pd.Series(0, range(10))
+    x = m.add_variables(lower, name="x")
+    y = m.add_variables(lower, name="y")
+
+    m.add_constraints(x + y >= 10)
+
+    m.add_objective(-2 * x + y + x * x)
+    return m
+
+
 @pytest.fixture
 def modified_model():
     m = Model()
@@ -203,6 +232,32 @@ def masked_constraint_model():
     return m
 
 
+def test_model_types(
+    model,
+    model_with_duplicated_variables,
+    milp_binary_model,
+    milp_binary_model_r,
+    milp_model,
+    milp_model_r,
+    quadratic_model,
+    quadratic_model_cross_terms,
+    modified_model,
+    masked_variable_model,
+    masked_constraint_model,
+):
+    assert model.type == "LP"
+    assert model_with_duplicated_variables.type == "LP"
+    assert milp_binary_model.type == "MILP"
+    assert milp_binary_model_r.type == "MILP"
+    assert milp_model.type == "MILP"
+    assert milp_model_r.type == "MILP"
+    assert quadratic_model.type == "QP"
+    assert quadratic_model_cross_terms.type == "QP"
+    assert modified_model.type == "MILP"
+    assert masked_variable_model.type == "LP"
+    assert masked_constraint_model.type == "LP"
+
+
 @pytest.mark.parametrize("solver,io_api", params)
 def test_default_setting(model, solver, io_api):
     status, condition = model.solve(solver, io_api=io_api)
@@ -296,7 +351,7 @@ def test_infeasible_model(model, solver, io_api):
     if solver == "gurobi":
         model.compute_set_of_infeasible_constraints()
     else:
-        with pytest.raises(NotImplementedError):
+        with pytest.raises((NotImplementedError, ImportError)):
             model.compute_set_of_infeasible_constraints()
 
 
@@ -345,6 +400,56 @@ def test_milp_model_r(milp_model_r, solver, io_api):
         assert ((milp_model_r.solution.x == 11) | (milp_model_r.solution.y == 0)).all()
 
 
+@pytest.mark.parametrize("solver,io_api", params)
+def test_quadratic_model(quadratic_model, solver, io_api):
+    if solver in quadratic_solvers:
+        status, condition = quadratic_model.solve(solver, io_api=io_api)
+        assert condition == "optimal"
+        assert (quadratic_model.solution.x.round(3) == 0).all()
+        assert (quadratic_model.solution.y.round(3) == 10).all()
+        assert round(quadratic_model.objective_value, 3) == 0
+    else:
+        with pytest.raises(ValueError):
+            quadratic_model.solve(solver, io_api=io_api)
+
+
+@pytest.mark.parametrize("solver,io_api", params)
+def test_quadratic_model_cross_terms(quadratic_model_cross_terms, solver, io_api):
+    if solver in quadratic_solvers:
+        status, condition = quadratic_model_cross_terms.solve(solver, io_api=io_api)
+        assert condition == "optimal"
+        assert (quadratic_model_cross_terms.solution.x.round(3) == 1.5).all()
+        assert (quadratic_model_cross_terms.solution.y.round(3) == 8.5).all()
+        assert round(quadratic_model_cross_terms.objective_value, 3) == 77.5
+    else:
+        with pytest.raises(ValueError):
+            quadratic_model_cross_terms.solve(solver, io_api=io_api)
+
+
+@pytest.mark.parametrize("solver,io_api", params)
+def test_quadratic_model_wo_constraint(quadratic_model, solver, io_api):
+    quadratic_model.constraints.remove("con0")
+    if solver in quadratic_solvers:
+        status, condition = quadratic_model.solve(solver, io_api=io_api)
+        assert condition == "optimal"
+        assert (quadratic_model.solution.x.round(3) == 0).all()
+        assert round(quadratic_model.objective_value, 3) == 0
+    else:
+        with pytest.raises(ValueError):
+            quadratic_model.solve(solver, io_api=io_api)
+
+
+@pytest.mark.parametrize("solver,io_api", params)
+def test_quadratic_model_unbounded(quadratic_model, solver, io_api):
+    quadratic_model.objective = -quadratic_model.objective
+    if solver in quadratic_solvers:
+        status, condition = quadratic_model.solve(solver, io_api=io_api)
+        assert condition in ["unbounded", "unknown", "infeasible_or_unbounded"]
+    else:
+        with pytest.raises(ValueError):
+            quadratic_model.solve(solver, io_api=io_api)
+
+
 @pytest.mark.parametrize("solver,io_api", params)
 def test_modified_model(modified_model, solver, io_api):
     status, condition = modified_model.solve(solver, io_api=io_api)
diff --git a/test/test_quadratic_expression.py b/test/test_quadratic_expression.py
new file mode 100644
index 0000000..56c2d87
--- /dev/null
+++ b/test/test_quadratic_expression.py
@@ -0,0 +1,151 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+
+import numpy as np
+import pandas as pd
+import pytest
+from scipy.sparse import csc_matrix
+
+from linopy import Model
+from linopy.constants import FACTOR_DIM
+from linopy.expressions import QuadraticExpression
+
+
+@pytest.fixture
+def model():
+    return Model()
+
+
+@pytest.fixture
+def x(model):
+    return model.add_variables(pd.Series([0, 0]), 1, name="x")
+
+
+@pytest.fixture
+def y(model):
+    return model.add_variables(4, pd.Series([8, 10]), name="y")
+
+
+@pytest.fixture
+def z(model):
+    return model.add_variables(4, pd.Series([8, 10]), name="z")
+
+
+def test_quadratic_expression_from_variables_multiplication(x, y):
+    expr = x * y
+    assert isinstance(expr, QuadraticExpression)
+    assert expr.data.dims[FACTOR_DIM] == 2
+
+
+def test_quadratic_expression_from_linexpr_multiplication(x, y):
+    expr = (10 * x + y) * y
+    assert isinstance(expr, QuadraticExpression)
+    assert expr.data.dims[FACTOR_DIM] == 2
+    assert expr.nterm == 2
+
+
+def test_quadratic_expression_wrong_assignment(x, y):
+    with pytest.raises(ValueError):
+        QuadraticExpression((x + y).data, x.model)
+
+    with pytest.raises(ValueError):
+        QuadraticExpression((x + y).data.expand_dims(FACTOR_DIM), x.model)
+
+
+def test_quadratic_expression_addition(model, x, y):
+    expr = x * y + x + 5
+    assert isinstance(expr, QuadraticExpression)
+    assert (expr.const == 5).all()
+    assert expr.nterm == 2
+
+
+def test_quadratic_expression_raddition(x, y):
+    expr = x + x * y + 5
+    assert isinstance(expr, QuadraticExpression)
+    assert (expr.const == 5).all()
+    assert expr.nterm == 2
+
+    with pytest.raises(TypeError):
+        5 + x * y + x
+
+
+def test_quadratic_expression_subtraction(x, y):
+    expr = x * y - x - 5
+    assert isinstance(expr, QuadraticExpression)
+    assert (expr.const == -5).all()
+    assert expr.nterm == 2
+
+
+def test_quadratic_expression_rsubtraction(x, y):
+    expr = x - x * y - 5
+    assert isinstance(expr, QuadraticExpression)
+    assert (expr.const == -5).all()
+    assert expr.nterm == 2
+
+
+def test_quadratic_expression_sum(x, y):
+    expr = x * y + x + 5
+
+    summed_expr = expr.sum(dims="dim_0")
+    assert isinstance(summed_expr, QuadraticExpression)
+    assert not summed_expr.coord_dims
+
+    summed_expr_all = expr.sum()
+    assert isinstance(summed_expr_all, QuadraticExpression)
+    assert not summed_expr_all.coord_dims
+
+
+def test_quadratic_expression_wrong_multiplication(x, y):
+    with pytest.raises(TypeError):
+        x * x * y
+
+
+def test_quadratic_expression_loc(x):
+    expr = x * x
+    assert expr.loc[0].size < expr.loc[:5].size
+
+
+def test_quadratic_expression_flat(x, y):
+    expr = x * y + x + 5
+    df = expr.flat
+    assert isinstance(df, pd.DataFrame)
+
+    expr = x * y + x * y
+    assert expr.nterm == 2
+    assert (expr.flat.coeffs == 2).all()
+    assert len(expr.flat) == 2
+
+
+def test_quadratic_expression_to_matrix(model, x, y):
+    expr = x * y + x + 5
+
+    Q = expr.to_matrix()
+    assert isinstance(Q, csc_matrix)
+    assert Q.shape == (model.nvars, model.nvars)
+
+
+def test_matrices_matrix(model, x, y):
+    expr = 10 * x * y
+    model.objective = expr
+
+    Q = model.matrices.Q
+    assert isinstance(Q, csc_matrix)
+    assert Q.shape == (model.nvars, model.nvars)
+
+
+def test_matrices_matrix_mixed_linear_and_quadratic(model, x, y):
+    expr = x * y + x
+    model.objective = expr + x
+
+    Q = model.matrices.Q
+    assert isinstance(Q, csc_matrix)
+    assert Q.shape == (model._xCounter, model._xCounter)
+
+    c = model.matrices.c
+    assert isinstance(c, np.ndarray)
+    assert c.shape == (model.nvars,)
+
+
+def test_quadratic_to_constraint(x, y):
+    with pytest.raises(NotImplementedError):
+        x * y <= 10
diff --git a/test/test_repr.py b/test/test_repr.py
index 2822a56..99530e6 100644
--- a/test/test_repr.py
+++ b/test/test_repr.py
@@ -44,7 +44,9 @@ lc = 1 * c
 ld = 1 * d
 lav = 1 * a + 1 * v
 luc = 1 * v + 10
-
+lq = x * x
+lq2 = x * x + 1 * x
+lq3 = x * x + 1 * x + 1 + 1 * y + 1 * z
 
 # create anonymous constraint for linear expression
 cu_ = lu >= 0
@@ -92,7 +94,7 @@ def test_single_variable_repr():
 
 
 def test_linear_expression_repr():
-    for expr in [lu, lv, lx, ly, lz, la, lb, lc, ld, lav, luc]:
+    for expr in [lu, lv, lx, ly, lz, la, lb, lc, ld, lav, luc, lq, lq2, lq3]:
         repr(expr)