diff --git a/.github/workflows/test_and_lint.yml b/.github/workflows/test_and_lint.yml index 655425d..5019dda 100644 --- a/.github/workflows/test_and_lint.yml +++ b/.github/workflows/test_and_lint.yml @@ -17,7 +17,7 @@ jobs: - name: Install Dependencies shell: bash -el {0} run: | - conda install ipopt + conda install ipopt pip python -m pip install --upgrade pip pip install .[test] - name: Test with pytest 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..18874a3 100644 --- a/eeco/costs.py +++ b/eeco/costs.py @@ -503,6 +503,56 @@ def get_charge_df( return charge_df +def get_prev_demand_dict( + charge_dict, + usage_data, + start_dt, + billing_period_starts, + prev_dict=None, +): + """Compute a new or updated previous max demand charge dict for use in costing. + + For each charge in charge_dict, resets the tracked maximum to zero at the + start of the relevant timeframe. Then accumulates the running maximum. + + Parameters + ---------- + charge_dict : dict + Maps charge key strings to charge arrays (numpy arrays or cp.Expression). + + usage_data : numpy.ndarray + Array of consumption values for the current timestep window. + + start_dt : datetime-like + Start datetime of the current optimization window. + + billing_period_starts : list + List of datetimes marking the start of each billing period. + + prev_dict : dict, optional + Existing previous-max dict to update. Defaults to empty dict. + + Returns + ------- + dict + Updated mapping of charge key to previous maximum demand value. + """ + prev_dict = dict(prev_dict or {}) + for charge_name, charge_array in charge_dict.items(): + prev_dict.setdefault(charge_name, 0) + if start_dt in billing_period_starts or ( + get_charge_array_duration(charge_name) == 1 + and pd.Timestamp(start_dt).hour == 0 + ): + prev_dict[charge_name] = 0 + else: + prev_dict[charge_name] = max( + prev_dict[charge_name], + np.max(usage_data * charge_array), + ) + return prev_dict + + def default_varstr_alias_func( utility, charge_type, name, start_date, end_date, charge_limit ): @@ -638,7 +688,9 @@ def calculate_demand_cost( in USD for the given `charge_array` and `consumption_data` and the second entry being the pyomo model object (or None) """ - if isinstance(consumption_estimate, (float, int)): + if isinstance(prev_demand, cp.Expression): # for cp.Parameter + consumption_max = None # evaluate later for parameterized prev_demand + elif isinstance(consumption_estimate, (float, int)): consumption_max = max(float(consumption_estimate), prev_demand) else: consumption_max = max(max(consumption_estimate), prev_demand) @@ -692,17 +744,18 @@ def const_rule(model, t): else: demand_charged = np.array([0]) elif isinstance(consumption_data, cp.Expression): - if consumption_max >= limit: - if consumption_max <= next_limit: + _use_param = consumption_max is None # True when prev_demand is a cp.Parameter + if _use_param or consumption_max >= limit: + if not _use_param and consumption_max > next_limit: demand_charged, model = ut.multiply( - consumption_data - limit, + next_limit - limit, charge_array, model=model, varstr=varstr + "_multiply", ) else: demand_charged, model = ut.multiply( - next_limit - limit, + consumption_data - limit, charge_array, model=model, varstr=varstr + "_multiply", @@ -795,7 +848,7 @@ def calculate_energy_cost( and the second entry being the pyomo model object (or None) """ cost = 0 - if model is None: + if hasattr(consumption_data, "shape"): n_steps = consumption_data.shape[0] else: # Pyomo does not support shape attribute n_steps = len(consumption_data) @@ -971,7 +1024,7 @@ def get_conversion_factors(electric_consumption_units, gas_consumption_units): def get_converted_consumption_data( - consumption_data_dict, conversion_factors, decomposition_type, model=None + consumption_data_dict, conversion_factors, decomposition_type, model=None, big_m=1e6 ): """Ensure consumption_data_dict has imports/exports structure with proper unit conversion. @@ -989,9 +1042,11 @@ def get_converted_consumption_data( Returns ------- - dict - Updated consumption_data_dict with imports/exports structure + tuple + - consumption_data_dict with imports/exports structure + - pyomo ConcreteModel, list of cvxpy constraints, or None """ + cvxpy_constraints = [] for utility in consumption_data_dict.keys(): conversion_factor = conversion_factors[utility] @@ -1030,7 +1085,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) @@ -1038,12 +1093,14 @@ def get_converted_consumption_data( # Decompose if not already done if model is None or not hasattr(model, pos_name): - imports, exports, model = ut.decompose_consumption( + imports, exports, model, new_constraints = ut.decompose_consumption( converted_consumption, model=model, varstr=utility, - decomposition_type="absolute_value", + decomposition_type=decomposition_type, + big_m=big_m, ) + cvxpy_constraints.extend(new_constraints) consumption_data_dict[utility] = { "imports": imports, @@ -1057,7 +1114,7 @@ def get_converted_consumption_data( else: raise NotImplementedError - return consumption_data_dict + return consumption_data_dict, model or cvxpy_constraints or None def get_charge_array_duration(key): @@ -1117,6 +1174,7 @@ def calculate_cost( demand_scale_factor=1, model=None, decomposition_type=None, + big_m=1e6, varstr_alias_func=default_varstr_alias_func, ): """Calculates the cost of given charges (demand or energy) for the given @@ -1201,8 +1259,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 or CVXPY. - None (default): No decomposition, treats all consumption as imports varstr_alias_func: function @@ -1232,11 +1292,11 @@ def calculate_cost( Returns ------- - (numpy.Array, cvxpy.Expression, or pyomo.environ.Var), pyomo.Model - tuple with the first entry being a float, - cvxpy Expression, or pyomo Var representing energy charge costs - in USD for the given `charge_array` and `consumption_data` - and the second entry being the pyomo model object (or None) + tuple + - First entry: float, cvxpy Expression, or pyomo Var representing energy + charge costs in USD for the given `charge_array` and `consumption_data` + - Second entry: pyomo ConcreteModel (pyomo path), list of cvxpy constraints + from decomposition (cvxpy + decomposition_type path), or None otherwise """ cost = 0 n_per_hour = int(60 / ut.get_freq_binsize_minutes(resolution)) @@ -1248,8 +1308,12 @@ def calculate_cost( electric_consumption_units, gas_consumption_units ) - consumption_data_dict = get_converted_consumption_data( - consumption_data_dict, conversion_factors, decomposition_type, model + consumption_data_dict, model_objects = get_converted_consumption_data( + consumption_data_dict, + conversion_factors, + decomposition_type, + model, + big_m=big_m, ) for key, charge_array in charge_dict.items(): @@ -1284,14 +1348,22 @@ def calculate_cost( if isinstance(consumption_estimate, (float, int)): # convert single kWh to the equivalent kW per timestep - demand_consumption_estimate = ( - consumption_estimate * n_per_hour / len(charge_array) + _n = ( + charge_array.size + if isinstance(charge_array, cp.Expression) + else len(charge_array) ) + demand_consumption_estimate = consumption_estimate * n_per_hour / _n elif isinstance(consumption_estimate, (dict)): demand_consumption_estimate = consumption_estimate[utility] if isinstance(demand_consumption_estimate, (float, int)): + _n = ( + charge_array.size + if isinstance(charge_array, cp.Expression) + else len(charge_array) + ) demand_consumption_estimate = ( - demand_consumption_estimate * n_per_hour / len(charge_array) + demand_consumption_estimate * n_per_hour / _n ) else: demand_consumption_estimate = consumption_estimate @@ -1355,7 +1427,7 @@ def calculate_cost( else: raise ValueError("Invalid charge_type: " + charge_type) - return cost, model + return cost, model_objects def build_pyomo_costing( @@ -1449,8 +1521,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 or CVXPY. - None (default): No decomposition, treats all consumption as imports varstr_alias_func: function @@ -1520,6 +1594,7 @@ def calculate_itemized_cost( demand_scale_factor=1, model=None, decomposition_type=None, + big_m=1e6, by_charge_key=False, varstr_alias_func=default_varstr_alias_func, ): @@ -1584,7 +1659,8 @@ def calculate_itemized_cost( decomposition_type : str or None Type of decomposition to use for consumption data. - "absolute_value": Linear problem using absolute value - - "binary_variable": `NotImplementedError` + - "binary_big_M": Uses binary indicator with Big-M constraints. + Creates a MILP (mixed-integer linear program) for Pyomo or CVXPY. - None (default): No decomposition, treats all consumption as imports by_charge_key : bool @@ -1627,23 +1703,20 @@ def calculate_itemized_cost( } - """ - # Check if decomposition_type is used with CVXPY objects - # (not yet supported because imports - exports creates non-DCP issues) - if decomposition_type is not None: - for utility in consumption_data_dict.keys(): - if isinstance(consumption_data_dict[utility], cp.Variable): - raise NotImplementedError( - "Decomposition types are not supported with CVXPY objects. " - "Use Pyomo instead for problems requiring decomposition_type." - ) + model_objects : pyomo ConcreteModel, list of cvxpy constraints, or None + Same as the second return value of `calculate_cost`. + """ conversion_factors = get_conversion_factors( electric_consumption_units, gas_consumption_units ) - consumption_data_dict = get_converted_consumption_data( - consumption_data_dict, conversion_factors, decomposition_type, model + consumption_data_dict, model_objects = get_converted_consumption_data( + consumption_data_dict, + conversion_factors, + decomposition_type, + model, + big_m=big_m, ) results_dict = {} @@ -1659,7 +1732,7 @@ def calculate_itemized_cost( for charge_key, charges in charge_items: charge_input = charges if charge_key is None else {charge_key: charges} - cost, model = calculate_cost( + cost, _ = calculate_cost( charge_input, consumption_data_dict, resolution=resolution, @@ -1670,7 +1743,6 @@ def calculate_itemized_cost( desired_charge_type=charge_type, demand_scale_factor=demand_scale_factor, model=model, - decomposition_type=decomposition_type, varstr_alias_func=varstr_alias_func, ) if by_charge_key: @@ -1702,7 +1774,7 @@ def calculate_itemized_cost( results_dict[utility]["total"] for utility in utilities ) - return results_dict, model + return results_dict, model_objects def detect_charge_periods( diff --git a/eeco/tests/test_costs.py b/eeco/tests/test_costs.py index bb15259..bc36900 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) @@ -1157,6 +1157,45 @@ def test_calculate_cost_np( None, pytest.approx(120.0), ), + ( + { + "electric_demand_peak-summer_2024-07-10_2024-07-10_0": np.concatenate( + [np.ones(48) * 0, np.ones(24) * 1, np.ones(24) * 0] + ), + "electric_demand_half-peak-summer_2024-07-10_2024-07-10_0": ( + np.concatenate( + [ + np.ones(34) * 0, + np.ones(14) * 2, + np.ones(24) * 0, + np.ones(14) * 2, + np.ones(10) * 0, + ] + ) + ), + "electric_demand_off-peak_2024-07-10_2024-07-10_0": np.ones(96) * 10, + }, + {ELECTRIC: np.arange(96), GAS: np.arange(96)}, + "15m", + { + "electric_demand_peak-summer_2024-07-10_2024-07-10_0": { + "demand": cp.Parameter(nonneg=True, value=0), + "cost": 0, + }, + "electric_demand_half-peak-summer_2024-07-10_2024-07-10_0": { + "demand": cp.Parameter(nonneg=True, value=0), + "cost": 0, + }, + "electric_demand_off-peak_2024-07-10_2024-07-10_0": { + "demand": cp.Parameter(nonneg=True, value=0), + "cost": 0, + }, + }, + 0, # consumption_estimate + None, # desired_utility + None, # desired_charge_type + pytest.approx(1191), # same expected as float-zero version + ), ], ) def test_calculate_cost_cvx( @@ -1480,6 +1519,64 @@ 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 + ), + ( + { + "electric_demand_off-peak_2024-07-10_2024-07-10_0": np.ones(96) * 10, + "electric_demand_off-peak_2024-07-10_2024-07-10_90": np.ones(96) * 5, + }, + {ELECTRIC: np.arange(96), GAS: np.arange(96)}, + "15m", + { + "electric_demand_off-peak_2024-07-10_2024-07-10_0": { + "demand": 100, + "cost": 0, + }, + "electric_demand_off-peak_2024-07-10_2024-07-10_90": { + "demand": 0, + "cost": 0, + }, + }, + 0, # consumption_estimate + None, # desired_utility + None, # desired_charge_type + None, # decomposition_type + pytest.approx(900), + ), ], ) def test_calculate_cost_pyo( @@ -3232,7 +3329,7 @@ def test_calculate_itemized_cost_np( expected_itemized, ): """Test calculate_itemized_cost with and without decomposition_type.""" - result, model = costs.calculate_itemized_cost( + result, _ = costs.calculate_itemized_cost( charge_dict, consumption_data_dict, resolution=resolution, @@ -3356,7 +3453,7 @@ def test_calculate_itemized_cost_cvx( by_charge_key=by_charge_key, ) else: - result, model = costs.calculate_itemized_cost( + result, _ = costs.calculate_itemized_cost( charge_dict, cvx_vars, resolution=resolution, @@ -3539,7 +3636,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 +3647,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 ( @@ -3615,39 +3725,33 @@ def test_calculate_itemized_cost_pyo( if gas_consumption_units is not None: kwargs["gas_consumption_units"] = gas_consumption_units - if decomposition_type == "binary_variable": - with pytest.raises(NotImplementedError): - result, model = costs.calculate_itemized_cost( - charge_dict, pyo_vars, **kwargs - ) - else: - result, model = costs.calculate_itemized_cost(charge_dict, pyo_vars, **kwargs) - total_expr = sum(result["total"].values()) if by_charge_key else result["total"] - solve_pyo_problem( - model, - total_expr, - decomposition_type, - charge_dict, - consumption_data_dict, - by_charge_key, - ) + result, model = costs.calculate_itemized_cost(charge_dict, pyo_vars, **kwargs) + total_expr = sum(result["total"].values()) if by_charge_key else result["total"] + solve_pyo_problem( + model, + total_expr, + decomposition_type, + charge_dict, + consumption_data_dict, + by_charge_key, + ) - total_value = ( - sum(pyo.value(v) for v in result["total"].values()) - if by_charge_key - else pyo.value(result["total"]) - ) - assert total_value == expected_cost - for utility in expected_itemized: - for charge_type in expected_itemized[utility]: - expected_value = expected_itemized[utility][charge_type] - actual_value = result[utility][charge_type] - if by_charge_key: - # Compare pyo.value for each key - actual_value = {k: pyo.value(v) for k, v in actual_value.items()} - assert actual_value == expected_value - else: - assert pyo.value(actual_value) == expected_value + total_value = ( + sum(pyo.value(v) for v in result["total"].values()) + if by_charge_key + else pyo.value(result["total"]) + ) + assert total_value == expected_cost + for utility in expected_itemized: + for charge_type in expected_itemized[utility]: + expected_value = expected_itemized[utility][charge_type] + actual_value = result[utility][charge_type] + if by_charge_key: + # Compare pyo.value for each key + actual_value = {k: pyo.value(v) for k, v in actual_value.items()} + assert actual_value == expected_value + else: + assert pyo.value(actual_value) == expected_value @pytest.mark.parametrize( @@ -3687,3 +3791,34 @@ def test_individual_charge(charge_list, expected_result): model=None, ) assert cost == expected_result + + +@pytest.mark.parametrize( + "start_dt, billing_period_starts, prev_dict, expected", + [ + # start_dt in billing_period_starts -> reset to 0 + ( + np.datetime64("2024-07-10"), + [np.datetime64("2024-07-10")], + {"electric_demand_peak_20240710_20240731_100": 5.0}, + {"electric_demand_peak_20240710_20240731_100": 0}, + ), + # start_dt not in billing_period_starts, multi-day key -> accumulate max + ( + np.datetime64("2024-07-11"), + [np.datetime64("2024-07-10")], + {"electric_demand_peak_20240710_20240731_100": 2.0}, + {"electric_demand_peak_20240710_20240731_100": 3.0}, + ), + ], +) +@pytest.mark.skipif(skip_all_tests, reason="Exclude all tests") +def test_get_prev_demand_dict(start_dt, billing_period_starts, prev_dict, expected): + charge_dict = { + "electric_demand_peak_20240710_20240731_100": np.ones(3), + } + usage_data = np.array([1.0, 2.0, 3.0]) + result = costs.get_prev_demand_dict( + charge_dict, usage_data, start_dt, billing_period_starts, prev_dict + ) + assert result == expected diff --git a/eeco/tests/test_utils.py b/eeco/tests/test_utils.py index c629fbb..5b981ab 100644 --- a/eeco/tests/test_utils.py +++ b/eeco/tests/test_utils.py @@ -243,31 +243,82 @@ def test_decompose_consumption_np( with pytest.raises(TypeError): ut.decompose_consumption(consumption_data) else: - positive_values, negative_values, model = ut.decompose_consumption( + positive_values, negative_values, model, constraints = ut.decompose_consumption( consumption_data ) assert np.array_equal(positive_values, expected_positive) assert np.array_equal(negative_values, expected_negative) assert model is None + assert constraints == [] assert np.array_equal(consumption_data, positive_values - negative_values) @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, model, 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 model is None + 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 +330,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( @@ -297,7 +351,7 @@ def test_decompose_consumption_pyo( if expect_warning: with pytest.warns(UserWarning): - positive_var, negative_var, model = ut.decompose_consumption( + positive_var, negative_var, model, constraints = ut.decompose_consumption( pyo_vars["electric"], model=model, varstr="electric", @@ -305,18 +359,25 @@ def test_decompose_consumption_pyo( ) assert positive_var is None assert negative_var is None + assert constraints == [] else: - positive_var, negative_var, model = ut.decompose_consumption( + positive_var, negative_var, model, constraints = ut.decompose_consumption( pyo_vars["electric"], model=model, varstr="electric", decomposition_type=decomposition_type, ) + assert constraints == [] # Check that variables exist and have the correct length 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..7d84c64 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,76 @@ 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, model, []) + - Pyomo: (positive_var, negative_var, model, []) - constraints added to model + - CVXPY: (positive_var, negative_var, None, constraints) - list of constraints + that must be added to the Problem; model is always None for CVXPY """ 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 + 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": + positive_var, negative_var, constraints = _decompose_binary_cvx( + expression, big_m + ) + return positive_var, negative_var, None, constraints 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 positive_var, negative_var, model + 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."