From c22c94b856e85bdd10ea2656cb770f72dd9edde5 Mon Sep 17 00:00:00 2001 From: Bobby Xiong Date: Tue, 17 Mar 2026 13:37:46 +0100 Subject: [PATCH 01/16] 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/16] 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/16] 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/16] 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/16] 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/16] 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/16] 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/16] 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/16] 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/16] 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/16] 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/16] 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/16] 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/16] 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/16] 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/16] 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