Use integer variables in appraisal for divisible assets

src/asset.rs

@@ -831,11 +831,42 @@ impl Asset {
         *mothballed_year
     }
 
+    /// Get the unit size for this asset's process (if any)
+    pub fn unit_size(&self) -> Option<Capacity> {
+        self.process.unit_size
+    }
+
     /// Checks if the asset corresponds to a process that has a `unit_size` and is therefore divisible.
     pub fn is_divisible(&self) -> bool {
         self.process.unit_size.is_some()
     }
 
+    /// Convert a capacity to number of units for a divisible asset.
+    ///
+    /// Divides the given capacity by the process unit size and rounds up. In other words, this is
+    /// the minimum number of units required to achieve at least the given capacity.
+    ///
+    /// Panics if the asset is not divisible.
+    pub fn capacity_to_units(&self, capacity: Capacity) -> u32 {
+        let unit_size = self.unit_size().expect("Asset must be divisible");
+        #[allow(clippy::cast_possible_truncation, clippy::cast_sign_loss)]
+        {
+            (capacity / unit_size).value().ceil() as u32
+        }
+    }
+
+    /// Round a capacity up to the nearest multiple of the unit size.
+    ///
+    /// For a divisible asset, returns the minimum capacity (as a multiple of `unit_size`)
+    /// that is at least as large as the given capacity.
+    ///
+    /// Panics if the asset is not divisible.
+    pub fn round_capacity_to_unit_size(&self, capacity: Capacity) -> Capacity {
+        let unit_size = self.unit_size().expect("Asset must be divisible");
+        let n_units = self.capacity_to_units(capacity);
+        Capacity(unit_size.value() * n_units as f64)
+    }
+
     /// Divides an asset if it is divisible and returns a vector of children
     ///
     /// The child assets are identical to the parent (including state) but with a capacity
@@ -860,10 +891,7 @@ impl Asset {
         );
 
         // Calculate the number of units corresponding to the asset's capacity
-        // Safe because capacity and unit_size are both positive finite numbers, so their ratio
-        // must also be positive and finite.
-        #[allow(clippy::cast_possible_truncation, clippy::cast_sign_loss)]
-        let n_units = (self.capacity / unit_size).value().ceil() as usize;
+        let n_units = self.capacity_to_units(self.capacity) as usize;
 
         // Divide the asset into `n_units` children of size `unit_size`
         let child_asset = Self {

src/simulation/investment.rs

@@ -617,12 +617,16 @@ fn get_candidate_assets<'a>(
             let mut asset =
                 Asset::new_candidate(process.clone(), region_id.clone(), Capacity(0.0), year)
                     .unwrap();
-            asset.set_capacity(get_demand_limiting_capacity(
-                time_slice_info,
-                &asset,
-                commodity,
-                demand,
-            ));
+
+            // Set capacity based on demand
+            // This will serve as the upper limit when appraising the asset
+            // If the asset is divisible, round capacity to the nearest multiple of the unit size
+            let mut capacity =
+                get_demand_limiting_capacity(time_slice_info, &asset, commodity, demand);
+            if asset.is_divisible() {
+                capacity = asset.round_capacity_to_unit_size(capacity);
+            }
+            asset.set_capacity(capacity);
 
             asset.into()
         })
@@ -708,7 +712,13 @@ fn select_best_assets(
         let mut outputs_for_opts = Vec::new();
         for asset in &opt_assets {
             let max_capacity = (!asset.is_commissioned()).then(|| {
-                let max_capacity = model.parameters.capacity_limit_factor * asset.capacity();
+                let mut max_capacity = model.parameters.capacity_limit_factor * asset.capacity();
+
+                // For divisible assets, round up to the nearest multiple of the process unit size
+                if asset.is_divisible() {
+                    max_capacity = asset.round_capacity_to_unit_size(max_capacity);
+                }
+
                 let remaining_capacity = remaining_candidate_capacity[asset];
                 max_capacity.min(remaining_capacity)
             });

src/simulation/investment/appraisal/constraints.rs

@@ -5,6 +5,7 @@ use crate::asset::{AssetRef, AssetState};
 use crate::commodity::Commodity;
 use crate::time_slice::{TimeSliceID, TimeSliceInfo};
 use crate::units::{Capacity, Flow};
+use float_cmp::approx_eq;
 use highs::RowProblem as Problem;
 use indexmap::IndexMap;
 
@@ -19,15 +20,26 @@ pub fn add_capacity_constraint(
     max_capacity: Option<Capacity>,
     capacity_var: Variable,
 ) {
-    let capacity = max_capacity.unwrap_or(asset.capacity());
+    let mut capacity_limit = max_capacity.unwrap_or(asset.capacity()).value();
+
+    // If asset is divisible, capacity_var represents number of units, so we must divide the
+    // capacity bounds by the unit size.
+    if let Some(unit_size) = asset.unit_size() {
+        capacity_limit /= unit_size.value();
+
+        // Sanity check: capacity_limit should be a whole number of units (i.e pre-adjusted
+        // capacity limit was a multiple of unit size)
+        assert!(approx_eq!(f64, capacity_limit, capacity_limit.round()));
+    }
+
     let bounds = match asset.state() {
         AssetState::Commissioned { .. } => {
             // Fixed capacity for commissioned assets
-            capacity.value()..=capacity.value()
+            capacity_limit..=capacity_limit
         }
         AssetState::Candidate => {
             // Variable capacity between 0 and max for candidate assets
-            0.0..=capacity.value()
+            0.0..=capacity_limit
         }
         _ => panic!(
             "add_capacity_constraint should only be called with Commissioned or Candidate assets"
@@ -100,8 +112,15 @@ fn add_activity_constraints_for_candidate(
     time_slice_info: &TimeSliceInfo,
 ) {
     for (ts_selection, limits) in asset.iter_activity_per_capacity_limits() {
-        let upper_limit = limits.end().value();
-        let lower_limit = limits.start().value();
+        let mut upper_limit = limits.end().value();
+        let mut lower_limit = limits.start().value();
+
+        // If the asset is divisible, the capacity variable represents number of units,
+        // so we need to multiply the per-capacity limits by the unit size.
+        if let Some(unit_size) = asset.unit_size() {
+            upper_limit *= unit_size.value();
+            lower_limit *= unit_size.value();
+        }
 
         // Collect capacity and activity terms
         // We have a single capacity term, and activity terms for all time slices in the selection

src/simulation/investment/appraisal/optimisation.rs

@@ -21,7 +21,10 @@ pub type Variable = highs::Col;
 
 /// Map storing variables for the optimisation problem
 struct VariableMap {
-    /// Capacity variable
+    /// Capacity variable.
+    ///
+    /// This represents absolute capacity for indivisible assets and number of units for
+    /// divisible assets.
     capacity_var: Variable,
     /// Activity variables in each time slice
     activity_vars: IndexMap<TimeSliceID, Variable>,
@@ -38,10 +41,23 @@ impl VariableMap {
     ///
     /// # Returns
     /// A new `VariableMap` containing all created decision variables
-    fn add_to_problem(problem: &mut Problem, cost_coefficients: &ObjectiveCoefficients) -> Self {
+    fn add_to_problem(
+        problem: &mut Problem,
+        cost_coefficients: &ObjectiveCoefficients,
+        capacity_unit_size: Option<Capacity>,
+    ) -> Self {
         // Create capacity variable with its associated cost
-        let capacity_var =
-            problem.add_column(cost_coefficients.capacity_coefficient.value(), 0.0..);
+        let capacity_coefficient = cost_coefficients.capacity_coefficient.value();
+        let capacity_var = match capacity_unit_size {
+            Some(unit_size) => {
+                // Divisible asset: capacity variable represents number of units
+                problem.add_integer_column(capacity_coefficient * unit_size.value(), 0.0..)
+            }
+            None => {
+                // Indivisible asset: capacity variable represents total capacity
+                problem.add_column(capacity_coefficient, 0.0..)
+            }
+        };
 
         // Create activity variables for each time slice
         let mut activity_vars = IndexMap::new();
@@ -118,7 +134,7 @@ pub fn perform_optimisation(
 ) -> Result<ResultsMap> {
     // Create problem and add variables
     let mut problem = Problem::default();
-    let variables = VariableMap::add_to_problem(&mut problem, coefficients);
+    let variables = VariableMap::add_to_problem(&mut problem, coefficients, asset.unit_size());
 
     // Add constraints
     add_constraints(
@@ -135,7 +151,12 @@ pub fn perform_optimisation(
     let solution = solve_optimal(problem.optimise(sense))?.get_solution();
     let solution_values = solution.columns();
     Ok(ResultsMap {
-        capacity: Capacity::new(solution_values[0]),
+        // If the asset is divisible, the capacity variable represents number of units, so convert
+        // to total capacity
+        capacity: match asset.unit_size() {
+            Some(unit_size) => Capacity::new(solution_values[0] * unit_size.value()),
+            None => Capacity::new(solution_values[0]),
+        },
         // The mapping below assumes the column ordering documented on `VariableMap::add_to_problem`:
         // index 0 = capacity, next `n` entries = activities (in the same key order as
         // `cost_coefficients.activity_coefficients`), remaining entries = unmet demand.

tests/data/simple_divisible/assets.csv

@@ -1,15 +1,15 @@
 asset_id,process_id,region_id,agent_id,group_id,commission_year,decommission_year,capacity
 0,GASDRV,GBR,A0_GEX,,2020,,4002.26
 1,GASPRC,GBR,A0_GPR,,2020,,3782.13
-2,WNDFRM,GBR,A0_ELC,,2020,2040,3.964844
-3,GASCGT,GBR,A0_ELC,,2020,2040,2.43
+2,WNDFRM,GBR,A0_ELC,,2020,,3.964844
+3,GASCGT,GBR,A0_ELC,,2020,,2.43
 4,RGASBR,GBR,A0_RES,,2020,2035,1000.0
 5,RGASBR,GBR,A0_RES,,2020,2035,1000.0
 6,RGASBR,GBR,A0_RES,,2020,2035,1000.0
 7,RELCHP,GBR,A0_RES,,2020,2035,399.98
-8,RGASBR,GBR,A0_RES,1,2030,,1000.0
-9,GASCGT,GBR,A0_ELC,,2030,2040,0.44245235762867363
-10,RGASBR,GBR,A0_RES,2,2040,,1000.0
-11,RGASBR,GBR,A0_RES,2,2040,,1000.0
-12,RGASBR,GBR,A0_RES,2,2040,,1000.0
-13,RGASBR,GBR,A0_RES,2,2040,,1000.0
+8,RELCHP,GBR,A0_RES,,2030,,255.83840587648046
+9,GASCGT,GBR,A0_ELC,,2030,2040,3.1192651014219064
+10,RGASBR,GBR,A0_RES,1,2040,,1000.0
+11,RGASBR,GBR,A0_RES,1,2040,,1000.0
+12,RGASBR,GBR,A0_RES,1,2040,,1000.0
+13,RGASBR,GBR,A0_RES,1,2040,,1000.0