Feature/manual train planner

altrios-core/src/consist/locomotive/powertrain/fuel_converter.rs

@@ -80,8 +80,8 @@ pub struct FuelConverter {
     pub(crate) pwr_for_peak_eff: si::Power,
     /// idle fuel power to overcome internal friction (not including aux load)
     pub pwr_idle_fuel: si::Power,
-    /// Interpolator for derating dynamic engine peak power based on altitude
-    /// and temperature. When interpolating, this returns fraction of normal
+    /// Interpolator for derating dynamic engine peak power based on altitude [m]
+    /// and temperature [*C]. When interpolating, this returns fraction of normal
     /// peak power, e.g. a value of 1 means no derating and a value of 0 means
     /// the engine is completely disabled.
     #[api(skip_get, skip_set)]
@@ -190,10 +190,12 @@ impl FuelConverter {
 
         let pwr_max_derated = match (&mut self.elev_and_temp_derate, elev_and_temp) {
             (Some(elev_and_temp_derate), Some(elev_and_temp)) => {
-                elev_and_temp_derate.interpolate(&[
+                let derate_factor = elev_and_temp_derate.interpolate(&[
                     elev_and_temp.0.get::<si::meter>(),
                     elev_and_temp.1.get::<si::degree_celsius>(),
-                ])? * self.pwr_out_max
+                ])?;
+                ensure!(derate_factor <= 1.0 && derate_factor >= 0.0, format_dbg!());
+                derate_factor * self.pwr_out_max
             }
             (None, Some(_)) => bail!(
                 "{}\nExpected (self.elev_and_temp_derate, elev_and_temp) to both be Some or None",
@@ -209,7 +211,6 @@ impl FuelConverter {
         self.pwr_out_max_init = self.pwr_out_max_init.max(self.pwr_out_max / 10.);
         self.state.pwr_out_max = (self.state.pwr_mech_out
             + (self.pwr_out_max / self.pwr_ramp_lag) * dt)
-            .min(self.pwr_out_max)
             .min(pwr_max_derated)
             .max(self.pwr_out_max_init);
         Ok(())
@@ -319,8 +320,9 @@ impl FuelConverter {
     fn set_default_elev_and_temp_derate(&mut self) {
         self.elev_and_temp_derate = Some(
             Interp2D::new(
-                array![0.0, 3_000.0, 6_000.0],
-                array![0.0, 35.0, 45.0, 50.0],
+                array![0.0, 1_000.0, 2_000.0], // elevation in meters
+                array![0.0, 35.0, 45.0, 50.0], // temperature in degrees Celsius
+                // Fraction of static peak power
                 array![
                     [1.0, 1.0, 0.95, 0.8],
                     [0.95, 0.95, 0.9025, 0.76],

python/altrios/demos/manual_train_planner_demo.py

@@ -0,0 +1,201 @@
+# %%
+from altrios import sim_manager_manual
+from altrios import utilities, defaults
+import altrios as alt
+from altrios.train_planner import planner_config, manual_train_planner
+import numpy as np
+import polars as pl
+import matplotlib.pyplot as plt
+import time
+import seaborn as sns
+from pathlib import Path
+
+sns.set_theme()
+
+t0_total = time.perf_counter()
+
+SHOW_PLOTS = alt.utils.show_plots()
+# %
+
+plot_dir = Path() / "plots"
+# make the dir if it doesn't exist
+plot_dir.mkdir(exist_ok=True)
+#%%
+
+consist_plan = pl.read_csv(alt.resources_root() / "demo_data/consist plan for manual planner demo.csv")
+
+location_map = alt.import_locations(
+    alt.resources_root() / "networks/default_locations.csv"
+)
+network = alt.Network.from_file(
+    alt.resources_root() / "networks/Taconite-NoBalloon.yaml"
+)
+
+loco_map = {'Diesel_Large' : 'Diesel_Large',
+            'BEL' : 'BEL'}
+
+rail_vehicles = [
+    alt.RailVehicle.from_file(vehicle_file, skip_init=False)
+    for vehicle_file in Path(alt.resources_root() / "rolling_stock/").glob("*.yaml")
+]
+
+t0_main = time.perf_counter()
+#not passing in a trainplanner config here because we do not need to specify train length or any other parameters like
+#sim_manager_demo.py.  This example is just replaying trains that have already been planned.
+(
+    train_consist_plan,
+    loco_pool,
+    refuel_facilities,
+    grid_emissions_factors,
+    nodal_energy_prices,
+    speed_limit_train_sims,
+    timed_paths,
+    train_consist_plan_untrimmed,
+) = sim_manager_manual.main(
+    network=network,
+    rail_vehicles=rail_vehicles, #double check this to see if it is actually need in sim_manager_manual.py
+    location_map=location_map,
+    consist_plan= consist_plan,
+    loco_map=loco_map,
+    debug=True,
+)
+
+
+# train_consist_plan, loco_pool, refuelers, speed_limit_train_sims, est_time_nets = (
+#     manual_train_planner(consist_plan, loco_map)
+# )
+
+#%%
+t1_main = time.perf_counter()
+print(f"Elapsed time to run `sim_manager.main()`: {t1_main-t0_main:.3g} s")
+
+# %%
+t0_train_sims = time.perf_counter()
+speed_limit_train_sims.set_save_interval(100)
+(sims, refuel_sessions) = alt.run_speed_limit_train_sims(
+    speed_limit_train_sims=speed_limit_train_sims,
+    network=network,
+    train_consist_plan_py=train_consist_plan,
+    loco_pool_py=loco_pool,
+    refuel_facilities_py=refuel_facilities,
+    timed_paths=timed_paths,
+)
+t1_train_sims = time.perf_counter()
+print(f"Elapsed time to run train sims: {t1_train_sims-t0_train_sims:.3g} s")
+t_train_time = sum([sim.state.time_seconds for sim in sims.tolist()])
+print(f"Total train-seconds simulated: {t_train_time} s")
+
+# %%
+t0_summary_sims = time.perf_counter()
+speed_limit_train_sims.set_save_interval(None)
+(summary_sims, summary_refuel_sessions) = alt.run_speed_limit_train_sims(
+    speed_limit_train_sims=speed_limit_train_sims,
+    network=network,
+    train_consist_plan_py=train_consist_plan,
+    loco_pool_py=loco_pool,
+    refuel_facilities_py=refuel_facilities,
+    timed_paths=timed_paths,
+)
+t1_summary_sims = time.perf_counter()
+print(
+    f"Elapsed time to build and run summary sims: {t1_summary_sims-t0_summary_sims:.3g} s"
+)
+
+# %%
+t0_tolist = time.perf_counter()
+sims_list = sims.tolist()
+t1_tolist = time.perf_counter()
+print(f"Elapsed time to run `tolist()`: {t1_tolist-t0_tolist:.3g} s")
+
+sim0 = sims_list[0]
+# sim0 = alt.SpeedLimitTrainSim.from_bincode(
+#     sim0.to_bincode())  # to support linting
+
+
+# %%
+
+t0_main = time.perf_counter()
+e_total_fuel_mj = summary_sims.get_energy_fuel_joules(annualize=False) / 1e9
+t1_main = time.perf_counter()
+
+print(f"Elapsed time to get total fuel energy: {t1_main-t0_main:.3g} s")
+print(f"Total fuel energy used: {e_total_fuel_mj:.3g} GJ")
+
+v_total_fuel_gal = (
+    summary_sims.get_energy_fuel_joules(annualize=False)
+    / 1e3
+    / defaults.LHV_DIESEL_KJ_PER_KG
+    / defaults.RHO_DIESEL_KG_PER_M3
+    * utilities.LITER_PER_M3
+    * utilities.GALLONS_PER_LITER
+)
+
+print(f"Total fuel used: {v_total_fuel_gal:.3g} gallons")
+print(f"Total elapsed time: {time.perf_counter() - t0_total} s")
+
+
+# %%
+
+if SHOW_PLOTS:
+    for idx, sim in enumerate(sims_list[:10]):
+        sim: alt.SpeedLimitTrainSim
+        fig, ax = plt.subplots(3, 1, sharex=True)
+
+        loco0 = sim.loco_con.loco_vec.tolist()[0]
+        # loco0 = alt.Locomotive.from_bincode(
+        #     loco0.to_bincode())  # to support linting
+        plt.suptitle(f"sim #: {idx}")
+
+        if loco0.fc is not None:
+            ax[0].plot(
+                np.array(sim.history.time_seconds) / 3_600,
+                np.array(loco0.fc.history.pwr_fuel_watts) / 1e6,
+                # label='fuel'
+            )
+        # ax[0].plot(
+        #     np.array(sim.history.time_seconds) / 3_600,
+        #     np.array(loco0.history.pwr_out_watts) / 1e6,
+        #     label='conv. loco. tractive'
+        # )
+        # ax[0].plot(
+        #     np.array(sim.history.time_seconds) / 3_600,
+        #     np.array(loco1.history.pwr_out_watts) / 1e6
+        #     label='BEL tractive'
+        # )
+        ax[0].set_ylabel("Single Loco.\nFuel Power [MW]")
+        # ax[0].legend()
+
+        # TODO: Figure out robust way to ensure one bel in demo consist
+        if len(sim.loco_con.loco_vec.tolist()) > 1:
+            loco1 = sim.loco_con.loco_vec.tolist()[1]
+            if loco1.res is not None:
+                ax[1].plot(
+                    np.array(sim.history.time_seconds) / 3_600, loco1.res.history.soc
+                )
+                ax[1].set_ylabel("SOC")
+
+        ax[-1].plot(
+            np.array(sim.history.time_seconds) / 3_600,
+            sim.history.speed_meters_per_second,
+            label="actual",
+        )
+        ax[-1].plot(
+            np.array(sim.history.time_seconds) / 3_600,
+            sim.history.speed_limit_meters_per_second,
+            label="limit",
+        )
+        # ax[-1].plot(
+        #     sim.history.time_seconds,
+        #     sim.history.speed_target_meters_per_second,
+        #     label='target',
+        #     linestyle="-."
+        # )
+        ax[-1].legend()
+        ax[-1].set_xlabel("Time [hr]")
+        ax[-1].set_ylabel("Speed [m/s]")
+
+        plt.tight_layout()
+        plt.savefig(plot_dir / f"sim num {idx}.png")
+        plt.savefig(plot_dir / f"sim num {idx}.svg")
+
+# %%

python/altrios/demos/sim_manager_demo.py

@@ -25,33 +25,40 @@
 t0_import = time.perf_counter()
 t0_total = time.perf_counter()
 
-rail_vehicles=[alt.RailVehicle.from_file(vehicle_file, skip_init=False) 
-               for vehicle_file in Path(alt.resources_root() / "rolling_stock/").glob('*.yaml')]
+rail_vehicles = [
+    alt.RailVehicle.from_file(vehicle_file, skip_init=False)
+    for vehicle_file in Path(alt.resources_root() / "rolling_stock/").glob("*.yaml")
+]
 
-location_map = alt.import_locations(alt.resources_root() / "networks/default_locations.csv")
-network = alt.Network.from_file(alt.resources_root() / "networks/Taconite-NoBalloon.yaml")
+location_map = alt.import_locations(
+    alt.resources_root() / "networks/default_locations.csv"
+)
+network = alt.Network.from_file(
+    alt.resources_root() / "networks/Taconite-NoBalloon.yaml"
+)
 
 t1_import = time.perf_counter()
 print(
     f"Elapsed time to import rail vehicles, locations, and network: {t1_import - t0_import:.3g} s"
 )
 
 train_planner_config = planner_config.TrainPlannerConfig(
-            cars_per_locomotive={"Default": 50},
-            target_cars_per_train={"Default": 90},
-            require_diesel=True)
+    cars_per_locomotive={"Default": 50},
+    target_cars_per_train={"Default": 90},
+    require_diesel=True,
+)
 
 t0_main = time.perf_counter()
 
 (
-    train_consist_plan, 
-    loco_pool, 
-    refuel_facilities, 
-    grid_emissions_factors, 
-    nodal_energy_prices, 
-    speed_limit_train_sims, 
+    train_consist_plan,
+    loco_pool,
+    refuel_facilities,
+    grid_emissions_factors,
+    nodal_energy_prices,
+    speed_limit_train_sims,
     timed_paths,
-    train_consist_plan_untrimmed
+    train_consist_plan_untrimmed,
 ) = sim_manager.main(
     network=network,
     rail_vehicles=rail_vehicles,
@@ -72,13 +79,11 @@
     train_consist_plan_py=train_consist_plan,
     loco_pool_py=loco_pool,
     refuel_facilities_py=refuel_facilities,
-    timed_paths=timed_paths
+    timed_paths=timed_paths,
 )
 t1_train_sims = time.perf_counter()
 print(f"Elapsed time to run train sims: {t1_train_sims-t0_train_sims:.3g} s")
-t_train_time = sum([
-    sim.state.time_seconds for sim in sims.tolist()
-])
+t_train_time = sum([sim.state.time_seconds for sim in sims.tolist()])
 print(f"Total train-seconds simulated: {t_train_time} s")
 
 # %%
@@ -117,8 +122,14 @@
 print(f"Elapsed time to get total fuel energy: {t1_main-t0_main:.3g} s")
 print(f"Total fuel energy used: {e_total_fuel_mj:.3g} GJ")
 
-v_total_fuel_gal = summary_sims.get_energy_fuel_joules(annualize=False) / 1e3 / defaults.LHV_DIESEL_KJ_PER_KG / \
-    defaults.RHO_DIESEL_KG_PER_M3 * utilities.LITER_PER_M3 * utilities.GALLONS_PER_LITER
+v_total_fuel_gal = (
+    summary_sims.get_energy_fuel_joules(annualize=False)
+    / 1e3
+    / defaults.LHV_DIESEL_KJ_PER_KG
+    / defaults.RHO_DIESEL_KG_PER_M3
+    * utilities.LITER_PER_M3
+    * utilities.GALLONS_PER_LITER
+)
 
 print(f"Total fuel used: {v_total_fuel_gal:.3g} gallons")
 print(f"Total elapsed time: {time.perf_counter() - t0_total} s")
@@ -160,10 +171,9 @@
             loco1 = sim.loco_con.loco_vec.tolist()[1]
             if loco1.res is not None:
                 ax[1].plot(
-                    np.array(sim.history.time_seconds) / 3_600,
-                    loco1.res.history.soc
+                    np.array(sim.history.time_seconds) / 3_600, loco1.res.history.soc
                 )
-                ax[1].set_ylabel('SOC')
+                ax[1].set_ylabel("SOC")
 
         ax[-1].plot(
             np.array(sim.history.time_seconds) / 3_600,