Adding binary (binary_big_M) decomposition type

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:

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

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
         (

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