From c22c94b856e85bdd10ea2656cb770f72dd9edde5 Mon Sep 17 00:00:00 2001 From: Bobby Xiong Date: Tue, 17 Mar 2026 13:37:46 +0100 Subject: [PATCH 01/30] Added m.copy() method. --- linopy/model.py | 72 +++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 72 insertions(+) diff --git a/linopy/model.py b/linopy/model.py index 06e814c6..bfa3ddd2 100644 --- a/linopy/model.py +++ b/linopy/model.py @@ -1877,6 +1877,78 @@ def reset_solution(self) -> None: self.variables.reset_solution() self.constraints.reset_dual() + def copy(self, include_solution: bool = False) -> Model: + """ + Return a deep copy of this model. + + Copies variables, constraints, objective, parameters, blocks, and all + scalar attributes (counters, flags). The copy is fully independent: + modifying one does not affect the other. + + Parameters + ---------- + include_solution : bool, optional + Whether to include the current solution and dual values in the copy. + If False (default), the copy is returned in an initialized state: + solution and dual data are excluded, objective value is set to None, + and status is set to 'initialized'. If True, solution, dual values, + solve status, and objective value are also copied. + + Returns + ------- + Model + A deep copy of the model. + """ + SOLVE_STATE_ATTRS = {"status", "termination_condition"} + + m = Model( + chunk=self._chunk, + force_dim_names=self._force_dim_names, + auto_mask=self._auto_mask, + solver_dir=self._solver_dir, + ) + + m._variables = Variables( + { + name: Variable( + var.data.copy() + if include_solution + else var.data[self.variables.dataset_attrs].copy(), + m, + name, + ) + for name, var in self.variables.items() + }, + m, + ) + + m._constraints = Constraints( + { + name: Constraint( + con.data.copy() + if include_solution + else con.data[self.constraints.dataset_attrs].copy(), + m, + name, + ) + for name, con in self.constraints.items() + }, + m, + ) + + obj_expr = LinearExpression(self.objective.expression.data.copy(), m) + m._objective = Objective(obj_expr, m, self.objective.sense) + m._objective._value = self.objective.value if include_solution else None + + m._parameters = self._parameters.copy(deep=True) + m._blocks = self._blocks.copy() if self._blocks is not None else None + + for attr in self.scalar_attrs: + if include_solution or attr not in SOLVE_STATE_ATTRS: + setattr(m, attr, getattr(self, attr)) + + return m + to_netcdf = to_netcdf to_file = to_file From 821850372acb4c436a1b6b5f5f51036efdb80629 Mon Sep 17 00:00:00 2001 From: Bobby Xiong Date: Tue, 17 Mar 2026 13:38:24 +0100 Subject: [PATCH 02/30] Added testing suite for m.copy(). --- test/test_model.py | 73 ++++++++++++++++++++++++++++++++++++++++++++-- 1 file changed, 71 insertions(+), 2 deletions(-) diff --git a/test/test_model.py b/test/test_model.py index c363fe4c..51241f87 100644 --- a/test/test_model.py +++ b/test/test_model.py @@ -12,8 +12,13 @@ import pytest import xarray as xr -from linopy import EQUAL, Model -from linopy.testing import assert_model_equal +from linopy import EQUAL, Model, available_solvers +from linopy.testing import ( + assert_conequal, + assert_equal, + assert_linequal, + assert_model_equal, +) target_shape: tuple[int, int] = (10, 10) @@ -163,3 +168,67 @@ def test_assert_model_equal() -> None: m.add_objective(obj) assert_model_equal(m, m) + + +def _build_model() -> Model: + """Small representative model used across copy tests.""" + m: Model = Model() + + lower: xr.DataArray = xr.DataArray( + np.zeros((10, 10)), coords=[range(10), range(10)] + ) + upper: xr.DataArray = xr.DataArray(np.ones((10, 10)), coords=[range(10), range(10)]) + x = m.add_variables(lower, upper, name="x") + y = m.add_variables(name="y") + + m.add_constraints(1 * x + 10 * y, EQUAL, 0) + m.add_objective((10 * x + 5 * y).sum()) + + return m + + +def test_model_copy_unsolved() -> None: + """Copy of unsolved model is structurally equal and independent.""" + m = _build_model() + c = m.copy(include_solution=False) + + assert_model_equal(m, c) + + # independence: mutating copy does not affect source + c.add_variables(name="z") + assert "z" not in m.variables + + +@pytest.mark.skipif(len(available_solvers) == 0, reason="No solver installed") +def test_model_copy_solved_with_solution() -> None: + """Copy with include_solution=True preserves solve state.""" + m = _build_model() + m.solve() + + c = m.copy(include_solution=True) + assert_model_equal(m, c) + + +@pytest.mark.skipif(len(available_solvers) == 0, reason="No solver installed") +def test_model_copy_solved_without_solution() -> None: + """Copy with include_solution=False (default) drops solve state but preserves problem structure.""" + m = _build_model() + m.solve() + + c = m.copy(include_solution=False) + + # solve state is dropped + assert c.status == "initialized" + assert c.termination_condition == "" + assert c.objective.value is None + + # problem structure is preserved — compare only dataset_attrs to exclude solution/dual + for v in m.variables: + assert_equal( + c.variables[v].data[c.variables.dataset_attrs], + m.variables[v].data[m.variables.dataset_attrs], + ) + for con in m.constraints: + assert_conequal(c.constraints[con], m.constraints[con], strict=False) + assert_linequal(c.objective.expression, m.objective.expression) + assert c.objective.sense == m.objective.sense From 2da2d3521ac655a33d0749dd1554530eeabe2523 Mon Sep 17 00:00:00 2001 From: Bobby Xiong Date: Tue, 17 Mar 2026 13:50:57 +0100 Subject: [PATCH 03/30] Fix solver_dir type annotation. --- linopy/model.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/linopy/model.py b/linopy/model.py index bfa3ddd2..2db7d7df 100644 --- a/linopy/model.py +++ b/linopy/model.py @@ -1905,7 +1905,7 @@ def copy(self, include_solution: bool = False) -> Model: chunk=self._chunk, force_dim_names=self._force_dim_names, auto_mask=self._auto_mask, - solver_dir=self._solver_dir, + solver_dir=str(self._solver_dir), ) m._variables = Variables( From bb6b8a4216f5d8a9a2467263c87360314cb3f173 Mon Sep 17 00:00:00 2001 From: Bobby Xiong Date: Wed, 18 Mar 2026 14:58:00 +0100 Subject: [PATCH 04/30] Bug fix: xarray copyies need to be . --- linopy/model.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/linopy/model.py b/linopy/model.py index 2db7d7df..55b9d5a5 100644 --- a/linopy/model.py +++ b/linopy/model.py @@ -1913,7 +1913,7 @@ def copy(self, include_solution: bool = False) -> Model: name: Variable( var.data.copy() if include_solution - else var.data[self.variables.dataset_attrs].copy(), + else var.data[self.variables.dataset_attrs].copy(deep=True), m, name, ) @@ -1927,7 +1927,7 @@ def copy(self, include_solution: bool = False) -> Model: name: Constraint( con.data.copy() if include_solution - else con.data[self.constraints.dataset_attrs].copy(), + else con.data[self.constraints.dataset_attrs].copy(deep=True), m, name, ) From 724b850c2c0aaf13840f49b40eaaf900aca0f347 Mon Sep 17 00:00:00 2001 From: Bobby Xiong Date: Wed, 18 Mar 2026 15:35:33 +0100 Subject: [PATCH 05/30] Intial commit: Working dualizer for LPs. --- linopy/dual.py | 548 ++++++++++++++++++++++++++++++++++++++++++++++++ linopy/model.py | 5 + 2 files changed, 553 insertions(+) create mode 100644 linopy/dual.py diff --git a/linopy/dual.py b/linopy/dual.py new file mode 100644 index 00000000..3542e1d9 --- /dev/null +++ b/linopy/dual.py @@ -0,0 +1,548 @@ +""" +Linopy dual module. + +This module contains implementations for constructing the dual of a linear optimization problem. +""" + +from __future__ import annotations + +import logging +from typing import TYPE_CHECKING + +import numpy as np +import pandas as pd +import xarray as xr + +from linopy.expressions import LinearExpression + +if TYPE_CHECKING: + from linopy.model import Model + +logger = logging.getLogger(__name__) + + +def _var_lookup(m: Model) -> dict: + """ + Build a flat label -> (var_name, coord_dict) lookup for all variables in m. + + Used to map entries in m.matrices.vlabels back to their variable name + and xarray coordinates for use in dual feasibility constraint construction. + + Skips masked entries (label == -1) and empty variables. + + Parameters + ---------- + m : Model + Primal linopy model. + + Returns + ------- + dict + Mapping from flat integer label to (var_name, coord_dict) tuple. + """ + var_lookup = {} + logger.debug("Building variable label lookup.") + for var_name, var in m.variables.items(): + labels = var.labels + flat_labels = labels.values.flatten() + + if len(flat_labels) == 0: + logger.debug(f"Skipping empty variable '{var_name}'.") + continue + if not (flat_labels != -1).any(): + logger.debug(f"Variable '{var_name}' is fully masked, skipping.") + continue + + logger.debug( + f"Creating label lookup for variable '{var_name}' with shape {labels.shape} and dims {labels.dims}." + ) + + coord_arrays = ( + np.meshgrid( + *[labels.coords[dim].values for dim in labels.dims], indexing="ij" + ) + if len(labels.dims) > 0 + else [] + ) + flat_coords = [arr.flatten() for arr in coord_arrays] + + for k, flat in enumerate(flat_labels): + if flat != -1: + var_lookup[int(flat)] = ( + var_name, + {dim: flat_coords[i][k] for i, dim in enumerate(labels.dims)}, + ) + + return var_lookup + + +def _con_lookup(m: Model) -> dict: + """ + Build a flat label -> (con_name, coord_dict) lookup for all constraints in m. + + Used to map entries in m.matrices.clabels back to their constraint name + and xarray coordinates for use in dual feasibility constraint construction. + + Skips masked entries (label == -1) and empty or fully-masked constraints. + + Parameters + ---------- + m : Model + Primal linopy model. + + Returns + ------- + dict + Mapping from flat integer label to (con_name, coord_dict) tuple. + """ + con_lookup = {} + logger.debug("Building constraint label lookup.") + for con_name, con in m.constraints.items(): + labels = con.labels + flat_labels = labels.values.flatten() + + if len(flat_labels) == 0: + logger.debug(f"Skipping empty constraint '{con_name}'.") + continue + if not (flat_labels != -1).any(): + logger.debug(f"Constraint '{con_name}' is fully masked, skipping.") + continue + + logger.debug( + f"Creating label lookup for constraint '{con_name}' with shape {labels.shape} and dims {labels.dims}." + ) + + coord_arrays = ( + np.meshgrid( + *[labels.coords[dim].values for dim in labels.dims], indexing="ij" + ) + if len(labels.dims) > 0 + else [] + ) + flat_coords = [arr.flatten() for arr in coord_arrays] + + for k, flat in enumerate(flat_labels): + if flat != -1: + con_lookup[int(flat)] = ( + con_name, + {dim: flat_coords[i][k] for i, dim in enumerate(labels.dims)}, + ) + + return con_lookup + + +def bounds_to_constraints(self) -> None: + """ + Add explicit bound constraints for variables with bounds set directly + in the variable rather than via explicit constraints. + + Adds constraints named '{var_name}-bound-lower' and '{var_name}-bound-upper' + to distinguish from PyPSA's automatic '-fix-*' constraints. + + Also resets variable bounds to [-inf, inf] after adding constraints, + to avoid double-counting in the dual. + """ + logger.debug("Converting variable bounds to explicit constraints.") + logger.debug("Relaxing variable bounds to [-inf, inf].") + for var_name, var in self.variables.items(): + mask = var.labels != -1 + lb = var.lower + ub = var.upper + + # lower bound + if f"{var_name}-bound-lower" not in self.constraints: + has_finite_lb = np.isfinite(lb.values[mask.values]).any() + if has_finite_lb: + self.add_constraints( + var >= lb, + name=f"{var_name}-bound-lower", + mask=mask, + ) + logger.debug(f"Added lower bound constraint for '{var_name}'.") + var.lower.values[mask.values] = -np.inf + # Remove bounds to avoid double-counting in the dual. Rely on the new constraints instead. + self.variables[var_name].lower.values[mask.values] = -np.inf + else: + logger.debug( + f"Variable '{var_name}' has no finite lower bound, skipping." + ) + + # upper bound + if f"{var_name}-bound-upper" not in self.constraints: + has_finite_ub = np.isfinite(ub.values[mask.values]).any() + if has_finite_ub: + self.add_constraints( + var <= ub, + name=f"{var_name}-bound-upper", + mask=mask, + ) + logger.debug(f"Added upper bound constraint for '{var_name}'.") + var.upper.values[mask.values] = np.inf + # Remove bounds to avoid double-counting in the dual. Rely on the new constraints instead. + self.variables[var_name].upper.values[mask.values] = np.inf + else: + logger.debug( + f"Variable '{var_name}' has no finite upper bound, skipping." + ) + + +def _add_dual_variables(m: Model, m2: Model) -> dict: + """ + Add dual variables to m2 corresponding to constraints in m. + + For each active constraint in m, adds a dual variable to m2 following + linopy's sign convention: + + - Equality constraints (=) -> free dual variable (lower=-inf, upper=inf) + - <= constraints -> non-positive dual variable (lower=-inf, upper=0) + - >= constraints -> non-negative dual variable (lower=0, upper=inf) + + This convention ensures that m2.variables[con_name].solution has the same + sign as m.constraints[con_name].dual after solving, allowing direct + comparison without sign adjustments. + + The sign encodes the direction of impact on the objective per unit RHS change: + - <= constraint dual (<=0): increasing RHS by 1 unit changes objective by + dual units (negative = cost decreases, i.e. relaxing the constraint). + - >= constraint dual (>=0): increasing RHS by 1 unit changes objective by + dual units (positive = cost increases, i.e. tightening the constraint). + + Skips constraints with no active rows (empty or fully masked). + + Parameters + ---------- + m : Model + Primal linopy model containing the constraints to dualize. + m2 : Model + Dual linopy model to which dual variables are added. + + Returns + ------- + dict + Mapping from constraint name (str) to the corresponding dual + variable (linopy.Variable) in m2. + """ + dual_vars = {} + for name, con in m.constraints.items(): + sign_vals = con.sign.values.flatten() + + if len(sign_vals) == 0: + logger.warning(f"Constraint '{name}' has no sign values, skipping.") + continue + + mask = con.labels != -1 + if not mask.any(): + logger.debug(f"Constraint '{name}' is fully masked, skipping.") + continue + + if sign_vals[0] == "=": + lower, upper = -np.inf, np.inf + var_type = "free" + elif sign_vals[0] == "<=": + lower, upper = -np.inf, 0 + var_type = "non-positive" + else: # >= + lower, upper = 0, np.inf + var_type = "non-negative" + + logger.debug( + f"Adding {var_type} dual variable for constraint '{name}' with shape {con.shape} and dims {con.labels.dims}." + ) + dual_vars[name] = m2.add_variables( + lower=lower, + upper=upper, + coords=list(con.coords.values()), + name=name, + mask=mask, + ) + + return dual_vars + + +def _build_dual_feas_terms( + m: Model, + dual_vars: dict, + var_lookup: dict, + con_lookup: dict, +) -> dict: + """ + Build dual feasibility terms for each primal variable in m. + + For each active primal variable x_j, collects the constraint matrix + entries A_ji and their corresponding constraint names and coordinates, + forming the terms of the stationarity condition: + sum_i (A_ji * lambda_i) = c_j + + Raw constraint matrix coefficients are used directly without sign + factors, as the sign convention is encoded in the dual variable bounds: + - <= constraints: lambda_i <= 0 + - >= constraints: lambda_i >= 0 + - = constraints: lambda_i free + + Parameters + ---------- + m : Model + Primal linopy model. + dual_vars : dict + Mapping from constraint name to dual variable in m2, + as returned by _add_dual_variables(). Used to skip constraints + that were not dualized (e.g. empty or fully masked). + var_lookup : dict + Mapping from flat variable label to (var_name, coord_dict), + as returned by _var_lookup(). + con_lookup : dict + Mapping from flat constraint label to (con_name, coord_dict), + as returned by _con_lookup(). + + Returns + ------- + dict + Nested dict: {var_name: {flat_label: (var_coords, terms, obj_coeff)}} + where terms is a list of (con_name, con_coords, coeff) tuples. + """ + A_csc = m.matrices.A.tocsc() + c = m.matrices.c + indptr = A_csc.indptr + indices = A_csc.indices + data = A_csc.data + vlabels = m.matrices.vlabels + clabels = m.matrices.clabels + + dual_feas_terms = {var_name: {} for var_name in m.variables} + + logger.debug("Building dual feasibility terms for each primal variable.") + + for i in range(A_csc.shape[1]): + flat_var = vlabels[i] + if flat_var == -1: + continue + if flat_var not in var_lookup: + continue + var_name, var_coords = var_lookup[flat_var] + terms = [] + for k in range(indptr[i], indptr[i + 1]): + j = indices[k] + flat_con = clabels[j] + if flat_con == -1: + continue + if flat_con not in con_lookup: + continue + con_name, con_coords = con_lookup[flat_con] + if con_name not in dual_vars: + continue + coeff = data[k] + terms.append((con_name, con_coords, coeff)) + dual_feas_terms[var_name][flat_var] = (var_coords, terms, c[i]) + + return dual_feas_terms + + +def _add_dual_feasibility_constraints( + m: Model, + m2: Model, + dual_vars: dict, + var_lookup: dict, + con_lookup: dict, +) -> None: + """ + Add dual feasibility constraints to m2. + + For each primal variable x_j in m, adds the stationarity constraint: + sum_i (A_ji * lambda_i) = c_j + where: + - A is the primal constraint matrix + - lambda_i are the dual variables in m2 + - c_j is the objective coefficient of x_j + + Raw constraint matrix coefficients are used directly without sign factors, + because the sign convention is encoded in the dual variable bounds: + - <= constraints: lambda_i <= 0 + - >= constraints: lambda_i >= 0 + - = constraints: lambda_i free + + Skips masked variable entries (label == -1) and variables not present + in var_lookup (e.g. from empty constraints). + + Parameters + ---------- + m : Model + Primal linopy model. + m2 : Model + Dual linopy model. + dual_vars : dict + Mapping from constraint name to dual variable in m2, + as returned by _add_dual_variables(). + var_lookup : dict + Mapping from flat variable label to (var_name, coord_dict), + as returned by _var_lookup(). + con_lookup : dict + Mapping from flat constraint label to (con_name, coord_dict), + as returned by _con_lookup(). + """ + + dual_feas_terms = _build_dual_feas_terms(m, dual_vars, var_lookup, con_lookup) + + c = m.matrices.c + vlabels = m.matrices.vlabels + + # build objective coefficient lookup by flat variable label + c_by_label = {vlabels[i]: c[i] for i in range(len(vlabels))} + + # add dual feasibility constraints to m2 + logger.debug("Adding dual feasibility constraints to model.") + for var_name, var in m.variables.items(): + coords = [ + pd.Index(var.labels.coords[dim].values, name=dim) for dim in var.labels.dims + ] + mask = var.labels != -1 + + c_vals = xr.DataArray( + np.vectorize(lambda flat: c_by_label.get(flat, 0.0))(var.labels.values), + coords=var.labels.coords, + ) + + def rule(m, *coord_vals, vname=var_name, vdims=var.labels.dims): + coord_dict = dict(zip(vdims, coord_vals)) + flat = var.labels.sel(**coord_dict).item() + if flat == -1: + return None + if flat not in dual_feas_terms[vname]: + return None + _, terms, _ = dual_feas_terms[vname][flat] + if not terms: + return None + return sum( + coeff * dual_vars[con_name].at[tuple(con_coords.values())] + for con_name, con_coords, coeff in terms + ) + + lhs = LinearExpression.from_rule(m2, rule, coords) + m2.add_constraints(lhs == c_vals, name=var_name, mask=mask) + + +def _add_dual_objective( + m: Model, + m2: Model, + dual_vars: dict, + add_objective_constant: float = 0.0, +) -> None: + """ + Construct and add the dual objective to m2. + + The dual objective is sum(rhs * dual) over all constraints, added uniformly + with a + sign. The sign convention is encoded in the dual variable bounds: + - <= constraints: dual <= 0, so rhs * dual contributes negatively + - >= constraints: dual >= 0, so rhs * dual contributes positively + - = constraints: dual free + + This matches linopy's and Gurobi's native dual sign convention, allowing + direct comparison between m2.variables[con_name].solution and + m.constraints[con_name].dual without sign adjustments. + + The dual objective sense is flipped relative to the primal: + - min primal -> max dual + - max primal -> min dual + + Parameters + ---------- + m : Model + Primal linopy model. + m2 : Model + Dual linopy model. + dual_vars : dict + Mapping from constraint name to dual variable in m2, + as returned by _add_dual_variables(). + add_objective_constant : float, optional + Constant term to add to the dual objective. Use this to pass through + a primal objective constant excluded via include_objective_constant=False + during model creation. Default is 0.0. + """ + dual_obj = 0 + sense = "max" if m.objective.sense == "min" else "min" + + for name, con in m.constraints.items(): + if name not in dual_vars: + continue + + mask = con.labels != -1 + rhs_masked = con.rhs.where(mask, 0) + dual_obj += (rhs_masked * dual_vars[name]).sum() + + if add_objective_constant != 0.0: + dual_obj += add_objective_constant + logger.debug(f"Added constant {add_objective_constant} to dual objective.") + + logger.debug(f"Constructed dual objective with {len(dual_obj.coeffs)} terms.") + logger.debug("Adding dual objective to model.") + m2.add_objective(dual_obj, sense=sense, overwrite=True) + + +def dualize( + self, + add_objective_constant: float = 0.0, +) -> Model: + """ + Construct the dual of a linopy LP model. + + Transforms the primal model into its dual equivalent m2 by: + 1. Converting variable bounds to explicit constraints + 2. Adding dual variables to m2 (one per active constraint) + 3. Adding dual feasibility constraints to m2 (one per primal variable) + 4. Adding the dual objective to m2 + + The dual is constructed following standard LP duality theory: + + Primal (min): Dual (max): + min c^T x max b_eq^T λ + b_leq^T μ + b_geq^T ν + s.t. A_eq x = b_eq : λ free s.t. A_eq^T λ + A_leq^T μ + A_geq^T ν = c + A_leq x <= b_leq : μ <= 0 λ free, μ <= 0, ν >= 0 + A_geq x >= b_geq : ν >= 0 + + Variable bounds are converted to explicit constraints before dualization + via bounds_to_constraints(), so that they appear in the constraint matrix + A and are correctly reflected in the dual. + + The dual variables in m2 are named identically to their corresponding + primal constraints and are accessible via m2.variables[con_name]. + + Strong duality guarantees that at optimality: + primal objective = dual objective + + Note: The standalone dual m2 may be unbounded if the primal is degenerate. + + Parameters + ---------- + add_objective_constant : float, optional + Constant term to add to the dual objective. Use this to pass through + a primal objective constant. Default is 0.0. + + Returns + ------- + Model + The dual linopy model. Dual variables are named after their + corresponding primal constraints. + + Examples + -------- + >>> m2 = m.dualize() + >>> m2.solve(solver_name="gurobi", Method=2, Crossover=0) + >>> gap = abs(m.objective.value - m2.objective.value) + """ + from linopy.model import Model + + m = self.copy() + m2 = Model() + + if not m.variables or not m.constraints: + logger.warning( + "Primal model has no variables or constraints. Returning empty dual model." + ) + return m2 + + m.bounds_to_constraints() + var_lup = _var_lookup(m) + con_lup = _con_lookup(m) + dual_vars = _add_dual_variables(m, m2) + _add_dual_feasibility_constraints(m, m2, dual_vars, var_lup, con_lup) + _add_dual_objective(m, m2, dual_vars, add_objective_constant=add_objective_constant) + return m2 diff --git a/linopy/model.py b/linopy/model.py index 55b9d5a5..ed4b84b6 100644 --- a/linopy/model.py +++ b/linopy/model.py @@ -47,6 +47,7 @@ TerminationCondition, ) from linopy.constraints import AnonymousScalarConstraint, Constraint, Constraints +from linopy.dual import bounds_to_constraints, dualize from linopy.expressions import ( LinearExpression, QuadraticExpression, @@ -1962,3 +1963,7 @@ def copy(self, include_solution: bool = False) -> Model: to_cupdlpx = to_cupdlpx to_block_files = to_block_files + + bounds_to_constraints = bounds_to_constraints + + dualize = dualize From c49a2b8a38f6d2a2e54be005a7aa98db95540022 Mon Sep 17 00:00:00 2001 From: Bobby Xiong Date: Wed, 18 Mar 2026 16:24:50 +0100 Subject: [PATCH 06/30] Fixed typing issues. Initialise dual objective with LinearExpression(None) instead of 0. --- linopy/dual.py | 19 ++++++++++++++----- 1 file changed, 14 insertions(+), 5 deletions(-) diff --git a/linopy/dual.py b/linopy/dual.py index 3542e1d9..7024d825 100644 --- a/linopy/dual.py +++ b/linopy/dual.py @@ -7,7 +7,7 @@ from __future__ import annotations import logging -from typing import TYPE_CHECKING +from typing import TYPE_CHECKING, Any import numpy as np import pandas as pd @@ -300,6 +300,9 @@ def _build_dual_feas_terms( Nested dict: {var_name: {flat_label: (var_coords, terms, obj_coeff)}} where terms is a list of (con_name, con_coords, coeff) tuples. """ + A = m.matrices.A + if A is None: + raise ValueError("Constraint matrix is None, model has no constraints.") A_csc = m.matrices.A.tocsc() c = m.matrices.c indptr = A_csc.indptr @@ -308,7 +311,9 @@ def _build_dual_feas_terms( vlabels = m.matrices.vlabels clabels = m.matrices.clabels - dual_feas_terms = {var_name: {} for var_name in m.variables} + dual_feas_terms: dict[str, dict[int, tuple]] = { + var_name: {} for var_name in m.variables + } logger.debug("Building dual feasibility terms for each primal variable.") @@ -379,7 +384,6 @@ def _add_dual_feasibility_constraints( Mapping from flat constraint label to (con_name, coord_dict), as returned by _con_lookup(). """ - dual_feas_terms = _build_dual_feas_terms(m, dual_vars, var_lookup, con_lookup) c = m.matrices.c @@ -401,7 +405,12 @@ def _add_dual_feasibility_constraints( coords=var.labels.coords, ) - def rule(m, *coord_vals, vname=var_name, vdims=var.labels.dims): + def rule( + m: Model, + *coord_vals: Any, + vname: str = var_name, + vdims: tuple = var.labels.dims, + ) -> LinearExpression | None: coord_dict = dict(zip(vdims, coord_vals)) flat = var.labels.sel(**coord_dict).item() if flat == -1: @@ -457,7 +466,7 @@ def _add_dual_objective( a primal objective constant excluded via include_objective_constant=False during model creation. Default is 0.0. """ - dual_obj = 0 + dual_obj = LinearExpression(None, m2) sense = "max" if m.objective.sense == "min" else "min" for name, con in m.constraints.items(): From b6c26febc9a23e07fe40874d4607191c9badc1bd Mon Sep 17 00:00:00 2001 From: Bobby Xiong Date: Wed, 18 Mar 2026 16:37:16 +0100 Subject: [PATCH 07/30] Fixed more typing issues. --- linopy/dual.py | 46 ++++++++++++++++++++++++++++------------------ 1 file changed, 28 insertions(+), 18 deletions(-) diff --git a/linopy/dual.py b/linopy/dual.py index 7024d825..47c12233 100644 --- a/linopy/dual.py +++ b/linopy/dual.py @@ -131,7 +131,7 @@ def _con_lookup(m: Model) -> dict: return con_lookup -def bounds_to_constraints(self) -> None: +def bounds_to_constraints(m) -> None: """ Add explicit bound constraints for variables with bounds set directly in the variable rather than via explicit constraints. @@ -141,19 +141,24 @@ def bounds_to_constraints(self) -> None: Also resets variable bounds to [-inf, inf] after adding constraints, to avoid double-counting in the dual. + + Parameters + ---------- + m : Model + Linopy model to convert variable bounds to constraints. Mutates the model in-place. """ logger.debug("Converting variable bounds to explicit constraints.") logger.debug("Relaxing variable bounds to [-inf, inf].") - for var_name, var in self.variables.items(): + for var_name, var in m.variables.items(): mask = var.labels != -1 lb = var.lower ub = var.upper # lower bound - if f"{var_name}-bound-lower" not in self.constraints: + if f"{var_name}-bound-lower" not in m.constraints: has_finite_lb = np.isfinite(lb.values[mask.values]).any() if has_finite_lb: - self.add_constraints( + m.add_constraints( var >= lb, name=f"{var_name}-bound-lower", mask=mask, @@ -161,17 +166,17 @@ def bounds_to_constraints(self) -> None: logger.debug(f"Added lower bound constraint for '{var_name}'.") var.lower.values[mask.values] = -np.inf # Remove bounds to avoid double-counting in the dual. Rely on the new constraints instead. - self.variables[var_name].lower.values[mask.values] = -np.inf + m.variables[var_name].lower.values[mask.values] = -np.inf else: logger.debug( f"Variable '{var_name}' has no finite lower bound, skipping." ) # upper bound - if f"{var_name}-bound-upper" not in self.constraints: + if f"{var_name}-bound-upper" not in m.constraints: has_finite_ub = np.isfinite(ub.values[mask.values]).any() if has_finite_ub: - self.add_constraints( + m.add_constraints( var <= ub, name=f"{var_name}-bound-upper", mask=mask, @@ -179,7 +184,7 @@ def bounds_to_constraints(self) -> None: logger.debug(f"Added upper bound constraint for '{var_name}'.") var.upper.values[mask.values] = np.inf # Remove bounds to avoid double-counting in the dual. Rely on the new constraints instead. - self.variables[var_name].upper.values[mask.values] = np.inf + m.variables[var_name].upper.values[mask.values] = np.inf else: logger.debug( f"Variable '{var_name}' has no finite upper bound, skipping." @@ -303,7 +308,7 @@ def _build_dual_feas_terms( A = m.matrices.A if A is None: raise ValueError("Constraint matrix is None, model has no constraints.") - A_csc = m.matrices.A.tocsc() + A_csc = A.tocsc() c = m.matrices.c indptr = A_csc.indptr indices = A_csc.indices @@ -466,7 +471,7 @@ def _add_dual_objective( a primal objective constant excluded via include_objective_constant=False during model creation. Default is 0.0. """ - dual_obj = LinearExpression(None, m2) + dual_obj: LinearExpression = LinearExpression(None, m2) sense = "max" if m.objective.sense == "min" else "min" for name, con in m.constraints.items(): @@ -487,7 +492,7 @@ def _add_dual_objective( def dualize( - self, + m, add_objective_constant: float = 0.0, ) -> Model: """ @@ -521,6 +526,9 @@ def dualize( Parameters ---------- + m : Model + Primal linopy model to dualize. Must have a linear objective and linear constraints. + add_objective_constant : float, optional Constant term to add to the dual objective. Use this to pass through a primal objective constant. Default is 0.0. @@ -539,7 +547,7 @@ def dualize( """ from linopy.model import Model - m = self.copy() + m1 = m.copy() m2 = Model() if not m.variables or not m.constraints: @@ -548,10 +556,12 @@ def dualize( ) return m2 - m.bounds_to_constraints() - var_lup = _var_lookup(m) - con_lup = _con_lookup(m) - dual_vars = _add_dual_variables(m, m2) - _add_dual_feasibility_constraints(m, m2, dual_vars, var_lup, con_lup) - _add_dual_objective(m, m2, dual_vars, add_objective_constant=add_objective_constant) + m1.bounds_to_constraints() + var_lup = _var_lookup(m1) + con_lup = _con_lookup(m1) + dual_vars = _add_dual_variables(m1, m2) + _add_dual_feasibility_constraints(m1, m2, dual_vars, var_lup, con_lup) + _add_dual_objective( + m1, m2, dual_vars, add_objective_constant=add_objective_constant + ) return m2 From 8b66ba0d16e90c1d9f586100d5961ca94f1226d4 Mon Sep 17 00:00:00 2001 From: Bobby Xiong Date: Wed, 18 Mar 2026 16:39:01 +0100 Subject: [PATCH 08/30] Added code-block to dualize() examples. --- linopy/dual.py | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/linopy/dual.py b/linopy/dual.py index 47c12233..59532e5a 100644 --- a/linopy/dual.py +++ b/linopy/dual.py @@ -541,9 +541,11 @@ def dualize( Examples -------- - >>> m2 = m.dualize() - >>> m2.solve(solver_name="gurobi", Method=2, Crossover=0) - >>> gap = abs(m.objective.value - m2.objective.value) + .. code-block:: python + + m2 = m.dualize() + m2.solve(solver_name="gurobi", Method=2, Crossover=1) + gap = abs(m.objective.value - m2.objective.value) """ from linopy.model import Model From 883d858c115178503b8f033196c9dba4b70c0d05 Mon Sep 17 00:00:00 2001 From: Bobby Xiong Date: Wed, 18 Mar 2026 16:45:34 +0100 Subject: [PATCH 09/30] Fix more typing issues. --- linopy/dual.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/linopy/dual.py b/linopy/dual.py index 59532e5a..e5483b06 100644 --- a/linopy/dual.py +++ b/linopy/dual.py @@ -131,7 +131,7 @@ def _con_lookup(m: Model) -> dict: return con_lookup -def bounds_to_constraints(m) -> None: +def bounds_to_constraints(m: Model) -> None: """ Add explicit bound constraints for variables with bounds set directly in the variable rather than via explicit constraints. @@ -492,7 +492,7 @@ def _add_dual_objective( def dualize( - m, + m: Model, add_objective_constant: float = 0.0, ) -> Model: """ @@ -558,7 +558,7 @@ def dualize( ) return m2 - m1.bounds_to_constraints() + bounds_to_constraints(m1) var_lup = _var_lookup(m1) con_lup = _con_lookup(m1) dual_vars = _add_dual_variables(m1, m2) From 1b6d09f740ddbd9de6b6c03c4c7e1ccb97d03cb8 Mon Sep 17 00:00:00 2001 From: Bobby Xiong Date: Wed, 18 Mar 2026 19:04:12 +0100 Subject: [PATCH 10/30] Fixed variable bounds for dualisation of max problems (primal) --> reverse sign. --- linopy/dual.py | 24 ++++++++++++++++++------ 1 file changed, 18 insertions(+), 6 deletions(-) diff --git a/linopy/dual.py b/linopy/dual.py index e5483b06..a0881433 100644 --- a/linopy/dual.py +++ b/linopy/dual.py @@ -227,6 +227,8 @@ def _add_dual_variables(m: Model, m2: Model) -> dict: Mapping from constraint name (str) to the corresponding dual variable (linopy.Variable) in m2. """ + primal_is_min = m.objective.sense == "min" + dual_vars = {} for name, con in m.constraints.items(): sign_vals = con.sign.values.flatten() @@ -244,19 +246,29 @@ def _add_dual_variables(m: Model, m2: Model) -> dict: lower, upper = -np.inf, np.inf var_type = "free" elif sign_vals[0] == "<=": - lower, upper = -np.inf, 0 - var_type = "non-positive" - else: # >= - lower, upper = 0, np.inf - var_type = "non-negative" + lower, upper = (-np.inf, 0) if primal_is_min else (0, np.inf) + var_type = "non-positive" if primal_is_min else "non-negative" + elif sign_vals[0] == ">=": + lower, upper = (0, np.inf) if primal_is_min else (-np.inf, 0) + var_type = "non-negative" if primal_is_min else "non-positive" + else: + logger.warning( + f"Constraint '{name}' has unrecognized sign '{sign_vals[0]}', skipping." + ) + continue logger.debug( f"Adding {var_type} dual variable for constraint '{name}' with shape {con.shape} and dims {con.labels.dims}." ) + coords = ( + [con.labels.coords[dim] for dim in con.labels.dims] + if con.labels.dims + else None + ) dual_vars[name] = m2.add_variables( lower=lower, upper=upper, - coords=list(con.coords.values()), + coords=coords, name=name, mask=mask, ) From f45c49cfe8b49edb5bb3d2e63cd65004f60c7ef1 Mon Sep 17 00:00:00 2001 From: Bobby Xiong Date: Wed, 18 Mar 2026 19:21:08 +0100 Subject: [PATCH 11/30] Updated docstrings for max -> min. --- linopy/dual.py | 44 ++++++++++++++++++++++---------------------- 1 file changed, 22 insertions(+), 22 deletions(-) diff --git a/linopy/dual.py b/linopy/dual.py index a0881433..0a787e68 100644 --- a/linopy/dual.py +++ b/linopy/dual.py @@ -193,24 +193,21 @@ def bounds_to_constraints(m: Model) -> None: def _add_dual_variables(m: Model, m2: Model) -> dict: """ - Add dual variables to m2 corresponding to constraints in m. + Add dual variables to m2 corresponding to constraints in m.. For each active constraint in m, adds a dual variable to m2 following - linopy's sign convention: + standard LP duality sign conventions. The sign of the dual variable bounds + depends on both the constraint type and the primal objective sense: + For a minimization primal: - Equality constraints (=) -> free dual variable (lower=-inf, upper=inf) - <= constraints -> non-positive dual variable (lower=-inf, upper=0) - >= constraints -> non-negative dual variable (lower=0, upper=inf) - This convention ensures that m2.variables[con_name].solution has the same - sign as m.constraints[con_name].dual after solving, allowing direct - comparison without sign adjustments. - - The sign encodes the direction of impact on the objective per unit RHS change: - - <= constraint dual (<=0): increasing RHS by 1 unit changes objective by - dual units (negative = cost decreases, i.e. relaxing the constraint). - - >= constraint dual (>=0): increasing RHS by 1 unit changes objective by - dual units (positive = cost increases, i.e. tightening the constraint). + For a maximization primal: + - Equality constraints (=) -> free dual variable (lower=-inf, upper=inf) + - <= constraints -> non-negative dual variable (lower=0, upper=inf) + - >= constraints -> non-positive dual variable (lower=-inf, upper=0) Skips constraints with no active rows (empty or fully masked). @@ -510,19 +507,21 @@ def dualize( """ Construct the dual of a linopy LP model. - Transforms the primal model into its dual equivalent m2 by: - 1. Converting variable bounds to explicit constraints - 2. Adding dual variables to m2 (one per active constraint) - 3. Adding dual feasibility constraints to m2 (one per primal variable) - 4. Adding the dual objective to m2 + Transforms the primal model into its dual equivalent m2 following + standard LP duality theory. The dual sense is flipped relative to the + primal (min -> max, max -> min), and dual variable bounds depend on + both constraint type and primal objective sense. + + For a minimization primal: - The dual is constructed following standard LP duality theory: + Primal (min): Dual (max): + min c^T x max b_eq^T λ + b_leq^T μ + b_geq^T ν + s.t. A_eq x = b_eq : λ free s.t. A_eq^T λ + A_leq^T μ + A_geq^T ν = c + A_leq x <= b_leq : μ <= 0 λ free, μ <= 0, ν >= 0 + A_geq x >= b_geq : ν >= 0 - Primal (min): Dual (max): - min c^T x max b_eq^T λ + b_leq^T μ + b_geq^T ν - s.t. A_eq x = b_eq : λ free s.t. A_eq^T λ + A_leq^T μ + A_geq^T ν = c - A_leq x <= b_leq : μ <= 0 λ free, μ <= 0, ν >= 0 - A_geq x >= b_geq : ν >= 0 + For a maximization primal the dual variable bounds are flipped: + μ >= 0 for <= constraints, ν <= 0 for >= constraints. Variable bounds are converted to explicit constraints before dualization via bounds_to_constraints(), so that they appear in the constraint matrix @@ -535,6 +534,7 @@ def dualize( primal objective = dual objective Note: The standalone dual m2 may be unbounded if the primal is degenerate. + Only linear programs (LP) are supported. Parameters ---------- From fe8d952aea94fb80e73d227719f10b61cbf7b36d Mon Sep 17 00:00:00 2001 From: Bobby Xiong Date: Wed, 18 Mar 2026 19:39:58 +0100 Subject: [PATCH 12/30] Moved copy to io.py, added deep-copy to all xarray operations. --- linopy/io.py | 85 +++++++++++++++++++++++++++++++++++++++++++++++++ linopy/model.py | 73 ++---------------------------------------- 2 files changed, 87 insertions(+), 71 deletions(-) diff --git a/linopy/io.py b/linopy/io.py index 2213cbb5..c380c931 100644 --- a/linopy/io.py +++ b/linopy/io.py @@ -1239,3 +1239,88 @@ def get_prefix(ds: xr.Dataset, prefix: str) -> xr.Dataset: setattr(m, k, ds.attrs.get(k)) return m + + +def copy(src: Model, include_solution: bool = False) -> Model: + """ + Return a deep copy of this model. + + Copies variables, constraints, objective, parameters, blocks, and all + scalar attributes (counters, flags). The copy is fully independent: + modifying one does not affect the other. + + Parameters + ---------- + src : Model + The model to copy. + include_solution : bool, optional + Whether to include the current solution and dual values in the copy. + If False (default), the copy is returned in an initialized state: + solution and dual data are excluded, objective value is set to None, + and status is set to 'initialized'. If True, solution, dual values, + solve status, and objective value are also copied. + + Returns + ------- + Model + A deep copy of the model. + """ + from linopy.model import ( + Constraint, + Constraints, + LinearExpression, + Model, + Objective, + Variable, + Variables, + ) + + SOLVE_STATE_ATTRS = {"status", "termination_condition"} + + m = Model( + chunk=src._chunk, + force_dim_names=src._force_dim_names, + auto_mask=src._auto_mask, + solver_dir=str(src._solver_dir), + ) + + m._variables = Variables( + { + name: Variable( + var.data.copy(deep=True) + if include_solution + else var.data[src.variables.dataset_attrs].copy(deep=True), + m, + name, + ) + for name, var in src.variables.items() + }, + m, + ) + + m._constraints = Constraints( + { + name: Constraint( + con.data.copy(deep=True) + if include_solution + else con.data[src.constraints.dataset_attrs].copy(deep=True), + m, + name, + ) + for name, con in src.constraints.items() + }, + m, + ) + + obj_expr = LinearExpression(src.objective.expression.data.copy(deep=True), m) + m._objective = Objective(obj_expr, m, src.objective.sense) + m._objective._value = src.objective.value if include_solution else None + + m._parameters = src._parameters.copy(deep=True) + m._blocks = src._blocks.copy(deep=True) if src._blocks is not None else None + + for attr in src.scalar_attrs: + if include_solution or attr not in SOLVE_STATE_ATTRS: + setattr(m, attr, getattr(src, attr)) + + return m diff --git a/linopy/model.py b/linopy/model.py index 55b9d5a5..48ba8b02 100644 --- a/linopy/model.py +++ b/linopy/model.py @@ -53,6 +53,7 @@ ScalarLinearExpression, ) from linopy.io import ( + copy, to_block_files, to_cupdlpx, to_file, @@ -1877,77 +1878,7 @@ def reset_solution(self) -> None: self.variables.reset_solution() self.constraints.reset_dual() - def copy(self, include_solution: bool = False) -> Model: - """ - Return a deep copy of this model. - - Copies variables, constraints, objective, parameters, blocks, and all - scalar attributes (counters, flags). The copy is fully independent: - modifying one does not affect the other. - - Parameters - ---------- - include_solution : bool, optional - Whether to include the current solution and dual values in the copy. - If False (default), the copy is returned in an initialized state: - solution and dual data are excluded, objective value is set to None, - and status is set to 'initialized'. If True, solution, dual values, - solve status, and objective value are also copied. - - Returns - ------- - Model - A deep copy of the model. - """ - SOLVE_STATE_ATTRS = {"status", "termination_condition"} - - m = Model( - chunk=self._chunk, - force_dim_names=self._force_dim_names, - auto_mask=self._auto_mask, - solver_dir=str(self._solver_dir), - ) - - m._variables = Variables( - { - name: Variable( - var.data.copy() - if include_solution - else var.data[self.variables.dataset_attrs].copy(deep=True), - m, - name, - ) - for name, var in self.variables.items() - }, - m, - ) - - m._constraints = Constraints( - { - name: Constraint( - con.data.copy() - if include_solution - else con.data[self.constraints.dataset_attrs].copy(deep=True), - m, - name, - ) - for name, con in self.constraints.items() - }, - m, - ) - - obj_expr = LinearExpression(self.objective.expression.data.copy(), m) - m._objective = Objective(obj_expr, m, self.objective.sense) - m._objective._value = self.objective.value if include_solution else None - - m._parameters = self._parameters.copy(deep=True) - m._blocks = self._blocks.copy() if self._blocks is not None else None - - for attr in self.scalar_attrs: - if include_solution or attr not in SOLVE_STATE_ATTRS: - setattr(m, attr, getattr(self, attr)) - - return m + copy = copy to_netcdf = to_netcdf From f4016a8811e196092f7983e7543d776c74f8c928 Mon Sep 17 00:00:00 2001 From: Bobby Xiong Date: Fri, 20 Mar 2026 15:09:44 +0100 Subject: [PATCH 13/30] Improved copy method: Strengtheninc copy protocol compatibility, check for deep copy independence. --- linopy/io.py | 99 ++++++++++++++++++--------- linopy/model.py | 6 ++ test/test_model.py | 163 ++++++++++++++++++++++++++++++++++++--------- 3 files changed, 204 insertions(+), 64 deletions(-) diff --git a/linopy/io.py b/linopy/io.py index c380c931..fd3c536a 100644 --- a/linopy/io.py +++ b/linopy/io.py @@ -1241,29 +1241,38 @@ def get_prefix(ds: xr.Dataset, prefix: str) -> xr.Dataset: return m -def copy(src: Model, include_solution: bool = False) -> Model: +def copy(m: Model, include_solution: bool = False, deep: bool = True) -> Model: """ - Return a deep copy of this model. + Return a copy of this model. - Copies variables, constraints, objective, parameters, blocks, and all - scalar attributes (counters, flags). The copy is fully independent: - modifying one does not affect the other. + With ``deep=True`` (default), variables, constraints, objective, + parameters, blocks, and scalar attributes are copied to a fully + independent model. With ``deep=False``, returns a shallow copy. + + Solver runtime metadata (for example, ``solver_name`` and + ``solver_model``) is intentionally not copied. Solver backend state + is recreated on ``solve()``. Parameters ---------- - src : Model + m : Model The model to copy. include_solution : bool, optional Whether to include the current solution and dual values in the copy. If False (default), the copy is returned in an initialized state: solution and dual data are excluded, objective value is set to None, and status is set to 'initialized'. If True, solution, dual values, - solve status, and objective value are also copied. + solve status, and objective value are also copied. If the model is + unsolved, this has no additional effect. + deep : bool, optional + Whether to return a deep copy (default) or shallow copy. If False, + the returned model uses independent wrapper objects that share + underlying data buffers with the source model. Returns ------- Model - A deep copy of the model. + A deep or shallow copy of the model. """ from linopy.model import ( Constraint, @@ -1277,50 +1286,74 @@ def copy(src: Model, include_solution: bool = False) -> Model: SOLVE_STATE_ATTRS = {"status", "termination_condition"} - m = Model( - chunk=src._chunk, - force_dim_names=src._force_dim_names, - auto_mask=src._auto_mask, - solver_dir=str(src._solver_dir), + new_model = Model( + chunk=m._chunk, + force_dim_names=m._force_dim_names, + auto_mask=m._auto_mask, + solver_dir=str(m._solver_dir), ) - m._variables = Variables( + new_model._variables = Variables( { name: Variable( - var.data.copy(deep=True) + var.data.copy(deep=deep) if include_solution - else var.data[src.variables.dataset_attrs].copy(deep=True), - m, + else var.data[m.variables.dataset_attrs].copy(deep=deep), + new_model, name, ) - for name, var in src.variables.items() + for name, var in m.variables.items() }, - m, + new_model, ) - m._constraints = Constraints( + new_model._constraints = Constraints( { name: Constraint( - con.data.copy(deep=True) + con.data.copy(deep=deep) if include_solution - else con.data[src.constraints.dataset_attrs].copy(deep=True), - m, + else con.data[m.constraints.dataset_attrs].copy(deep=deep), + new_model, name, ) - for name, con in src.constraints.items() + for name, con in m.constraints.items() }, - m, + new_model, ) - obj_expr = LinearExpression(src.objective.expression.data.copy(deep=True), m) - m._objective = Objective(obj_expr, m, src.objective.sense) - m._objective._value = src.objective.value if include_solution else None + obj_expr = LinearExpression(m.objective.expression.data.copy(deep=deep), new_model) + new_model._objective = Objective(obj_expr, new_model, m.objective.sense) + new_model._objective._value = m.objective.value if include_solution else None - m._parameters = src._parameters.copy(deep=True) - m._blocks = src._blocks.copy(deep=True) if src._blocks is not None else None + new_model._parameters = m._parameters.copy(deep=deep) + new_model._blocks = m._blocks.copy(deep=deep) if m._blocks is not None else None - for attr in src.scalar_attrs: + for attr in m.scalar_attrs: if include_solution or attr not in SOLVE_STATE_ATTRS: - setattr(m, attr, getattr(src, attr)) + setattr(new_model, attr, getattr(m, attr)) - return m + return new_model + + +def shallowcopy(m: Model) -> Model: + """ + Support Python's ``copy.copy`` protocol for ``Model``. + + Returns a shallow copy with independent wrapper objects that share + underlying array buffers with ``m``. Solve artifacts are excluded, + matching :meth:`Model.copy` defaults. + """ + return copy(m, include_solution=False, deep=False) + + +def deepcopy(m: Model, memo: dict[int, Any]) -> Model: + """ + Support Python's ``copy.deepcopy`` protocol for ``Model``. + + Returns a deep, structurally independent copy and records it in ``memo`` + as required by Python's copy protocol. Solve artifacts are excluded, + matching :meth:`Model.copy` defaults. + """ + new_model = copy(m, include_solution=False, deep=True) + memo[id(m)] = new_model + return new_model diff --git a/linopy/model.py b/linopy/model.py index 48ba8b02..84275c8b 100644 --- a/linopy/model.py +++ b/linopy/model.py @@ -54,6 +54,8 @@ ) from linopy.io import ( copy, + deepcopy, + shallowcopy, to_block_files, to_cupdlpx, to_file, @@ -1880,6 +1882,10 @@ def reset_solution(self) -> None: copy = copy + __copy__ = shallowcopy + + __deepcopy__ = deepcopy + to_netcdf = to_netcdf to_file = to_file diff --git a/test/test_model.py b/test/test_model.py index 51241f87..c0988c26 100644 --- a/test/test_model.py +++ b/test/test_model.py @@ -5,6 +5,7 @@ from __future__ import annotations +import copy as pycopy from pathlib import Path from tempfile import gettempdir @@ -170,7 +171,8 @@ def test_assert_model_equal() -> None: assert_model_equal(m, m) -def _build_model() -> Model: +@pytest.fixture(scope="module") +def copy_test_model() -> Model: """Small representative model used across copy tests.""" m: Model = Model() @@ -187,9 +189,17 @@ def _build_model() -> Model: return m -def test_model_copy_unsolved() -> None: +@pytest.fixture(scope="module") +def solved_copy_test_model(copy_test_model: Model) -> Model: + """Solved representative model used across solved-copy tests.""" + m = copy_test_model.copy(deep=True) + m.solve() + return m + + +def test_model_copy_unsolved(copy_test_model: Model) -> None: """Copy of unsolved model is structurally equal and independent.""" - m = _build_model() + m = copy_test_model.copy(deep=True) c = m.copy(include_solution=False) assert_model_equal(m, c) @@ -199,36 +209,127 @@ def test_model_copy_unsolved() -> None: assert "z" not in m.variables -@pytest.mark.skipif(len(available_solvers) == 0, reason="No solver installed") -def test_model_copy_solved_with_solution() -> None: - """Copy with include_solution=True preserves solve state.""" - m = _build_model() - m.solve() +def test_model_copy_unsolved_with_solution_flag(copy_test_model: Model) -> None: + """Unsolved model with include_solution=True has no extra solve artifacts.""" + m = copy_test_model.copy(deep=True) - c = m.copy(include_solution=True) - assert_model_equal(m, c) + c_include_solution = m.copy(include_solution=True) + c_exclude_solution = m.copy(include_solution=False) + assert_model_equal(c_include_solution, c_exclude_solution) + assert c_include_solution.status == "initialized" + assert c_include_solution.termination_condition == "" + assert c_include_solution.objective.value is None -@pytest.mark.skipif(len(available_solvers) == 0, reason="No solver installed") -def test_model_copy_solved_without_solution() -> None: - """Copy with include_solution=False (default) drops solve state but preserves problem structure.""" - m = _build_model() - m.solve() - c = m.copy(include_solution=False) +def test_model_copy_shallow(copy_test_model: Model) -> None: + """Shallow copy has independent wrappers sharing underlying data buffers.""" + m = copy_test_model.copy(deep=True) + c = m.copy(deep=False) + + assert c is not m + assert c.variables is not m.variables + assert c.constraints is not m.constraints + assert c.objective is not m.objective + + # wrappers are distinct, but shallow copy shares payload buffers + c.variables["x"].lower.values[0, 0] = 123.0 + assert m.variables["x"].lower.values[0, 0] == 123.0 + + +def test_model_deepcopy_protocol(copy_test_model: Model) -> None: + """copy.deepcopy(model) dispatches to Model.__deepcopy__ and stays independent.""" + m = copy_test_model.copy(deep=True) + c = pycopy.deepcopy(m) + + assert_model_equal(m, c) + + # Test independence: mutations to copy do not affect source + # 1. Variable mutation: add new variable + c.add_variables(name="z") + assert "z" not in m.variables - # solve state is dropped - assert c.status == "initialized" - assert c.termination_condition == "" - assert c.objective.value is None - - # problem structure is preserved — compare only dataset_attrs to exclude solution/dual - for v in m.variables: - assert_equal( - c.variables[v].data[c.variables.dataset_attrs], - m.variables[v].data[m.variables.dataset_attrs], - ) - for con in m.constraints: - assert_conequal(c.constraints[con], m.constraints[con], strict=False) - assert_linequal(c.objective.expression, m.objective.expression) - assert c.objective.sense == m.objective.sense + # 2. Variable data mutation (bounds): verify buffers are independent + original_lower = m.variables["x"].lower.values[0, 0].item() + new_lower = 999 + c.variables["x"].lower.values[0, 0] = new_lower + assert c.variables["x"].lower.values[0, 0] == new_lower + assert m.variables["x"].lower.values[0, 0] == original_lower + + # 3. Constraint coefficient mutation: deep copy must not leak back + original_con_coeff = m.constraints["con0"].coeffs.values.flat[0].item() + new_con_coeff = original_con_coeff + 42 + c.constraints["con0"].coeffs.values.flat[0] = new_con_coeff + assert c.constraints["con0"].coeffs.values.flat[0] == new_con_coeff + assert m.constraints["con0"].coeffs.values.flat[0] == original_con_coeff + + # 4. Objective expression coefficient mutation: deep copy must not leak back + original_obj_coeff = m.objective.expression.coeffs.values.flat[0].item() + new_obj_coeff = original_obj_coeff + 20 + c.objective.expression.coeffs.values.flat[0] = new_obj_coeff + assert c.objective.expression.coeffs.values.flat[0] == new_obj_coeff + assert m.objective.expression.coeffs.values.flat[0] == original_obj_coeff + + # 5. Objective sense mutation + original_sense = m.objective.sense + c.objective.sense = "max" + assert c.objective.sense == "max" + assert m.objective.sense == original_sense + + +@pytest.mark.skipif(not available_solvers, reason="No solver installed") +class TestModelCopySolved: + def test_model_deepcopy_protocol_excludes_solution( + self, solved_copy_test_model: Model + ) -> None: + """copy.deepcopy on solved model drops solve state by default.""" + m = solved_copy_test_model + + c = pycopy.deepcopy(m) + + assert c.status == "initialized" + assert c.termination_condition == "" + assert c.objective.value is None + + for v in m.variables: + assert_equal( + c.variables[v].data[c.variables.dataset_attrs], + m.variables[v].data[m.variables.dataset_attrs], + ) + for con in m.constraints: + assert_conequal(c.constraints[con], m.constraints[con], strict=False) + assert_linequal(c.objective.expression, m.objective.expression) + assert c.objective.sense == m.objective.sense + + def test_model_copy_solved_with_solution( + self, solved_copy_test_model: Model + ) -> None: + """Copy with include_solution=True preserves solve state.""" + m = solved_copy_test_model + + c = m.copy(include_solution=True) + assert_model_equal(m, c) + + def test_model_copy_solved_without_solution( + self, solved_copy_test_model: Model + ) -> None: + """Copy with include_solution=False (default) drops solve state but preserves problem structure.""" + m = solved_copy_test_model + + c = m.copy(include_solution=False) + + # solve state is dropped + assert c.status == "initialized" + assert c.termination_condition == "" + assert c.objective.value is None + + # problem structure is preserved — compare only dataset_attrs to exclude solution/dual + for v in m.variables: + assert_equal( + c.variables[v].data[c.variables.dataset_attrs], + m.variables[v].data[m.variables.dataset_attrs], + ) + for con in m.constraints: + assert_conequal(c.constraints[con], m.constraints[con], strict=False) + assert_linequal(c.objective.expression, m.objective.expression) + assert c.objective.sense == m.objective.sense From 49e9246371dba823a6ab28b29054226ee419441b Mon Sep 17 00:00:00 2001 From: Bobby Xiong Date: Fri, 20 Mar 2026 15:36:01 +0100 Subject: [PATCH 14/30] Added release notes. --- doc/release_notes.rst | 1 + 1 file changed, 1 insertion(+) diff --git a/doc/release_notes.rst b/doc/release_notes.rst index 35b21c67..a138fbbb 100644 --- a/doc/release_notes.rst +++ b/doc/release_notes.rst @@ -4,6 +4,7 @@ Release Notes Upcoming Version ---------------- +* Add ``Model.copy()`` method with ``deep`` and ``include_solution`` options; support Python ``copy.copy`` and ``copy.deepcopy`` protocols via ``__copy__`` and ``__deepcopy__``. * Harmonize coordinate alignment for operations with subset/superset objects: - Multiplication and division fill missing coords with 0 (variable doesn't participate) - Addition and subtraction of constants fill missing coords with 0 (identity element) and pin result to LHS coords From 4089a85226638caaee8974388cfcf7a351104ee8 Mon Sep 17 00:00:00 2001 From: Bobby Xiong Date: Fri, 20 Mar 2026 15:50:59 +0100 Subject: [PATCH 15/30] Made Model.copy defaulting to deep copy more explicit. --- linopy/io.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/linopy/io.py b/linopy/io.py index fd3c536a..f083b685 100644 --- a/linopy/io.py +++ b/linopy/io.py @@ -1249,6 +1249,9 @@ def copy(m: Model, include_solution: bool = False, deep: bool = True) -> Model: parameters, blocks, and scalar attributes are copied to a fully independent model. With ``deep=False``, returns a shallow copy. + :meth:`Model.copy` defaults to deep copy for workflow safety, while + Python's ``copy.copy(model)`` performs a shallow copy via ``__copy__``. + Solver runtime metadata (for example, ``solver_name`` and ``solver_model``) is intentionally not copied. Solver backend state is recreated on ``solve()``. From 0199d67406b6bb911eaf7d13c9e313e73b317a31 Mon Sep 17 00:00:00 2001 From: Bobby Xiong Date: Fri, 20 Mar 2026 16:00:15 +0100 Subject: [PATCH 16/30] Fine-tuned docs and added to read the docs api.rst. --- doc/api.rst | 1 + doc/release_notes.rst | 2 +- linopy/io.py | 16 ++++++++-------- 3 files changed, 10 insertions(+), 9 deletions(-) diff --git a/doc/api.rst b/doc/api.rst index 20958857..1554ce60 100644 --- a/doc/api.rst +++ b/doc/api.rst @@ -24,6 +24,7 @@ Creating a model piecewise.segments model.Model.linexpr model.Model.remove_constraints + model.Model.copy Classes under the hook diff --git a/doc/release_notes.rst b/doc/release_notes.rst index a138fbbb..26007b5c 100644 --- a/doc/release_notes.rst +++ b/doc/release_notes.rst @@ -4,7 +4,7 @@ Release Notes Upcoming Version ---------------- -* Add ``Model.copy()`` method with ``deep`` and ``include_solution`` options; support Python ``copy.copy`` and ``copy.deepcopy`` protocols via ``__copy__`` and ``__deepcopy__``. +* Add ``Model.copy()`` (default deep copy) with ``deep`` and ``include_solution`` options; support Python ``copy.copy`` and ``copy.deepcopy`` protocols via ``__copy__`` and ``__deepcopy__``. * Harmonize coordinate alignment for operations with subset/superset objects: - Multiplication and division fill missing coords with 0 (variable doesn't participate) - Addition and subtraction of constants fill missing coords with 0 (identity element) and pin result to LHS coords diff --git a/linopy/io.py b/linopy/io.py index f083b685..e7353b60 100644 --- a/linopy/io.py +++ b/linopy/io.py @@ -1249,8 +1249,9 @@ def copy(m: Model, include_solution: bool = False, deep: bool = True) -> Model: parameters, blocks, and scalar attributes are copied to a fully independent model. With ``deep=False``, returns a shallow copy. - :meth:`Model.copy` defaults to deep copy for workflow safety, while - Python's ``copy.copy(model)`` performs a shallow copy via ``__copy__``. + :meth:`Model.copy` defaults to deep copy for workflow safety. + In contrast, ``copy.copy(model)`` is shallow via ``__copy__``, and + ``copy.deepcopy(model)`` is deep via ``__deepcopy__``. Solver runtime metadata (for example, ``solver_name`` and ``solver_model``) is intentionally not copied. Solver backend state @@ -1261,12 +1262,11 @@ def copy(m: Model, include_solution: bool = False, deep: bool = True) -> Model: m : Model The model to copy. include_solution : bool, optional - Whether to include the current solution and dual values in the copy. - If False (default), the copy is returned in an initialized state: - solution and dual data are excluded, objective value is set to None, - and status is set to 'initialized'. If True, solution, dual values, - solve status, and objective value are also copied. If the model is - unsolved, this has no additional effect. + Whether to include solution and dual values in the copy. + If False (default), solve artifacts are excluded: solution/dual data, + objective value, and solve status are reset to initialized state. + If True, these values are copied when present. For unsolved models, + this has no additional effect. deep : bool, optional Whether to return a deep copy (default) or shallow copy. If False, the returned model uses independent wrapper objects that share From 3e0f5124994f208538ccb067149632795a8cacc1 Mon Sep 17 00:00:00 2001 From: Bryn Pickering <17178478+brynpickering@users.noreply.github.com> Date: Wed, 8 Apr 2026 12:53:19 +0100 Subject: [PATCH 17/30] Clean up some methods in `dual` module (#629) * Clean up some methods in `dual` module * Typing fix (L434). * Ignore type assessment in expressions.py --------- Co-authored-by: Bobby Xiong Co-authored-by: Lukas Trippe --- linopy/dual.py | 219 +++++++++++++++++++++--------------------- linopy/expressions.py | 2 +- 2 files changed, 113 insertions(+), 108 deletions(-) diff --git a/linopy/dual.py b/linopy/dual.py index 0a787e68..1a81ee22 100644 --- a/linopy/dual.py +++ b/linopy/dual.py @@ -7,10 +7,9 @@ from __future__ import annotations import logging -from typing import TYPE_CHECKING, Any +from typing import TYPE_CHECKING, Any, Literal import numpy as np -import pandas as pd import xarray as xr from linopy.expressions import LinearExpression @@ -21,7 +20,86 @@ logger = logging.getLogger(__name__) -def _var_lookup(m: Model) -> dict: +def _skip( + da: xr.DataArray, component_type: Literal["variable", "constraint"], name: str +) -> bool: + """ + Determine whether to skip processing a variable or constraint based on its labels. + + Parameters + ---------- + da : xr.DataArray + The labels DataArray of the variable or constraint. + component_type : Literal["variable", "constraint"] + The type of component being checked, used for logging. + name : str + The name of the variable or constraint, used for logging. + + Returns + ------- + bool + True if the component should be skipped (empty or fully masked), False otherwise. + """ + if da.size == 0: + logger.debug(f"Skipping empty {component_type} '{name}'.") + return True + + if (da == -1).all(): + logger.debug(f"{component_type} '{name}' is fully masked, skipping.") + return True + return False + + +def _lookup( + labels: xr.DataArray, name: str, component_type: Literal["variable", "constraint"] +) -> dict[int, tuple[str, dict]]: + """ + Create a lookup dictionary mapping labels to their corresponding names and coordinates. + + Parameters + ---------- + labels : xr.DataArray + Array of labels. + name : str + Name of the component. + component_type : Literal["variable", "constraint"] + Type of the component. + + Returns + ------- + dict[int, tuple[str, dict]] + Mapping from flat integer label to (name, coord_dict) tuple. + """ + lookup: dict[int, tuple[str, dict]] = {} + + vals = labels.values + if _skip(labels, component_type, name): + return lookup + + logger.debug( + f"Creating label lookup for {component_type} '{name}' with shape {labels.shape} and dims {labels.dims}." + ) + + if labels.ndim == 0: + lookup[int(vals.item())] = (name, {}) + return lookup + + coord_values = [labels.coords[d].values for d in labels.dims] + + # Choosing np.ndindex over np.argwhere or da.to_series for memory efficiency on large n-dimensional arrays + for idx in np.ndindex(vals.shape): + label = int(vals[idx]) + if label == -1: + continue + lookup[label] = ( + name, + {dim: coord_values[i][idx[i]] for i, dim in enumerate(labels.dims)}, + ) + + return lookup + + +def _var_lookup(m: Model) -> dict[int, tuple[str, dict]]: """ Build a flat label -> (var_name, coord_dict) lookup for all variables in m. @@ -43,40 +121,12 @@ def _var_lookup(m: Model) -> dict: var_lookup = {} logger.debug("Building variable label lookup.") for var_name, var in m.variables.items(): - labels = var.labels - flat_labels = labels.values.flatten() - - if len(flat_labels) == 0: - logger.debug(f"Skipping empty variable '{var_name}'.") - continue - if not (flat_labels != -1).any(): - logger.debug(f"Variable '{var_name}' is fully masked, skipping.") - continue - - logger.debug( - f"Creating label lookup for variable '{var_name}' with shape {labels.shape} and dims {labels.dims}." - ) - - coord_arrays = ( - np.meshgrid( - *[labels.coords[dim].values for dim in labels.dims], indexing="ij" - ) - if len(labels.dims) > 0 - else [] - ) - flat_coords = [arr.flatten() for arr in coord_arrays] - - for k, flat in enumerate(flat_labels): - if flat != -1: - var_lookup[int(flat)] = ( - var_name, - {dim: flat_coords[i][k] for i, dim in enumerate(labels.dims)}, - ) - + lookup = _lookup(var.labels, var_name, "variable") + var_lookup.update(lookup) return var_lookup -def _con_lookup(m: Model) -> dict: +def _con_lookup(m: Model) -> dict[int, tuple[str, dict]]: """ Build a flat label -> (con_name, coord_dict) lookup for all constraints in m. @@ -98,36 +148,8 @@ def _con_lookup(m: Model) -> dict: con_lookup = {} logger.debug("Building constraint label lookup.") for con_name, con in m.constraints.items(): - labels = con.labels - flat_labels = labels.values.flatten() - - if len(flat_labels) == 0: - logger.debug(f"Skipping empty constraint '{con_name}'.") - continue - if not (flat_labels != -1).any(): - logger.debug(f"Constraint '{con_name}' is fully masked, skipping.") - continue - - logger.debug( - f"Creating label lookup for constraint '{con_name}' with shape {labels.shape} and dims {labels.dims}." - ) - - coord_arrays = ( - np.meshgrid( - *[labels.coords[dim].values for dim in labels.dims], indexing="ij" - ) - if len(labels.dims) > 0 - else [] - ) - flat_coords = [arr.flatten() for arr in coord_arrays] - - for k, flat in enumerate(flat_labels): - if flat != -1: - con_lookup[int(flat)] = ( - con_name, - {dim: flat_coords[i][k] for i, dim in enumerate(labels.dims)}, - ) - + lookup = _lookup(con.labels, con_name, "constraint") + con_lookup.update(lookup) return con_lookup @@ -228,44 +250,35 @@ def _add_dual_variables(m: Model, m2: Model) -> dict: dual_vars = {} for name, con in m.constraints.items(): - sign_vals = con.sign.values.flatten() - - if len(sign_vals) == 0: - logger.warning(f"Constraint '{name}' has no sign values, skipping.") + if _skip(con.labels, "constraint", name): continue mask = con.labels != -1 - if not mask.any(): - logger.debug(f"Constraint '{name}' is fully masked, skipping.") - continue - - if sign_vals[0] == "=": - lower, upper = -np.inf, np.inf - var_type = "free" - elif sign_vals[0] == "<=": - lower, upper = (-np.inf, 0) if primal_is_min else (0, np.inf) - var_type = "non-positive" if primal_is_min else "non-negative" - elif sign_vals[0] == ">=": - lower, upper = (0, np.inf) if primal_is_min else (-np.inf, 0) - var_type = "non-negative" if primal_is_min else "non-positive" - else: - logger.warning( - f"Constraint '{name}' has unrecognized sign '{sign_vals[0]}', skipping." - ) - continue + sign = con.sign.isel({d: 0 for d in con.sign.dims}).item() + + match sign: + case "=": + lower, upper = -np.inf, np.inf + var_type = "free" + case "<=": + lower, upper = (-np.inf, 0) if primal_is_min else (0, np.inf) + var_type = "non-positive" if primal_is_min else "non-negative" + case ">=": + lower, upper = (0, np.inf) if primal_is_min else (-np.inf, 0) + var_type = "non-negative" if primal_is_min else "non-positive" + case _: + logger.warning( + f"Constraint '{name}' has unrecognized sign '{sign}', skipping." + ) + continue logger.debug( f"Adding {var_type} dual variable for constraint '{name}' with shape {con.shape} and dims {con.labels.dims}." ) - coords = ( - [con.labels.coords[dim] for dim in con.labels.dims] - if con.labels.dims - else None - ) dual_vars[name] = m2.add_variables( lower=lower, upper=upper, - coords=coords, + coords=con.labels.coords, name=name, mask=mask, ) @@ -409,9 +422,6 @@ def _add_dual_feasibility_constraints( # add dual feasibility constraints to m2 logger.debug("Adding dual feasibility constraints to model.") for var_name, var in m.variables.items(): - coords = [ - pd.Index(var.labels.coords[dim].values, name=dim) for dim in var.labels.dims - ] mask = var.labels != -1 c_vals = xr.DataArray( @@ -419,19 +429,14 @@ def _add_dual_feasibility_constraints( coords=var.labels.coords, ) - def rule( - m: Model, - *coord_vals: Any, - vname: str = var_name, - vdims: tuple = var.labels.dims, - ) -> LinearExpression | None: - coord_dict = dict(zip(vdims, coord_vals)) + def __rule(m: Model, *coord_vals: Any) -> LinearExpression | None: + coord_dict = { + str(dim): val for dim, val in zip(var.labels.dims, coord_vals) + } flat = var.labels.sel(**coord_dict).item() - if flat == -1: - return None - if flat not in dual_feas_terms[vname]: + if flat == -1 or flat not in (term_dict := dual_feas_terms[var_name]): return None - _, terms, _ = dual_feas_terms[vname][flat] + _, terms, _ = term_dict[flat] if not terms: return None return sum( @@ -439,7 +444,7 @@ def rule( for con_name, con_coords, coeff in terms ) - lhs = LinearExpression.from_rule(m2, rule, coords) + lhs = LinearExpression.from_rule(m2, __rule, var.labels.coords) m2.add_constraints(lhs == c_vals, name=var_name, mask=mask) diff --git a/linopy/expressions.py b/linopy/expressions.py index ca491c3e..382315a1 100644 --- a/linopy/expressions.py +++ b/linopy/expressions.py @@ -2384,7 +2384,7 @@ def merge( has_quad_expression = any(type(e) is QuadraticExpression for e in exprs) has_linear_expression = any(type(e) is LinearExpression for e in exprs) if cls is None: - cls = QuadraticExpression if has_quad_expression else LinearExpression + cls = QuadraticExpression if has_quad_expression else LinearExpression # type: ignore[assignment] if cls is QuadraticExpression and dim == TERM_DIM and has_linear_expression: raise ValueError( From 70d0cabcaba4494564a55e1e6b5901181c6e70dc Mon Sep 17 00:00:00 2001 From: Bobby Xiong Date: Wed, 8 Apr 2026 14:03:58 +0200 Subject: [PATCH 18/30] fix: use dimension coords only when adding dual variables to avoid scalar coord conflicts. --- linopy/dual.py | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/linopy/dual.py b/linopy/dual.py index 1a81ee22..dd970d0d 100644 --- a/linopy/dual.py +++ b/linopy/dual.py @@ -275,14 +275,18 @@ def _add_dual_variables(m: Model, m2: Model) -> dict: logger.debug( f"Adding {var_type} dual variable for constraint '{name}' with shape {con.shape} and dims {con.labels.dims}." ) + coords = ( + [con.labels.coords[dim] for dim in con.labels.dims] + if con.labels.dims + else None + ) dual_vars[name] = m2.add_variables( lower=lower, upper=upper, - coords=con.labels.coords, + coords=coords, name=name, mask=mask, ) - return dual_vars From 066de31f1275ed2b9a45a2059fe8d786d1003932 Mon Sep 17 00:00:00 2001 From: bobbyxng Date: Mon, 1 Jun 2026 13:04:07 +0200 Subject: [PATCH 19/30] fix CI testing issues --- .github/workflows/test-models.yml | 4 ++-- linopy/expressions.py | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/.github/workflows/test-models.yml b/.github/workflows/test-models.yml index 14eedae5..5444455f 100644 --- a/.github/workflows/test-models.yml +++ b/.github/workflows/test-models.yml @@ -77,12 +77,12 @@ jobs: uses: actions/cache@v5 with: path: ${{ env.CONDA }}/envs - key: conda-pypsa-eur-${{ env.week }}-${{ hashFiles('envs/linux-64.lock.yaml') }} + key: conda-pypsa-eur-${{ env.week }}-${{ hashFiles('envs/environment.yaml') }} id: cache-env - name: Update environment if: env.pinned == 'false' && steps.cache-env.outputs.cache-hit != 'true' - run: conda env update -n pypsa-eur -f envs/linux-64.lock.yaml + run: conda env update -n pypsa-eur -f envs/environment.yaml - name: Install package from ref if: env.pinned == 'false' diff --git a/linopy/expressions.py b/linopy/expressions.py index c96d0816..31868234 100644 --- a/linopy/expressions.py +++ b/linopy/expressions.py @@ -2386,7 +2386,7 @@ def merge( has_quad_expression = any(isinstance(e, QuadraticExpression) for e in exprs) has_linear_expression = any(isinstance(e, LinearExpression) for e in exprs) if cls is None: - cls = QuadraticExpression if has_quad_expression else LinearExpression # type: ignore[assignment] + cls = QuadraticExpression if has_quad_expression else LinearExpression if ( issubclass(cls, QuadraticExpression) From 9e3c770af4a11f1b0c35bd6396ea432cbdd58bc8 Mon Sep 17 00:00:00 2001 From: bobbyxng Date: Mon, 1 Jun 2026 22:45:06 +0200 Subject: [PATCH 20/30] refactor dualisation, update tests, modularise functions and vectorise code --- linopy/dual.py | 570 ++++++++++++++++++++++------------------------ test/test_dual.py | 240 +++++++++++++++++++ 2 files changed, 517 insertions(+), 293 deletions(-) create mode 100644 test/test_dual.py diff --git a/linopy/dual.py b/linopy/dual.py index dd970d0d..fadcc671 100644 --- a/linopy/dual.py +++ b/linopy/dual.py @@ -23,23 +23,7 @@ def _skip( da: xr.DataArray, component_type: Literal["variable", "constraint"], name: str ) -> bool: - """ - Determine whether to skip processing a variable or constraint based on its labels. - - Parameters - ---------- - da : xr.DataArray - The labels DataArray of the variable or constraint. - component_type : Literal["variable", "constraint"] - The type of component being checked, used for logging. - name : str - The name of the variable or constraint, used for logging. - - Returns - ------- - bool - True if the component should be skipped (empty or fully masked), False otherwise. - """ + """Return True if the label array is empty or entirely masked (all -1).""" if da.size == 0: logger.debug(f"Skipping empty {component_type} '{name}'.") return True @@ -50,124 +34,19 @@ def _skip( return False -def _lookup( - labels: xr.DataArray, name: str, component_type: Literal["variable", "constraint"] -) -> dict[int, tuple[str, dict]]: - """ - Create a lookup dictionary mapping labels to their corresponding names and coordinates. - - Parameters - ---------- - labels : xr.DataArray - Array of labels. - name : str - Name of the component. - component_type : Literal["variable", "constraint"] - Type of the component. - - Returns - ------- - dict[int, tuple[str, dict]] - Mapping from flat integer label to (name, coord_dict) tuple. - """ - lookup: dict[int, tuple[str, dict]] = {} - - vals = labels.values - if _skip(labels, component_type, name): - return lookup - - logger.debug( - f"Creating label lookup for {component_type} '{name}' with shape {labels.shape} and dims {labels.dims}." - ) - - if labels.ndim == 0: - lookup[int(vals.item())] = (name, {}) - return lookup - - coord_values = [labels.coords[d].values for d in labels.dims] - - # Choosing np.ndindex over np.argwhere or da.to_series for memory efficiency on large n-dimensional arrays - for idx in np.ndindex(vals.shape): - label = int(vals[idx]) - if label == -1: - continue - lookup[label] = ( - name, - {dim: coord_values[i][idx[i]] for i, dim in enumerate(labels.dims)}, - ) - - return lookup - - -def _var_lookup(m: Model) -> dict[int, tuple[str, dict]]: - """ - Build a flat label -> (var_name, coord_dict) lookup for all variables in m. - - Used to map entries in m.matrices.vlabels back to their variable name - and xarray coordinates for use in dual feasibility constraint construction. - - Skips masked entries (label == -1) and empty variables. - - Parameters - ---------- - m : Model - Primal linopy model. - - Returns - ------- - dict - Mapping from flat integer label to (var_name, coord_dict) tuple. - """ - var_lookup = {} - logger.debug("Building variable label lookup.") - for var_name, var in m.variables.items(): - lookup = _lookup(var.labels, var_name, "variable") - var_lookup.update(lookup) - return var_lookup - - -def _con_lookup(m: Model) -> dict[int, tuple[str, dict]]: - """ - Build a flat label -> (con_name, coord_dict) lookup for all constraints in m. - - Used to map entries in m.matrices.clabels back to their constraint name - and xarray coordinates for use in dual feasibility constraint construction. - - Skips masked entries (label == -1) and empty or fully-masked constraints. - - Parameters - ---------- - m : Model - Primal linopy model. - - Returns - ------- - dict - Mapping from flat integer label to (con_name, coord_dict) tuple. - """ - con_lookup = {} - logger.debug("Building constraint label lookup.") - for con_name, con in m.constraints.items(): - lookup = _lookup(con.labels, con_name, "constraint") - con_lookup.update(lookup) - return con_lookup - - def bounds_to_constraints(m: Model) -> None: """ - Add explicit bound constraints for variables with bounds set directly - in the variable rather than via explicit constraints. + Convert finite variable bounds to explicit ``>=`` / ``<=`` constraints. - Adds constraints named '{var_name}-bound-lower' and '{var_name}-bound-upper' - to distinguish from PyPSA's automatic '-fix-*' constraints. - - Also resets variable bounds to [-inf, inf] after adding constraints, - to avoid double-counting in the dual. + Each finite lower bound becomes a ``'{var_name}-bound-lower'`` constraint and + each finite upper bound becomes a ``'{var_name}-bound-upper'`` constraint. + The variable bounds are then relaxed to ``[-inf, inf]`` to avoid + double-counting when the dual is formed. Parameters ---------- m : Model - Linopy model to convert variable bounds to constraints. Mutates the model in-place. + Model to mutate in-place. """ logger.debug("Converting variable bounds to explicit constraints.") logger.debug("Relaxing variable bounds to [-inf, inf].") @@ -215,36 +94,34 @@ def bounds_to_constraints(m: Model) -> None: def _add_dual_variables(m: Model, m2: Model) -> dict: """ - Add dual variables to m2 corresponding to constraints in m.. + Add one dual variable to m2 for each active constraint in m. - For each active constraint in m, adds a dual variable to m2 following - standard LP duality sign conventions. The sign of the dual variable bounds - depends on both the constraint type and the primal objective sense: + Dual variable bounds encode the sign convention for each constraint type + and primal objective sense: - For a minimization primal: - - Equality constraints (=) -> free dual variable (lower=-inf, upper=inf) - - <= constraints -> non-positive dual variable (lower=-inf, upper=0) - - >= constraints -> non-negative dual variable (lower=0, upper=inf) + ============ =========== ================ ================ + Constraint Primal sense lower upper + ============ =========== ================ ================ + = min / max -inf +inf (free) + <= min -inf 0 + <= max 0 +inf + >= min 0 +inf + >= max -inf 0 + ============ =========== ================ ================ - For a maximization primal: - - Equality constraints (=) -> free dual variable (lower=-inf, upper=inf) - - <= constraints -> non-negative dual variable (lower=0, upper=inf) - - >= constraints -> non-positive dual variable (lower=-inf, upper=0) - - Skips constraints with no active rows (empty or fully masked). + Fully masked or empty constraints are skipped. Parameters ---------- m : Model - Primal linopy model containing the constraints to dualize. + Primal model. m2 : Model - Dual linopy model to which dual variables are added. + Dual model to populate. Returns ------- dict - Mapping from constraint name (str) to the corresponding dual - variable (linopy.Variable) in m2. + ``{constraint_name: dual_variable}`` for every dualized constraint. """ primal_is_min = m.objective.sense == "min" @@ -290,165 +167,282 @@ def _add_dual_variables(m: Model, m2: Model) -> dict: return dual_vars -def _build_dual_feas_terms( - m: Model, - dual_vars: dict, - var_lookup: dict, - con_lookup: dict, -) -> dict: +def _build_con_to_dual_label(m: Model, dual_vars: dict) -> np.ndarray: """ - Build dual feasibility terms for each primal variable in m. + Build a flat constraint-label → flat dual-variable-label lookup array. - For each active primal variable x_j, collects the constraint matrix - entries A_ji and their corresponding constraint names and coordinates, - forming the terms of the stationarity condition: - sum_i (A_ji * lambda_i) = c_j - - Raw constraint matrix coefficients are used directly without sign - factors, as the sign convention is encoded in the dual variable bounds: - - <= constraints: lambda_i <= 0 - - >= constraints: lambda_i >= 0 - - = constraints: lambda_i free + ``result[flat_con_label]`` gives the flat dual variable label in m2, + or -1 for constraints not in *dual_vars* or with masked labels. Parameters ---------- m : Model - Primal linopy model. + Primal model. dual_vars : dict - Mapping from constraint name to dual variable in m2, - as returned by _add_dual_variables(). Used to skip constraints - that were not dualized (e.g. empty or fully masked). - var_lookup : dict - Mapping from flat variable label to (var_name, coord_dict), - as returned by _var_lookup(). - con_lookup : dict - Mapping from flat constraint label to (con_name, coord_dict), - as returned by _con_lookup(). + ``{constraint_name: dual_variable}`` as returned by _add_dual_variables(). Returns ------- - dict - Nested dict: {var_name: {flat_label: (var_coords, terms, obj_coeff)}} - where terms is a list of (con_name, con_coords, coeff) tuples. + np.ndarray of int64 + Lookup array of length ``max_flat_con_label + 1``; empty if no valid + constraint labels exist. """ - A = m.matrices.A - if A is None: - raise ValueError("Constraint matrix is None, model has no constraints.") - A_csc = A.tocsc() - c = m.matrices.c - indptr = A_csc.indptr - indices = A_csc.indices - data = A_csc.data - vlabels = m.matrices.vlabels - clabels = m.matrices.clabels - - dual_feas_terms: dict[str, dict[int, tuple]] = { - var_name: {} for var_name in m.variables - } - - logger.debug("Building dual feasibility terms for each primal variable.") - - for i in range(A_csc.shape[1]): - flat_var = vlabels[i] - if flat_var == -1: - continue - if flat_var not in var_lookup: - continue - var_name, var_coords = var_lookup[flat_var] - terms = [] - for k in range(indptr[i], indptr[i + 1]): - j = indices[k] - flat_con = clabels[j] - if flat_con == -1: - continue - if flat_con not in con_lookup: - continue - con_name, con_coords = con_lookup[flat_con] - if con_name not in dual_vars: - continue - coeff = data[k] - terms.append((con_name, con_coords, coeff)) - dual_feas_terms[var_name][flat_var] = (var_coords, terms, c[i]) + max_flat_con = -1 + for con_name in dual_vars: + flat = m.constraints[con_name].labels.values.ravel() + valid = flat[flat != -1] + if len(valid): + max_flat_con = max(max_flat_con, int(valid.max())) + + if max_flat_con < 0: + return np.array([], dtype=np.int64) + + lookup = np.full(max_flat_con + 1, -1, dtype=np.int64) + for con_name, dv in dual_vars.items(): + con_flat = m.constraints[con_name].labels.values.ravel().astype(np.int64) + dv_flat = dv.labels.values.ravel().astype(np.int64) + valid = (con_flat != -1) & (dv_flat != -1) + if valid.any(): + lookup[con_flat[valid]] = dv_flat[valid] + + return lookup + + +def _extract_dual_feas_entries( + A: Any, + vlabels: np.ndarray, + clabels: np.ndarray, + flat_con_to_dual: np.ndarray, +) -> tuple[np.ndarray, np.ndarray, np.ndarray]: + """ + Return ``(flat_var_label, flat_dual_label, coeff)`` for each active NNZ entry of A. + + Converts A to COO, maps rows/columns to flat labels via *clabels* / *vlabels*, + and drops entries that are masked or belong to a constraint without a dual variable. + + Parameters + ---------- + A : scipy sparse matrix + Primal constraint matrix (any scipy sparse format). + vlabels : np.ndarray of int64 + Flat variable label per column of A; -1 for masked columns. + clabels : np.ndarray of int64 + Flat constraint label per row of A; -1 for masked rows. + flat_con_to_dual : np.ndarray of int64 + Lookup array as returned by _build_con_to_dual_label(). + + Returns + ------- + flat_v : np.ndarray of int64 + Flat primal variable label for each retained entry. + flat_d : np.ndarray of int64 + Corresponding flat dual variable label. + coeffs : np.ndarray of float64 + Corresponding A matrix coefficient. + """ + A_coo = A.tocoo() + flat_v = vlabels[A_coo.col].astype(np.int64) + flat_c = clabels[A_coo.row].astype(np.int64) + coeffs = A_coo.data + + # Drop entries where either label is masked. + active = (flat_v != -1) & (flat_c != -1) + flat_v, flat_c, coeffs = flat_v[active], flat_c[active], coeffs[active] + + # Map primal constraint labels to flat dual variable labels. + n = len(flat_con_to_dual) + in_range = flat_c < n + flat_d = np.full(len(flat_c), -1, dtype=np.int64) + flat_d[in_range] = flat_con_to_dual[flat_c[in_range]] + + # Drop entries with no corresponding dual variable. + has_dual = flat_d != -1 + return flat_v[has_dual], flat_d[has_dual], coeffs[has_dual] + + +def _build_obj_coeff_lookup(vlabels: np.ndarray, c_vec: np.ndarray) -> np.ndarray: + """ + Build a flat variable-label → objective-coefficient lookup array. + + ``result[flat_var_label]`` gives the objective coefficient; unlabelled + positions are 0.0. + + Parameters + ---------- + vlabels : np.ndarray of int64 + Flat variable label per column of A; -1 for masked entries. + c_vec : np.ndarray of float64 + Objective coefficient in the same column ordering. + + Returns + ------- + np.ndarray of float64 + Lookup array of length ``max_flat_var_label + 1``; empty if no valid + variable labels exist. + """ + valid_mask = vlabels != -1 + valid_vlabels = vlabels[valid_mask] + if not len(valid_vlabels): + return np.array([], dtype=np.float64) + + lookup = np.zeros(int(valid_vlabels.max()) + 1, dtype=np.float64) + lookup[valid_vlabels.astype(np.int64)] = c_vec[valid_mask] + return lookup - return dual_feas_terms + +def _make_feas_lhs( + var: Any, + flat_v: np.ndarray, + flat_d: np.ndarray, + nnz_data: np.ndarray, + m2: Model, +) -> LinearExpression: + """ + Build ``sum_i(A_ji * lambda_i)`` as a LinearExpression over *var*'s coordinates. + + Parameters + ---------- + var : linopy.Variable + Primal variable whose stationarity LHS is being constructed. + flat_v : np.ndarray of int64 + Flat primal variable labels for all active NNZ entries, as returned + by _extract_dual_feas_entries(). + flat_d : np.ndarray of int64 + Corresponding flat dual variable labels. + nnz_data : np.ndarray of float64 + Corresponding A matrix coefficients. + m2 : Model + Dual model owning the dual variables. + + Returns + ------- + LinearExpression + Stationarity LHS over *var*'s coordinate space. Returns an empty + expression (constant zero) when the variable has no constraint connections. + """ + var_flat = var.labels.values.ravel().astype(np.int64) + n_elements = len(var_flat) + valid_mask = var_flat != -1 + + if not valid_mask.any(): + return LinearExpression(None, m2) + + # Map flat variable labels to ravelled positions (unique per variable type). + max_fv = int(var_flat[valid_mask].max()) + var_to_idx = np.full(max_fv + 1, -1, dtype=np.int64) + var_to_idx[var_flat[valid_mask]] = np.where(valid_mask)[0].astype(np.int64) + + # Locate NNZ entries that belong to this variable via bounded lookup. + lin_idx = np.full(len(flat_v), -1, dtype=np.int64) + in_range = flat_v <= max_fv + if in_range.any(): + lin_idx[in_range] = var_to_idx[flat_v[in_range]] + + keep = lin_idx != -1 + lin_idx_f, flat_d_f, nnz_f = lin_idx[keep], flat_d[keep], nnz_data[keep] + + if not len(lin_idx_f): + return LinearExpression(None, m2) + + # Sort by ravelled position so entries for the same element are contiguous. + order = np.argsort(lin_idx_f, kind="stable") + lin_idx_s = lin_idx_f[order] + flat_d_s = flat_d_f[order] + nnz_s = nnz_f[order] + + # Compute within-group term-slot positions (fully vectorised, no Python loop). + # group_start[i] is True at the first entry of each new element group. + group_start = np.empty(len(lin_idx_s), dtype=bool) + group_start[0] = True + group_start[1:] = lin_idx_s[1:] != lin_idx_s[:-1] + group_ids = np.cumsum(group_start) - 1 + group_start_pos = np.where(group_start)[0] + col_idx = np.arange(len(lin_idx_s), dtype=np.int64) - group_start_pos[group_ids] + + max_terms = int(col_idx.max()) + 1 + + # Populate (n_elements, max_terms) arrays via advanced indexing (C-level). + dual_labels_2d = np.full((n_elements, max_terms), -1, dtype=np.int64) + dual_coeffs_2d = np.zeros((n_elements, max_terms), dtype=np.float64) + dual_labels_2d[lin_idx_s, col_idx] = flat_d_s + dual_coeffs_2d[lin_idx_s, col_idx] = nnz_s + + # Wrap in a LinearExpression, reshaping to (*var_dims, _term). + target_shape = var.labels.shape + (max_terms,) + dims = list(var.labels.dims) + ["_term"] + ds = xr.Dataset( + { + "vars": xr.DataArray(dual_labels_2d.reshape(target_shape), dims=dims), + "coeffs": xr.DataArray(dual_coeffs_2d.reshape(target_shape), dims=dims), + }, + coords={dim: var.labels.coords[dim] for dim in var.labels.dims}, + ) + return LinearExpression(ds, m2) def _add_dual_feasibility_constraints( m: Model, m2: Model, dual_vars: dict, - var_lookup: dict, - con_lookup: dict, ) -> None: """ - Add dual feasibility constraints to m2. + Add the stationarity constraint ``sum_i(A_ji * lambda_i) = c_j`` for each primal variable. - For each primal variable x_j in m, adds the stationarity constraint: - sum_i (A_ji * lambda_i) = c_j - where: - - A is the primal constraint matrix - - lambda_i are the dual variables in m2 - - c_j is the objective coefficient of x_j - - Raw constraint matrix coefficients are used directly without sign factors, - because the sign convention is encoded in the dual variable bounds: - - <= constraints: lambda_i <= 0 - - >= constraints: lambda_i >= 0 - - = constraints: lambda_i free - - Skips masked variable entries (label == -1) and variables not present - in var_lookup (e.g. from empty constraints). + Sign conventions are already encoded in the dual variable bounds produced by + _add_dual_variables(), so raw A coefficients are used without adjustment. + Variables with no constraint connections (e.g. unconstrained free variables) + are skipped; the dual is infeasible for such problems if ``c_j != 0``. Parameters ---------- m : Model - Primal linopy model. + Primal model. m2 : Model - Dual linopy model. + Dual model to populate. dual_vars : dict - Mapping from constraint name to dual variable in m2, - as returned by _add_dual_variables(). - var_lookup : dict - Mapping from flat variable label to (var_name, coord_dict), - as returned by _var_lookup(). - con_lookup : dict - Mapping from flat constraint label to (con_name, coord_dict), - as returned by _con_lookup(). + ``{constraint_name: dual_variable}`` as returned by _add_dual_variables(). """ - dual_feas_terms = _build_dual_feas_terms(m, dual_vars, var_lookup, con_lookup) + A = m.matrices.A + if A is None: + raise ValueError("Constraint matrix is None, model has no constraints.") + + vlabels = np.asarray(m.matrices.vlabels, dtype=np.int64) + clabels = np.asarray(m.matrices.clabels, dtype=np.int64) - c = m.matrices.c - vlabels = m.matrices.vlabels + flat_con_to_dual = _build_con_to_dual_label(m, dual_vars) + if not len(flat_con_to_dual): + logger.warning("No valid constraint labels found, skipping dual feasibility constraints.") + return - # build objective coefficient lookup by flat variable label - c_by_label = {vlabels[i]: c[i] for i in range(len(vlabels))} + flat_v, flat_d, nnz_data = _extract_dual_feas_entries(A, vlabels, clabels, flat_con_to_dual) + c_lookup = _build_obj_coeff_lookup(vlabels, m.matrices.c) - # add dual feasibility constraints to m2 - logger.debug("Adding dual feasibility constraints to model.") + logger.debug("Building dual feasibility constraints for each primal variable.") for var_name, var in m.variables.items(): + if _skip(var.labels, "variable", var_name): + continue + mask = var.labels != -1 + var_flat = var.labels.values.ravel().astype(np.int64) + # RHS: objective coefficient for each element of this variable. + in_c_range = (var_flat != -1) & (var_flat < len(c_lookup)) + safe_flat = np.where(in_c_range, var_flat, 0) c_vals = xr.DataArray( - np.vectorize(lambda flat: c_by_label.get(flat, 0.0))(var.labels.values), + np.where(in_c_range, c_lookup[safe_flat], 0.0).reshape(var.labels.shape), coords=var.labels.coords, ) - def __rule(m: Model, *coord_vals: Any) -> LinearExpression | None: - coord_dict = { - str(dim): val for dim, val in zip(var.labels.dims, coord_vals) - } - flat = var.labels.sel(**coord_dict).item() - if flat == -1 or flat not in (term_dict := dual_feas_terms[var_name]): - return None - _, terms, _ = term_dict[flat] - if not terms: - return None - return sum( - coeff * dual_vars[con_name].at[tuple(con_coords.values())] - for con_name, con_coords, coeff in terms + lhs = _make_feas_lhs(var, flat_v, flat_d, nnz_data, m2) + if lhs.is_constant: + # Variable has no constraint connections (free, unconstrained variable). + # The stationarity condition 0 = c_j cannot be expressed as a linopy + # constraint; the dual is infeasible if c_j != 0, trivial if c_j == 0. + logger.debug( + f"Variable '{var_name}' has no constraint connections; " + "skipping dual feasibility constraint." ) - - lhs = LinearExpression.from_rule(m2, __rule, var.labels.coords) + continue m2.add_constraints(lhs == c_vals, name=var_name, mask=mask) @@ -459,35 +453,28 @@ def _add_dual_objective( add_objective_constant: float = 0.0, ) -> None: """ - Construct and add the dual objective to m2. - - The dual objective is sum(rhs * dual) over all constraints, added uniformly - with a + sign. The sign convention is encoded in the dual variable bounds: - - <= constraints: dual <= 0, so rhs * dual contributes negatively - - >= constraints: dual >= 0, so rhs * dual contributes positively - - = constraints: dual free + Construct and add ``sum(rhs_i * lambda_i)`` as the dual objective of m2. - This matches linopy's and Gurobi's native dual sign convention, allowing - direct comparison between m2.variables[con_name].solution and - m.constraints[con_name].dual without sign adjustments. + The uniform ``+`` sign is correct because sign conventions are already + encoded in the dual variable bounds: a ``<=`` constraint has ``lambda <= 0``, + so ``rhs * lambda`` contributes negatively without an explicit sign flip. + This aligns with linopy's native dual convention, so + ``m2.variables[con_name].solution`` can be compared directly with + ``m.constraints[con_name].dual`` after solving. - The dual objective sense is flipped relative to the primal: - - min primal -> max dual - - max primal -> min dual + The objective sense is flipped: ``min`` primal → ``max`` dual, and vice versa. Parameters ---------- m : Model - Primal linopy model. + Primal model. m2 : Model - Dual linopy model. + Dual model to populate. dual_vars : dict - Mapping from constraint name to dual variable in m2, - as returned by _add_dual_variables(). + ``{constraint_name: dual_variable}`` as returned by _add_dual_variables(). add_objective_constant : float, optional - Constant term to add to the dual objective. Use this to pass through - a primal objective constant excluded via include_objective_constant=False - during model creation. Default is 0.0. + Constant added to the dual objective, e.g. to pass through a primal + objective constant that was excluded from the model. """ dual_obj: LinearExpression = LinearExpression(None, m2) sense = "max" if m.objective.sense == "min" else "min" @@ -551,14 +538,13 @@ def dualize( Primal linopy model to dualize. Must have a linear objective and linear constraints. add_objective_constant : float, optional - Constant term to add to the dual objective. Use this to pass through - a primal objective constant. Default is 0.0. + Constant added to the dual objective, e.g. to pass through a primal + objective constant that was excluded from the model. Returns ------- Model - The dual linopy model. Dual variables are named after their - corresponding primal constraints. + Dual model whose variables are named after the primal constraints. Examples -------- @@ -580,10 +566,8 @@ def dualize( return m2 bounds_to_constraints(m1) - var_lup = _var_lookup(m1) - con_lup = _con_lookup(m1) dual_vars = _add_dual_variables(m1, m2) - _add_dual_feasibility_constraints(m1, m2, dual_vars, var_lup, con_lup) + _add_dual_feasibility_constraints(m1, m2, dual_vars) _add_dual_objective( m1, m2, dual_vars, add_objective_constant=add_objective_constant ) diff --git a/test/test_dual.py b/test/test_dual.py new file mode 100644 index 00000000..5a84fff5 --- /dev/null +++ b/test/test_dual.py @@ -0,0 +1,240 @@ +"""Tests for linopy/dual.py - LP dualization.""" + +from __future__ import annotations + +import numpy as np +import pandas as pd +import pytest +import xarray as xr + +from linopy import Model +from linopy.dual import bounds_to_constraints, dualize +from linopy.solvers import licensed_solvers + +_lp_solver = next((s for s in ("highs", "glpk", "scip") if s in licensed_solvers), None) +needs_solver = pytest.mark.skipif(_lp_solver is None, reason="No LP solver available") + + +# Structural tests (no solver required) +def test_dualize_empty_model(): + m = Model() + m2 = dualize(m) + assert len(m2.variables) == 0 + assert len(m2.constraints) == 0 + + +def test_dual_variables_named_after_primal_constraints(): + m = Model() + x = m.add_variables(lower=0, coords=[pd.RangeIndex(3)], name="x") + m.add_constraints(x >= 1.0, name="lb_con") + m.add_objective(2.0 * x) + + m2 = dualize(m) + assert "lb_con" in m2.variables + assert "x-bound-lower" in m2.variables + + +def test_dual_feasibility_constraints_named_after_primal_variables(): + m = Model() + x = m.add_variables(lower=0, coords=[pd.RangeIndex(3)], name="x") + m.add_constraints(x >= 1.0, name="lb_con") + m.add_objective(2.0 * x) + + m2 = dualize(m) + assert "x" in m2.constraints + + +def test_dual_objective_sense_flipped_min(): + m = Model() + x = m.add_variables(lower=0, name="x") + m.add_constraints(x >= 1, name="c") + m.add_objective(x) # min + + m2 = dualize(m) + assert m2.objective.sense == "max" + + +def test_dual_objective_sense_flipped_max(): + m = Model() + x = m.add_variables(upper=10, name="x") + m.add_constraints(x <= 5, name="c") + m.add_objective(x, sense="max") + + m2 = dualize(m) + assert m2.objective.sense == "min" + + +def test_dual_sign_conventions_min(): + """For a min primal: <= -> dual <= 0, >= -> dual >= 0, = -> dual free.""" + m = Model() + x = m.add_variables(lower=-np.inf, upper=np.inf, name="x") + m.add_constraints(x == 5, name="eq") + m.add_constraints(x <= 10, name="leq") + m.add_constraints(x >= -10, name="geq") + m.add_objective(x) + + m2 = dualize(m) + assert m2.variables["eq"].lower.item() == -np.inf + assert m2.variables["eq"].upper.item() == np.inf + assert m2.variables["leq"].upper.item() == 0 + assert m2.variables["geq"].lower.item() == 0 + + +def test_dual_sign_conventions_max(): + """For a max primal: <= -> dual >= 0, >= -> dual <= 0.""" + m = Model() + x = m.add_variables(lower=-np.inf, upper=np.inf, name="x") + m.add_constraints(x <= 10, name="leq") + m.add_constraints(x >= -10, name="geq") + m.add_objective(x, sense="max") + + m2 = dualize(m) + assert m2.variables["leq"].lower.item() == 0 + assert m2.variables["geq"].upper.item() == 0 + + +def test_dual_feasibility_rhs_equals_objective_coefficients(): + """The dual feasibility constraint RHS must equal the primal objective coefficient.""" + m = Model() + x = m.add_variables(lower=-np.inf, upper=np.inf, coords=[pd.RangeIndex(4)], name="x") + m.add_constraints(x == np.array([1.0, 2.0, 3.0, 4.0]), name="eq") + c = np.array([10.0, 20.0, 30.0, 40.0]) + m.add_objective(c * x) + + m2 = dualize(m) + np.testing.assert_allclose(m2.constraints["x"].rhs.values, c) + + +def test_dual_multi_constraint_per_variable(): + """Each primal variable appearing in k constraints produces k terms in its dual feas constraint.""" + m = Model() + n = 3 + x = m.add_variables(lower=-np.inf, upper=np.inf, coords=[pd.RangeIndex(n)], name="x") + # Two equality constraints, both using x + m.add_constraints(x == 1.0, name="c1") + m.add_constraints(2.0 * x == 2.0, name="c2") + m.add_objective(5.0 * x) + + m2 = dualize(m) + # dual feas constraint for x: lambda_c1 + 2*lambda_c2 = 5 + con_x = m2.constraints["x"] + # Should have 2 terms (one from c1, one from c2) + n_terms = (con_x.vars != -1).sum(dim="_term").max().item() + assert n_terms == 2 + + +def test_bounds_to_constraints_removes_bounds(): + """bounds_to_constraints should add constraints and relax variable bounds.""" + m = Model() + m.add_variables(lower=1.0, upper=5.0, coords=[pd.RangeIndex(3)], name="x") + + bounds_to_constraints(m) + + assert "x-bound-lower" in m.constraints + assert "x-bound-upper" in m.constraints + assert np.all(np.isinf(m.variables["x"].lower.values)) + assert np.all(np.isinf(m.variables["x"].upper.values)) + + +def test_dual_with_masked_variable(): + """Partially masked variable: active elements produce a constraint, masked ones do not.""" + m = Model() + mask = xr.DataArray([True, False, True], dims=["dim_0"]) + x = m.add_variables(lower=0, coords=[pd.RangeIndex(3)], name="x", mask=mask) + m.add_constraints(x >= 1.0, name="c", mask=mask) + m.add_objective(x) + + m2 = dualize(m) + assert "x" in m2.constraints + assert (m2.constraints["x"].labels != -1).sum().item() == 2 # only 2 active rows + + +def test_dual_free_unconstrained_variable_no_error(): + """A free variable with no constraint connections is silently skipped.""" + m = Model() + m.add_variables(lower=-np.inf, upper=np.inf, name="x") # no bounds → no bound constraints + y = m.add_variables(lower=-np.inf, upper=np.inf, name="y") + m.add_constraints(y == 5, name="eq") + m.add_objective(y) # x has no connections at all + + m2 = dualize(m) # must not raise + assert "y" in m2.constraints + assert "x" not in m2.constraints # no constraint connections → silently skipped + + +# Numerical tests (require solver) +def _solve(model, **kwargs): + model.solve(solver_name=_lp_solver, io_api="lp", **kwargs) + return model.objective.value + + +@needs_solver +def test_strong_duality_simple(): + """Strong duality: primal obj == dual obj at optimality.""" + m = Model() + x = m.add_variables(lower=0, name="x") + y = m.add_variables(lower=0, name="y") + m.add_constraints(x + y >= 3, name="c1") + m.add_constraints(x + 2 * y >= 4, name="c2") + m.add_objective(5 * x + 4 * y) + + primal_obj = _solve(m) + dual_obj = _solve(m.dualize()) + assert abs(primal_obj - dual_obj) < 1e-5 + + +@needs_solver +def test_strong_duality_array_variable(): + """Strong duality with multi-dimensional variables.""" + rng = np.random.default_rng(0) + n = 6 + A = rng.random((4, n)) - 0.5 + b = rng.random(4) + 0.5 + c_obj = rng.random(n) + 0.1 + + m = Model() + x = m.add_variables(lower=0, coords=[pd.RangeIndex(n)], name="x") + for i in range(4): + m.add_constraints( + sum(float(A[i, j]) * x[j] for j in range(n)) <= float(b[i]), + name=f"c{i}", + ) + m.add_objective(sum(float(c_obj[j]) * x[j] for j in range(n))) + + primal_obj = _solve(m) + dual_obj = _solve(m.dualize()) + assert abs(primal_obj - dual_obj) < 1e-4 + + +@needs_solver +def test_strong_duality_mixed_constraint_types(): + """Strong duality with =, <=, >= constraints.""" + m = Model() + x = m.add_variables(lower=-np.inf, upper=np.inf, name="x") + y = m.add_variables(lower=-np.inf, upper=np.inf, name="y") + m.add_constraints(x + y == 5, name="eq") + m.add_constraints(x - y <= 2, name="leq") + m.add_constraints(x >= 0, name="geq") + m.add_objective(2 * x + 3 * y) + + primal_obj = _solve(m) + dual_obj = _solve(m.dualize()) + assert abs(primal_obj - dual_obj) < 1e-5 + + +@needs_solver +def test_strong_duality_maximization(): + """Strong duality for a maximization primal.""" + m = Model() + x = m.add_variables(lower=0, name="x") + y = m.add_variables(lower=0, name="y") + m.add_constraints(x + y <= 10, name="c1") + m.add_constraints(x <= 6, name="c2") + m.add_constraints(y <= 8, name="c3") + m.add_objective(3 * x + 5 * y, sense="max") + + primal_obj = _solve(m) + dual_obj = _solve(m.dualize()) + assert abs(primal_obj - dual_obj) < 1e-5 + + From 9cd954339a09c7dba870a8189fdf9b22f31a9e7e Mon Sep 17 00:00:00 2001 From: bobbyxng Date: Tue, 2 Jun 2026 10:24:19 +0200 Subject: [PATCH 21/30] rename to _lift_bounds_to_constraints and make private --- linopy/dual.py | 16 ++++++++++------ linopy/model.py | 4 +--- test/test_dual.py | 29 ++++++++++------------------- 3 files changed, 21 insertions(+), 28 deletions(-) diff --git a/linopy/dual.py b/linopy/dual.py index fadcc671..420d62fd 100644 --- a/linopy/dual.py +++ b/linopy/dual.py @@ -34,7 +34,7 @@ def _skip( return False -def bounds_to_constraints(m: Model) -> None: +def _lift_bounds_to_constraints(m: Model) -> None: """ Convert finite variable bounds to explicit ``>=`` / ``<=`` constraints. @@ -217,7 +217,7 @@ def _extract_dual_feas_entries( """ Return ``(flat_var_label, flat_dual_label, coeff)`` for each active NNZ entry of A. - Converts A to COO, maps rows/columns to flat labels via *clabels* / *vlabels*, + Converts sparse matrix A to coordinate format, maps rows/columns to flat labels via *clabels* / *vlabels*, and drops entries that are masked or belong to a constraint without a dual variable. Parameters @@ -411,10 +411,14 @@ def _add_dual_feasibility_constraints( flat_con_to_dual = _build_con_to_dual_label(m, dual_vars) if not len(flat_con_to_dual): - logger.warning("No valid constraint labels found, skipping dual feasibility constraints.") + logger.warning( + "No valid constraint labels found, skipping dual feasibility constraints." + ) return - flat_v, flat_d, nnz_data = _extract_dual_feas_entries(A, vlabels, clabels, flat_con_to_dual) + flat_v, flat_d, nnz_data = _extract_dual_feas_entries( + A, vlabels, clabels, flat_con_to_dual + ) c_lookup = _build_obj_coeff_lookup(vlabels, m.matrices.c) logger.debug("Building dual feasibility constraints for each primal variable.") @@ -520,7 +524,7 @@ def dualize( μ >= 0 for <= constraints, ν <= 0 for >= constraints. Variable bounds are converted to explicit constraints before dualization - via bounds_to_constraints(), so that they appear in the constraint matrix + via _lift_bounds_to_constraints(), so that they appear in the constraint matrix A and are correctly reflected in the dual. The dual variables in m2 are named identically to their corresponding @@ -565,7 +569,7 @@ def dualize( ) return m2 - bounds_to_constraints(m1) + _lift_bounds_to_constraints(m1) dual_vars = _add_dual_variables(m1, m2) _add_dual_feasibility_constraints(m1, m2, dual_vars) _add_dual_objective( diff --git a/linopy/model.py b/linopy/model.py index f3483eb7..ee78a114 100644 --- a/linopy/model.py +++ b/linopy/model.py @@ -55,7 +55,7 @@ Constraints, CSRConstraint, ) -from linopy.dual import bounds_to_constraints, dualize +from linopy.dual import dualize from linopy.expressions import ( LinearExpression, QuadraticExpression, @@ -2194,6 +2194,4 @@ def reset_solution(self) -> None: to_block_files = to_block_files - bounds_to_constraints = bounds_to_constraints - dualize = dualize diff --git a/test/test_dual.py b/test/test_dual.py index 5a84fff5..3dbf04e5 100644 --- a/test/test_dual.py +++ b/test/test_dual.py @@ -8,7 +8,7 @@ import xarray as xr from linopy import Model -from linopy.dual import bounds_to_constraints, dualize +from linopy.dual import dualize from linopy.solvers import licensed_solvers _lp_solver = next((s for s in ("highs", "glpk", "scip") if s in licensed_solvers), None) @@ -96,7 +96,9 @@ def test_dual_sign_conventions_max(): def test_dual_feasibility_rhs_equals_objective_coefficients(): """The dual feasibility constraint RHS must equal the primal objective coefficient.""" m = Model() - x = m.add_variables(lower=-np.inf, upper=np.inf, coords=[pd.RangeIndex(4)], name="x") + x = m.add_variables( + lower=-np.inf, upper=np.inf, coords=[pd.RangeIndex(4)], name="x" + ) m.add_constraints(x == np.array([1.0, 2.0, 3.0, 4.0]), name="eq") c = np.array([10.0, 20.0, 30.0, 40.0]) m.add_objective(c * x) @@ -109,7 +111,9 @@ def test_dual_multi_constraint_per_variable(): """Each primal variable appearing in k constraints produces k terms in its dual feas constraint.""" m = Model() n = 3 - x = m.add_variables(lower=-np.inf, upper=np.inf, coords=[pd.RangeIndex(n)], name="x") + x = m.add_variables( + lower=-np.inf, upper=np.inf, coords=[pd.RangeIndex(n)], name="x" + ) # Two equality constraints, both using x m.add_constraints(x == 1.0, name="c1") m.add_constraints(2.0 * x == 2.0, name="c2") @@ -123,19 +127,6 @@ def test_dual_multi_constraint_per_variable(): assert n_terms == 2 -def test_bounds_to_constraints_removes_bounds(): - """bounds_to_constraints should add constraints and relax variable bounds.""" - m = Model() - m.add_variables(lower=1.0, upper=5.0, coords=[pd.RangeIndex(3)], name="x") - - bounds_to_constraints(m) - - assert "x-bound-lower" in m.constraints - assert "x-bound-upper" in m.constraints - assert np.all(np.isinf(m.variables["x"].lower.values)) - assert np.all(np.isinf(m.variables["x"].upper.values)) - - def test_dual_with_masked_variable(): """Partially masked variable: active elements produce a constraint, masked ones do not.""" m = Model() @@ -152,7 +143,9 @@ def test_dual_with_masked_variable(): def test_dual_free_unconstrained_variable_no_error(): """A free variable with no constraint connections is silently skipped.""" m = Model() - m.add_variables(lower=-np.inf, upper=np.inf, name="x") # no bounds → no bound constraints + m.add_variables( + lower=-np.inf, upper=np.inf, name="x" + ) # no bounds → no bound constraints y = m.add_variables(lower=-np.inf, upper=np.inf, name="y") m.add_constraints(y == 5, name="eq") m.add_objective(y) # x has no connections at all @@ -236,5 +229,3 @@ def test_strong_duality_maximization(): primal_obj = _solve(m) dual_obj = _solve(m.dualize()) assert abs(primal_obj - dual_obj) < 1e-5 - - From ecf661aa17dbeb5cc961cc75647fe462cb6d9634 Mon Sep 17 00:00:00 2001 From: bobbyxng Date: Tue, 2 Jun 2026 11:03:01 +0200 Subject: [PATCH 22/30] refactor, improve docstrings --- linopy/dual.py | 95 +++++++++++++++++++++++++---------------------- test/test_dual.py | 62 ++++++++++++++++--------------- 2 files changed, 83 insertions(+), 74 deletions(-) diff --git a/linopy/dual.py b/linopy/dual.py index 420d62fd..96b825ba 100644 --- a/linopy/dual.py +++ b/linopy/dual.py @@ -94,7 +94,7 @@ def _lift_bounds_to_constraints(m: Model) -> None: def _add_dual_variables(m: Model, m2: Model) -> dict: """ - Add one dual variable to m2 for each active constraint in m. + Add one dual variable to m2 for each included constraint in m. Dual variable bounds encode the sign convention for each constraint type and primal objective sense: @@ -167,12 +167,13 @@ def _add_dual_variables(m: Model, m2: Model) -> dict: return dual_vars -def _build_con_to_dual_label(m: Model, dual_vars: dict) -> np.ndarray: +def _build_flat_con_to_dual_label_lookup(m: Model, dual_vars: dict) -> np.ndarray: """ - Build a flat constraint-label → flat dual-variable-label lookup array. + Build a lookup from flat primal constraint labels to flat dual variable labels. - ``result[flat_con_label]`` gives the flat dual variable label in m2, - or -1 for constraints not in *dual_vars* or with masked labels. + The returned array maps each flat constraint label to the corresponding + flat dual variable label in ``m2``. Entries are -1 for constraints not in + ``dual_vars`` or with masked labels. Parameters ---------- @@ -215,30 +216,31 @@ def _extract_dual_feas_entries( flat_con_to_dual: np.ndarray, ) -> tuple[np.ndarray, np.ndarray, np.ndarray]: """ - Return ``(flat_var_label, flat_dual_label, coeff)`` for each active NNZ entry of A. + Return ``(flat_var_label, flat_dual_label, coeff)`` for each included nonzero entry of ``A``. - Converts sparse matrix A to coordinate format, maps rows/columns to flat labels via *clabels* / *vlabels*, - and drops entries that are masked or belong to a constraint without a dual variable. + Converts sparse matrix ``A`` to coordinate format, maps rows/columns to + flat labels via ``clabels`` / ``vlabels``, and drops entries that are masked + or belong to a constraint without a dual variable. Parameters ---------- A : scipy sparse matrix - Primal constraint matrix (any scipy sparse format). + Primal constraint matrix in any scipy sparse format. vlabels : np.ndarray of int64 - Flat variable label per column of A; -1 for masked columns. + Flat variable label per column of ``A``; -1 for masked columns. clabels : np.ndarray of int64 - Flat constraint label per row of A; -1 for masked rows. + Flat constraint label per row of ``A``; -1 for masked rows. flat_con_to_dual : np.ndarray of int64 - Lookup array as returned by _build_con_to_dual_label(). + Lookup array as returned by ``_build_flat_con_to_dual_label_lookup()``. Returns ------- flat_v : np.ndarray of int64 - Flat primal variable label for each retained entry. + Flat primal variable label for each retained nonzero entry. flat_d : np.ndarray of int64 Corresponding flat dual variable label. coeffs : np.ndarray of float64 - Corresponding A matrix coefficient. + Corresponding coefficient from ``A``. """ A_coo = A.tocoo() flat_v = vlabels[A_coo.col].astype(np.int64) @@ -246,8 +248,8 @@ def _extract_dual_feas_entries( coeffs = A_coo.data # Drop entries where either label is masked. - active = (flat_v != -1) & (flat_c != -1) - flat_v, flat_c, coeffs = flat_v[active], flat_c[active], coeffs[active] + has_labels = (flat_v != -1) & (flat_c != -1) + flat_v, flat_c, coeffs = flat_v[has_labels], flat_c[has_labels], coeffs[has_labels] # Map primal constraint labels to flat dual variable labels. n = len(flat_con_to_dual) @@ -262,17 +264,18 @@ def _extract_dual_feas_entries( def _build_obj_coeff_lookup(vlabels: np.ndarray, c_vec: np.ndarray) -> np.ndarray: """ - Build a flat variable-label → objective-coefficient lookup array. + Build a lookup from flat variable labels to objective coefficients. - ``result[flat_var_label]`` gives the objective coefficient; unlabelled - positions are 0.0. + The returned array maps each valid flat variable label to its objective + coefficient. Labels that do not occur in ``vlabels`` but fall within the + lookup range have value 0.0. Parameters ---------- vlabels : np.ndarray of int64 - Flat variable label per column of A; -1 for masked entries. + Flat variable label per column of ``A``; -1 for masked columns. c_vec : np.ndarray of float64 - Objective coefficient in the same column ordering. + Objective coefficient per column, in the same ordering as ``vlabels``. Returns ------- @@ -290,67 +293,71 @@ def _build_obj_coeff_lookup(vlabels: np.ndarray, c_vec: np.ndarray) -> np.ndarra return lookup -def _make_feas_lhs( +def _build_dual_feas_lhs( var: Any, flat_v: np.ndarray, flat_d: np.ndarray, nnz_data: np.ndarray, - m2: Model, + m_dual: Model, ) -> LinearExpression: """ - Build ``sum_i(A_ji * lambda_i)`` as a LinearExpression over *var*'s coordinates. + Build the dual-feasibility LHS for one primal variable. + + For each coordinate of ``var``, constructs the linear expression + ``sum_j A[j, var] * lambda[j]`` over the corresponding dual variables. Parameters ---------- var : linopy.Variable - Primal variable whose stationarity LHS is being constructed. + Primal variable whose dual-feasibility LHS is being constructed. flat_v : np.ndarray of int64 - Flat primal variable labels for all active NNZ entries, as returned - by _extract_dual_feas_entries(). + Flat primal variable labels for all included nonzero entries, as returned + by ``_extract_dual_feas_entries()``. flat_d : np.ndarray of int64 Corresponding flat dual variable labels. nnz_data : np.ndarray of float64 - Corresponding A matrix coefficients. - m2 : Model + Corresponding nonzero coefficients from the primal constraint matrix. + m_dual : Model Dual model owning the dual variables. Returns ------- LinearExpression - Stationarity LHS over *var*'s coordinate space. Returns an empty - expression (constant zero) when the variable has no constraint connections. + Dual-feasibility LHS over ``var``'s coordinate space. Returns an empty + expression (constant zero) when ``var`` has no included constraint-matrix + entries. """ var_flat = var.labels.values.ravel().astype(np.int64) n_elements = len(var_flat) valid_mask = var_flat != -1 if not valid_mask.any(): - return LinearExpression(None, m2) + return LinearExpression(None, m_dual) - # Map flat variable labels to ravelled positions (unique per variable type). + # Map flat variable labels to ravelled positions in this variable. max_fv = int(var_flat[valid_mask].max()) var_to_idx = np.full(max_fv + 1, -1, dtype=np.int64) var_to_idx[var_flat[valid_mask]] = np.where(valid_mask)[0].astype(np.int64) - # Locate NNZ entries that belong to this variable via bounded lookup. + # Locate entries that belong to this variable via bounded lookup. lin_idx = np.full(len(flat_v), -1, dtype=np.int64) - in_range = flat_v <= max_fv + in_range = (flat_v >= 0) & (flat_v <= max_fv) if in_range.any(): lin_idx[in_range] = var_to_idx[flat_v[in_range]] keep = lin_idx != -1 - lin_idx_f, flat_d_f, nnz_f = lin_idx[keep], flat_d[keep], nnz_data[keep] + lin_idx_f, flat_d_f, coeffs_f = lin_idx[keep], flat_d[keep], nnz_data[keep] if not len(lin_idx_f): - return LinearExpression(None, m2) + return LinearExpression(None, m_dual) # Sort by ravelled position so entries for the same element are contiguous. order = np.argsort(lin_idx_f, kind="stable") lin_idx_s = lin_idx_f[order] flat_d_s = flat_d_f[order] - nnz_s = nnz_f[order] + coeffs_s = coeffs_f[order] - # Compute within-group term-slot positions (fully vectorised, no Python loop). + # Compute within-group term-slot positions. # group_start[i] is True at the first entry of each new element group. group_start = np.empty(len(lin_idx_s), dtype=bool) group_start[0] = True @@ -361,11 +368,11 @@ def _make_feas_lhs( max_terms = int(col_idx.max()) + 1 - # Populate (n_elements, max_terms) arrays via advanced indexing (C-level). + # Populate (n_elements, max_terms) arrays via advanced indexing. dual_labels_2d = np.full((n_elements, max_terms), -1, dtype=np.int64) dual_coeffs_2d = np.zeros((n_elements, max_terms), dtype=np.float64) dual_labels_2d[lin_idx_s, col_idx] = flat_d_s - dual_coeffs_2d[lin_idx_s, col_idx] = nnz_s + dual_coeffs_2d[lin_idx_s, col_idx] = coeffs_s # Wrap in a LinearExpression, reshaping to (*var_dims, _term). target_shape = var.labels.shape + (max_terms,) @@ -377,7 +384,7 @@ def _make_feas_lhs( }, coords={dim: var.labels.coords[dim] for dim in var.labels.dims}, ) - return LinearExpression(ds, m2) + return LinearExpression(ds, m_dual) def _add_dual_feasibility_constraints( @@ -409,7 +416,7 @@ def _add_dual_feasibility_constraints( vlabels = np.asarray(m.matrices.vlabels, dtype=np.int64) clabels = np.asarray(m.matrices.clabels, dtype=np.int64) - flat_con_to_dual = _build_con_to_dual_label(m, dual_vars) + flat_con_to_dual = _build_flat_con_to_dual_label_lookup(m, dual_vars) if not len(flat_con_to_dual): logger.warning( "No valid constraint labels found, skipping dual feasibility constraints." @@ -437,7 +444,7 @@ def _add_dual_feasibility_constraints( coords=var.labels.coords, ) - lhs = _make_feas_lhs(var, flat_v, flat_d, nnz_data, m2) + lhs = _build_dual_feas_lhs(var, flat_v, flat_d, nnz_data, m2) if lhs.is_constant: # Variable has no constraint connections (free, unconstrained variable). # The stationarity condition 0 = c_j cannot be expressed as a linopy diff --git a/test/test_dual.py b/test/test_dual.py index 3dbf04e5..1b093202 100644 --- a/test/test_dual.py +++ b/test/test_dual.py @@ -18,9 +18,9 @@ # Structural tests (no solver required) def test_dualize_empty_model(): m = Model() - m2 = dualize(m) - assert len(m2.variables) == 0 - assert len(m2.constraints) == 0 + m_dual = dualize(m) + assert len(m_dual.variables) == 0 + assert len(m_dual.constraints) == 0 def test_dual_variables_named_after_primal_constraints(): @@ -29,9 +29,9 @@ def test_dual_variables_named_after_primal_constraints(): m.add_constraints(x >= 1.0, name="lb_con") m.add_objective(2.0 * x) - m2 = dualize(m) - assert "lb_con" in m2.variables - assert "x-bound-lower" in m2.variables + m_dual = dualize(m) + assert "lb_con" in m_dual.variables + assert "x-bound-lower" in m_dual.variables def test_dual_feasibility_constraints_named_after_primal_variables(): @@ -40,8 +40,8 @@ def test_dual_feasibility_constraints_named_after_primal_variables(): m.add_constraints(x >= 1.0, name="lb_con") m.add_objective(2.0 * x) - m2 = dualize(m) - assert "x" in m2.constraints + m_dual = dualize(m) + assert "x" in m_dual.constraints def test_dual_objective_sense_flipped_min(): @@ -50,8 +50,8 @@ def test_dual_objective_sense_flipped_min(): m.add_constraints(x >= 1, name="c") m.add_objective(x) # min - m2 = dualize(m) - assert m2.objective.sense == "max" + m_dual = dualize(m) + assert m_dual.objective.sense == "max" def test_dual_objective_sense_flipped_max(): @@ -60,8 +60,8 @@ def test_dual_objective_sense_flipped_max(): m.add_constraints(x <= 5, name="c") m.add_objective(x, sense="max") - m2 = dualize(m) - assert m2.objective.sense == "min" + m_dual = dualize(m) + assert m_dual.objective.sense == "min" def test_dual_sign_conventions_min(): @@ -73,11 +73,11 @@ def test_dual_sign_conventions_min(): m.add_constraints(x >= -10, name="geq") m.add_objective(x) - m2 = dualize(m) - assert m2.variables["eq"].lower.item() == -np.inf - assert m2.variables["eq"].upper.item() == np.inf - assert m2.variables["leq"].upper.item() == 0 - assert m2.variables["geq"].lower.item() == 0 + m_dual = dualize(m) + assert m_dual.variables["eq"].lower.item() == -np.inf + assert m_dual.variables["eq"].upper.item() == np.inf + assert m_dual.variables["leq"].upper.item() == 0 + assert m_dual.variables["geq"].lower.item() == 0 def test_dual_sign_conventions_max(): @@ -88,9 +88,9 @@ def test_dual_sign_conventions_max(): m.add_constraints(x >= -10, name="geq") m.add_objective(x, sense="max") - m2 = dualize(m) - assert m2.variables["leq"].lower.item() == 0 - assert m2.variables["geq"].upper.item() == 0 + m_dual = dualize(m) + assert m_dual.variables["leq"].lower.item() == 0 + assert m_dual.variables["geq"].upper.item() == 0 def test_dual_feasibility_rhs_equals_objective_coefficients(): @@ -103,8 +103,8 @@ def test_dual_feasibility_rhs_equals_objective_coefficients(): c = np.array([10.0, 20.0, 30.0, 40.0]) m.add_objective(c * x) - m2 = dualize(m) - np.testing.assert_allclose(m2.constraints["x"].rhs.values, c) + m_dual = dualize(m) + np.testing.assert_allclose(m_dual.constraints["x"].rhs.values, c) def test_dual_multi_constraint_per_variable(): @@ -119,9 +119,9 @@ def test_dual_multi_constraint_per_variable(): m.add_constraints(2.0 * x == 2.0, name="c2") m.add_objective(5.0 * x) - m2 = dualize(m) + m_dual = dualize(m) # dual feas constraint for x: lambda_c1 + 2*lambda_c2 = 5 - con_x = m2.constraints["x"] + con_x = m_dual.constraints["x"] # Should have 2 terms (one from c1, one from c2) n_terms = (con_x.vars != -1).sum(dim="_term").max().item() assert n_terms == 2 @@ -135,9 +135,11 @@ def test_dual_with_masked_variable(): m.add_constraints(x >= 1.0, name="c", mask=mask) m.add_objective(x) - m2 = dualize(m) - assert "x" in m2.constraints - assert (m2.constraints["x"].labels != -1).sum().item() == 2 # only 2 active rows + m_dual = dualize(m) + assert "x" in m_dual.constraints + assert ( + m_dual.constraints["x"].labels != -1 + ).sum().item() == 2 # only 2 active rows def test_dual_free_unconstrained_variable_no_error(): @@ -150,9 +152,9 @@ def test_dual_free_unconstrained_variable_no_error(): m.add_constraints(y == 5, name="eq") m.add_objective(y) # x has no connections at all - m2 = dualize(m) # must not raise - assert "y" in m2.constraints - assert "x" not in m2.constraints # no constraint connections → silently skipped + m_dual = dualize(m) # must not raise + assert "y" in m_dual.constraints + assert "x" not in m_dual.constraints # no constraint connections → silently skipped # Numerical tests (require solver) From 9babff9835e42cb20dae13ccf962243f0fb8cd79 Mon Sep 17 00:00:00 2001 From: bobbyxng Date: Tue, 2 Jun 2026 11:58:05 +0200 Subject: [PATCH 23/30] improve readability, docstring formats, testing suite --- linopy/dual.py | 253 ++++++++++++++++++++++++++++++---------------- test/test_dual.py | 136 ++++++++++++++++++++++--- 2 files changed, 290 insertions(+), 99 deletions(-) diff --git a/linopy/dual.py b/linopy/dual.py index 96b825ba..aa6749d1 100644 --- a/linopy/dual.py +++ b/linopy/dual.py @@ -65,9 +65,9 @@ def _lift_bounds_to_constraints(m: Model) -> None: mask=mask, ) logger.debug(f"Added lower bound constraint for '{var_name}'.") + # Remove bounds to avoid double-counting in the dual. + # Rely on the new constraints instead. var.lower.values[mask.values] = -np.inf - # Remove bounds to avoid double-counting in the dual. Rely on the new constraints instead. - m.variables[var_name].lower.values[mask.values] = -np.inf else: logger.debug( f"Variable '{var_name}' has no finite lower bound, skipping." @@ -83,18 +83,18 @@ def _lift_bounds_to_constraints(m: Model) -> None: mask=mask, ) logger.debug(f"Added upper bound constraint for '{var_name}'.") + # Remove bounds to avoid double-counting in the dual. + # Rely on the new constraints instead. var.upper.values[mask.values] = np.inf - # Remove bounds to avoid double-counting in the dual. Rely on the new constraints instead. - m.variables[var_name].upper.values[mask.values] = np.inf else: logger.debug( f"Variable '{var_name}' has no finite upper bound, skipping." ) -def _add_dual_variables(m: Model, m2: Model) -> dict: +def _add_dual_variables(m: Model, m_dual: Model) -> dict: """ - Add one dual variable to m2 for each included constraint in m. + Add one dual variable to m_dual for each included constraint in m. Dual variable bounds encode the sign convention for each constraint type and primal objective sense: @@ -115,7 +115,7 @@ def _add_dual_variables(m: Model, m2: Model) -> dict: ---------- m : Model Primal model. - m2 : Model + m_dual : Model Dual model to populate. Returns @@ -150,14 +150,15 @@ def _add_dual_variables(m: Model, m2: Model) -> dict: continue logger.debug( - f"Adding {var_type} dual variable for constraint '{name}' with shape {con.shape} and dims {con.labels.dims}." + f"Adding {var_type} dual variable for constraint '{name}' " + f"with shape {con.shape} and dims {con.labels.dims}." ) coords = ( [con.labels.coords[dim] for dim in con.labels.dims] if con.labels.dims else None ) - dual_vars[name] = m2.add_variables( + dual_vars[name] = m_dual.add_variables( lower=lower, upper=upper, coords=coords, @@ -172,15 +173,15 @@ def _build_flat_con_to_dual_label_lookup(m: Model, dual_vars: dict) -> np.ndarra Build a lookup from flat primal constraint labels to flat dual variable labels. The returned array maps each flat constraint label to the corresponding - flat dual variable label in ``m2``. Entries are -1 for constraints not in - ``dual_vars`` or with masked labels. + flat dual variable label from ``dual_vars``. Entries are -1 for constraints + not in ``dual_vars`` or with masked labels. Parameters ---------- m : Model Primal model. dual_vars : dict - ``{constraint_name: dual_variable}`` as returned by _add_dual_variables(). + ``{constraint_name: dual_variable}`` as returned by ``_add_dual_variables()``. Returns ------- @@ -216,7 +217,8 @@ def _extract_dual_feas_entries( flat_con_to_dual: np.ndarray, ) -> tuple[np.ndarray, np.ndarray, np.ndarray]: """ - Return ``(flat_var_label, flat_dual_label, coeff)`` for each included nonzero entry of ``A``. + Return ``(flat_var_label, flat_dual_label, coeff)`` for each included + nonzero entry of ``A``. Converts sparse matrix ``A`` to coordinate format, maps rows/columns to flat labels via ``clabels`` / ``vlabels``, and drops entries that are masked @@ -253,7 +255,7 @@ def _extract_dual_feas_entries( # Map primal constraint labels to flat dual variable labels. n = len(flat_con_to_dual) - in_range = flat_c < n + in_range = (flat_c >= 0) & (flat_c < n) flat_d = np.full(len(flat_c), -1, dtype=np.int64) flat_d[in_range] = flat_con_to_dual[flat_c[in_range]] @@ -293,6 +295,68 @@ def _build_obj_coeff_lookup(vlabels: np.ndarray, c_vec: np.ndarray) -> np.ndarra return lookup +def _build_label_to_flat_index_lookup(labels: np.ndarray) -> np.ndarray: + """ + Build a lookup from flat labels to flat indices. + + ``labels`` is expected to contain nonnegative integer labels, with -1 used + as the masked-label sentinel. + + ``lookup[label]`` gives the flat index of ``label`` in ``labels``. + Labels that do not occur in ``labels`` map to -1. If no valid labels exist, + returns an empty int64 array. + """ + valid = labels != -1 + if not valid.any(): + return np.array([], dtype=np.int64) + + lookup = np.full(int(labels[valid].max()) + 1, -1, dtype=np.int64) + lookup[labels[valid].astype(np.int64)] = np.where(valid)[0].astype(np.int64) + return lookup + + +def _lookup_flat_indices(labels: np.ndarray, lookup: np.ndarray) -> np.ndarray: + """ + Look up flat indices for flat labels. + + ``labels`` is expected to contain nonnegative integer labels, with -1 used + as the masked-label sentinel. ``lookup`` is expected to be an array as + returned by ``_build_label_to_flat_index_lookup()``. + + Labels outside the lookup range, including masked labels, map to -1. + """ + flat_indices = np.full(len(labels), -1, dtype=np.int64) + + in_range = (labels >= 0) & (labels < len(lookup)) + if in_range.any(): + flat_indices[in_range] = lookup[labels[in_range].astype(np.int64)] + + return flat_indices + + +def _term_slots_for_sorted_flat_indices(sorted_flat_indices: np.ndarray) -> np.ndarray: + """ + Return the term slot within each run of equal sorted flat indices. + + ``sorted_flat_indices`` must be non-empty and sorted in nondecreasing order. + Each run of equal values is treated as one group. + + Example: + ------- + ``[2, 2, 2, 5, 5, 9]`` becomes ``[0, 1, 2, 0, 1, 0]``. + """ + group_start = np.empty(len(sorted_flat_indices), dtype=bool) + group_start[0] = True + group_start[1:] = sorted_flat_indices[1:] != sorted_flat_indices[:-1] + + group_ids = np.cumsum(group_start) - 1 + group_start_pos = np.where(group_start)[0] + + return ( + np.arange(len(sorted_flat_indices), dtype=np.int64) - group_start_pos[group_ids] + ) + + def _build_dual_feas_lhs( var: Any, flat_v: np.ndarray, @@ -324,55 +388,45 @@ def _build_dual_feas_lhs( ------- LinearExpression Dual-feasibility LHS over ``var``'s coordinate space. Returns an empty - expression (constant zero) when ``var`` has no included constraint-matrix - entries. + expression, equivalent to constant zero, when ``var`` has no included + constraint-matrix entries. """ var_flat = var.labels.values.ravel().astype(np.int64) n_elements = len(var_flat) - valid_mask = var_flat != -1 - if not valid_mask.any(): + # Map flat variable labels to flat indices in this variable. + var_to_idx = _build_label_to_flat_index_lookup(var_flat) + if not len(var_to_idx): return LinearExpression(None, m_dual) - # Map flat variable labels to ravelled positions in this variable. - max_fv = int(var_flat[valid_mask].max()) - var_to_idx = np.full(max_fv + 1, -1, dtype=np.int64) - var_to_idx[var_flat[valid_mask]] = np.where(valid_mask)[0].astype(np.int64) - # Locate entries that belong to this variable via bounded lookup. - lin_idx = np.full(len(flat_v), -1, dtype=np.int64) - in_range = (flat_v >= 0) & (flat_v <= max_fv) - if in_range.any(): - lin_idx[in_range] = var_to_idx[flat_v[in_range]] + lin_idx = _lookup_flat_indices(flat_v, var_to_idx) + # Keep only entries that belong to this variable. keep = lin_idx != -1 - lin_idx_f, flat_d_f, coeffs_f = lin_idx[keep], flat_d[keep], nnz_data[keep] - - if not len(lin_idx_f): + if not keep.any(): return LinearExpression(None, m_dual) - # Sort by ravelled position so entries for the same element are contiguous. - order = np.argsort(lin_idx_f, kind="stable") - lin_idx_s = lin_idx_f[order] - flat_d_s = flat_d_f[order] - coeffs_s = coeffs_f[order] + lin_idx = lin_idx[keep] + dual_labels = flat_d[keep] + coeffs = nnz_data[keep] - # Compute within-group term-slot positions. - # group_start[i] is True at the first entry of each new element group. - group_start = np.empty(len(lin_idx_s), dtype=bool) - group_start[0] = True - group_start[1:] = lin_idx_s[1:] != lin_idx_s[:-1] - group_ids = np.cumsum(group_start) - 1 - group_start_pos = np.where(group_start)[0] - col_idx = np.arange(len(lin_idx_s), dtype=np.int64) - group_start_pos[group_ids] + # Sort by flat index so entries for the same element are contiguous. + order = np.argsort(lin_idx, kind="stable") + lin_idx = lin_idx[order] + dual_labels = dual_labels[order] + coeffs = coeffs[order] + # Assign each term to a column slot within its variable-element group. + # Each run of equal ``lin_idx`` values corresponds to one variable element. + col_idx = _term_slots_for_sorted_flat_indices(lin_idx) max_terms = int(col_idx.max()) + 1 # Populate (n_elements, max_terms) arrays via advanced indexing. dual_labels_2d = np.full((n_elements, max_terms), -1, dtype=np.int64) dual_coeffs_2d = np.zeros((n_elements, max_terms), dtype=np.float64) - dual_labels_2d[lin_idx_s, col_idx] = flat_d_s - dual_coeffs_2d[lin_idx_s, col_idx] = coeffs_s + dual_labels_2d[lin_idx, col_idx] = dual_labels + dual_coeffs_2d[lin_idx, col_idx] = coeffs # Wrap in a LinearExpression, reshaping to (*var_dims, _term). target_shape = var.labels.shape + (max_terms,) @@ -389,14 +443,15 @@ def _build_dual_feas_lhs( def _add_dual_feasibility_constraints( m: Model, - m2: Model, + m_dual: Model, dual_vars: dict, ) -> None: """ - Add the stationarity constraint ``sum_i(A_ji * lambda_i) = c_j`` for each primal variable. + Add the stationarity constraint ``sum_i(A_ji * lambda_i) = c_j`` for each + primal variable. Sign conventions are already encoded in the dual variable bounds produced by - _add_dual_variables(), so raw A coefficients are used without adjustment. + ``_add_dual_variables()``, so raw A coefficients are used without adjustment. Variables with no constraint connections (e.g. unconstrained free variables) are skipped; the dual is infeasible for such problems if ``c_j != 0``. @@ -404,10 +459,10 @@ def _add_dual_feasibility_constraints( ---------- m : Model Primal model. - m2 : Model + m_dual : Model Dual model to populate. dual_vars : dict - ``{constraint_name: dual_variable}`` as returned by _add_dual_variables(). + ``{constraint_name: dual_variable}`` as returned by ``_add_dual_variables()``. """ A = m.matrices.A if A is None: @@ -437,40 +492,55 @@ def _add_dual_feasibility_constraints( var_flat = var.labels.values.ravel().astype(np.int64) # RHS: objective coefficient for each element of this variable. - in_c_range = (var_flat != -1) & (var_flat < len(c_lookup)) - safe_flat = np.where(in_c_range, var_flat, 0) + c_arr = np.zeros(len(var_flat), dtype=np.float64) + in_c_range = (var_flat >= 0) & (var_flat < len(c_lookup)) + + if in_c_range.any(): + c_arr[in_c_range] = c_lookup[var_flat[in_c_range]] + c_vals = xr.DataArray( - np.where(in_c_range, c_lookup[safe_flat], 0.0).reshape(var.labels.shape), + c_arr.reshape(var.labels.shape), coords=var.labels.coords, + dims=var.labels.dims, ) - lhs = _build_dual_feas_lhs(var, flat_v, flat_d, nnz_data, m2) + lhs = _build_dual_feas_lhs(var, flat_v, flat_d, nnz_data, m_dual) + if lhs.is_constant: - # Variable has no constraint connections (free, unconstrained variable). - # The stationarity condition 0 = c_j cannot be expressed as a linopy - # constraint; the dual is infeasible if c_j != 0, trivial if c_j == 0. - logger.debug( - f"Variable '{var_name}' has no constraint connections; " - "skipping dual feasibility constraint." - ) + # If c_j is zero, the condition is redundant. If c_j is nonzero, the + # condition is impossible because it reduces to 0 == c_j. Since this + # constant-only condition has no variables to add to the dual model, we + # report it and skip adding the constraint. + if np.any(c_arr[mask.values.ravel()] != 0): + logger.warning( + f"Variable '{var_name}' has no constraint connections but has " + "nonzero objective coefficients; the corresponding dual-feasibility " + "condition is infeasible." + ) + else: + logger.debug( + f"Variable '{var_name}' has no constraint connections and zero " + "objective coefficients; skipping redundant dual feasibility constraint." + ) continue - m2.add_constraints(lhs == c_vals, name=var_name, mask=mask) + + m_dual.add_constraints(lhs == c_vals, name=var_name, mask=mask) def _add_dual_objective( m: Model, - m2: Model, + m_dual: Model, dual_vars: dict, add_objective_constant: float = 0.0, ) -> None: """ - Construct and add ``sum(rhs_i * lambda_i)`` as the dual objective of m2. + Construct and add ``sum(rhs_i * lambda_i)`` as the dual objective of m_dual. The uniform ``+`` sign is correct because sign conventions are already encoded in the dual variable bounds: a ``<=`` constraint has ``lambda <= 0``, so ``rhs * lambda`` contributes negatively without an explicit sign flip. This aligns with linopy's native dual convention, so - ``m2.variables[con_name].solution`` can be compared directly with + ``m_dual.variables[con_name].solution`` can be compared directly with ``m.constraints[con_name].dual`` after solving. The objective sense is flipped: ``min`` primal → ``max`` dual, and vice versa. @@ -479,15 +549,15 @@ def _add_dual_objective( ---------- m : Model Primal model. - m2 : Model + m_dual : Model Dual model to populate. dual_vars : dict - ``{constraint_name: dual_variable}`` as returned by _add_dual_variables(). + ``{constraint_name: dual_variable}`` as returned by ``_add_dual_variables()``. add_objective_constant : float, optional Constant added to the dual objective, e.g. to pass through a primal objective constant that was excluded from the model. """ - dual_obj: LinearExpression = LinearExpression(None, m2) + dual_obj: LinearExpression = LinearExpression(None, m_dual) sense = "max" if m.objective.sense == "min" else "min" for name, con in m.constraints.items(): @@ -504,7 +574,7 @@ def _add_dual_objective( logger.debug(f"Constructed dual objective with {len(dual_obj.coeffs)} terms.") logger.debug("Adding dual objective to model.") - m2.add_objective(dual_obj, sense=sense, overwrite=True) + m_dual.add_objective(dual_obj, sense=sense, overwrite=True) def dualize( @@ -514,7 +584,7 @@ def dualize( """ Construct the dual of a linopy LP model. - Transforms the primal model into its dual equivalent m2 following + Transforms the primal model into its dual equivalent m_dual following standard LP duality theory. The dual sense is flipped relative to the primal (min -> max, max -> min), and dual variable bounds depend on both constraint type and primal objective sense. @@ -534,19 +604,22 @@ def dualize( via _lift_bounds_to_constraints(), so that they appear in the constraint matrix A and are correctly reflected in the dual. - The dual variables in m2 are named identically to their corresponding - primal constraints and are accessible via m2.variables[con_name]. + The dual variables in m_dual are named identically to their corresponding + primal constraints and are accessible via m_dual.variables[con_name]. Strong duality guarantees that at optimality: primal objective = dual objective - Note: The standalone dual m2 may be unbounded if the primal is degenerate. - Only linear programs (LP) are supported. + Note: This constructs a standalone dual model. Pathological or unsupported + primal formulations may lead to infeasible or unbounded dual models or may + cause this function to raise an error. Only linear primal models with linear + objectives and constraints are supported. Parameters ---------- m : Model - Primal linopy model to dualize. Must have a linear objective and linear constraints. + Primal linopy model to dualize. Must have a linear objective and either + linear constraints or finite variable bounds. add_objective_constant : float, optional Constant added to the dual objective, e.g. to pass through a primal @@ -561,25 +634,31 @@ def dualize( -------- .. code-block:: python - m2 = m.dualize() - m2.solve(solver_name="gurobi", Method=2, Crossover=1) - gap = abs(m.objective.value - m2.objective.value) + m_dual = m.dualize() + m_dual.solve(solver_name="gurobi", Method=2, Crossover=1) + gap = abs(m.objective.value - m_dual.objective.value) """ from linopy.model import Model m1 = m.copy() - m2 = Model() + m_dual = Model() + + if not m1.variables: + logger.warning("Primal model has no variables. Returning empty dual model.") + return m_dual - if not m.variables or not m.constraints: + _lift_bounds_to_constraints(m1) + + if not m1.constraints: logger.warning( - "Primal model has no variables or constraints. Returning empty dual model." + "Primal model has no constraints after lifting variable bounds. " + "Returning empty dual model." ) - return m2 + return m_dual - _lift_bounds_to_constraints(m1) - dual_vars = _add_dual_variables(m1, m2) - _add_dual_feasibility_constraints(m1, m2, dual_vars) + dual_vars = _add_dual_variables(m1, m_dual) + _add_dual_feasibility_constraints(m1, m_dual, dual_vars) _add_dual_objective( - m1, m2, dual_vars, add_objective_constant=add_objective_constant + m1, m_dual, dual_vars, add_objective_constant=add_objective_constant ) - return m2 + return m_dual diff --git a/test/test_dual.py b/test/test_dual.py index 1b093202..07904b48 100644 --- a/test/test_dual.py +++ b/test/test_dual.py @@ -8,22 +8,89 @@ import xarray as xr from linopy import Model -from linopy.dual import dualize +from linopy.dual import ( + _build_label_to_flat_index_lookup, + _lookup_flat_indices, + _term_slots_for_sorted_flat_indices, + dualize, +) from linopy.solvers import licensed_solvers _lp_solver = next((s for s in ("highs", "glpk", "scip") if s in licensed_solvers), None) needs_solver = pytest.mark.skipif(_lp_solver is None, reason="No LP solver available") +# Structural tests for important internal functions +def test_build_label_to_flat_index_lookup(): + """Flat labels are mapped to their positions in the flattened label array.""" + labels = np.array([10, -1, 12, 99], dtype=np.int64) + + lookup = _build_label_to_flat_index_lookup(labels) + + assert lookup[10] == 0 + assert lookup[12] == 2 + assert lookup[99] == 3 + assert lookup[11] == -1 + + +def test_build_label_to_flat_index_lookup_all_masked(): + """An all-masked label array produces an empty lookup.""" + labels = np.array([-1, -1], dtype=np.int64) + + lookup = _build_label_to_flat_index_lookup(labels) + + assert lookup.dtype == np.int64 + assert len(lookup) == 0 + + +def test_lookup_flat_indices_bounds_safe(): + """Out-of-range and masked labels safely map to -1.""" + lookup = np.array([-1, -1, 5, -1, 7], dtype=np.int64) + labels = np.array([-1, 0, 2, 4, 99], dtype=np.int64) + + flat_indices = _lookup_flat_indices(labels, lookup) + + np.testing.assert_array_equal( + flat_indices, + np.array([-1, -1, 5, 7, -1], dtype=np.int64), + ) + + +def test_term_slots_for_sorted_flat_indices(): + """Repeated sorted flat indices are assigned increasing term slots.""" + sorted_flat_indices = np.array([2, 2, 2, 5, 5, 9], dtype=np.int64) + + slots = _term_slots_for_sorted_flat_indices(sorted_flat_indices) + + np.testing.assert_array_equal( + slots, + np.array([0, 1, 2, 0, 1, 0], dtype=np.int64), + ) + + # Structural tests (no solver required) def test_dualize_empty_model(): + """Dualizing an empty model returns an empty dual model.""" m = Model() m_dual = dualize(m) assert len(m_dual.variables) == 0 assert len(m_dual.constraints) == 0 +def test_variable_bounds_lifted_to_dual_variables(): + """Finite variable bounds are converted into dual variables.""" + m = Model() + x = m.add_variables(lower=1, upper=3, name="x") + m.add_objective(x) + + m_dual = dualize(m) + + assert "x-bound-lower" in m_dual.variables + assert "x-bound-upper" in m_dual.variables + + def test_dual_variables_named_after_primal_constraints(): + """Dual variables use the names of their corresponding primal constraints.""" m = Model() x = m.add_variables(lower=0, coords=[pd.RangeIndex(3)], name="x") m.add_constraints(x >= 1.0, name="lb_con") @@ -35,6 +102,7 @@ def test_dual_variables_named_after_primal_constraints(): def test_dual_feasibility_constraints_named_after_primal_variables(): + """Dual-feasibility constraints use the names of the primal variables.""" m = Model() x = m.add_variables(lower=0, coords=[pd.RangeIndex(3)], name="x") m.add_constraints(x >= 1.0, name="lb_con") @@ -45,6 +113,7 @@ def test_dual_feasibility_constraints_named_after_primal_variables(): def test_dual_objective_sense_flipped_min(): + """A minimization primal produces a maximization dual.""" m = Model() x = m.add_variables(lower=0, name="x") m.add_constraints(x >= 1, name="c") @@ -55,6 +124,7 @@ def test_dual_objective_sense_flipped_min(): def test_dual_objective_sense_flipped_max(): + """A maximization primal produces a minimization dual.""" m = Model() x = m.add_variables(upper=10, name="x") m.add_constraints(x <= 5, name="c") @@ -94,7 +164,7 @@ def test_dual_sign_conventions_max(): def test_dual_feasibility_rhs_equals_objective_coefficients(): - """The dual feasibility constraint RHS must equal the primal objective coefficient.""" + """The dual feasibility RHS equals the primal objective coefficient.""" m = Model() x = m.add_variables( lower=-np.inf, upper=np.inf, coords=[pd.RangeIndex(4)], name="x" @@ -108,7 +178,7 @@ def test_dual_feasibility_rhs_equals_objective_coefficients(): def test_dual_multi_constraint_per_variable(): - """Each primal variable appearing in k constraints produces k terms in its dual feas constraint.""" + """A variable appearing in k constraints gets k dual-feasibility terms.""" m = Model() n = 3 x = m.add_variables( @@ -128,7 +198,7 @@ def test_dual_multi_constraint_per_variable(): def test_dual_with_masked_variable(): - """Partially masked variable: active elements produce a constraint, masked ones do not.""" + """Partially masked variables only produce constraints for unmasked elements.""" m = Model() mask = xr.DataArray([True, False, True], dims=["dim_0"]) x = m.add_variables(lower=0, coords=[pd.RangeIndex(3)], name="x", mask=mask) @@ -139,11 +209,11 @@ def test_dual_with_masked_variable(): assert "x" in m_dual.constraints assert ( m_dual.constraints["x"].labels != -1 - ).sum().item() == 2 # only 2 active rows + ).sum().item() == 2 # only 2 unmasked elements -def test_dual_free_unconstrained_variable_no_error(): - """A free variable with no constraint connections is silently skipped.""" +def test_dual_free_unconstrained_zero_cost_variable_no_error(): + """A disconnected zero-cost free variable is skipped.""" m = Model() m.add_variables( lower=-np.inf, upper=np.inf, name="x" @@ -154,11 +224,49 @@ def test_dual_free_unconstrained_variable_no_error(): m_dual = dualize(m) # must not raise assert "y" in m_dual.constraints - assert "x" not in m_dual.constraints # no constraint connections → silently skipped + assert "x" not in m_dual.constraints # no constraint connections + + +def test_dual_free_unconstrained_variable_with_objective_warns(caplog): + """A disconnected variable with nonzero objective coefficient is reported.""" + m = Model() + x = m.add_variables(lower=-np.inf, upper=np.inf, name="x") + y = m.add_variables(lower=-np.inf, upper=np.inf, name="y") + + m.add_constraints(y == 5, name="eq") + m.add_objective(x + y) + + with caplog.at_level("WARNING"): + m_dual = dualize(m) + + assert "y" in m_dual.constraints + assert "x" not in m_dual.constraints + assert any( + "corresponding dual-feasibility condition is infeasible" in rec.message + for rec in caplog.records + ) + + +def test_dual_multi_constraint_per_variable_coefficients(): + """Dual-feasibility terms keep the correct A-matrix coefficients.""" + m = Model() + x = m.add_variables(lower=-np.inf, upper=np.inf, name="x") + + m.add_constraints(x == 1.0, name="c1") + m.add_constraints(2.0 * x == 2.0, name="c2") + m.add_objective(5.0 * x) + + m_dual = dualize(m) + coeffs = m_dual.constraints["x"].coeffs.values.ravel() + coeffs = coeffs[m_dual.constraints["x"].vars.values.ravel() != -1] + + np.testing.assert_allclose(np.sort(coeffs), np.array([1.0, 2.0])) # Numerical tests (require solver) def _solve(model, **kwargs): + """Solve a model with the available LP solver and return its objective value.""" + assert _lp_solver is not None model.solve(solver_name=_lp_solver, io_api="lp", **kwargs) return model.objective.value @@ -180,21 +288,25 @@ def test_strong_duality_simple(): @needs_solver def test_strong_duality_array_variable(): - """Strong duality with multi-dimensional variables.""" + """Strong duality with array variables.""" rng = np.random.default_rng(0) n = 6 - A = rng.random((4, n)) - 0.5 - b = rng.random(4) + 0.5 + + A = rng.random((4, n)) + 0.1 + b = rng.random(4) + 1.0 c_obj = rng.random(n) + 0.1 m = Model() x = m.add_variables(lower=0, coords=[pd.RangeIndex(n)], name="x") + for i in range(4): m.add_constraints( sum(float(A[i, j]) * x[j] for j in range(n)) <= float(b[i]), name=f"c{i}", ) - m.add_objective(sum(float(c_obj[j]) * x[j] for j in range(n))) + + # Bounded because x >= 0 and A x <= b with A > 0. + m.add_objective(sum(float(c_obj[j]) * x[j] for j in range(n)), sense="max") primal_obj = _solve(m) dual_obj = _solve(m.dualize()) From ae3cc7c9a1835645301c604c75e286a53a69b3be Mon Sep 17 00:00:00 2001 From: bobbyxng Date: Tue, 2 Jun 2026 12:15:25 +0200 Subject: [PATCH 24/30] fix typing in pytests --- test/test_dual.py | 53 ++++++++++++++++++++++++++--------------------- 1 file changed, 29 insertions(+), 24 deletions(-) diff --git a/test/test_dual.py b/test/test_dual.py index 07904b48..764f7a94 100644 --- a/test/test_dual.py +++ b/test/test_dual.py @@ -2,10 +2,13 @@ from __future__ import annotations +from typing import Any + import numpy as np import pandas as pd import pytest import xarray as xr +from _pytest.logging import LogCaptureFixture from linopy import Model from linopy.dual import ( @@ -21,7 +24,7 @@ # Structural tests for important internal functions -def test_build_label_to_flat_index_lookup(): +def test_build_label_to_flat_index_lookup() -> None: """Flat labels are mapped to their positions in the flattened label array.""" labels = np.array([10, -1, 12, 99], dtype=np.int64) @@ -33,7 +36,7 @@ def test_build_label_to_flat_index_lookup(): assert lookup[11] == -1 -def test_build_label_to_flat_index_lookup_all_masked(): +def test_build_label_to_flat_index_lookup_all_masked() -> None: """An all-masked label array produces an empty lookup.""" labels = np.array([-1, -1], dtype=np.int64) @@ -43,7 +46,7 @@ def test_build_label_to_flat_index_lookup_all_masked(): assert len(lookup) == 0 -def test_lookup_flat_indices_bounds_safe(): +def test_lookup_flat_indices_bounds_safe() -> None: """Out-of-range and masked labels safely map to -1.""" lookup = np.array([-1, -1, 5, -1, 7], dtype=np.int64) labels = np.array([-1, 0, 2, 4, 99], dtype=np.int64) @@ -56,7 +59,7 @@ def test_lookup_flat_indices_bounds_safe(): ) -def test_term_slots_for_sorted_flat_indices(): +def test_term_slots_for_sorted_flat_indices() -> None: """Repeated sorted flat indices are assigned increasing term slots.""" sorted_flat_indices = np.array([2, 2, 2, 5, 5, 9], dtype=np.int64) @@ -69,7 +72,7 @@ def test_term_slots_for_sorted_flat_indices(): # Structural tests (no solver required) -def test_dualize_empty_model(): +def test_dualize_empty_model() -> None: """Dualizing an empty model returns an empty dual model.""" m = Model() m_dual = dualize(m) @@ -77,7 +80,7 @@ def test_dualize_empty_model(): assert len(m_dual.constraints) == 0 -def test_variable_bounds_lifted_to_dual_variables(): +def test_variable_bounds_lifted_to_dual_variables() -> None: """Finite variable bounds are converted into dual variables.""" m = Model() x = m.add_variables(lower=1, upper=3, name="x") @@ -89,7 +92,7 @@ def test_variable_bounds_lifted_to_dual_variables(): assert "x-bound-upper" in m_dual.variables -def test_dual_variables_named_after_primal_constraints(): +def test_dual_variables_named_after_primal_constraints() -> None: """Dual variables use the names of their corresponding primal constraints.""" m = Model() x = m.add_variables(lower=0, coords=[pd.RangeIndex(3)], name="x") @@ -101,7 +104,7 @@ def test_dual_variables_named_after_primal_constraints(): assert "x-bound-lower" in m_dual.variables -def test_dual_feasibility_constraints_named_after_primal_variables(): +def test_dual_feasibility_constraints_named_after_primal_variables() -> None: """Dual-feasibility constraints use the names of the primal variables.""" m = Model() x = m.add_variables(lower=0, coords=[pd.RangeIndex(3)], name="x") @@ -112,7 +115,7 @@ def test_dual_feasibility_constraints_named_after_primal_variables(): assert "x" in m_dual.constraints -def test_dual_objective_sense_flipped_min(): +def test_dual_objective_sense_flipped_min() -> None: """A minimization primal produces a maximization dual.""" m = Model() x = m.add_variables(lower=0, name="x") @@ -123,7 +126,7 @@ def test_dual_objective_sense_flipped_min(): assert m_dual.objective.sense == "max" -def test_dual_objective_sense_flipped_max(): +def test_dual_objective_sense_flipped_max() -> None: """A maximization primal produces a minimization dual.""" m = Model() x = m.add_variables(upper=10, name="x") @@ -134,7 +137,7 @@ def test_dual_objective_sense_flipped_max(): assert m_dual.objective.sense == "min" -def test_dual_sign_conventions_min(): +def test_dual_sign_conventions_min() -> None: """For a min primal: <= -> dual <= 0, >= -> dual >= 0, = -> dual free.""" m = Model() x = m.add_variables(lower=-np.inf, upper=np.inf, name="x") @@ -150,7 +153,7 @@ def test_dual_sign_conventions_min(): assert m_dual.variables["geq"].lower.item() == 0 -def test_dual_sign_conventions_max(): +def test_dual_sign_conventions_max() -> None: """For a max primal: <= -> dual >= 0, >= -> dual <= 0.""" m = Model() x = m.add_variables(lower=-np.inf, upper=np.inf, name="x") @@ -163,7 +166,7 @@ def test_dual_sign_conventions_max(): assert m_dual.variables["geq"].upper.item() == 0 -def test_dual_feasibility_rhs_equals_objective_coefficients(): +def test_dual_feasibility_rhs_equals_objective_coefficients() -> None: """The dual feasibility RHS equals the primal objective coefficient.""" m = Model() x = m.add_variables( @@ -177,7 +180,7 @@ def test_dual_feasibility_rhs_equals_objective_coefficients(): np.testing.assert_allclose(m_dual.constraints["x"].rhs.values, c) -def test_dual_multi_constraint_per_variable(): +def test_dual_multi_constraint_per_variable() -> None: """A variable appearing in k constraints gets k dual-feasibility terms.""" m = Model() n = 3 @@ -197,7 +200,7 @@ def test_dual_multi_constraint_per_variable(): assert n_terms == 2 -def test_dual_with_masked_variable(): +def test_dual_with_masked_variable() -> None: """Partially masked variables only produce constraints for unmasked elements.""" m = Model() mask = xr.DataArray([True, False, True], dims=["dim_0"]) @@ -212,7 +215,7 @@ def test_dual_with_masked_variable(): ).sum().item() == 2 # only 2 unmasked elements -def test_dual_free_unconstrained_zero_cost_variable_no_error(): +def test_dual_free_unconstrained_zero_cost_variable_no_error() -> None: """A disconnected zero-cost free variable is skipped.""" m = Model() m.add_variables( @@ -227,7 +230,9 @@ def test_dual_free_unconstrained_zero_cost_variable_no_error(): assert "x" not in m_dual.constraints # no constraint connections -def test_dual_free_unconstrained_variable_with_objective_warns(caplog): +def test_dual_free_unconstrained_variable_with_objective_warns( + caplog: LogCaptureFixture, +) -> None: """A disconnected variable with nonzero objective coefficient is reported.""" m = Model() x = m.add_variables(lower=-np.inf, upper=np.inf, name="x") @@ -247,7 +252,7 @@ def test_dual_free_unconstrained_variable_with_objective_warns(caplog): ) -def test_dual_multi_constraint_per_variable_coefficients(): +def test_dual_multi_constraint_per_variable_coefficients() -> None: """Dual-feasibility terms keep the correct A-matrix coefficients.""" m = Model() x = m.add_variables(lower=-np.inf, upper=np.inf, name="x") @@ -264,15 +269,15 @@ def test_dual_multi_constraint_per_variable_coefficients(): # Numerical tests (require solver) -def _solve(model, **kwargs): +def _solve(model: Model, **kwargs: Any) -> float: """Solve a model with the available LP solver and return its objective value.""" assert _lp_solver is not None model.solve(solver_name=_lp_solver, io_api="lp", **kwargs) - return model.objective.value + return float(model.objective.value) @needs_solver -def test_strong_duality_simple(): +def test_strong_duality_simple() -> None: """Strong duality: primal obj == dual obj at optimality.""" m = Model() x = m.add_variables(lower=0, name="x") @@ -287,7 +292,7 @@ def test_strong_duality_simple(): @needs_solver -def test_strong_duality_array_variable(): +def test_strong_duality_array_variable() -> None: """Strong duality with array variables.""" rng = np.random.default_rng(0) n = 6 @@ -314,7 +319,7 @@ def test_strong_duality_array_variable(): @needs_solver -def test_strong_duality_mixed_constraint_types(): +def test_strong_duality_mixed_constraint_types() -> None: """Strong duality with =, <=, >= constraints.""" m = Model() x = m.add_variables(lower=-np.inf, upper=np.inf, name="x") @@ -330,7 +335,7 @@ def test_strong_duality_mixed_constraint_types(): @needs_solver -def test_strong_duality_maximization(): +def test_strong_duality_maximization() -> None: """Strong duality for a maximization primal.""" m = Model() x = m.add_variables(lower=0, name="x") From fe1e3893bc9abbeb955bf05be5af510318ce99e9 Mon Sep 17 00:00:00 2001 From: bobbyxng Date: Tue, 2 Jun 2026 12:19:25 +0200 Subject: [PATCH 25/30] fix typing in pytests --- test/test_dual.py | 14 ++++++++------ 1 file changed, 8 insertions(+), 6 deletions(-) diff --git a/test/test_dual.py b/test/test_dual.py index 764f7a94..7ea31862 100644 --- a/test/test_dual.py +++ b/test/test_dual.py @@ -273,7 +273,10 @@ def _solve(model: Model, **kwargs: Any) -> float: """Solve a model with the available LP solver and return its objective value.""" assert _lp_solver is not None model.solve(solver_name=_lp_solver, io_api="lp", **kwargs) - return float(model.objective.value) + + value = model.objective.value + assert value is not None + return float(value) @needs_solver @@ -305,13 +308,12 @@ def test_strong_duality_array_variable() -> None: x = m.add_variables(lower=0, coords=[pd.RangeIndex(n)], name="x") for i in range(4): - m.add_constraints( - sum(float(A[i, j]) * x[j] for j in range(n)) <= float(b[i]), - name=f"c{i}", - ) + lhs: Any = sum(float(A[i, j]) * x[j] for j in range(n)) + m.add_constraints(lhs <= float(b[i]), name=f"c{i}") # Bounded because x >= 0 and A x <= b with A > 0. - m.add_objective(sum(float(c_obj[j]) * x[j] for j in range(n)), sense="max") + obj: Any = sum(float(c_obj[j]) * x[j] for j in range(n)) + m.add_objective(obj, sense="max") primal_obj = _solve(m) dual_obj = _solve(m.dualize()) From bac70c31eb6a4f973c22a22439b2bc4fbdd1cd7c Mon Sep 17 00:00:00 2001 From: bobbyxng Date: Tue, 2 Jun 2026 13:18:09 +0200 Subject: [PATCH 26/30] improve test coverage and remove non-supported objective constant --- linopy/dual.py | 21 ++----------- test/test_dual.py | 76 +++++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 79 insertions(+), 18 deletions(-) diff --git a/linopy/dual.py b/linopy/dual.py index aa6749d1..0e09ebf9 100644 --- a/linopy/dual.py +++ b/linopy/dual.py @@ -531,7 +531,6 @@ def _add_dual_objective( m: Model, m_dual: Model, dual_vars: dict, - add_objective_constant: float = 0.0, ) -> None: """ Construct and add ``sum(rhs_i * lambda_i)`` as the dual objective of m_dual. @@ -553,9 +552,6 @@ def _add_dual_objective( Dual model to populate. dual_vars : dict ``{constraint_name: dual_variable}`` as returned by ``_add_dual_variables()``. - add_objective_constant : float, optional - Constant added to the dual objective, e.g. to pass through a primal - objective constant that was excluded from the model. """ dual_obj: LinearExpression = LinearExpression(None, m_dual) sense = "max" if m.objective.sense == "min" else "min" @@ -568,10 +564,6 @@ def _add_dual_objective( rhs_masked = con.rhs.where(mask, 0) dual_obj += (rhs_masked * dual_vars[name]).sum() - if add_objective_constant != 0.0: - dual_obj += add_objective_constant - logger.debug(f"Added constant {add_objective_constant} to dual objective.") - logger.debug(f"Constructed dual objective with {len(dual_obj.coeffs)} terms.") logger.debug("Adding dual objective to model.") m_dual.add_objective(dual_obj, sense=sense, overwrite=True) @@ -579,7 +571,6 @@ def _add_dual_objective( def dualize( m: Model, - add_objective_constant: float = 0.0, ) -> Model: """ Construct the dual of a linopy LP model. @@ -621,17 +612,13 @@ def dualize( Primal linopy model to dualize. Must have a linear objective and either linear constraints or finite variable bounds. - add_objective_constant : float, optional - Constant added to the dual objective, e.g. to pass through a primal - objective constant that was excluded from the model. - Returns ------- Model Dual model whose variables are named after the primal constraints. - Examples - -------- + Example + ------- .. code-block:: python m_dual = m.dualize() @@ -658,7 +645,5 @@ def dualize( dual_vars = _add_dual_variables(m1, m_dual) _add_dual_feasibility_constraints(m1, m_dual, dual_vars) - _add_dual_objective( - m1, m_dual, dual_vars, add_objective_constant=add_objective_constant - ) + _add_dual_objective(m1, m_dual, dual_vars) return m_dual diff --git a/test/test_dual.py b/test/test_dual.py index 7ea31862..0836c88e 100644 --- a/test/test_dual.py +++ b/test/test_dual.py @@ -13,7 +13,9 @@ from linopy import Model from linopy.dual import ( _build_label_to_flat_index_lookup, + _build_obj_coeff_lookup, _lookup_flat_indices, + _skip, _term_slots_for_sorted_flat_indices, dualize, ) @@ -24,6 +26,20 @@ # Structural tests for important internal functions +def test_skip_empty_label_array() -> None: + """Empty label arrays are skipped.""" + labels = xr.DataArray(np.array([], dtype=np.int64), dims=["dim_0"]) + + assert _skip(labels, "variable", "x") + + +def test_skip_fully_masked_label_array() -> None: + """Fully masked label arrays are skipped.""" + labels = xr.DataArray(np.array([-1, -1], dtype=np.int64), dims=["dim_0"]) + + assert _skip(labels, "constraint", "c") + + def test_build_label_to_flat_index_lookup() -> None: """Flat labels are mapped to their positions in the flattened label array.""" labels = np.array([10, -1, 12, 99], dtype=np.int64) @@ -71,6 +87,17 @@ def test_term_slots_for_sorted_flat_indices() -> None: ) +def test_build_obj_coeff_lookup_all_masked() -> None: + """All-masked variable labels produce an empty objective-coefficient lookup.""" + lookup = _build_obj_coeff_lookup( + np.array([-1, -1], dtype=np.int64), + np.array([1.0, 2.0], dtype=np.float64), + ) + + assert lookup.dtype == np.float64 + assert len(lookup) == 0 + + # Structural tests (no solver required) def test_dualize_empty_model() -> None: """Dualizing an empty model returns an empty dual model.""" @@ -80,6 +107,55 @@ def test_dualize_empty_model() -> None: assert len(m_dual.constraints) == 0 +def test_only_lower_bound_lifted_to_dual_variable() -> None: + """A finite lower bound is lifted even when the upper bound is infinite.""" + m = Model() + x = m.add_variables(lower=1, upper=np.inf, name="x") + m.add_objective(x) + + m_dual = dualize(m) + + assert "x-bound-lower" in m_dual.variables + assert "x-bound-upper" not in m_dual.variables + + +def test_only_upper_bound_lifted_to_dual_variable() -> None: + """A finite upper bound is lifted even when the lower bound is infinite.""" + m = Model() + x = m.add_variables(lower=-np.inf, upper=3, name="x") + m.add_objective(x) + + m_dual = dualize(m) + + assert "x-bound-lower" not in m_dual.variables + assert "x-bound-upper" in m_dual.variables + + +def test_unbounded_variable_bounds_do_not_create_dual_variables() -> None: + """Infinite variable bounds are not lifted into dual variables.""" + m = Model() + x = m.add_variables(lower=-np.inf, upper=np.inf, name="x") + m.add_constraints(x == 1, name="c") + m.add_objective(x) + + m_dual = dualize(m) + + assert "x-bound-lower" not in m_dual.variables + assert "x-bound-upper" not in m_dual.variables + + +def test_dualize_model_with_variables_but_no_constraints_or_finite_bounds() -> None: + """A model with variables but no constraints or finite bounds returns an empty dual.""" + m = Model() + x = m.add_variables(lower=-np.inf, upper=np.inf, name="x") + m.add_objective(x) + + m_dual = dualize(m) + + assert len(m_dual.variables) == 0 + assert len(m_dual.constraints) == 0 + + def test_variable_bounds_lifted_to_dual_variables() -> None: """Finite variable bounds are converted into dual variables.""" m = Model() From ea778b0ea96264964e52673c15f32dd2b660bf7a Mon Sep 17 00:00:00 2001 From: bobbyxng Date: Tue, 2 Jun 2026 13:51:49 +0200 Subject: [PATCH 27/30] Fix cases where constraints and variables of the same block/name can have different signs and/or bounds. Add tests --- linopy/dual.py | 141 ++++++++++++++++++++++++++--------- test/test_dual.py | 182 +++++++++++++++++++++++++++++++++++++++++++++- 2 files changed, 286 insertions(+), 37 deletions(-) diff --git a/linopy/dual.py b/linopy/dual.py index 0e09ebf9..6a9b7c17 100644 --- a/linopy/dual.py +++ b/linopy/dual.py @@ -43,55 +43,120 @@ def _lift_bounds_to_constraints(m: Model) -> None: The variable bounds are then relaxed to ``[-inf, inf]`` to avoid double-counting when the dual is formed. + Bounds may vary elementwise within one variable block. Infinite bounds are + not converted into constraints. + Parameters ---------- m : Model Model to mutate in-place. """ - logger.debug("Converting variable bounds to explicit constraints.") - logger.debug("Relaxing variable bounds to [-inf, inf].") + logger.debug("Converting finite variable bounds to explicit constraints.") + logger.debug("Relaxing converted variable bounds to [-inf, inf].") + for var_name, var in m.variables.items(): - mask = var.labels != -1 + label_mask = var.labels != -1 lb = var.lower ub = var.upper - # lower bound + finite_lb = xr.DataArray( + np.isfinite(lb.values), + coords=lb.coords, + dims=lb.dims, + ) + finite_ub = xr.DataArray( + np.isfinite(ub.values), + coords=ub.coords, + dims=ub.dims, + ) + + lower_mask = label_mask & finite_lb + upper_mask = label_mask & finite_ub + + # Lower bound. if f"{var_name}-bound-lower" not in m.constraints: - has_finite_lb = np.isfinite(lb.values[mask.values]).any() - if has_finite_lb: + if bool(lower_mask.any()): m.add_constraints( var >= lb, name=f"{var_name}-bound-lower", - mask=mask, + mask=lower_mask, ) logger.debug(f"Added lower bound constraint for '{var_name}'.") - # Remove bounds to avoid double-counting in the dual. + + # Remove converted bounds to avoid double-counting in the dual. # Rely on the new constraints instead. - var.lower.values[mask.values] = -np.inf + var.lower.values[lower_mask.values] = -np.inf else: logger.debug( f"Variable '{var_name}' has no finite lower bound, skipping." ) - # upper bound + # Upper bound. if f"{var_name}-bound-upper" not in m.constraints: - has_finite_ub = np.isfinite(ub.values[mask.values]).any() - if has_finite_ub: + if bool(upper_mask.any()): m.add_constraints( var <= ub, name=f"{var_name}-bound-upper", - mask=mask, + mask=upper_mask, ) logger.debug(f"Added upper bound constraint for '{var_name}'.") - # Remove bounds to avoid double-counting in the dual. + + # Remove converted bounds to avoid double-counting in the dual. # Rely on the new constraints instead. - var.upper.values[mask.values] = np.inf + var.upper.values[upper_mask.values] = np.inf else: logger.debug( f"Variable '{var_name}' has no finite upper bound, skipping." ) +def _dual_bounds_from_constraint_signs( + signs: xr.DataArray, + labels: xr.DataArray, + primal_is_min: bool, +) -> tuple[xr.DataArray, xr.DataArray, xr.DataArray]: + """ + Return elementwise dual-variable bounds for constraint signs. + + ``signs`` is broadcast to ``labels``. Valid signs are ``=``, ``<=``, and + ``>=``. ``valid_sign`` is True where the broadcast sign is valid. + """ + signs = signs.broadcast_like(labels) + + is_eq = signs == "=" + is_le = signs == "<=" + is_ge = signs == ">=" + valid_sign = is_eq | is_le | is_ge + + # Bounds for invalid signs are irrelevant because invalid entries are + # masked out by the caller. Use finite placeholders to avoid NaNs. + lower = xr.zeros_like(labels, dtype=float) + upper = xr.zeros_like(labels, dtype=float) + + # Equality constraints: dual variable is free. + lower = lower.where(~is_eq, -np.inf) + upper = upper.where(~is_eq, np.inf) + + if primal_is_min: + # <= constraints in a min primal: dual <= 0. + lower = lower.where(~is_le, -np.inf) + upper = upper.where(~is_le, 0.0) + + # >= constraints in a min primal: dual >= 0. + lower = lower.where(~is_ge, 0.0) + upper = upper.where(~is_ge, np.inf) + else: + # <= constraints in a max primal: dual >= 0. + lower = lower.where(~is_le, 0.0) + upper = upper.where(~is_le, np.inf) + + # >= constraints in a max primal: dual <= 0. + lower = lower.where(~is_ge, -np.inf) + upper = upper.where(~is_ge, 0.0) + + return lower, upper, valid_sign + + def _add_dual_variables(m: Model, m_dual: Model) -> dict: """ Add one dual variable to m_dual for each included constraint in m. @@ -109,7 +174,8 @@ def _add_dual_variables(m: Model, m_dual: Model) -> dict: >= max -inf 0 ============ =========== ================ ================ - Fully masked or empty constraints are skipped. + Fully masked or empty constraints are skipped. Constraint arrays may contain + mixed signs; dual variable bounds are assigned elementwise. Parameters ---------- @@ -130,27 +196,29 @@ def _add_dual_variables(m: Model, m_dual: Model) -> dict: if _skip(con.labels, "constraint", name): continue - mask = con.labels != -1 - sign = con.sign.isel({d: 0 for d in con.sign.dims}).item() - - match sign: - case "=": - lower, upper = -np.inf, np.inf - var_type = "free" - case "<=": - lower, upper = (-np.inf, 0) if primal_is_min else (0, np.inf) - var_type = "non-positive" if primal_is_min else "non-negative" - case ">=": - lower, upper = (0, np.inf) if primal_is_min else (-np.inf, 0) - var_type = "non-negative" if primal_is_min else "non-positive" - case _: - logger.warning( - f"Constraint '{name}' has unrecognized sign '{sign}', skipping." - ) - continue + lower, upper, valid_sign = _dual_bounds_from_constraint_signs( + con.sign, + con.labels, + primal_is_min, + ) + + unmasked = con.labels != -1 + invalid_unmasked = unmasked & ~valid_sign + if bool(invalid_unmasked.any()): + logger.warning( + f"Constraint '{name}' has unrecognized signs; invalid entries " + "will be skipped." + ) + + mask = unmasked & valid_sign + if not bool(mask.any()): + logger.warning( + f"Constraint '{name}' has no entries with valid signs, skipping." + ) + continue logger.debug( - f"Adding {var_type} dual variable for constraint '{name}' " + f"Adding dual variable for constraint '{name}' " f"with shape {con.shape} and dims {con.labels.dims}." ) coords = ( @@ -165,6 +233,7 @@ def _add_dual_variables(m: Model, m_dual: Model) -> dict: name=name, mask=mask, ) + return dual_vars @@ -604,7 +673,7 @@ def dualize( Note: This constructs a standalone dual model. Pathological or unsupported primal formulations may lead to infeasible or unbounded dual models or may cause this function to raise an error. Only linear primal models with linear - objectives and constraints are supported. + objectives and either linear constraints or finite variable bounds are supported. Parameters ---------- diff --git a/test/test_dual.py b/test/test_dual.py index 0836c88e..1aab2ede 100644 --- a/test/test_dual.py +++ b/test/test_dual.py @@ -14,6 +14,7 @@ from linopy.dual import ( _build_label_to_flat_index_lookup, _build_obj_coeff_lookup, + _dual_bounds_from_constraint_signs, _lookup_flat_indices, _skip, _term_slots_for_sorted_flat_indices, @@ -98,6 +99,40 @@ def test_build_obj_coeff_lookup_all_masked() -> None: assert len(lookup) == 0 +def test_dual_bounds_from_mixed_constraint_signs_min() -> None: + """Mixed signs produce elementwise dual bounds for a minimization primal.""" + idx = pd.RangeIndex(3, name="i") + labels = xr.DataArray([0, 1, 2], dims=["i"], coords={"i": idx}) + signs = xr.DataArray(["<=", ">=", "="], dims=["i"], coords={"i": idx}) + + lower, upper, valid_sign = _dual_bounds_from_constraint_signs( + signs, + labels, + primal_is_min=True, + ) + + np.testing.assert_allclose(lower.values, np.array([-np.inf, 0.0, -np.inf])) + np.testing.assert_allclose(upper.values, np.array([0.0, np.inf, np.inf])) + np.testing.assert_array_equal(valid_sign.values, np.array([True, True, True])) + + +def test_dual_bounds_from_mixed_constraint_signs_max() -> None: + """Mixed signs produce elementwise dual bounds for a maximization primal.""" + idx = pd.RangeIndex(3, name="i") + labels = xr.DataArray([0, 1, 2], dims=["i"], coords={"i": idx}) + signs = xr.DataArray(["<=", ">=", "="], dims=["i"], coords={"i": idx}) + + lower, upper, valid_sign = _dual_bounds_from_constraint_signs( + signs, + labels, + primal_is_min=False, + ) + + np.testing.assert_allclose(lower.values, np.array([0.0, -np.inf, -np.inf])) + np.testing.assert_allclose(upper.values, np.array([np.inf, 0.0, np.inf])) + np.testing.assert_array_equal(valid_sign.values, np.array([True, True, True])) + + # Structural tests (no solver required) def test_dualize_empty_model() -> None: """Dualizing an empty model returns an empty dual model.""" @@ -157,7 +192,7 @@ def test_dualize_model_with_variables_but_no_constraints_or_finite_bounds() -> N def test_variable_bounds_lifted_to_dual_variables() -> None: - """Finite variable bounds are converted into dual variables.""" + """Finite variable bounds are lifted and dualized.""" m = Model() x = m.add_variables(lower=1, upper=3, name="x") m.add_objective(x) @@ -344,6 +379,79 @@ def test_dual_multi_constraint_per_variable_coefficients() -> None: np.testing.assert_allclose(np.sort(coeffs), np.array([1.0, 2.0])) +def test_dual_mixed_sign_constraint_array_min() -> None: + """Mixed elementwise constraint signs produce elementwise dual bounds.""" + m = Model() + idx = pd.RangeIndex(3, name="i") + x = m.add_variables(lower=-np.inf, upper=np.inf, coords=[idx], name="x") + + signs = xr.DataArray(["<=", ">=", "="], dims=["i"], coords={"i": idx}) + rhs = xr.DataArray([1.0, 2.0, 3.0], dims=["i"], coords={"i": idx}) + + m.add_constraints(x, signs, rhs, name="mixed") + m.add_objective(x.sum()) + + m_dual = dualize(m) + dv = m_dual.variables["mixed"] + + np.testing.assert_allclose( + dv.lower.values, + np.array([-np.inf, 0.0, -np.inf]), + ) + np.testing.assert_allclose( + dv.upper.values, + np.array([0.0, np.inf, np.inf]), + ) + + +def test_dual_mixed_sign_constraint_array_max() -> None: + """Mixed elementwise constraint signs follow maximization dual bounds.""" + m = Model() + idx = pd.RangeIndex(3, name="i") + x = m.add_variables(lower=-np.inf, upper=np.inf, coords=[idx], name="x") + + signs = xr.DataArray(["<=", ">=", "="], dims=["i"], coords={"i": idx}) + rhs = xr.DataArray([1.0, 2.0, 3.0], dims=["i"], coords={"i": idx}) + + m.add_constraints(x, signs, rhs, name="mixed") + m.add_objective(x.sum(), sense="max") + + m_dual = dualize(m) + dv = m_dual.variables["mixed"] + + np.testing.assert_allclose( + dv.lower.values, + np.array([0.0, -np.inf, -np.inf]), + ) + np.testing.assert_allclose( + dv.upper.values, + np.array([np.inf, 0.0, np.inf]), + ) + + +def test_mixed_variable_bounds_lift_only_finite_entries() -> None: + """Mixed finite/infinite bounds are lifted only for finite entries.""" + m = Model() + idx = pd.RangeIndex(3, name="i") + + lower = xr.DataArray([0.0, -np.inf, -np.inf], dims=["i"], coords={"i": idx}) + upper = xr.DataArray([np.inf, 5.0, np.inf], dims=["i"], coords={"i": idx}) + + x = m.add_variables(lower=lower, upper=upper, coords=[idx], name="x") + m.add_objective(x.sum()) + + m_dual = dualize(m) + + lower_labels = m_dual.variables["x-bound-lower"].labels + upper_labels = m_dual.variables["x-bound-upper"].labels + + assert (lower_labels != -1).sum().item() == 1 + assert (upper_labels != -1).sum().item() == 1 + + assert lower_labels.sel(i=0).item() != -1 + assert upper_labels.sel(i=1).item() != -1 + + # Numerical tests (require solver) def _solve(model: Model, **kwargs: Any) -> float: """Solve a model with the available LP solver and return its objective value.""" @@ -426,3 +534,75 @@ def test_strong_duality_maximization() -> None: primal_obj = _solve(m) dual_obj = _solve(m.dualize()) assert abs(primal_obj - dual_obj) < 1e-5 + + +@needs_solver +def test_dualize_mixed_signs_and_mixed_variable_bounds() -> None: + """Dualization handles mixed constraint signs and mixed variable bounds.""" + m = Model() + idx = pd.RangeIndex(3, name="i") + + lower = xr.DataArray([0.0, -np.inf, -np.inf], dims=["i"], coords={"i": idx}) + upper = xr.DataArray([np.inf, 5.0, np.inf], dims=["i"], coords={"i": idx}) + + x = m.add_variables( + lower=lower, + upper=upper, + coords=[idx], + name="x", + ) + + signs = xr.DataArray(["<=", ">=", "="], dims=["i"], coords={"i": idx}) + rhs = xr.DataArray([4.0, 2.0, 3.0], dims=["i"], coords={"i": idx}) + + m.add_constraints(x, signs, rhs, name="mixed") + m.add_objective(x.sum()) + + primal_obj = _solve(m) + + m_dual = dualize(m) + dual_obj = _solve(m_dual) + + assert abs(primal_obj - dual_obj) < 1e-5 + + mixed_dual = m_dual.variables["mixed"] + + np.testing.assert_allclose( + mixed_dual.lower.values, + np.array([-np.inf, 0.0, -np.inf]), + ) + np.testing.assert_allclose( + mixed_dual.upper.values, + np.array([0.0, np.inf, np.inf]), + ) + + assert "x-bound-lower" in m_dual.variables + assert "x-bound-upper" in m_dual.variables + + lower_bound_labels = m_dual.variables["x-bound-lower"].labels + upper_bound_labels = m_dual.variables["x-bound-upper"].labels + + assert (lower_bound_labels != -1).sum().item() == 1 + assert (upper_bound_labels != -1).sum().item() == 1 + + assert lower_bound_labels.sel(i=0).item() != -1 + assert upper_bound_labels.sel(i=1).item() != -1 + + con_x = m_dual.constraints["x"] + + # x[0]: mixed[0] + x-bound-lower[0] = 1 + # x[1]: mixed[1] + x-bound-upper[1] = 1 + # x[2]: mixed[2] = 1 + n_terms = (con_x.vars != -1).sum(dim="_term") + np.testing.assert_array_equal( + n_terms.values, + np.array([2, 2, 1]), + ) + + # The dual values from the standalone dual model should match Linopy's + # primal constraint duals for the original mixed constraint block. + np.testing.assert_allclose( + m.constraints["mixed"].dual.values, + m_dual.variables["mixed"].solution.values, + atol=1e-6, + ) From 215febdce68030daab02d6d9f4df72f00a1f14c0 Mon Sep 17 00:00:00 2001 From: bobbyxng Date: Tue, 2 Jun 2026 14:16:54 +0200 Subject: [PATCH 28/30] revert test-models.yml environments to upstream, see PR 522 --- .github/workflows/test-models.yml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/.github/workflows/test-models.yml b/.github/workflows/test-models.yml index 5444455f..14eedae5 100644 --- a/.github/workflows/test-models.yml +++ b/.github/workflows/test-models.yml @@ -77,12 +77,12 @@ jobs: uses: actions/cache@v5 with: path: ${{ env.CONDA }}/envs - key: conda-pypsa-eur-${{ env.week }}-${{ hashFiles('envs/environment.yaml') }} + key: conda-pypsa-eur-${{ env.week }}-${{ hashFiles('envs/linux-64.lock.yaml') }} id: cache-env - name: Update environment if: env.pinned == 'false' && steps.cache-env.outputs.cache-hit != 'true' - run: conda env update -n pypsa-eur -f envs/environment.yaml + run: conda env update -n pypsa-eur -f envs/linux-64.lock.yaml - name: Install package from ref if: env.pinned == 'false' From 66410ad293abae6d9c6baa7c69bb68751da95fae Mon Sep 17 00:00:00 2001 From: bobbyxng Date: Tue, 2 Jun 2026 14:23:48 +0200 Subject: [PATCH 29/30] add release notes and improve docstrings of m.dualize() --- doc/release_notes.rst | 4 ++++ linopy/dual.py | 28 +++++++++++++++++++++++----- 2 files changed, 27 insertions(+), 5 deletions(-) diff --git a/doc/release_notes.rst b/doc/release_notes.rst index 0268e08b..0bac59be 100644 --- a/doc/release_notes.rst +++ b/doc/release_notes.rst @@ -13,6 +13,10 @@ Upcoming Version * After ``model.solve()``, the solver object stays available on ``model.solver``. You can inspect it, reuse it, or release the underlying solver (and its license) by calling ``model.solver.close()`` or assigning ``model.solver = None``. It is also released automatically when the model is garbage-collected. * New ``SolverReport`` on the result (``result.report``) reports runtime, MIP gap, dual (best) bound, and iteration counts. It is shown in ``repr(result)`` and currently populated by CBC, HiGHS, Gurobi, Knitro, and cuPDLPx. +*Dualize LP models* + +* New ``Model.dualize()`` constructs the LP dual as a standalone model, lifting finite variable bounds into explicit constraints so they are reflected in the dual. Dual variables are named after the primal constraints. Works for linear problems with linear objective and constraints. (https://github.com/PyPSA/linopy/pull/626) + *A new way to call a solver (advanced)* Most users should keep calling ``model.solve(...)``. If you want more control, you can now build the solver yourself and run it in two steps: diff --git a/linopy/dual.py b/linopy/dual.py index 6a9b7c17..db616df3 100644 --- a/linopy/dual.py +++ b/linopy/dual.py @@ -651,15 +651,32 @@ def dualize( For a minimization primal: - Primal (min): Dual (max): - min c^T x max b_eq^T λ + b_leq^T μ + b_geq^T ν - s.t. A_eq x = b_eq : λ free s.t. A_eq^T λ + A_leq^T μ + A_geq^T ν = c - A_leq x <= b_leq : μ <= 0 λ free, μ <= 0, ν >= 0 + Primal (min): + min c^T x + s.t. A_eq x = b_eq : λ free + A_leq x <= b_leq : μ <= 0 A_geq x >= b_geq : ν >= 0 + Dual (max): + max b_eq^T λ + b_leq^T μ + b_geq^T ν + s.t. A_eq^T λ + A_leq^T μ + A_geq^T ν = c + λ free, μ <= 0, ν >= 0 + For a maximization primal the dual variable bounds are flipped: μ >= 0 for <= constraints, ν <= 0 for >= constraints. + The corresponding bound conventions are: + + ============ =========== ================ ================ + Constraint Primal sense lower upper + ============ =========== ================ ================ + = min / max -inf +inf (free) + <= min -inf 0 + <= max 0 +inf + >= min 0 +inf + >= max -inf 0 + ============ =========== ================ ================ + Variable bounds are converted to explicit constraints before dualization via _lift_bounds_to_constraints(), so that they appear in the constraint matrix A and are correctly reflected in the dual. @@ -673,7 +690,8 @@ def dualize( Note: This constructs a standalone dual model. Pathological or unsupported primal formulations may lead to infeasible or unbounded dual models or may cause this function to raise an error. Only linear primal models with linear - objectives and either linear constraints or finite variable bounds are supported. + objectives and linear constraints are supported; finite variable bounds are + lifted into explicit constraints before dualization. Parameters ---------- From 14f0c245457935c45c25a6cd1c0e82c8d9604153 Mon Sep 17 00:00:00 2001 From: bobbyxng Date: Tue, 2 Jun 2026 14:27:03 +0200 Subject: [PATCH 30/30] docstring improvements --- linopy/dual.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/linopy/dual.py b/linopy/dual.py index db616df3..9bfea082 100644 --- a/linopy/dual.py +++ b/linopy/dual.py @@ -702,7 +702,8 @@ def dualize( Returns ------- Model - Dual model whose variables are named after the primal constraints. + Dual model whose variables are named after the primal constraints and + constraints are named after the primal variables. Example -------