Add commodity balance epsilon to force candidate assets to activate

src/simulation/optimisation.rs

@@ -732,6 +732,7 @@ impl<'model, 'run> DispatchRun<'model, 'run> {
             &all_assets,
             markets_to_balance,
             self.year,
+            self.candidate_assets,
         );
 
         // Create model and apply any user-supplied HiGHS options to it

src/simulation/optimisation/constraints.rs

@@ -5,10 +5,13 @@ use crate::commodity::{CommodityID, CommodityType};
 use crate::model::Model;
 use crate::region::RegionID;
 use crate::time_slice::{TimeSliceInfo, TimeSliceSelection};
-use crate::units::UnitType;
+use crate::units::{Activity, UnitType};
 use highs::RowProblem as Problem;
 use indexmap::IndexMap;
 
+/// Small epsilon to add to the commodity balance constraints for candidate assets
+const COMMODITY_BALANCE_EPSILON_FOR_CANDIDATES: f64 = 1e-6;
+
 /// Corresponding variables for a constraint along with the row offset in the solution
 pub struct KeysWithOffset<T> {
     /// Row offset in the solver's row ordering corresponding to the first key in `keys`.
@@ -80,6 +83,7 @@ pub fn add_model_constraints<'a, I>(
     assets: &I,
     markets_to_balance: &'a [(CommodityID, RegionID)],
     year: u32,
+    candidate_assets: &'a [AssetRef],
 ) -> ConstraintKeys
 where
     I: Iterator<Item = &'a AssetRef> + Clone + 'a,
@@ -91,6 +95,7 @@ where
         assets,
         markets_to_balance,
         year,
+        candidate_assets,
     );
 
     let activity_keys =
@@ -121,6 +126,7 @@ fn add_commodity_balance_constraints<'a, I>(
     assets: &I,
     markets_to_balance: &'a [(CommodityID, RegionID)],
     year: u32,
+    candidate_assets: &'a [AssetRef],
 ) -> CommodityBalanceKeys
 where
     I: Iterator<Item = &'a AssetRef> + Clone + 'a,
@@ -173,13 +179,30 @@ where
                 }
             }
 
-            // For SED commodities, the LHS must be >=0 and for SVD commodities, it must be >=
-            // the exogenous demand supplied by the user
-            let min = if commodity.kind == CommodityType::ServiceDemand {
-                commodity.demand[&(region_id.clone(), year, ts_selection.clone())].value()
+            // If any candidate asset in this region produces this commodity, and has a
+            // nonzero activity limit for this time slice selection, then we add a small epsilon
+            // to the *lower bound* of the balance constraints to force some dispatch.
+            let epsilon = if candidate_assets.iter().any(|a| {
+                a.region_id() == region_id
+                    && a.iter_output_flows()
+                        .any(|flow| &flow.commodity.id == commodity_id)
+                    && a.get_activity_limits_for_selection(&ts_selection).end() > &Activity(0.0)
+            }) {
+                COMMODITY_BALANCE_EPSILON_FOR_CANDIDATES
             } else {
                 0.0
             };
+
+            // For SED commodities, enforce net production >= epsilon, and for SVD commodities
+            // enforce net production >= exogenous demand (+ epsilon).
+            let min = match commodity.kind {
+                CommodityType::ServiceDemand => {
+                    commodity.demand[&(region_id.clone(), year, ts_selection.clone())].value()
+                        + epsilon
+                }
+                _ => epsilon,
+            };
+
             // Consume collected terms into a row. `terms.drain(..)` ensures the vector is
             // emptied for the next selection.
             problem.add_row(min.., terms.drain(..));

tests/data/circularity/asset_capacities.csv

@@ -14,6 +14,7 @@ milestone_year,asset_id,group_id,capacity,num_units
 2020,12,,15.0,
 2020,13,,2900.0,
 2020,14,,399.98,
+2030,0,,4010.2,
 2030,1,,1738.05,
 2030,2,,3789.64,
 2030,5,,3.964844,
@@ -24,20 +25,19 @@ milestone_year,asset_id,group_id,capacity,num_units
 2030,16,,1827.9092029382402,
 2030,17,,4076.5105773701275,
 2030,18,,1643.7237844403953,
-2030,19,,558.9731652805061,
-2030,20,,2649.731875290583,
-2030,21,,2761.4559689056887,
+2030,19,,2649.731875290583,
+2030,20,,1625.1150136624153,
 2040,1,,1738.05,
 2040,5,,3.964844,
 2040,16,,1827.9092029382402,
 2040,17,,4076.5105773701275,
 2040,18,,1643.7237844403953,
-2040,20,,2649.731875290583,
-2040,21,,2761.4559689056887,
-2040,22,,912.8939641298439,
-2040,23,,2005.8286877382402,
-2040,24,,177.9194848000003,
-2040,25,,0.26061763467850396,
-2040,26,,1.3918603724966105,
-2040,27,,1905.1678522886418,
-2040,28,,3111.954097466583,
+2040,19,,2649.731875290583,
+2040,20,,1625.1150136624153,
+2040,21,,912.8939641298439,
+2040,22,,2005.8286877382402,
+2040,23,,177.9194848000003,
+2040,24,,1.3918603724966105,
+2040,25,,1905.1678522886418,
+2040,26,,3111.954097466583,
+2040,27,,2000.426244903074,

tests/data/circularity/assets.csv

@@ -18,13 +18,12 @@ asset_id,group_id,process_id,region_id,agent_id,commission_year
 16,,RBIOBL,GBR,A0_RES,2030
 17,,OILREF,GBR,A0_REF,2030
 18,,BIOPLL,GBR,A0_BPL,2030
-19,,GASDRV,GBR,A0_OAG,2030
-20,,OAGRSV,GBR,A0_OAG,2030
-21,,BIOPRO,GBR,A0_BPD,2030
-22,,TDIECR,GBR,A0_TRA,2040
-23,,RBIOBL,GBR,A0_RES,2040
-24,,RELCHP,GBR,A0_RES,2040
-25,,H2YGEN,GBR,A0_ELC,2040
-26,,WNDFRM,GBR,A0_ELC,2040
-27,,BIOPLL,GBR,A0_BPL,2040
-28,,OAGRSV,GBR,A0_OAG,2040
+19,,OAGRSV,GBR,A0_OAG,2030
+20,,BIOPRO,GBR,A0_BPD,2030
+21,,TDIECR,GBR,A0_TRA,2040
+22,,RBIOBL,GBR,A0_RES,2040
+23,,RELCHP,GBR,A0_RES,2040
+24,,WNDFRM,GBR,A0_ELC,2040
+25,,BIOPLL,GBR,A0_BPL,2040
+26,,OAGRSV,GBR,A0_OAG,2040
+27,,BIOPRO,GBR,A0_BPD,2040