diff --git a/docs/how_to_advanced.rst b/docs/how_to_advanced.rst index 78d64b9..f703882 100644 --- a/docs/how_to_advanced.rst +++ b/docs/how_to_advanced.rst @@ -145,7 +145,7 @@ When `demand_scale_factor < 1.0`, demand charges are proportionally reduced to r .. code-block:: python - from electric_emission_cost import costs + from eeco import costs # E.g. solving for 3 days out of a 30-day billing period demand_scale_factor = 3 / 30 @@ -190,34 +190,67 @@ How to Use `decomposition_type` The `decomposition_type` parameter allows you to decompose consumption data into positive (imports) and negative (exports) components. This is useful when you have export charges or credits in your rate structure. Options include: -- Default `None` -- `"binary_variable"`: To be implemented -- `"absolute_value"` +- Default `None`: No decomposition, consumption treated as imports only. +- `"binary_big_M"`: Uses binary variable with Big-M constraints. Creates a MILP requiring a MIP solver (e.g., Gurobi). Supported for both CVXPY and Pyomo. +- `"absolute_value"`: Uses absolute value constraints. Creates a nonlinear problem for Pyomo. **Not supported for CVXPY** (not DCP-compliant). + +For numpy arrays, `decomposition_type` is ignored since decomposition is a direct calculation. + +**Using Pyomo** .. code-block:: python - from electric_emission_cost import costs + from eeco import costs + import pyomo.environ as pyo - # Example with export charges + # Example with export charges using Pyomo charge_dict = { + "electric_energy_0_2024-07-10_2024-07-10_0": np.ones(96) * 0.05, "electric_export_0_2024-07-10_2024-07-10_0": np.ones(96) * 0.025, } + # Create Pyomo model with consumption variable + model = pyo.ConcreteModel() + model.t = pyo.RangeSet(0, 95) + model.consumption = pyo.Var(model.t, bounds=(None, None)) + consumption_data = { - "electric": np.concatenate([np.ones(48) * 10, -np.ones(48) * 5]), - "gas": np.ones(96), + "electric": model.consumption, + "gas": pyo.Param(model.t, initialize=1), } # Decompose consumption into imports and exports result, model = costs.calculate_cost( charge_dict, consumption_data, - decomposition_type="absolute_value" + decomposition_type="binary_big_M", # or "absolute_value" for Pyomo + model=model, + ) + +**Using CVXPY** + +.. code-block:: python + + from eeco import costs, utils + import cvxpy as cp + + # Create CVXPY variable for consumption + consumption = cp.Variable(96) + + # Decompose into imports/exports (returns constraints for CVXPY) + imports, exports, decomp_constraints = utils.decompose_consumption( + consumption, decomposition_type="binary_big_M" ) + + # Build your optimization problem with decomposition constraints + objective = cp.Minimize(...) # your objective + constraints = decomp_constraints + [...] # add other constraints + + prob = cp.Problem(objective, constraints) + prob.solve(solver=cp.GUROBI) # requires MIP solver -When decomposition_type is not None the function creates separate variables for positive consumption (imports) and negative consumption (exports) -and applies export charges only to the export component. -For Pyomo models, decomposition_type adds a constraint total_consumption = imports - exports +When `decomposition_type` is not `None`, the function creates separate variables for positive consumption (imports) and negative consumption (exports), applying export charges only to the export component. +A constraint `total_consumption = imports - exports` is added to balance the decomposition. .. _varstr-alias: diff --git a/eeco/costs.py b/eeco/costs.py index ef483fa..8682c11 100644 --- a/eeco/costs.py +++ b/eeco/costs.py @@ -1030,7 +1030,7 @@ def get_converted_consumption_data( varstr=converted_varstr, ) - if decomposition_type == "absolute_value": + if decomposition_type in ("absolute_value", "binary_big_M"): # Decompose consumption data into positive and negative components # with constraint that total = positive - negative # (where negative is stored as positive magnitude) @@ -1042,7 +1042,7 @@ def get_converted_consumption_data( converted_consumption, model=model, varstr=utility, - decomposition_type="absolute_value", + decomposition_type=decomposition_type, ) consumption_data_dict[utility] = { @@ -1201,8 +1201,10 @@ def calculate_cost( decomposition_type : str or None Type of decomposition to use for consumption data. - - "absolute_value": Linear problem using absolute value - - "binary_variable": `NotImplementedError` + - "absolute_value": Uses max(x, 0) constraints. Creates nonlinear problem + for Pyomo due to abs() constraint. + - "binary_big_M": Uses binary indicator with Big-M constraints. + Creates a MILP (mixed-integer linear program) for Pyomo. - None (default): No decomposition, treats all consumption as imports varstr_alias_func: function @@ -1449,8 +1451,10 @@ def build_pyomo_costing( decomposition_type : str or None Type of decomposition to use for consumption data. - - "absolute_value": Linear problem using absolute value - - "binary_variable": `NotImplementedError` + - "absolute_value": Uses max(x, 0) constraints. Creates nonlinear problem + for Pyomo due to abs() constraint. + - "binary_big_M": Uses binary indicator with Big-M constraints. + Creates a MILP (mixed-integer linear program) for Pyomo. - None (default): No decomposition, treats all consumption as imports varstr_alias_func: function diff --git a/eeco/tests/test_costs.py b/eeco/tests/test_costs.py index bb15259..707a6e6 100644 --- a/eeco/tests/test_costs.py +++ b/eeco/tests/test_costs.py @@ -109,9 +109,9 @@ def solve_pyo_problem( """Helper function to solve Pyomo optimization problem.""" model.obj = pyo.Objective(expr=objective) - if decomposition_type is not None: # Nonlinear when decomposition_type used + if decomposition_type == "absolute_value": # Nonlinear due to abs() constraint solver = pyo.SolverFactory("ipopt") - else: # Gurobi otherwise + else: # MILP or LP can use Gurobi solver = pyo.SolverFactory("gurobi") solver.solve(model) @@ -1480,6 +1480,41 @@ def test_calculate_cost_cvx( None, pytest.approx(9.0), ), + # binary_big_M decomposition: export charges only + ( + { + "electric_export_0_2024-07-10_2024-07-10_0": np.ones(96) * 0.025, + }, + { + ELECTRIC: np.concatenate([np.ones(48), -np.ones(48)]), + GAS: np.ones(96), + }, + "15m", + None, + 0, + None, + None, + "binary_big_M", + pytest.approx(-0.3), + ), + # binary_big_M decomposition: energy and export charges + ( + { + "electric_energy_0_2024-07-10_2024-07-10_0": np.ones(96) * 0.05, + "electric_export_0_2024-07-10_2024-07-10_0": np.ones(96) * 0.025, + }, + { + ELECTRIC: np.concatenate([np.ones(48), -np.ones(48)]), + GAS: np.ones(96), + }, + "15m", + None, + 0, + None, + None, + "binary_big_M", + pytest.approx(0.6 - 0.3), # 48*1*0.05/4 - 48*1*0.025/4 = 0.6 - 0.3 = 0.3 + ), ], ) def test_calculate_cost_pyo( @@ -3539,7 +3574,7 @@ def test_calculate_itemized_cost_cvx( }, }, ), - # `binary_variable` for `decomposition_type` should raise `NotImplementedError` + # binary_big_M decomposition (MILP) ( { "electric_energy_0_2024-07-10_2024-07-10_0": np.ones(96) * 0.05, @@ -3550,13 +3585,26 @@ def test_calculate_itemized_cost_cvx( GAS: np.ones(96), }, "15m", - "binary_variable", + "binary_big_M", 240, None, None, False, - None, # No expected cost - should raise NotImplementedError - None, # No expected itemized - should raise NotImplementedError + pytest.approx(6.0 - 1.5), # 48*10*0.05/4 - 48*5*0.025/4 = 6.0 - 1.5 = 4.5 + { + "electric": { + "energy": pytest.approx(6.0), # 48*10*0.05/4 = 6.0 + "export": pytest.approx(-1.5), # -48*5*0.025/4 = 1.5 + "customer": 0.0, + "demand": 0.0, + }, + "gas": { + "energy": 0.0, + "export": 0.0, + "customer": 0.0, + "demand": 0.0, + }, + }, ), # by_charge_key=True ( diff --git a/eeco/tests/test_utils.py b/eeco/tests/test_utils.py index c629fbb..b4f201f 100644 --- a/eeco/tests/test_utils.py +++ b/eeco/tests/test_utils.py @@ -254,20 +254,69 @@ def test_decompose_consumption_np( @pytest.mark.skipif(skip_all_tests, reason="Exclude all tests") -def test_decompose_consumption_cvx(): - """Test decompose_consumption with cvxpy expressions.""" - x = cp.Variable(5) - positive_values, negative_values, model = ut.decompose_consumption(x) - assert isinstance(positive_values, cp.Expression) - assert isinstance(negative_values, cp.Expression) - - # Test warning for unimplemented decomposition_type - with pytest.warns(UserWarning): - positive_values, negative_values, model = ut.decompose_consumption( - x, decomposition_type="unimplemented" +@pytest.mark.parametrize( + "consumption_data, expected_positive, expected_negative, " + "decomposition_type, expect_error", + [ + # absolute_value not DCP-compliant, raises NotImplementedError + (np.array([1, -2, 3]), None, None, "absolute_value", True), + (np.array([1, -2, 3]), None, None, None, True), # default is absolute_value + # binary_big_M works with MIP solver + ( + np.array([10, -5, 3, -2, 0]), + np.array([10, 0, 3, 0, 0]), + np.array([0, 5, 0, 2, 0]), + "binary_big_M", + False, + ), + ( + np.array([0, 0, 0]), + np.array([0, 0, 0]), + np.array([0, 0, 0]), + "binary_big_M", + False, + ), + ( + np.array([-10, -5, -1]), + np.array([0, 0, 0]), + np.array([10, 5, 1]), + "binary_big_M", + False, + ), + ], +) +def test_decompose_consumption_cvx( + consumption_data, + expected_positive, + expected_negative, + decomposition_type, + expect_error, +): + """Test decompose_consumption with CVXPY expressions.""" + x = cp.Variable(len(consumption_data)) + + if expect_error: + with pytest.raises(NotImplementedError): + if decomposition_type is None: + ut.decompose_consumption(x) + else: + ut.decompose_consumption(x, decomposition_type=decomposition_type) + else: + positive_var, negative_var, constraints = ut.decompose_consumption( + x, decomposition_type=decomposition_type ) - assert positive_values is None - assert negative_values is None + assert isinstance(positive_var, cp.Variable) + assert isinstance(negative_var, cp.Variable) + assert isinstance(constraints, list) + assert len(constraints) == 3 # decomposition + 2 Big-M constraints + + # Solve with Gurobi and verify values + constraints.append(x == consumption_data) + prob = cp.Problem(cp.Minimize(0), constraints) + prob.solve(solver=cp.GUROBI) + + np.testing.assert_array_almost_equal(positive_var.value, expected_positive) + np.testing.assert_array_almost_equal(negative_var.value, expected_negative) @pytest.mark.skipif(skip_all_tests, reason="Exclude all tests") @@ -279,7 +328,10 @@ def test_decompose_consumption_cvx(): (np.array([0, 0, 0]), 0, 0, "absolute_value", False), (np.array([-10, -5, -1]), 0, 16, "absolute_value", False), (np.array([10, 5, 1]), 16, 0, "absolute_value", False), - (np.array([1, -2, 3]), 4, 2, "binary_variable", True), + (np.array([1, -2, 3]), 4, 2, "binary_big_M", False), + (np.array([-10, -5, -1]), 0, 16, "binary_big_M", False), + (np.array([10, 5, 1]), 16, 0, "binary_big_M", False), + (np.array([1, -2, 3]), 4, 2, "unsupported_type", True), ], ) def test_decompose_consumption_pyo( @@ -316,7 +368,12 @@ def test_decompose_consumption_pyo( assert hasattr(model, "electric_positive") assert hasattr(model, "electric_negative") assert hasattr(model, "electric_decomposition_constraint") - assert hasattr(model, "electric_magnitude_constraint") + if decomposition_type == "absolute_value": + assert hasattr(model, "electric_magnitude_constraint") + elif decomposition_type == "binary_big_M": + assert hasattr(model, "electric_is_importing") + assert hasattr(model, "electric_import_bigm_constraint") + assert hasattr(model, "electric_export_bigm_constraint") assert len(positive_var) == len(consumption_data) assert len(negative_var) == len(consumption_data) # Testing of values handled after solving problem in test_costs.py diff --git a/eeco/utils.py b/eeco/utils.py index b9c0182..4ace39d 100644 --- a/eeco/utils.py +++ b/eeco/utils.py @@ -448,11 +448,168 @@ def get_decomposed_var_names(utility): return f"{utility}_positive", f"{utility}_negative" +def _decompose_binary_cvx(expression, big_m=1e6): + """Decompose CVXPY expression using binary variable method with Big-M. + + Creates a mixed-integer program (MIP) where a binary variable indicates + whether we are importing (1) or exporting (0) at each timestep. + Requires a MIP solver (e.g., Gurobi, MOSEK). + + Parameters + ---------- + expression : cvxpy.Expression + CVXPY expression representing net consumption + big_m : float, optional + Big-M value for constraints. Default is 1e6. + + Returns + ------- + tuple + (positive_var, negative_var, constraints) where constraints is a list + of CVXPY constraints that must be added to the problem + """ + n = expression.shape[0] if hasattr(expression, "shape") else 1 + + # Binary variable: 1 = importing, 0 = exporting + is_importing = cp.Variable(n, boolean=True) + + # Import/export variables (non-negative) + positive_var = cp.Variable(n, nonneg=True) + negative_var = cp.Variable(n, nonneg=True) + + # Constraints for binary decomposition + constraints = [ + # Decomposition balance: expression = imports - exports + expression == positive_var - negative_var, + # Big-M: imports <= big_m * is_importing + positive_var <= big_m * is_importing, + # Big-M: exports <= big_m * (1 - is_importing) + negative_var <= big_m * (1 - is_importing), + ] + + return positive_var, negative_var, constraints + + +def _decompose_absolute_value_pyo(expression, model, varstr): + """Create Pyomo vars and add absolute value specific constraints. + + Uses max_pos constraints and magnitude constraint with abs(). + Creates a nonlinear problem. + + Parameters + ---------- + expression : pyomo.environ.Var or pyomo.environ.Param + Pyomo variable representing net consumption + model : pyomo.environ.Model + The Pyomo model object + varstr : str + Name prefix for created variables + + Returns + ------- + tuple + (positive_var, negative_var, model) + """ + pos_name, neg_name = get_decomposed_var_names(varstr) + positive_var, model = max_pos(expression, model, pos_name) + + # Create negative expression since pyomo won't take -expression directly + def negative_rule(model, t): + return -expression[t] + + negative_expr = pyo.Expression(model.t, rule=negative_rule) + model.add_component(f"{varstr}_negative_expr", negative_expr) + negative_var, model = max_pos(negative_expr, model, neg_name) + + # Add constraint to ensure positive_var + negative_var = |expression| + # This prevents both variables becoming larger due to artificial arbitrage + def magnitude_rule(model, t): + return positive_var[t] + negative_var[t] == abs(expression[t]) + + model.add_component( + f"{varstr}_magnitude_constraint", + pyo.Constraint(model.t, rule=magnitude_rule), + ) + + return positive_var, negative_var, model + + +def _decompose_binary_pyo(expression, model, varstr, big_m=1e6): + """Create Pyomo vars and add binary/Big-M specific constraints. + + Creates a MILP where a binary variable indicates import (1) or export (0). + + Parameters + ---------- + expression : pyomo.environ.Var or pyomo.environ.Param + Pyomo variable representing net consumption + model : pyomo.environ.Model + The Pyomo model object + varstr : str + Name prefix for created variables + big_m : float, optional + Big-M value for constraints. Default is 1e6. + + Returns + ------- + tuple + (positive_var, negative_var, model) + """ + pos_name, neg_name = get_decomposed_var_names(varstr) + + # Binary variable: 1 = importing, 0 = exporting + binary_name = f"{varstr}_is_importing" + model.add_component( + binary_name, + pyo.Var(model.t, within=pyo.Binary, initialize=1), + ) + binary_var = model.find_component(binary_name) + + # Import variable (positive consumption) + model.add_component( + pos_name, + pyo.Var(model.t, bounds=(0, None), initialize=0), + ) + positive_var = model.find_component(pos_name) + + # Export variable (magnitude of negative consumption, stored as positive) + model.add_component( + neg_name, + pyo.Var(model.t, bounds=(0, None), initialize=0), + ) + negative_var = model.find_component(neg_name) + + # Big-M constraints to enforce mutual exclusivity: + # If binary=1: imports can be positive, exports must be 0 + # If binary=0: imports must be 0, exports can be positive + + # Constraint: imports <= big_m * binary (imports=0 when binary=0) + def import_bigm_rule(model, t): + return positive_var[t] <= big_m * binary_var[t] + + model.add_component( + f"{varstr}_import_bigm_constraint", + pyo.Constraint(model.t, rule=import_bigm_rule), + ) + + # Constraint: exports <= big_m * (1 - binary) (exports=0 when binary=1) + def export_bigm_rule(model, t): + return negative_var[t] <= big_m * (1 - binary_var[t]) + + model.add_component( + f"{varstr}_export_bigm_constraint", + pyo.Constraint(model.t, rule=export_bigm_rule), + ) + + return positive_var, negative_var, model + + def decompose_consumption( - expression, model=None, varstr=None, decomposition_type="absolute_value" + expression, model=None, varstr=None, decomposition_type="absolute_value", big_m=1e6 ): - """Decomposes consumption data into positive and negative components - And adds constraint such that total consumption equals + """Decomposes consumption data into positive and negative components. + + Adds constraint such that total consumption equals positive values minus negative values (where negative values are stored as positive magnitudes). @@ -477,77 +634,73 @@ def decompose_consumption( decomposition_type : str Type of decomposition to use. - - "binary_variable": To be implemented - - "absolute_value": Creates nonlinear problem + - "absolute_value": Uses max(x, 0) constraints. Creates nonlinear problem + for Pyomo due to abs() constraint. Not supported for CVXPY. + - "binary_big_M": Uses binary indicator with Big-M constraints. + Creates a MILP (mixed-integer linear program). + Supported for both Pyomo and CVXPY (requires MIP solver). + + Note: For numpy.ndarray inputs, decomposition_type is ignored since + the decomposition is a direct calculation, not an optimization variable. + + big_m : float, optional + Big-M value for binary decomposition. Should be larger than maximum + possible consumption magnitude. Default is 1e6. Only used when + decomposition_type="binary_big_M". Returns ------- tuple - (positive_values, negative_values, model) where - positive_values and negative_values are both positive - with the constraint that total = positive - negative + - numpy: (positive_values, negative_values, None) + - Pyomo: (positive_var, negative_var, model) - constraints added to model + - CVXPY: (positive_var, negative_var, constraints) - list of constraints + that must be added to the Problem """ if isinstance(expression, np.ndarray): positive_values = np.maximum(expression, 0) negative_values = np.maximum(-expression, 0) # magnitude as positive return positive_values, negative_values, model + elif isinstance(expression, cp.Expression): - if decomposition_type == "absolute_value": - positive_values, _ = max_pos(expression) - negative_values, _ = max_pos(-expression) # magnitude as positive + if decomposition_type == "binary_big_M": + return _decompose_binary_cvx(expression, big_m) else: - warnings.warn( - f"Decomposition type '{decomposition_type}' is not implemented yet. " - "Please use available type. Skipping decomposition.", - UserWarning, + raise NotImplementedError( + f"Decomposition type '{decomposition_type}' not supported for CVXPY. " + "Only 'binary_big_M' is available (requires MIP solver like Gurobi). " + "Use Pyomo for 'absolute_value' decomposition." ) - positive_values = None - negative_values = None - return positive_values, negative_values, model + elif isinstance(expression, (pyo.Var, pyo.Param)): + # Call mode-specific function to create vars and add mode-specific constraints if decomposition_type == "absolute_value": - - # Use max_pos to create positive_var and negative_var - pos_name, neg_name = get_decomposed_var_names(varstr) - positive_var, model = max_pos(expression, model, pos_name) - - # Create negative expression since pyomo won't take -expression directly - def negative_rule(model, t): - return -expression[t] - - negative_expr = pyo.Expression(model.t, rule=negative_rule) - model.add_component(f"{varstr}_negative_expr", negative_expr) - negative_var, model = max_pos(negative_expr, model, neg_name) - - # Add constraint to balance import and export decomposed values - def decomposition_rule(model, t): - return expression[t] == positive_var[t] - negative_var[t] - - model.add_component( - f"{varstr}_decomposition_constraint", - pyo.Constraint(model.t, rule=decomposition_rule), + positive_var, negative_var, model = _decompose_absolute_value_pyo( + expression, model, varstr ) - - # Add constraint to ensure positive_var + negative_var = |expression| - # This both variables becoming larger due to artificial arbitrage - def magnitude_rule(model, t): - return positive_var[t] + negative_var[t] == abs(expression[t]) - - model.add_component( - f"{varstr}_magnitude_constraint", - pyo.Constraint(model.t, rule=magnitude_rule), + elif decomposition_type == "binary_big_M": + positive_var, negative_var, model = _decompose_binary_pyo( + expression, model, varstr, big_m ) - - return positive_var, negative_var, model else: warnings.warn( - f"Decomposition type '{decomposition_type}' is not implemented yet. " - "Please use available type. Skipping decomposition.", + f"Decomposition type '{decomposition_type}' is not implemented. " + "Available types: 'absolute_value', 'binary_big_M'. " + "Skipping decomposition.", UserWarning, ) - positive_var = None - negative_var = None + return None, None, model + + # Add common decomposition constraint: expression = imports - exports + def decomposition_rule(model, t): + return expression[t] == positive_var[t] - negative_var[t] + + model.add_component( + f"{varstr}_decomposition_constraint", + pyo.Constraint(model.t, rule=decomposition_rule), + ) + return positive_var, negative_var, model + else: raise TypeError( "Only CVXPY or Pyomo variables and NumPy arrays are currently supported."