Allow passing custom options to HiGHS via model.toml

clippy.toml

@@ -4,3 +4,5 @@
 # The AssetRef type uses Rc<Asset> internally for shared ownership, but the hash
 # implementation is carefully designed to be consistent regardless of interior mutability
 ignore-interior-mutability = ["muse2::asset::AssetRef"]
+# Things that should not be treated as identifiers. The ".." represents the default options.
+doc-valid-idents = ["HiGHS", ".."]

docs/SUMMARY.md

@@ -18,6 +18,7 @@
   - [Setting up your development environment](developer_guide/setup.md)
   - [Building and developing MUSE2](developer_guide/coding.md)
   - [Architecture and coding style](developer_guide/architecture_quickstart.md)
+  - [Applying custom HiGHS options](developer_guide/custom_highs_options.md)
   - [Developing the documentation](developer_guide/docs.md)
 - [API documentation](./api/muse2/README.md)
 - [Release notes](release_notes/README.md)

docs/developer_guide/custom_highs_options.md

@@ -0,0 +1,36 @@
+# Applying custom HiGHS options
+
+As part of development, you may wish to directly set custom options for the HiGHS solver. Note that
+while some of these options will not affect results of simulations (e.g. to enable console logging
+for HiGHS), as we cannot guarantee this for all options, in order to use this feature, you have to
+set `please_give_me_broken_results = true` in your [`model.toml` file][model.toml].
+
+You can change any of the options exposed by the HiGHS solver; for more information, see [the HiGHS
+documentation][highs-opts].
+
+You can set options to be applied to all optimisations, just dispatch or just appraisal.
+
+Here is an example:
+
+```toml
+please_give_me_broken_results = true
+milestone_years = [2020, 2030, 2040]
+
+# These options are applied to all optimisations
+[highs.global_options]
+# These two options are required to be enabled to log to console
+log_to_console = true
+output_flag = true
+
+# These ones are just applied to dispatch
+[highs.dispatch_options]
+# Increase to higher than default
+primal_feasibility_tolerance = 10e-6
+
+# These ones are just applied to appraisal
+[highs.appraisal_options]
+optimality_tolerance = 10e-6
+```
+
+[model.toml]: https://energysystemsmodellinglab.github.io/MUSE2/file_formats/input_files.html#model-parameters-modeltoml
+[highs-opts]: https://ergo-code.github.io/HiGHS/stable/options/definitions/

docs/release_notes/upcoming.md

@@ -16,7 +16,7 @@ ready to be released, carry out the following steps:
 
 ## New features
 
-<!-- TODO -->
+- Users can now optionally pass [custom options][highs-options] to the HiGHS solver [#1276]
 
 ## Breaking changes
 
@@ -26,4 +26,6 @@ ready to be released, carry out the following steps:
 
 - Fix misleading warning message for assets decommissioned before simulation start ([#1259])
 
+[highs-options]: https://ergo-code.github.io/HiGHS/stable/options/definitions/
 [#1259]: https://github.com/EnergySystemsModellingLab/MUSE2/pull/1259
+[#1276]: https://github.com/EnergySystemsModellingLab/MUSE2/pull/1276

schemas/input/model.yaml

@@ -59,5 +59,26 @@ properties:
       investment cycle. Changing the value of this parameter is potentially dangerous,
       so it requires setting `please_give_me_broken_results` to true.
     default: 1e-12
+  highs:
+    type: object
+    description: |
+      Used for setting custom HiGHS options. As this is unsafe, it requires setting
+      `please_give_me_broken_results` to true. For more information, see the relevant [section of
+      the developer documentation][highs-opts-docs].
+
+      [highs-opts-docs]: https://energysystemsmodellinglab.github.io/MUSE2/developer_guide/custom_highs_options.html
+    properties:
+      global_options:
+        type: object
+        description: HiGHS options applied to all optimisations
+        properties: {}
+      dispatch_options:
+        type: object
+        description: HiGHS options applied only to dispatch optimisations
+        properties: {}
+      appraisal_options:
+        type: object
+        description: HiGHS options applied only to appraisal optimisations
+        properties: {}
 
 required: [milestone_years]

src/model/parameters.rs

@@ -13,6 +13,7 @@ use log::warn;
 use serde::Deserialize;
 use std::path::Path;
 use std::sync::OnceLock;
+use toml::Table;
 
 const MODEL_PARAMETERS_FILE_NAME: &str = "model.toml";
 
@@ -63,7 +64,7 @@ fn set_dangerous_model_options_flag(enabled: bool) {
 ///
 /// NOTE: If you add or change a field in this struct, you must also update the schema in
 /// `schemas/input/model.yaml`.
-#[derive(Debug, Deserialize, PartialEq)]
+#[derive(Deserialize)]
 #[serde(default)]
 pub struct ModelParameters {
     /// Milestone years
@@ -98,6 +99,12 @@ pub struct ModelParameters {
     pub mothball_years: u32,
     /// Absolute tolerance when checking if remaining demand is close enough to zero
     pub remaining_demand_absolute_tolerance: Flow,
+    /// Options for the HiGHS solver.
+    ///
+    /// For a full list of options, see [the HiGHS documentation].
+    ///
+    /// [the HiGHS documentation]: https://ergo-code.github.io/HiGHS/stable/options/definitions/
+    pub highs: HighsOptions,
 }
 
 impl Default for ModelParameters {
@@ -117,10 +124,47 @@ impl Default for ModelParameters {
             capacity_margin: Dimensionless(0.2),
             mothball_years: 0,
             remaining_demand_absolute_tolerance: DEFAULT_REMAINING_DEMAND_ABSOLUTE_TOLERANCE,
+            highs: HighsOptions::default(),
         }
     }
 }
 
+/// Defines the TOML table holding the sub-tables to define HiGHS options
+#[derive(Default, Deserialize)]
+#[serde(default)]
+pub struct HighsOptions {
+    /// HiGHS options applied to both dispatch and appraisal
+    global_options: Table,
+    /// HiGHS options applied only to dispatch
+    pub dispatch_options: Table,
+    /// HiGHS options applied only to appraisal
+    pub appraisal_options: Table,
+}
+
+impl HighsOptions {
+    /// Copy the options specified in `global_options` to `dispatch_options` and `appraisal_options`
+    fn apply_global_options(&mut self) {
+        let append_global_opts_to = |opts: &mut Table| {
+            for (option, value) in &self.global_options {
+                opts.entry(option).or_insert(value.clone());
+            }
+        };
+
+        append_global_opts_to(&mut self.dispatch_options);
+        append_global_opts_to(&mut self.appraisal_options);
+
+        // Now that we're finished with global options we can discard
+        self.global_options.clear();
+    }
+
+    /// Check whether any options have been set
+    pub fn is_empty(&self) -> bool {
+        self.global_options.is_empty()
+            && self.dispatch_options.is_empty()
+            && self.appraisal_options.is_empty()
+    }
+}
+
 /// Check that the `milestone_years` parameter is valid
 fn check_milestone_years(years: &[u32]) -> Result<()> {
     ensure!(
@@ -195,6 +239,20 @@ fn check_capacity_margin(value: Dimensionless) -> Result<()> {
     Ok(())
 }
 
+/// Check the custom HiGHS options are valid.
+///
+/// Note that we cannot know whether the options specified exist and are of the correct type until
+/// we attempt to use them. We could check for types that are never valid (e.g. an array), but as
+/// we're checking later anyway, we don't bother.
+fn check_highs_options(dangerous_options_enabled: bool, highs: &HighsOptions) -> Result<()> {
+    ensure!(
+        dangerous_options_enabled || highs.is_empty(),
+        "Cannot set custom HiGHS options without enabling {ALLOW_DANGEROUS_OPTION_NAME}"
+    );
+
+    Ok(())
+}
+
 impl ModelParameters {
     /// Read a model file from the specified directory.
     ///
@@ -207,14 +265,17 @@ impl ModelParameters {
     /// The model file contents as a [`ModelParameters`] struct or an error if the file is invalid
     pub fn from_path<P: AsRef<Path>>(model_dir: P) -> Result<ModelParameters> {
         let file_path = model_dir.as_ref().join(MODEL_PARAMETERS_FILE_NAME);
-        let model_params: ModelParameters = read_toml(&file_path)?;
+        let mut model_params: ModelParameters = read_toml(&file_path)?;
 
         set_dangerous_model_options_flag(model_params.allow_dangerous_options);
 
         model_params
             .validate()
             .with_context(|| input_err_msg(file_path))?;
 
+        // Copy global options to other tables
+        model_params.highs.apply_global_options();
+
         Ok(model_params)
     }
 
@@ -256,6 +317,8 @@ impl ModelParameters {
             self.remaining_demand_absolute_tolerance,
         )?;
 
+        check_highs_options(self.allow_dangerous_options, &self.highs)?;
+
         Ok(())
     }
 }

src/simulation/investment/appraisal.rs

@@ -258,12 +258,12 @@ fn calculate_lcox(
     demand: &DemandMap,
 ) -> Result<AppraisalOutput> {
     let results = perform_optimisation(
+        model,
         asset,
         max_capacity,
         commodity,
         coefficients,
         demand,
-        &model.time_slice_info,
         highs::Sense::Minimise,
     )?;
 
@@ -296,12 +296,12 @@ fn calculate_npv(
     demand: &DemandMap,
 ) -> Result<AppraisalOutput> {
     let results = perform_optimisation(
+        model,
         asset,
         max_capacity,
         commodity,
         coefficients,
         demand,
-        &model.time_slice_info,
         highs::Sense::Maximise,
     )?;
 

src/simulation/investment/appraisal/optimisation.rs

@@ -6,10 +6,13 @@ use super::constraints::{
 };
 use crate::asset::{AssetCapacity, AssetRef};
 use crate::commodity::Commodity;
+use crate::model::Model;
+use crate::simulation::optimisation::ModelError;
+use crate::simulation::optimisation::apply_highs_options_from_toml;
 use crate::simulation::optimisation::solve_optimal;
 use crate::time_slice::{TimeSliceID, TimeSliceInfo};
 use crate::units::{Activity, Capacity, Flow};
-use anyhow::Result;
+use anyhow::{Context, Result};
 use highs::{RowProblem as Problem, Sense};
 use indexmap::IndexMap;
 
@@ -126,12 +129,12 @@ fn add_constraints(
 /// The optimisation will use continuous or integer capacity variables depending on whether the
 /// asset has a defined unit size.
 pub fn perform_optimisation(
+    model: &Model,
     asset: &AssetRef,
     max_capacity: Option<AssetCapacity>,
     commodity: &Commodity,
     coefficients: &ObjectiveCoefficients,
     demand: &DemandMap,
-    time_slice_info: &TimeSliceInfo,
     sense: Sense,
 ) -> Result<ResultsMap> {
     // Create problem and add variables
@@ -146,11 +149,16 @@ pub fn perform_optimisation(
         commodity,
         &variables,
         demand,
-        time_slice_info,
+        &model.time_slice_info,
     );
 
     // Solve model
-    let solution = solve_optimal(problem.optimise(sense))?.get_solution();
+    let mut highs_model = problem.optimise(sense);
+    apply_highs_options_from_toml(&mut highs_model, &model.parameters.highs.appraisal_options)
+        .context("Failed to apply custom HiGHS options to appraisal optimisation")?;
+    let solution = solve_optimal(highs_model)
+        .map_err(ModelError::into_anyhow)?
+        .get_solution();
     let solution_values = solution.columns();
     Ok(ResultsMap {
         // If the asset has a defined unit size, the capacity variable represents number of units,