[ENG-1694]: Fix storage memory (#147)

* perf: attempt fix for storage memory explosion

* perf: make where mask drop explicit and set some true

* perf: only drop masked coefficients in storage rate of (dis)charge expressions

tests/test_integration/test_linopy_model.py

@@ -24,7 +24,7 @@ def test_linopy_model():
 
             model = Model.from_otoole_csv(sample_path)
 
-            model.solve()
+            model.solve(solver="highs")
 
             ref_results_df = pd.read_csv(Path(results_path) / "TotalDiscountedCost.csv")
 

tests/test_solve/test_growth_rates.py

@@ -80,7 +80,7 @@ def test_growth_rate_floor():
         technologies=technologies,
     )
 
-    model.solve()
+    model.solve(solver="highs")
 
     assert model.solution.NewCapacity.sel(YEAR=2022, TECHNOLOGY="generator") == 0.0606
     assert model.solution.NewCapacity.sel(YEAR=2020, TECHNOLOGY="unmet-demand") == 99.0
@@ -134,7 +134,7 @@ def test_growth_rate_ceil():
         technologies=technologies,
     )
 
-    model.solve()
+    model.solve(solver="highs")
 
     assert model.solution.NewCapacity.sel(YEAR=2021, TECHNOLOGY="generator") == 15.0
     assert model.solution.NewCapacity.sel(YEAR=2024, TECHNOLOGY="generator") == 20.0
@@ -201,7 +201,7 @@ def test_growth_rate_min():
         technologies=technologies,
     )
 
-    model.solve()
+    model.solve(solver="highs")
 
     assert model.solution.NewCapacity.sel(YEAR=2021, TECHNOLOGY="bad-generator") == 1.0
     assert model.solution.NewCapacity.sel(YEAR=2022, TECHNOLOGY="bad-generator") == 1.1

tests/test_solve/test_min_capacity_factors.py

@@ -47,7 +47,7 @@ def test_min_capacity_factors():
         technologies=technologies,
     )
 
-    model.solve()
+    model.solve(solver="highs")
 
     assert (
         model.solution["TotalTechnologyAnnualActivity"].sel(YEAR=2025, TECHNOLOGY="unmet-demand")

tests/test_solve/test_salvage.py

@@ -33,7 +33,7 @@ def test_salvage_value_straight_line():
         technologies=technologies,
     )
 
-    model.solve()
+    model.solve(solver="highs")
 
     assert np.round(model.solution.DiscountedSalvageValue.sum().values) == 1203.0
 
@@ -68,6 +68,6 @@ def test_salvage_value_sinking_fund():
         technologies=technologies,
     )
 
-    model.solve()
+    model.solve(solver="highs")
 
     assert np.round(model.solution.DiscountedSalvageValue.sum().values) == 1349.0

tests/test_solve/test_solve.py

@@ -14,7 +14,7 @@ def test_model_construction_from_yaml():
 
     model._build()
 
-    model._m.solve()
+    model._m.solve(solver_name="highs")
 
     assert model._m.termination_condition == "optimal"
     assert np.round(model._m.objective.value) == 29044.0
@@ -28,7 +28,7 @@ def test_model_solve_from_otoole_csv():
     path = "examples/otoole_compat/input_csv/otoole-simple-hydro"
 
     model = Model.from_otoole_csv(path)
-    model.solve()
+    model.solve(solver="highs")
 
     assert model._m.termination_condition == "optimal"
     assert np.round(model._m.objective.value) == 5591653.0
@@ -59,7 +59,7 @@ def test_most_simple():
 
     model._build()
 
-    model._m.solve()
+    model._m.solve(solver_name="highs")
 
     assert model._m.termination_condition == "optimal"
     assert np.round(model._m.objective.value) == 45736.0
@@ -182,7 +182,7 @@ def test_simple_storage():
         storage=storage,
         technologies=technologies,
     )
-    model.solve()
+    model.solve(solver="highs")
 
     assert model.solution.NewStorageCapacity.values[0][0][0] == 12.5
     assert model.solution.NetCharge[0][1][0][0] == 75  # bat-storage 2020 Day charge
@@ -235,7 +235,7 @@ def test_simple_trade():
         ],
     )
 
-    model.solve()
+    model.solve(solver="highs")
 
     assert model.solution["NetTrade"].values[0][2][0] == 15
     assert np.round(model._m.objective.value) == 28200.0
@@ -288,7 +288,7 @@ def test_simple_re_target():
 
     model._build()
 
-    model._m.solve()
+    model._m.solve(solver_name="highs")
 
     assert model._m.termination_condition == "optimal"
     assert np.round(model._m.objective.value) == 54671.0

tz/osemosys/model/constraints/re_targets.py

@@ -53,12 +53,7 @@ def add_re_targets_constraints(ds: xr.Dataset, m: Model, lex: Dict[str, LinearEx
         UseByTechnologyAnnual[r,t,f,y];
     ```
     """
-
-    con = (
-        lex["ProductionAnnualRE"]
-        >= lex["ProductionAnnual"].assign_coords({"FUEL": ds["RETagFuel"] == 1})
-        * ds["REMinProductionTarget"]
-    )
+    con = lex["ProductionAnnualRE"] >= lex["ProductionAnnual"] * ds["REMinProductionTarget"]
     mask = ds["RETagFuel"] == 1
     m.add_constraints(con, name="RE1_RenewableProduction_MinConstraint", mask=mask)
 

tz/osemosys/model/linear_expressions/emissions.py

@@ -8,18 +8,19 @@ def add_lex_emissions(ds: xr.Dataset, m: Model, lex: Dict[str, LinearExpression]
     AnnualTechnologyEmissionByMode = (
         (ds["EmissionActivityRatio"] * ds["YearSplit"] * m["RateOfActivity"])
         .sum("TIMESLICE")
-        .where(ds["EmissionActivityRatio"].notnull())
+        .where(ds["EmissionActivityRatio"].notnull(), drop=False)
     )
 
     AnnualTechnologyEmission = AnnualTechnologyEmissionByMode.sum(dims="MODE_OF_OPERATION").where(
-        ds["EmissionActivityRatio"].sum("MODE_OF_OPERATION") != 0
+        ds["EmissionActivityRatio"].sum("MODE_OF_OPERATION") != 0, drop=False
     )
 
     AnnualTechnologyEmissionPenaltyByEmission = (
         AnnualTechnologyEmission * ds["EmissionsPenalty"]
     ).where(
         ds["EmissionsPenalty"].notnull()
-        & (ds["EmissionActivityRatio"].sum("MODE_OF_OPERATION") != 0)
+        & (ds["EmissionActivityRatio"].sum("MODE_OF_OPERATION") != 0),
+        drop=False,
     )
 
     AnnualTechnologyEmissionsPenalty = AnnualTechnologyEmissionPenaltyByEmission.sum(

tz/osemosys/model/linear_expressions/financials.py

@@ -5,7 +5,6 @@
 
 
 def add_lex_financials(ds: xr.Dataset, m: Model, lex: Dict[str, LinearExpression]):
-
     CapitalInvestment = (
         ds["CapitalCost"].fillna(0)
         * m["NewCapacity"]
@@ -21,7 +20,8 @@ def add_lex_financials(ds: xr.Dataset, m: Model, lex: Dict[str, LinearExpression
         .sum(dims="MODE_OF_OPERATION")
         .where(
             (ds["VariableCost"].sum(dim="MODE_OF_OPERATION") != 0)
-            & (~ds["VariableCost"].sum(dim="MODE_OF_OPERATION").isnull())
+            & (~ds["VariableCost"].sum(dim="MODE_OF_OPERATION").isnull()),
+            drop=False,
         )
     )
     AnnualFixedOperatingCost = lex["GrossCapacity"] * ds["FixedCost"].fillna(0)
@@ -38,8 +38,8 @@ def add_lex_financials(ds: xr.Dataset, m: Model, lex: Dict[str, LinearExpression
     DiscountedCapitalInvestment = CapitalInvestment / lex["DiscountFactor"]
 
     SalvageValue = (
-        m["NewCapacity"] * SV1Cost.where(lex["sv1_mask"])
-        + m["NewCapacity"] * SV2Cost.where(lex["sv2_mask"])
+        m["NewCapacity"] * SV1Cost.where(lex["sv1_mask"], drop=False)
+        + m["NewCapacity"] * SV2Cost.where(lex["sv2_mask"], drop=False)
     ).fillna(0)
 
     DiscountedSalvageValue = SalvageValue / lex["DiscountFactorSalvage"]

tz/osemosys/model/linear_expressions/production.py

@@ -7,21 +7,21 @@
 def add_lex_quantities(ds: xr.Dataset, m: Model, lex: Dict[str, LinearExpression]):
     # Production
     RateOfProductionByTechnologyByMode = m["RateOfActivity"] * ds["OutputActivityRatio"].where(
-        ds["OutputActivityRatio"].notnull()
+        ds["OutputActivityRatio"].notnull(), drop=False
     )
     RateOfProductionByTechnology = RateOfProductionByTechnologyByMode.where(
-        ds["OutputActivityRatio"].sum("MODE_OF_OPERATION") != 0
+        ds["OutputActivityRatio"].sum("MODE_OF_OPERATION") != 0, drop=False
     ).sum(dims="MODE_OF_OPERATION")
     RateOfProduction = RateOfProductionByTechnology.sum(dims="TECHNOLOGY")
     ProductionByTechnology = RateOfProductionByTechnology * ds["YearSplit"]
     Production = RateOfProduction * ds["YearSplit"]
     ProductionAnnual = Production.sum(dims="TIMESLICE")
 
     RateOfUseByTechnologyByMode = m["RateOfActivity"] * ds["InputActivityRatio"].where(
-        ds["InputActivityRatio"].notnull()
+        ds["InputActivityRatio"].notnull(), drop=False
     )
     RateOfUseByTechnology = RateOfUseByTechnologyByMode.where(
-        ds["InputActivityRatio"].sum("MODE_OF_OPERATION") != 0
+        ds["InputActivityRatio"].sum("MODE_OF_OPERATION") != 0, drop=False
     ).sum(dims="MODE_OF_OPERATION")
     RateOfUse = RateOfUseByTechnology.sum(dims="TECHNOLOGY")
     Use = RateOfUse * ds["YearSplit"]

tz/osemosys/model/linear_expressions/re_production.py

@@ -5,12 +5,11 @@
 
 
 def add_lex_re_production(ds: xr.Dataset, m: Model, lex: Dict[str, LinearExpression]):
-
     RateOfProductionByTechnologyByModeRE = m["RateOfActivity"] * ds["OutputActivityRatio"].where(
-        ds["OutputActivityRatio"].notnull() & (ds["RETagTechnology"] == 1)
+        ds["OutputActivityRatio"].notnull() & (ds["RETagTechnology"] == 1), drop=False
     )
     RateOfProductionByTechnologyRE = RateOfProductionByTechnologyByModeRE.where(
-        ds["OutputActivityRatio"].sum("MODE_OF_OPERATION") != 0
+        ds["OutputActivityRatio"].sum("MODE_OF_OPERATION") != 0, drop=False
     ).sum(dims="MODE_OF_OPERATION")
     RateOfProductionRE = RateOfProductionByTechnologyRE.sum(dims="TECHNOLOGY")
     ProductionByTechnologyRE = RateOfProductionByTechnologyRE * ds["YearSplit"]

tz/osemosys/model/linear_expressions/reserve_margin.py

@@ -11,7 +11,8 @@ def add_lex_reserve_margin(ds: xr.Dataset, m: Model, lex: Dict[str, LinearExpres
         ).where(
             (ds["ReserveMargin"] > 0)
             & (ds["ReserveMarginTagTechnology"] == 1)
-            & (ds["ReserveMarginTagTechnology"] * ds["CapacityToActivityUnit"]).notnull()
+            & (ds["ReserveMarginTagTechnology"] * ds["CapacityToActivityUnit"]).notnull(),
+            drop=False,
         )
     ).sum("TECHNOLOGY")
 
@@ -21,28 +22,31 @@ def add_lex_reserve_margin(ds: xr.Dataset, m: Model, lex: Dict[str, LinearExpres
         (ds["OutputActivityRatio"].notnull())
         & (ds["ReserveMargin"] > 0)
         & (ds["ReserveMarginTagFuel"] == 1)
-        & (ds["ReserveMarginTagTechnology"] == 1)
+        & (ds["ReserveMarginTagTechnology"] == 1),
+        drop=False,
     )
 
     RateOfProductionByTechnologyWithReserveMargin = (
         RateOfProductionByTechnologyByModeWithReserveMargin.where(
             (ds["OutputActivityRatio"].notnull())
             & (ds["ReserveMargin"] > 0)
             & (ds["ReserveMarginTagFuel"] == 1)
-            & (ds["ReserveMarginTagTechnology"] == 1)
+            & (ds["ReserveMarginTagTechnology"] == 1),
+            drop=False,
         ).sum(dims="MODE_OF_OPERATION")
     )
 
     RateOfProductionWithReserveMargin = RateOfProductionByTechnologyWithReserveMargin.where(
         (ds["OutputActivityRatio"].notnull())
         & (ds["ReserveMargin"] > 0)
         & (ds["ReserveMarginTagFuel"] == 1)
-        & (ds["ReserveMarginTagTechnology"] == 1)
+        & (ds["ReserveMarginTagTechnology"] == 1),
+        drop=False,
     ).sum(dims="TECHNOLOGY")
 
     DemandNeedingReserveMargin = (
         (lex["RateOfProduction"] * ds["ReserveMarginTagFuel"])
-        .where((ds["ReserveMargin"] > 0) & (ds["ReserveMarginTagFuel"] == 1))
+        .where((ds["ReserveMargin"] > 0) & (ds["ReserveMarginTagFuel"] == 1), drop=False)
         .sum("FUEL")
     )
 

tz/osemosys/model/linear_expressions/storage.py

@@ -5,14 +5,13 @@
 
 
 def add_lex_storage(ds: xr.Dataset, m: Model, lex: Dict[str, LinearExpression]):
-
     DiscountFactorStorage = (1 + ds["DiscountRateStorage"]) ** (
         1 + ds.coords["YEAR"][-1] - ds.coords["YEAR"][0]
     )
 
     RateOfStorageCharge = (
         (ds["TechnologyToStorage"] * m["RateOfActivity"]).where(
-            (ds["TechnologyToStorage"].notnull()) & (ds["TechnologyToStorage"] != 0)
+            (ds["TechnologyToStorage"].notnull()) & (ds["TechnologyToStorage"] != 0), drop=True
         )
     ).sum(["TECHNOLOGY", "MODE_OF_OPERATION"])
 
@@ -29,13 +28,14 @@ def add_lex_storage(ds: xr.Dataset, m: Model, lex: Dict[str, LinearExpression]):
         ).where(
             (ds["TechnologyToStorage"].notnull())
             & (ds["StorageBalanceDay"] != 0)
-            & (ds["Conversionls"] != 0)
+            & (ds["Conversionls"] != 0),
+            drop=False,
         )
     ).sum(["TECHNOLOGY", "MODE_OF_OPERATION", "TIMESLICE"])
 
     RateOfStorageDischarge = (
         (ds["TechnologyFromStorage"] * m["RateOfActivity"]).where(
-            (ds["TechnologyFromStorage"].notnull()) & (ds["TechnologyFromStorage"] != 0)
+            (ds["TechnologyFromStorage"].notnull()) & (ds["TechnologyFromStorage"] != 0), drop=True
         )
     ).sum(["TECHNOLOGY", "MODE_OF_OPERATION"])
 
@@ -52,7 +52,8 @@ def add_lex_storage(ds: xr.Dataset, m: Model, lex: Dict[str, LinearExpression]):
         ).where(
             (ds["TechnologyFromStorage"].notnull())
             & (ds["StorageBalanceDay"] != 0)
-            & (ds["Conversionls"] != 0)
+            & (ds["Conversionls"] != 0),
+            drop=False,
         )
     ).sum(["TECHNOLOGY", "MODE_OF_OPERATION", "TIMESLICE"])
 
@@ -85,8 +86,8 @@ def add_lex_storage(ds: xr.Dataset, m: Model, lex: Dict[str, LinearExpression]):
     )
 
     SalvageValueStorage = (
-        m["NewStorageCapacity"] * SV1CostStorage.where(lex["sv1_mask"])
-        + m["NewStorageCapacity"] * SV2CostStorage.where(lex["sv2_mask"])
+        m["NewStorageCapacity"] * SV1CostStorage.where(lex["sv1_mask"], drop=False)
+        + m["NewStorageCapacity"] * SV2CostStorage.where(lex["sv2_mask"], drop=False)
     ).fillna(0)
 
     DiscountedSalvageValueStorage = SalvageValueStorage / DiscountFactorStorage

tz/osemosys/model/linear_expressions/trade.py

@@ -5,7 +5,6 @@
 
 
 def add_lex_trade(ds: xr.Dataset, m: Model, lex: Dict[str, LinearExpression]):
-
     # Capacity #
     NewTradeCapacity = m["NewTradeCapacity"].rename(YEAR="BUILDYEAR")
     mask = (ds.YEAR - NewTradeCapacity.data.BUILDYEAR >= 0) & (
@@ -17,7 +16,7 @@ def add_lex_trade(ds: xr.Dataset, m: Model, lex: Dict[str, LinearExpression]):
     # Activity #
     NetTrade = (
         ((m["Export"] / (1 - ds["TradeLossBetweenRegions"])) - m["Import"])
-        .where(ds["TradeRoute"].notnull())
+        .where(ds["TradeRoute"].notnull(), drop=False)
         .sum("_REGION")
         .fillna(0)
     )
@@ -87,8 +86,8 @@ def add_lex_trade(ds: xr.Dataset, m: Model, lex: Dict[str, LinearExpression]):
 
     # salvage value (trade)
     SalvageValueTrade = (
-        m["NewTradeCapacity"] * SV1CostTrade.where(sv1_trade_mask)
-        + m["NewTradeCapacity"] * SV2CostTrade.where(sv2_trade_mask)
+        m["NewTradeCapacity"] * SV1CostTrade.where(sv1_trade_mask, drop=False)
+        + m["NewTradeCapacity"] * SV2CostTrade.where(sv2_trade_mask, drop=False)
     ).fillna(0)
 
     DiscountedSalvageValueTrade = SalvageValueTrade / DiscountFactorSalvageTrade