Repository: quaquel/EMAworkbench (language: python)

Issue:
Make sure EMAworkbench outcomes and parameter names are valid Python identifiers
Make sure EMAworkbench outcomes and parameter names are valid Python identifiers
Make sure EMAworkbench outcomes and parameter names are valid Python identifiers
Make sure EMAworkbench outcomes and parameter names are valid Python identifiers


Discussion / hints:
@quaquel Since you seem in quite a productive mood, would you like to also resolve this issue before releasing 2.4.1?
The fix is straightforward. There is a [method for it](https://docs.python.org/3/library/stdtypes.html#str.isidentifier).
It needs to be done in Variable (root of uncertainties, outcomes, constants, etc.), AbstractModel (root for all model classes), and Point (root for Scenario and Policy) I guess. 


Main issue is that it will break a lot of existing code (including the final assignment for epa 1361). 

Would be better off in a 2.5/3.0 release. 


@quaquel Since you seem in quite a productive mood, would you like to also resolve this issue before releasing 2.4.1?
The fix is straightforward. There is a [method for it](https://docs.python.org/3/library/stdtypes.html#str.isidentifier).
It needs to be done in Variable (root of uncertainties, outcomes, constants, etc.), AbstractModel (root for all model classes), and Point (root for Scenario and Policy) I guess. 


Main issue is that it will break a lot of existing code (including the final assignment for epa 1361). 

Would be better off in a 2.5/3.0 release. 



Fix patch:
diff --git a/ema_workbench/em_framework/parameters.py b/ema_workbench/em_framework/parameters.py
index ed3108fdf..a924dd3fc 100644
--- a/ema_workbench/em_framework/parameters.py
+++ b/ema_workbench/em_framework/parameters.py
@@ -42,6 +42,10 @@ def __set_name__(self, cls, name):
         self.name = name
         self.internal_name = "_" + name
 
+    @abc.abstractmethod
+    def get_bound(self, instance):
+        ...
+
 
 class UpperBound(Bound):
     def get_bound(self, instance):
@@ -63,7 +67,7 @@ def get_bound(self, owner):
         return bound
 
 
-class Constant(NamedObject):
+class Constant(Variable):
     """Constant class,
 
     can be used for any parameter that has to be set to a fixed value
@@ -78,8 +82,10 @@ def __repr__(self, *args, **kwargs):
         return f"{self.__class__.__name__}('{self.name}', {self.value})"
 
 
-class Category(Constant):
-    pass
+class Category(NamedObject):
+    def __init__(self, name, value):
+        super().__init__(name)
+        self.value = value
 
 
 def create_category(cat):
diff --git a/ema_workbench/em_framework/util.py b/ema_workbench/em_framework/util.py
index d7b493a7a..acb8e1003 100644
--- a/ema_workbench/em_framework/util.py
+++ b/ema_workbench/em_framework/util.py
@@ -62,6 +62,13 @@ def variable_name(self, name):
             name = [name]
         self._variable_name = name
 
+    def __init__(self, name):
+        if not name.isidentifier():
+            DeprecationWarning(
+                f"'{name}' is not a valid Python identifier. Starting from version 3.0 of the EMAworkbench, names must be valid python identifiers"
+            )
+        super().__init__(name)
+
 
 class NamedObjectMap:
     def __init__(self, kind):  # @ReservedAssignment
diff --git a/ema_workbench/examples/data/1000 flu cases no policy.tar.gz b/ema_workbench/examples/data/1000 flu cases no policy.tar.gz
index b8bd17b6d..aa5a1cf82 100644
Binary files a/ema_workbench/examples/data/1000 flu cases no policy.tar.gz and b/ema_workbench/examples/data/1000 flu cases no policy.tar.gz differ
diff --git a/ema_workbench/examples/data/1000 flu cases with policies.tar.gz b/ema_workbench/examples/data/1000 flu cases with policies.tar.gz
index 078441d35..0781f2821 100644
Binary files a/ema_workbench/examples/data/1000 flu cases with policies.tar.gz and b/ema_workbench/examples/data/1000 flu cases with policies.tar.gz differ
diff --git a/ema_workbench/examples/example_pysd_teacup.py b/ema_workbench/examples/example_pysd_teacup.py
index bc3240836..e39d7bb4f 100644
--- a/ema_workbench/examples/example_pysd_teacup.py
+++ b/ema_workbench/examples/example_pysd_teacup.py
@@ -18,7 +18,9 @@
 
     model = PysdModel(mdl_file=mdl_file)
 
-    model.uncertainties = [RealParameter("Room Temperature", 33, 120)]
-    model.outcomes = [TimeSeriesOutcome("Teacup Temperature")]
+    model.uncertainties = [
+        RealParameter("room_temperature", 33, 120, variable_name="Room Temperature")
+    ]
+    model.outcomes = [TimeSeriesOutcome("teacup_temperature", variable_name="Teacup Temperature")]
 
     perform_experiments(model, 100)
diff --git a/ema_workbench/examples/example_vensim_advanced_flu.py b/ema_workbench/examples/example_vensim_advanced_flu.py
index 63f074309..ab4cfca57 100644
--- a/ema_workbench/examples/example_vensim_advanced_flu.py
+++ b/ema_workbench/examples/example_vensim_advanced_flu.py
@@ -35,17 +35,20 @@ def time_of_max(infected_fraction, time):
 
     # outcomes
     model.outcomes = [
-        TimeSeriesOutcome("deceased population region 1"),
-        TimeSeriesOutcome("infected fraction R1"),
+        TimeSeriesOutcome(
+            "deceased_population_region_1", variable_name="deceased population region 1"
+        ),
+        TimeSeriesOutcome("infected_fraction_R1", variable_name="infected fraction R1"),
         ScalarOutcome(
-            "max infection fraction", variable_name="infected fraction R1", function=np.max
+            "max_infection_fraction", variable_name="infected fraction R1", function=np.max
         ),
         ScalarOutcome(
-            "time of max", variable_name=["infected fraction R1", "TIME"], function=time_of_max
+            "time_of_max", variable_name=["infected fraction R1", "TIME"], function=time_of_max
         ),
     ]
 
     # create uncertainties based on csv
+    # FIXME csv is missing
     model.uncertainties = parameters_from_csv("./models/flu/flu_uncertainties.csv")
 
     # add policies
diff --git a/ema_workbench/examples/example_vensim_energy.py b/ema_workbench/examples/example_vensim_energy.py
new file mode 100644
index 000000000..b432d7354
--- /dev/null
+++ b/ema_workbench/examples/example_vensim_energy.py
@@ -0,0 +1,217 @@
+"""
+Created on 27 Jan 2014
+
+@author: jhkwakkel
+"""
+from ema_workbench import (
+    RealParameter,
+    TimeSeriesOutcome,
+    ema_logging,
+    MultiprocessingEvaluator,
+    ScalarOutcome,
+    perform_experiments,
+    CategoricalParameter,
+    save_results,
+    Policy,
+)
+from ema_workbench.connectors.vensim import VensimModel
+from ema_workbench.em_framework.evaluators import SequentialEvaluator
+
+
+def get_energy_model():
+    model = VensimModel(
+        "energyTransition",
+        wd="./models",
+        model_file="RB_V25_ets_1_policy_modified_adaptive_extended_outcomes.vpm",
+    )
+
+    model.outcomes = [
+        TimeSeriesOutcome(
+            "cumulative_carbon_emissions", variable_name="cumulative carbon emissions"
+        ),
+        TimeSeriesOutcome(
+            "carbon_emissions_reduction_fraction",
+            variable_name="carbon emissions reduction fraction",
+        ),
+        TimeSeriesOutcome("fraction_renewables", variable_name="fraction renewables"),
+        TimeSeriesOutcome("average_total_costs", variable_name="average total costs"),
+        TimeSeriesOutcome("total_costs_of_electricity", variable_name="total costs of electricity"),
+    ]
+
+    model.uncertainties = [
+        RealParameter(
+            "demand_fuel_price_elasticity_factor",
+            0,
+            0.5,
+            variable_name="demand fuel price elasticity factor",
+        ),
+        RealParameter(
+            "economic_lifetime_biomass", 30, 50, variable_name="economic lifetime biomass"
+        ),
+        RealParameter("economic_lifetime_coal", 30, 50, variable_name="economic lifetime coal"),
+        RealParameter("economic_lifetime_gas", 25, 40, variable_name="economic lifetime gas"),
+        RealParameter("economic_lifetime_igcc", 30, 50, variable_name="economic lifetime igcc"),
+        RealParameter("economic_lifetime_ngcc", 25, 40, variable_name="economic lifetime ngcc"),
+        RealParameter(
+            "economic_lifetime_nuclear", 50, 70, variable_name="economic lifetime nuclear"
+        ),
+        RealParameter("economic_lifetime_pv", 20, 30, variable_name="economic lifetime pv"),
+        RealParameter("economic_lifetime_wind", 20, 30, variable_name="economic lifetime wind"),
+        RealParameter("economic_lifetime_hydro", 50, 70, variable_name="economic lifetime hydro"),
+        RealParameter(
+            "uncertainty_initial_gross_fuel_costs",
+            0.5,
+            1.5,
+            variable_name="uncertainty initial gross fuel costs",
+        ),
+        RealParameter(
+            "investment_proportionality_constant",
+            0.5,
+            4,
+            variable_name="investment proportionality constant",
+        ),
+        RealParameter(
+            "investors_desired_excess_capacity_investment",
+            0.2,
+            2,
+            variable_name="investors desired excess capacity investment",
+        ),
+        RealParameter(
+            "price_demand_elasticity_factor",
+            -0.07,
+            -0.001,
+            variable_name="price demand elasticity factor",
+        ),
+        RealParameter(
+            "price_volatility_global_resource_markets",
+            0.1,
+            0.2,
+            variable_name="price volatility global resource markets",
+        ),
+        RealParameter("progress_ratio_biomass", 0.85, 1, variable_name="progress ratio biomass"),
+        RealParameter("progress_ratio_coal", 0.9, 1.05, variable_name="progress ratio coal"),
+        RealParameter("progress_ratio_gas", 0.85, 1, variable_name="progress ratio gas"),
+        RealParameter("progress_ratio_igcc", 0.9, 1.05, variable_name="progress ratio igcc"),
+        RealParameter("progress_ratio_ngcc", 0.85, 1, variable_name="progress ratio ngcc"),
+        RealParameter("progress_ratio_nuclear", 0.9, 1.05, variable_name="progress ratio nuclear"),
+        RealParameter("progress_ratio_pv", 0.75, 0.9, variable_name="progress ratio pv"),
+        RealParameter("progress_ratio_wind", 0.85, 1, variable_name="progress ratio wind"),
+        RealParameter("progress_ratio_hydro", 0.9, 1.05, variable_name="progress ratio hydro"),
+        RealParameter(
+            "starting_construction_time", 0.1, 3, variable_name="starting construction time"
+        ),
+        RealParameter(
+            "time_of_nuclear_power_plant_ban",
+            2013,
+            2100,
+            variable_name="time of nuclear power plant ban",
+        ),
+        RealParameter(
+            "weight_factor_carbon_abatement", 1, 10, variable_name="weight factor carbon abatement"
+        ),
+        RealParameter(
+            "weight_factor_marginal_investment_costs",
+            1,
+            10,
+            variable_name="weight factor marginal investment costs",
+        ),
+        RealParameter(
+            "weight_factor_technological_familiarity",
+            1,
+            10,
+            variable_name="weight factor technological familiarity",
+        ),
+        RealParameter(
+            "weight_factor_technological_growth_potential",
+            1,
+            10,
+            variable_name="weight factor technological growth potential",
+        ),
+        RealParameter(
+            "maximum_battery_storage_uncertainty_constant",
+            0.2,
+            3,
+            variable_name="maximum battery storage uncertainty constant",
+        ),
+        RealParameter(
+            "maximum_no_storage_penetration_rate_wind",
+            0.2,
+            0.6,
+            variable_name="maximum no storage penetration rate wind",
+        ),
+        RealParameter(
+            "maximum_no_storage_penetration_rate_pv",
+            0.1,
+            0.4,
+            variable_name="maximum no storage penetration rate pv",
+        ),
+        CategoricalParameter(
+            "SWITCH_lookup_curve_TGC", (1, 2, 3, 4), variable_name="SWITCH lookup curve TGC"
+        ),
+        CategoricalParameter(
+            "SWTICH_preference_carbon_curve", (1, 2), variable_name="SWTICH preference carbon curve"
+        ),
+        CategoricalParameter(
+            "SWITCH_economic_growth", (1, 2, 3, 4, 5, 6), variable_name="SWITCH economic growth"
+        ),
+        CategoricalParameter(
+            "SWITCH_electrification_rate",
+            (1, 2, 3, 4, 5, 6),
+            variable_name="SWITCH electrification rate",
+        ),
+        CategoricalParameter(
+            "SWITCH_Market_price_determination",
+            (1, 2),
+            variable_name="SWITCH Market price determination",
+        ),
+        CategoricalParameter(
+            "SWITCH_physical_limits", (1, 2), variable_name="SWITCH physical limits"
+        ),
+        CategoricalParameter(
+            "SWITCH_low_reserve_margin_price_markup",
+            (1, 2, 3, 4),
+            variable_name="SWITCH low reserve margin price markup",
+        ),
+        CategoricalParameter(
+            "SWITCH_interconnection_capacity_expansion",
+            (1, 2, 3, 4),
+            variable_name="SWITCH interconnection capacity expansion",
+        ),
+        CategoricalParameter(
+            "SWITCH_storage_for_intermittent_supply",
+            (1, 2, 3, 4, 5, 6, 7),
+            variable_name="SWITCH storage for intermittent supply",
+        ),
+        CategoricalParameter("SWITCH_carbon_cap", (1, 2, 3), variable_name="SWITCH carbon cap"),
+        CategoricalParameter(
+            "SWITCH_TGC_obligation_curve", (1, 2, 3), variable_name="SWITCH TGC obligation curve"
+        ),
+        CategoricalParameter(
+            "SWITCH_carbon_price_determination",
+            (1, 2, 3),
+            variable_name="SWITCH carbon price determination",
+        ),
+    ]
+    return model
+
+
+if __name__ == "__main__":
+    ema_logging.log_to_stderr(ema_logging.INFO)
+    model = get_energy_model()
+
+    policies = [
+        Policy("no policy", model_file="RB_V25_ets_1_extended_outcomes.vpm"),
+        Policy("static policy", model_file="ETSPolicy.vpm"),
+        Policy(
+            "adaptive policy",
+            model_file="RB_V25_ets_1_policy_modified_adaptive_extended_outcomes.vpm",
+        ),
+    ]
+
+    n = 100000
+    with MultiprocessingEvaluator(model) as evaluator:
+        experiments, outcomes = evaluator.perform_experiments(n, policies=policies)
+    #
+    # outcomes.pop("TIME")
+    # results = experiments, outcomes
+    # save_results(results, f"./data/{n}_lhs.tar.gz")
diff --git a/ema_workbench/examples/example_vensim_flu.py b/ema_workbench/examples/example_vensim_flu.py
index 035fab3c1..a12ed632c 100644
--- a/ema_workbench/examples/example_vensim_flu.py
+++ b/ema_workbench/examples/example_vensim_flu.py
@@ -16,9 +16,10 @@
     ema_logging,
     ScalarOutcome,
     perform_experiments,
+    Policy,
+    save_results,
 )
 from ema_workbench.connectors.vensim import VensimModel
-from ema_workbench.em_framework.parameters import Policy
 
 if __name__ == "__main__":
     ema_logging.log_to_stderr(ema_logging.INFO)
@@ -27,34 +28,95 @@
 
     # outcomes
     model.outcomes = [
-        TimeSeriesOutcome("deceased population region 1"),
-        TimeSeriesOutcome("infected fraction R1"),
+        TimeSeriesOutcome(
+            "deceased_population_region_1", variable_name="deceased population region 1"
+        ),
+        TimeSeriesOutcome("infected_fraction_R1", variable_name="infected fraction R1"),
         ScalarOutcome(
-            "max infection fraction", variable_name="infected fraction R1", function=np.max
+            "max_infection_fraction", variable_name="infected fraction R1", function=np.max
         ),
     ]
 
     # Plain Parametric Uncertainties
     model.uncertainties = [
-        RealParameter("additional seasonal immune population fraction R1", 0, 0.5),
-        RealParameter("additional seasonal immune population fraction R2", 0, 0.5),
-        RealParameter("fatality ratio region 1", 0.0001, 0.1),
-        RealParameter("fatality rate region 2", 0.0001, 0.1),
-        RealParameter("initial immune fraction of the population of region 1", 0, 0.5),
-        RealParameter("initial immune fraction of the population of region 2", 0, 0.5),
-        RealParameter("normal interregional contact rate", 0, 0.9),
-        RealParameter("permanent immune population fraction R1", 0, 0.5),
-        RealParameter("permanent immune population fraction R2", 0, 0.5),
-        RealParameter("recovery time region 1", 0.1, 0.75),
-        RealParameter("recovery time region 2", 0.1, 0.75),
-        RealParameter("susceptible to immune population delay time region 1", 0.5, 2),
-        RealParameter("susceptible to immune population delay time region 2", 0.5, 2),
-        RealParameter("root contact rate region 1", 0.01, 5),
-        RealParameter("root contact ratio region 2", 0.01, 5),
-        RealParameter("infection ratio region 1", 0, 0.15),
-        RealParameter("infection rate region 2", 0, 0.15),
-        RealParameter("normal contact rate region 1", 10, 100),
-        RealParameter("normal contact rate region 2", 10, 200),
+        RealParameter(
+            "additional_seasonal_immune_population_fraction_R1",
+            0,
+            0.5,
+            variable_name="additional seasonal immune population fraction R1",
+        ),
+        RealParameter(
+            "additional_seasonal_immune_population_fraction_R2",
+            0,
+            0.5,
+            variable_name="additional seasonal immune population fraction R2",
+        ),
+        RealParameter(
+            "fatality_ratio_region_1", 0.0001, 0.1, variable_name="fatality ratio region 1"
+        ),
+        RealParameter(
+            "fatality_rate_region_2", 0.0001, 0.1, variable_name="fatality rate region 2"
+        ),
+        RealParameter(
+            "initial_immune_fraction_of_the_population_of_region_1",
+            0,
+            0.5,
+            variable_name="initial immune fraction of the population of region 1",
+        ),
+        RealParameter(
+            "initial_immune_fraction_of_the_population_of_region_2",
+            0,
+            0.5,
+            variable_name="initial immune fraction of the population of region 2",
+        ),
+        RealParameter(
+            "normal_interregional_contact_rate",
+            0,
+            0.9,
+            variable_name="normal interregional contact rate",
+        ),
+        RealParameter(
+            "permanent_immune_population_fraction_R1",
+            0,
+            0.5,
+            variable_name="permanent immune population fraction R1",
+        ),
+        RealParameter(
+            "permanent_immune_population_fraction_R2",
+            0,
+            0.5,
+            variable_name="permanent immune population fraction R2",
+        ),
+        RealParameter("recovery_time_region_1", 0.1, 0.75, variable_name="recovery time region 1"),
+        RealParameter("recovery_time_region_2", 0.1, 0.75, variable_name="recovery time region 2"),
+        RealParameter(
+            "susceptible_to_immune_population_delay_time_region_1",
+            0.5,
+            2,
+            variable_name="susceptible to immune population delay time region 1",
+        ),
+        RealParameter(
+            "susceptible_to_immune_population_delay_time_region_2",
+            0.5,
+            2,
+            variable_name="susceptible to immune population delay time region 2",
+        ),
+        RealParameter(
+            "root_contact_rate_region_1", 0.01, 5, variable_name="root contact rate region 1"
+        ),
+        RealParameter(
+            "root_contact_ratio_region_2", 0.01, 5, variable_name="root contact ratio region 2"
+        ),
+        RealParameter(
+            "infection_ratio_region_1", 0, 0.15, variable_name="infection ratio region 1"
+        ),
+        RealParameter("infection_rate_region_2", 0, 0.15, variable_name="infection rate region 2"),
+        RealParameter(
+            "normal_contact_rate_region_1", 10, 100, variable_name="normal contact rate region 1"
+        ),
+        RealParameter(
+            "normal_contact_rate_region_2", 10, 200, variable_name="normal contact rate region 2"
+        ),
     ]
 
     # add policies
@@ -65,3 +127,4 @@
     ]
 
     results = perform_experiments(model, 1000, policies=policies)
+    save_results(results, "./data/1000 flu cases with policies.tar.gz")
diff --git a/ema_workbench/examples/example_vensim_lookup.py b/ema_workbench/examples/example_vensim_lookup.py
index 741f3a421..c9fafb6e4 100644
--- a/ema_workbench/examples/example_vensim_lookup.py
+++ b/ema_workbench/examples/example_vensim_lookup.py
@@ -18,10 +18,10 @@
 class Burnout(VensimModel):
     model_file = r"\BURNOUT.vpm"
     outcomes = [
-        TimeSeriesOutcome("Accomplishments to Date"),
-        TimeSeriesOutcome("Energy Level"),
-        TimeSeriesOutcome("Hours Worked Per Week"),
-        TimeSeriesOutcome("accomplishments per hour"),
+        TimeSeriesOutcome("accomplishments_to_date", variable_name="Accomplishments to Date"),
+        TimeSeriesOutcome("energy_level", variable_name="Energy Level"),
+        TimeSeriesOutcome("hours_worked_per_week", variable_name="Hours Worked Per Week"),
+        TimeSeriesOutcome("accomplishments_per_hour", variable_name="accomplishments per hour"),
     ]
 
     def __init__(self, working_directory, name):
diff --git a/ema_workbench/examples/example_vensim_no_policy_flu.py b/ema_workbench/examples/example_vensim_no_policy_flu.py
index ab9ad70d3..bf3c32851 100644
--- a/ema_workbench/examples/example_vensim_no_policy_flu.py
+++ b/ema_workbench/examples/example_vensim_no_policy_flu.py
@@ -14,6 +14,7 @@
     ema_logging,
     perform_experiments,
     MultiprocessingEvaluator,
+    save_results,
 )
 
 from ema_workbench.connectors.vensim import VensimModel
@@ -21,37 +22,100 @@
 if __name__ == "__main__":
     ema_logging.log_to_stderr(ema_logging.INFO)
 
-    model = VensimModel("fluCase", wd="./models/flu", model_file="FLUvensimV1basecase.vpm")
+    model = VensimModel("fluCase", wd=r"./models/flu", model_file=r"FLUvensimV1basecase.vpm")
 
     # outcomes
     model.outcomes = [
-        TimeSeriesOutcome("deceased population region 1"),
-        TimeSeriesOutcome("infected fraction R1"),
+        TimeSeriesOutcome(
+            "deceased_population_region_1", variable_name="deceased population region 1"
+        ),
+        TimeSeriesOutcome("infected_fraction_R1", variable_name="infected fraction R1"),
     ]
 
     # Plain Parametric Uncertainties
     model.uncertainties = [
-        RealParameter("additional seasonal immune population fraction R1", 0, 0.5),
-        RealParameter("additional seasonal immune population fraction R2", 0, 0.5),
-        RealParameter("fatality ratio region 1", 0.0001, 0.1),
-        RealParameter("fatality rate region 2", 0.0001, 0.1),
-        RealParameter("initial immune fraction of the population of region 1", 0, 0.5),
-        RealParameter("initial immune fraction of the population of region 2", 0, 0.5),
-        RealParameter("normal interregional contact rate", 0, 0.9),
-        RealParameter("permanent immune population fraction R1", 0, 0.5),
-        RealParameter("permanent immune population fraction R2", 0, 0.5),
-        RealParameter("recovery time region 1", 0.1, 0.75),
-        RealParameter("recovery time region 2", 0.1, 0.75),
-        RealParameter("susceptible to immune population delay time region 1", 0.5, 2),
-        RealParameter("susceptible to immune population delay time region 2", 0.5, 2),
-        RealParameter("root contact rate region 1", 0.01, 5),
-        RealParameter("root contact ratio region 2", 0.01, 5),
-        RealParameter("infection ratio region 1", 0, 0.15),
-        RealParameter("infection rate region 2", 0, 0.15),
-        RealParameter("normal contact rate region 1", 10, 100),
-        RealParameter("normal contact rate region 2", 10, 200),
+        RealParameter(
+            "additional_seasonal_immune_population_fraction_R1",
+            0,
+            0.5,
+            variable_name="additional seasonal immune population fraction R1",
+        ),
+        RealParameter(
+            "additional_seasonal_immune_population_fraction_R2",
+            0,
+            0.5,
+            variable_name="additional seasonal immune population fraction R2",
+        ),
+        RealParameter(
+            "fatality_ratio_region_1", 0.0001, 0.1, variable_name="fatality ratio region 1"
+        ),
+        RealParameter(
+            "fatality_rate_region_2", 0.0001, 0.1, variable_name="fatality rate region 2"
+        ),
+        RealParameter(
+            "initial_immune_fraction_of_the_population_of_region_1",
+            0,
+            0.5,
+            variable_name="initial immune fraction of the population of region 1",
+        ),
+        RealParameter(
+            "initial_immune_fraction_of_the_population_of_region_2",
+            0,
+            0.5,
+            variable_name="initial immune fraction of the population of region 2",
+        ),
+        RealParameter(
+            "normal_interregional_contact_rate",
+            0,
+            0.9,
+            variable_name="normal interregional contact rate",
+        ),
+        RealParameter(
+            "permanent_immune_population_fraction_R1",
+            0,
+            0.5,
+            variable_name="permanent immune population fraction R1",
+        ),
+        RealParameter(
+            "permanent_immune_population_fraction_R2",
+            0,
+            0.5,
+            variable_name="permanent immune population fraction R2",
+        ),
+        RealParameter("recovery_time_region_1", 0.1, 0.75, variable_name="recovery time region 1"),
+        RealParameter("recovery_time_region_2", 0.1, 0.75, variable_name="recovery time region 2"),
+        RealParameter(
+            "susceptible_to_immune_population_delay_time_region_1",
+            0.5,
+            2,
+            variable_name="susceptible to immune population delay time region 1",
+        ),
+        RealParameter(
+            "susceptible_to_immune_population_delay_time_region_2",
+            0.5,
+            2,
+            variable_name="susceptible to immune population delay time region 2",
+        ),
+        RealParameter(
+            "root_contact_rate_region_1", 0.01, 5, variable_name="root contact rate region 1"
+        ),
+        RealParameter(
+            "root_contact_ratio_region_2", 0.01, 5, variable_name="root contact ratio region 2"
+        ),
+        RealParameter(
+            "infection_ratio_region_1", 0, 0.15, variable_name="infection ratio region 1"
+        ),
+        RealParameter("infection_rate_region_2", 0, 0.15, variable_name="infection rate region 2"),
+        RealParameter(
+            "normal_contact_rate_region_1", 10, 100, variable_name="normal contact rate region 1"
+        ),
+        RealParameter(
+            "normal_contact_rate_region_2", 10, 200, variable_name="normal contact rate region 2"
+        ),
     ]
 
-    nr_experiments = 10
+    nr_experiments = 1000
     with MultiprocessingEvaluator(model) as evaluator:
         results = perform_experiments(model, nr_experiments, evaluator=evaluator)
+
+    save_results(results, "./data/1000 flu cases no policy.tar.gz")
diff --git a/ema_workbench/examples/example_vensim_scarcity.py b/ema_workbench/examples/example_vensim_scarcity.py
index 2a326dd0f..cefdd1c23 100644
--- a/ema_workbench/examples/example_vensim_scarcity.py
+++ b/ema_workbench/examples/example_vensim_scarcity.py
@@ -38,25 +38,25 @@ def priceSubstite(self, x, speed, begin, end):
     def run_model(self, scenario, policy):
         """Method for running an instantiated model structure"""
         kwargs = scenario
-        loc = kwargs.pop("lookup shortage loc")
-        speed = kwargs.pop("lookup shortage speed")
+        loc = kwargs.pop("lookup_shortage_loc")
+        speed = kwargs.pop("lookup_shortage_speed")
         lookup = [self.f(x / 10, speed, loc) for x in range(0, 100)]
         kwargs["shortage price effect lookup"] = lookup
 
-        speed = kwargs.pop("lookup price substitute speed")
-        begin = kwargs.pop("lookup price substitute begin")
-        end = kwargs.pop("lookup price substitute end")
+        speed = kwargs.pop("lookup_price_substitute_speed")
+        begin = kwargs.pop("lookup_price_substitute_begin")
+        end = kwargs.pop("lookup_price_substitute_end")
         lookup = [self.priceSubstite(x, speed, begin, end) for x in range(0, 100, 10)]
         kwargs["relative price substitute lookup"] = lookup
 
-        scale = kwargs.pop("lookup returns to scale speed")
-        speed = kwargs.pop("lookup returns to scale scale")
+        scale = kwargs.pop("lookup_returns_to_scale_speed")
+        speed = kwargs.pop("lookup_returns_to_scale_scale")
         lookup = [self.returnsToScale(x, speed, scale) for x in range(0, 101, 10)]
         kwargs["returns to scale lookup"] = lookup
 
-        scale = kwargs.pop("lookup approximated learning speed")
-        speed = kwargs.pop("lookup approximated learning scale")
-        start = kwargs.pop("lookup approximated learning start")
+        scale = kwargs.pop("lookup_approximated_learning_speed")
+        speed = kwargs.pop("lookup_approximated_learning_scale")
+        start = kwargs.pop("lookup_approximated_learning_start")
         lookup = [self.approxLearning(x, speed, scale, start) for x in range(0, 101, 10)]
         kwargs["approximated learning effect lookup"] = lookup
 
@@ -66,48 +66,143 @@ def run_model(self, scenario, policy):
 if __name__ == "__main__":
     ema_logging.log_to_stderr(ema_logging.DEBUG)
 
-    model = ScarcityModel("scarcity", wd=r"./models/scarcity", model_file=r"\MetalsEMA.vpm")
+    model = ScarcityModel("scarcity", wd="./models/scarcity", model_file="MetalsEMA.vpm")
 
     model.outcomes = [
-        TimeSeriesOutcome("relative market price"),
-        TimeSeriesOutcome("supply demand ratio"),
-        TimeSeriesOutcome("real annual demand"),
-        TimeSeriesOutcome("produced of intrinsically demanded"),
-        TimeSeriesOutcome("supply"),
-        TimeSeriesOutcome("Installed Recycling Capacity"),
-        TimeSeriesOutcome("Installed Extraction Capacity"),
+        TimeSeriesOutcome("relative_market_price", variable_name="relative market price"),
+        TimeSeriesOutcome("supply_demand_ratio", variable_name="supply demand ratio"),
+        TimeSeriesOutcome("real_annual_demand", variable_name="real annual demand"),
+        TimeSeriesOutcome(
+            "produced_of_intrinsically_demanded", variable_name="produced of intrinsically demanded"
+        ),
+        TimeSeriesOutcome("supply", variable_name="supply"),
+        TimeSeriesOutcome(
+            "Installed_Recycling_Capacity", variable_name="Installed Recycling Capacity"
+        ),
+        TimeSeriesOutcome(
+            "Installed_Extraction_Capacity", variable_name="Installed Extraction Capacity"
+        ),
     ]
 
     model.uncertainties = [
-        RealParameter("price elasticity of demand", 0, 0.5),
-        RealParameter("fraction of maximum extraction capacity used", 0.6, 1.2),
-        RealParameter("initial average recycling cost", 1, 4),
-        RealParameter("exogenously planned extraction capacity", 0, 15000),
-        RealParameter("absolute recycling loss fraction", 0.1, 0.5),
-        RealParameter("normal profit margin", 0, 0.4),
-        RealParameter("initial annual supply", 100000, 120000),
-        RealParameter("initial in goods", 1500000, 2500000),
-        RealParameter("average construction time extraction capacity", 1, 10),
-        RealParameter("average lifetime extraction capacity", 20, 40),
-        RealParameter("average lifetime recycling capacity", 20, 40),
-        RealParameter("initial extraction capacity under construction", 5000, 20000),
-        RealParameter("initial recycling capacity under construction", 5000, 20000),
-        RealParameter("initial recycling infrastructure", 5000, 20000),
+        RealParameter(
+            "price_elasticity_of_demand", 0, 0.5, variable_name="price elasticity of demand"
+        ),
+        RealParameter(
+            "fraction_of_maximum_extraction_capacity_used",
+            0.6,
+            1.2,
+            variable_name="fraction of maximum extraction capacity used",
+        ),
+        RealParameter(
+            "initial_average_recycling_cost", 1, 4, variable_name="initial average recycling cost"
+        ),
+        RealParameter(
+            "exogenously_planned_extraction_capacity",
+            0,
+            15000,
+            variable_name="exogenously planned extraction capacity",
+        ),
+        RealParameter(
+            "absolute_recycling_loss_fraction",
+            0.1,
+            0.5,
+            variable_name="absolute recycling loss fraction",
+        ),
+        RealParameter("normal_profit_margin", 0, 0.4, variable_name="normal profit margin"),
+        RealParameter(
+            "initial_annual_supply", 100000, 120000, variable_name="initial annual supply"
+        ),
+        RealParameter("initial_in_goods", 1500000, 2500000, variable_name="initial in goods"),
+        RealParameter(
+            "average_construction_time_extraction_capacity",
+            1,
+            10,
+            variable_name="average construction time extraction capacity",
+        ),
+        RealParameter(
+            "average_lifetime_extraction_capacity",
+            20,
+            40,
+            variable_name="average lifetime extraction capacity",
+        ),
+        RealParameter(
+            "average_lifetime_recycling_capacity",
+            20,
+            40,
+            variable_name="average lifetime recycling capacity",
+        ),
+        RealParameter(
+            "initial_extraction_capacity_under_construction",
+            5000,
+            20000,
+            variable_name="initial extraction capacity under construction",
+        ),
+        RealParameter(
+            "initial_recycling_capacity_under_construction",
+            5000,
+            20000,
+            variable_name="initial recycling capacity under construction",
+        ),
+        RealParameter(
+            "initial_recycling_infrastructure",
+            5000,
+            20000,
+            variable_name="initial recycling infrastructure",
+        ),
         # order of delay
-        CategoricalParameter("order in goods delay", (1, 4, 10, 1000)),
-        CategoricalParameter("order recycling capacity delay", (1, 4, 10)),
-        CategoricalParameter("order extraction capacity delay", (1, 4, 10)),
+        CategoricalParameter(
+            "order_in_goods_delay", (1, 4, 10, 1000), variable_name="order in goods delay"
+        ),
+        CategoricalParameter(
+            "order_recycling_capacity_delay",
+            (1, 4, 10),
+            variable_name="order recycling capacity delay",
+        ),
+        CategoricalParameter(
+            "order_extraction_capacity_delay",
+            (1, 4, 10),
+            variable_name="order extraction capacity delay",
+        ),
         # uncertainties associated with lookups
-        RealParameter("lookup shortage loc", 20, 50),
-        RealParameter("lookup shortage speed", 1, 5),
-        RealParameter("lookup price substitute speed", 0.1, 0.5),
-        RealParameter("lookup price substitute begin", 3, 7),
-        RealParameter("lookup price substitute end", 15, 25),
-        RealParameter("lookup returns to scale speed", 0.01, 0.2),
-        RealParameter("lookup returns to scale scale", 0.3, 0.7),
-        RealParameter("lookup approximated learning speed", 0.01, 0.2),
-        RealParameter("lookup approximated learning scale", 0.3, 0.6),
-        RealParameter("lookup approximated learning start", 30, 60),
+        RealParameter("lookup_shortage_loc", 20, 50, variable_name="lookup shortage loc"),
+        RealParameter("lookup_shortage_speed", 1, 5, variable_name="lookup shortage speed"),
+        RealParameter(
+            "lookup_price_substitute_speed", 0.1, 0.5, variable_name="lookup price substitute speed"
+        ),
+        RealParameter(
+            "lookup_price_substitute_begin", 3, 7, variable_name="lookup price substitute begin"
+        ),
+        RealParameter(
+            "lookup_price_substitute_end", 15, 25, variable_name="lookup price substitute end"
+        ),
+        RealParameter(
+            "lookup_returns_to_scale_speed",
+            0.01,
+            0.2,
+            variable_name="lookup returns to scale speed",
+        ),
+        RealParameter(
+            "lookup_returns_to_scale_scale", 0.3, 0.7, variable_name="lookup returns to scale scale"
+        ),
+        RealParameter(
+            "lookup_approximated_learning_speed",
+            0.01,
+            0.2,
+            variable_name="lookup approximated learning speed",
+        ),
+        RealParameter(
+            "lookup_approximated_learning_scale",
+            0.3,
+            0.6,
+            variable_name="lookup approximated learning scale",
+        ),
+        RealParameter(
+            "lookup_approximated_learning_start",
+            30,
+            60,
+            variable_name="lookup approximated learning start",
+        ),
     ]
 
-    results = perform_experiments(model, 50)
+    results = perform_experiments(model, 1000)
diff --git a/ema_workbench/examples/feature_scoring_flu.py b/ema_workbench/examples/feature_scoring_flu.py
index 57777902f..16e442619 100644
--- a/ema_workbench/examples/feature_scoring_flu.py
+++ b/ema_workbench/examples/feature_scoring_flu.py
@@ -18,8 +18,8 @@
 
 # we have timeseries so we need scalars
 y = {
-    "deceased population": outcomes["deceased population region 1"][:, -1],
-    "max. infected fraction": np.max(outcomes["infected fraction R1"], axis=1),
+    "deceased population": outcomes["deceased_population_region_1"][:, -1],
+    "max. infected fraction": np.max(outcomes["infected_fraction_R1"], axis=1),
 }
 
 scores = feature_scoring.get_feature_scores_all(x, y)
diff --git a/ema_workbench/examples/feature_scoring_flu_confidence.py b/ema_workbench/examples/feature_scoring_flu_confidence.py
index ea3a947d8..697472b1a 100644
--- a/ema_workbench/examples/feature_scoring_flu_confidence.py
+++ b/ema_workbench/examples/feature_scoring_flu_confidence.py
@@ -18,7 +18,7 @@
 x, outcomes = load_results(fn)
 
 x = x.drop(["model", "policy"], axis=1)
-y = np.max(outcomes["infected fraction R1"], axis=1)
+y = np.max(outcomes["infected_fraction_R1"], axis=1)
 
 all_scores = []
 for i in range(100):
diff --git a/ema_workbench/examples/feature_scoring_flu_overtime.py b/ema_workbench/examples/feature_scoring_flu_overtime.py
index 41fd1e0b9..5879d3ccf 100644
--- a/ema_workbench/examples/feature_scoring_flu_overtime.py
+++ b/ema_workbench/examples/feature_scoring_flu_overtime.py
@@ -18,7 +18,7 @@
 x, outcomes = load_results(fn)
 x = x.drop(["model", "policy"], axis=1)
 
-y = outcomes["deceased population region 1"]
+y = outcomes["deceased_population_region_1"]
 
 all_scores = []
 
diff --git a/ema_workbench/examples/models/flu/Current.vdfx b/ema_workbench/examples/models/flu/Current.vdfx
index 58349bce3..7398e772f 100644
Binary files a/ema_workbench/examples/models/flu/Current.vdfx and b/ema_workbench/examples/models/flu/Current.vdfx differ
diff --git a/ema_workbench/examples/models/hamarat/ETSPolicy.vpm b/ema_workbench/examples/models/hamarat/ETSPolicy.vpm
new file mode 100644
index 000000000..ef6d5f814
Binary files /dev/null and b/ema_workbench/examples/models/hamarat/ETSPolicy.vpm differ
diff --git a/ema_workbench/examples/models/hamarat/RB_V25_ets_1_extended_outcomes.vpm b/ema_workbench/examples/models/hamarat/RB_V25_ets_1_extended_outcomes.vpm
new file mode 100644
index 000000000..f3b3e09f4
Binary files /dev/null and b/ema_workbench/examples/models/hamarat/RB_V25_ets_1_extended_outcomes.vpm differ
diff --git a/ema_workbench/examples/models/hamarat/RB_V25_ets_1_policy_modified_adaptive_extended_outcomes.vpm b/ema_workbench/examples/models/hamarat/RB_V25_ets_1_policy_modified_adaptive_extended_outcomes.vpm
new file mode 100644
index 000000000..14bf77ae0
Binary files /dev/null and b/ema_workbench/examples/models/hamarat/RB_V25_ets_1_policy_modified_adaptive_extended_outcomes.vpm differ
diff --git a/ema_workbench/examples/models/scarcity/Current.vdfx b/ema_workbench/examples/models/scarcity/Current.vdfx
new file mode 100644
index 000000000..c83246d58
Binary files /dev/null and b/ema_workbench/examples/models/scarcity/Current.vdfx differ
diff --git a/ema_workbench/examples/plotting_pairsplot_flu.py b/ema_workbench/examples/plotting_pairsplot_flu.py
index 9e2ff92a7..4f994f5fe 100644
--- a/ema_workbench/examples/plotting_pairsplot_flu.py
+++ b/ema_workbench/examples/plotting_pairsplot_flu.py
@@ -24,7 +24,7 @@
 time = outcomes.pop("TIME")
 
 for key, value in outcomes.items():
-    if key == "deceased population region 1":
+    if key == "deceased_population_region_1":
         tr[key] = value[:, -1]  # we want the end value
     else:
         # we want the maximum value of the peak
diff --git a/ema_workbench/examples/sd_boostedtrees_flu.py b/ema_workbench/examples/sd_boostedtrees_flu.py
index 073624d68..30d2c62e0 100644
--- a/ema_workbench/examples/sd_boostedtrees_flu.py
+++ b/ema_workbench/examples/sd_boostedtrees_flu.py
@@ -49,7 +49,7 @@ def plot_diag(x1, ax):
 
 # transform to numpy array with proper recoding of cateogorical variables
 x, columns = feature_scoring._prepare_experiments(experiments)
-y = outcomes["deceased population region 1"][:, -1] > 1000000
+y = outcomes["deceased_population_region 1"][:, -1] > 1000000
 
 # establish mean case for factor maps
 # this is questionable in particular in case of categorical dimensions
diff --git a/ema_workbench/examples/sd_cart_flu.py b/ema_workbench/examples/sd_cart_flu.py
index 399a9686d..248f00b5e 100644
--- a/ema_workbench/examples/sd_cart_flu.py
+++ b/ema_workbench/examples/sd_cart_flu.py
@@ -13,7 +13,7 @@
 
 def classify(data):
     # get the output for deceased population
-    result = data["deceased population region 1"]
+    result = data["deceased_population_region_1"]
 
     # if deceased population is higher then 1.000.000 people,
     # classify as 1
diff --git a/ema_workbench/examples/sd_dimensional_stacking_flu.py b/ema_workbench/examples/sd_dimensional_stacking_flu.py
index 6ddbdd302..c99883467 100644
--- a/ema_workbench/examples/sd_dimensional_stacking_flu.py
+++ b/ema_workbench/examples/sd_dimensional_stacking_flu.py
@@ -18,7 +18,7 @@
 fn = "./data/1000 flu cases no policy.tar.gz"
 x, outcomes = load_results(fn)
 
-y = outcomes["deceased population region 1"][:, -1] > 1000000
+y = outcomes["deceased_population_region_1"][:, -1] > 1000000
 
 fig = dimensional_stacking.create_pivot_plot(x, y, 2, bin_labels=True)
 
diff --git a/ema_workbench/examples/sd_logit_flu_example.py b/ema_workbench/examples/sd_logit_flu_example.py
index f81d3b4a0..b3e8c999b 100644
--- a/ema_workbench/examples/sd_logit_flu_example.py
+++ b/ema_workbench/examples/sd_logit_flu_example.py
@@ -15,7 +15,7 @@
 experiments, outcomes = load_results("./data/1000 flu cases no policy.tar.gz")
 
 x = experiments.drop(["model", "policy"], axis=1)
-y = outcomes["deceased population region 1"][:, -1] > 1000000
+y = outcomes["deceased_population_region_1"][:, -1] > 1000000
 
 logit = logistic_regression.Logit(x, y)
 logit.run()
diff --git a/ema_workbench/examples/sd_prim_PCA_flu.py b/ema_workbench/examples/sd_prim_PCA_flu.py
index 1982500fa..c6d151c15 100644
--- a/ema_workbench/examples/sd_prim_PCA_flu.py
+++ b/ema_workbench/examples/sd_prim_PCA_flu.py
@@ -23,7 +23,7 @@
 x, outcomes = load_results(fn)
 
 # specify y
-y = outcomes["deceased population region 1"][:, -1] > 1000000
+y = outcomes["deceased_population_region_1"][:, -1] > 1000000
 
 rotated_experiments, rotation_matrix = prim.pca_preprocess(x, y, exclude=["model", "policy"])
 
diff --git a/ema_workbench/examples/sd_prim_flu.py b/ema_workbench/examples/sd_prim_flu.py
index 706d1b47e..abb159103 100644
--- a/ema_workbench/examples/sd_prim_flu.py
+++ b/ema_workbench/examples/sd_prim_flu.py
@@ -21,7 +21,7 @@
 
 def classify(data):
     # get the output for deceased population
-    ooi = data["deceased population region 1"]
+    ooi = data["deceased_population_region_1"]
     return ooi[:, -1] > 1000000
 
 
diff --git a/ema_workbench/examples/timeseries_clustering_flu.py b/ema_workbench/examples/timeseries_clustering_flu.py
index 38c892c9b..b1dd8c024 100644
--- a/ema_workbench/examples/timeseries_clustering_flu.py
+++ b/ema_workbench/examples/timeseries_clustering_flu.py
@@ -10,7 +10,7 @@
 from ema_workbench.analysis import clusterer, plotting, Density
 
 experiments, outcomes = load_results("./data/1000 flu cases no policy.tar.gz")
-data = outcomes["infected fraction R1"]
+data = outcomes["infected_fraction_R1"]
 
 # calculate distances
 distances = clusterer.calculate_cid(data)


Test patch:
diff --git a/test/data/1000 flu cases no policy.tar.gz b/test/data/1000 flu cases no policy.tar.gz
index c55082d9c..aa5a1cf82 100644
Binary files a/test/data/1000 flu cases no policy.tar.gz and b/test/data/1000 flu cases no policy.tar.gz differ
diff --git a/test/data/1000 runs scarcity.tar.gz b/test/data/1000 runs scarcity.tar.gz
index 905389040..a571a5629 100644
Binary files a/test/data/1000 runs scarcity.tar.gz and b/test/data/1000 runs scarcity.tar.gz differ
diff --git a/test/data/eng_trans.tar.gz b/test/data/eng_trans.tar.gz
index d86b7b62c..82e77f1c7 100644
Binary files a/test/data/eng_trans.tar.gz and b/test/data/eng_trans.tar.gz differ
diff --git a/test/test_analysis/test_cart.py b/test/test_analysis/test_cart.py
index 878af1a19..f4006a06b 100644
--- a/test/test_analysis/test_cart.py
+++ b/test/test_analysis/test_cart.py
@@ -16,7 +16,7 @@
 
 def flu_classify(data):
     # get the output for deceased population
-    result = data["deceased population region 1"]
+    result = data["deceased_population_region_1"]
 
     # make an empty array of length equal to number of cases
     classes = np.zeros(result.shape[0])
@@ -28,7 +28,7 @@ def flu_classify(data):
 
 
 def scarcity_classify(outcomes):
-    outcome = outcomes["relative market price"]
+    outcome = outcomes["relative_market_price"]
     change = np.abs(outcome[:, 1::] - outcome[:, 0:-1])
 
     neg_change = np.min(change, axis=1)
@@ -57,7 +57,7 @@ def test_setup_cart(self):
             y[k] = v[:, -1]
 
         temp_results = (x, y)
-        alg = cart.setup_cart(temp_results, "deceased population region 1", mass_min=0.05)
+        alg = cart.setup_cart(temp_results, "deceased_population_region_1", mass_min=0.05)
         self.assertTrue(alg.mode == RuleInductionType.REGRESSION)
 
         n_cols = 5
@@ -137,7 +137,7 @@ def test_stats_to_dataframe(self):
         alg.build_tree()
         stats = alg.stats_to_dataframe()
 
-        y = outcomes["deceased population region 1"][:, -1]
+        y = outcomes["deceased_population_region_1"][:, -1]
         alg = cart.CART(x, y, mode=RuleInductionType.REGRESSION)
         alg.build_tree()
         stats = alg.stats_to_dataframe()
@@ -163,7 +163,7 @@ def test_build_tree(self):
             y[k] = v
 
         temp_results = (x, y)
-        alg = cart.setup_cart(temp_results, "deceased population region 1", mass_min=0.05)
+        alg = cart.setup_cart(temp_results, "deceased_population_region_1", mass_min=0.05)
         alg.build_tree()
         self.assertTrue(isinstance(alg.clf, cart.tree.DecisionTreeRegressor))
 
diff --git a/test/test_analysis/test_clusterer.py b/test/test_analysis/test_clusterer.py
index 6d5a2acbd..40415045e 100644
--- a/test/test_analysis/test_clusterer.py
+++ b/test/test_analysis/test_clusterer.py
@@ -10,7 +10,7 @@ class ClusterTestCase(unittest.TestCase):
     def test_cluster(self):
         n = 10
         experiments, outcomes = utilities.load_flu_data()
-        data = outcomes["infected fraction R1"][0:n, :]
+        data = outcomes["infected_fraction_R1"][0:n, :]
 
         distances = clusterer.calculate_cid(data)
         self.assertEqual(distances.shape, (n, n))
diff --git a/test/test_analysis/test_dimensional_stacking.py b/test/test_analysis/test_dimensional_stacking.py
index 3c81a7c64..a247a2281 100644
--- a/test/test_analysis/test_dimensional_stacking.py
+++ b/test/test_analysis/test_dimensional_stacking.py
@@ -35,7 +35,7 @@ def test_discretize(self):
 
     def test_create_pivot_plot(self):
         x, outcomes = utilities.load_flu_data()
-        y = outcomes["deceased population region 1"][:, -1] > 1000000
+        y = outcomes["deceased_population_region_1"][:, -1] > 1000000
 
         dimensional_stacking.create_pivot_plot(x, y, 2)
         dimensional_stacking.create_pivot_plot(x, y, 2, labels=False, bin_labels=True)
diff --git a/test/test_analysis/test_feature_scoring.py b/test/test_analysis/test_feature_scoring.py
index e819d0796..f602348d1 100644
--- a/test/test_analysis/test_feature_scoring.py
+++ b/test/test_analysis/test_feature_scoring.py
@@ -64,7 +64,7 @@ def test_prepare_outcomes(self):
 
         # string type correct
         ooi = "nr deaths"
-        outcomes[ooi] = outcomes["deceased population region 1"][:, -1]
+        outcomes[ooi] = outcomes["deceased_population_region_1"][:, -1]
         y, categorical = fs._prepare_outcomes(outcomes, ooi)
 
         self.assertFalse(categorical)
@@ -76,7 +76,7 @@ def test_prepare_outcomes(self):
 
         # classify function correct
         def classify(data):
-            result = data["deceased population region 1"]
+            result = data["deceased_population_region_1"]
             classes = np.zeros(result.shape[0])
             classes[result[:, -1] > 1000000] = 1
             return classes
@@ -95,7 +95,7 @@ def test_get_univariate_feature_scores(self):
 
         def classify(data):
             # get the output for deceased population
-            result = data["deceased population region 1"]
+            result = data["deceased_population_region_1"]
 
             # make an empty array of length equal to number of cases
             classes = np.zeros(result.shape[0])
@@ -116,7 +116,7 @@ def classify(data):
         self.assertEqual(len(scores), len(x.columns) - 3)
 
         # f regression
-        y = outcomes["deceased population region 1"][:, -1]
+        y = outcomes["deceased_population_region_1"][:, -1]
         scores = fs.get_univariate_feature_scores(x, y, score_func=F_REGRESSION)
         self.assertEqual(len(scores), len(x.columns) - 3)
 
@@ -125,7 +125,7 @@ def test_get_rf_feature_scores(self):
 
         def classify(data):
             # get the output for deceased population
-            result = data["deceased population region 1"]
+            result = data["deceased_population_region_1"]
 
             # make an empty array of length equal to number of cases
             classes = np.zeros(result.shape[0])
@@ -147,7 +147,7 @@ def classify(data):
 
         self.assertRaises(ValueError, fs.get_rf_feature_scores, x, y, mode="illegal argument")
 
-        y = outcomes["deceased population region 1"][:, -1]
+        y = outcomes["deceased_population_region_1"][:, -1]
         scores, forest = fs.get_rf_feature_scores(
             x, y, mode=RuleInductionType.REGRESSION, random_state=10
         )
@@ -157,7 +157,7 @@ def classify(data):
 
     def test_get_ex_feature_scores(self):
         x, outcomes = utilities.load_flu_data()
-        y = outcomes["deceased population region 1"][:, -1] > 1000000
+        y = outcomes["deceased_population_region_1"][:, -1] > 1000000
 
         scores, forest = fs.get_ex_feature_scores(
             x, y, mode=RuleInductionType.CLASSIFICATION, random_state=10
@@ -168,7 +168,7 @@ def test_get_ex_feature_scores(self):
 
         self.assertRaises(ValueError, fs.get_ex_feature_scores, x, y, mode="illegal argument")
 
-        y = outcomes["deceased population region 1"][:, -1]
+        y = outcomes["deceased_population_region_1"][:, -1]
         scores, forest = fs.get_ex_feature_scores(
             x, y, mode=RuleInductionType.REGRESSION, random_state=10
         )
@@ -181,8 +181,8 @@ def test_get_feature_scores_all(self):
 
         # we have timeseries so we need scalars
         y = {
-            "deceased population": outcomes["deceased population region 1"][:, -1],
-            "max. infected fraction": np.max(outcomes["infected fraction R1"], axis=1),
+            "deceased population": outcomes["deceased_population_region_1"][:, -1],
+            "max. infected fraction": np.max(outcomes["infected_fraction_R1"], axis=1),
         }
         scores = fs.get_feature_scores_all(x, y)
 
diff --git a/test/test_analysis/test_logistic_regression.py b/test/test_analysis/test_logistic_regression.py
index 941a9a926..3578c5604 100644
--- a/test/test_analysis/test_logistic_regression.py
+++ b/test/test_analysis/test_logistic_regression.py
@@ -15,7 +15,7 @@
 
 def flu_classify(data):
     # get the output for deceased population
-    result = data["deceased population region 1"][:, -1]
+    result = data["deceased_population_region_1"][:, -1]
     return result > 1000000
 
 
diff --git a/test/test_analysis/test_plotting.py b/test/test_analysis/test_plotting.py
index e24a91797..f1fba9451 100644
--- a/test/test_analysis/test_plotting.py
+++ b/test/test_analysis/test_plotting.py
@@ -73,12 +73,12 @@ def test_group_results(self):
         print(experiments.shape[0], total_data)
 
         # test continuous parameter type
-        array = experiments["average planning and construction period T1"]
+        array = experiments["average_planning_and_construction_period_T1"]
         grouping_specifiers = make_continuous_grouping_specifiers(array, nr_of_groups=5)
         groups = group_results(
             experiments,
             outcomes,
-            group_by="average planning and construction period T1",
+            group_by="average_planning_and_construction_period_T1",
             grouping_specifiers=grouping_specifiers,
             grouping_labels=[str(entry) for entry in grouping_specifiers],
         )
@@ -88,12 +88,12 @@ def test_group_results(self):
         print(experiments.shape[0], total_data)
 
         # test integer type
-        array = experiments["seed PR T1"]
+        array = experiments["seed_PR_T1"]
         grouping_specifiers = make_continuous_grouping_specifiers(array, nr_of_groups=10)
         groups = group_results(
             experiments,
             outcomes,
-            group_by="seed PR T1",
+            group_by="seed_PR_T1",
             grouping_specifiers=grouping_specifiers,
             grouping_labels=[str(entry) for entry in grouping_specifiers],
         )
@@ -122,7 +122,7 @@ def test_lines(self):
         lines(
             experiments,
             outcomes,
-            outcomes_to_show="total fraction new technologies",
+            outcomes_to_show="total_fraction_new_technologies",
             experiments_to_show=np.arange(0, 600, 20),
             group_by="policy",
             grouping_specifiers="basic policy",
@@ -292,7 +292,7 @@ def test_envelopes(self):
         envelopes(experiments, outcomes, density=None, titles=None)
         envelopes(experiments, outcomes, density=None, titles={})
         envelopes(
-            experiments, outcomes, density=None, titles={"total fraction new technologies": "a"}
+            experiments, outcomes, density=None, titles={"total_fraction_new_technologies": "a"}
         )
 
         plt.draw()
@@ -302,7 +302,7 @@ def test_envelopes(self):
         envelopes(experiments, outcomes, density=None, ylabels=None)
         envelopes(experiments, outcomes, density=None, ylabels={})
         envelopes(
-            experiments, outcomes, density=None, ylabels={"total fraction new technologies": "a"}
+            experiments, outcomes, density=None, ylabels={"total_fraction_new_technologies": "a"}
         )
 
         plt.draw()
@@ -401,7 +401,7 @@ def test_kde_over_time(self):
 
     def test_multiple_densities(self):
         experiments, outcomes = utilities.load_eng_trans_data()
-        ooi = "total fraction new technologies"
+        ooi = "total_fraction_new_technologies"
 
         multiple_densities(experiments, outcomes, group_by="policy", points_in_time=[2010])
         multiple_densities(
diff --git a/test/test_analysis/test_prim.py b/test/test_analysis/test_prim.py
index bf047491c..4897c77c7 100644
--- a/test/test_analysis/test_prim.py
+++ b/test/test_analysis/test_prim.py
@@ -19,7 +19,7 @@
 
 def flu_classify(data):
     # get the output for deceased population
-    result = data["deceased population region 1"]
+    result = data["deceased_population_region_1"]
 
     # make an empty array of length equal to number of cases
     classes = np.zeros(result.shape[0])
@@ -178,7 +178,7 @@ def test_setup_prim(self):
 
         # test initialization, including t_coi calculation in case of searching
         # for results equal to or higher than the threshold
-        outcomes["death toll"] = outcomes["deceased population region 1"][:, -1]
+        outcomes["death toll"] = outcomes["deceased_population_region_1"][:, -1]
         results = experiments, outcomes
         threshold = 10000
         prim_obj = prim.setup_prim(
@@ -228,7 +228,7 @@ def test_prim_init_select(self):
 
         # test initialization, including t_coi calculation in case of searching
         # for results equal to or higher than the threshold
-        outcomes["death toll"] = outcomes["deceased population region 1"][:, -1]
+        outcomes["death toll"] = outcomes["deceased_population_region_1"][:, -1]
         results = experiments, outcomes
         threshold = 10000
         prim_obj = prim.setup_prim(
@@ -323,7 +323,7 @@ def test_box_init(self):
     def test_prim_exceptions(self):
         results = utilities.load_flu_data()
         x, outcomes = results
-        y = outcomes["deceased population region 1"]
+        y = outcomes["deceased_population_region_1"]
 
         self.assertRaises(
             prim.PrimException, prim.Prim, x, y, threshold=0.8, mode=RuleInductionType.REGRESSION
diff --git a/test/test_analysis/test_regional_sa.py b/test/test_analysis/test_regional_sa.py
index b95bec5b8..81501b06b 100644
--- a/test/test_analysis/test_regional_sa.py
+++ b/test/test_analysis/test_regional_sa.py
@@ -14,7 +14,7 @@
 class Test(unittest.TestCase):
     def test_plot_cdfs(self):
         x, outcomes = utilities.load_flu_data()
-        y = outcomes["deceased population region 1"][:, -1] > 1000000
+        y = outcomes["deceased_population_region_1"][:, -1] > 1000000
 
         regional_sa.plot_cdfs(x, y)
         regional_sa.plot_cdfs(x, y, ccdf=True)
@@ -24,10 +24,10 @@ def test_plot_cdfs(self):
 
     def test_plot__individual_cdf(self):
         x, outcomes = utilities.load_flu_data()
-        y = outcomes["deceased population region 1"][:, -1] > 1000000
+        y = outcomes["deceased_population_region_1"][:, -1] > 1000000
 
         fig, ax = plt.subplots()
-        unc = "fatality ratio region 1"
+        unc = "fatality_ratio_region_1"
 
         regional_sa.plot_individual_cdf(
             ax,
diff --git a/test/test_em_framework/test_optimization.py b/test/test_em_framework/test_optimization.py
index 65714f06d..a028835c4 100644
--- a/test/test_em_framework/test_optimization.py
+++ b/test/test_em_framework/test_optimization.py
@@ -170,8 +170,8 @@ def test_to_robust_problem(self, mocked_platypus):
 
         scenarios = 5
         robustness_functions = [
-            ScalarOutcome("mean x", variable_name="x", function=mock.Mock(), kind="maximize"),
-            ScalarOutcome("mean y", variable_name="y", function=mock.Mock(), kind="maximize"),
+            ScalarOutcome("mean_x", variable_name="x", function=mock.Mock(), kind="maximize"),
+            ScalarOutcome("mean_y", variable_name="y", function=mock.Mock(), kind="maximize"),
         ]
 
         problem = to_robust_problem(mocked_model, scenarios, robustness_functions)
@@ -179,7 +179,7 @@ def test_to_robust_problem(self, mocked_platypus):
         self.assertEqual("robust", problem.searchover)
         for entry in problem.parameters:
             self.assertIn(entry.name, mocked_model.levers.keys())
-        self.assertEqual(["mean x", "mean y"], problem.outcome_names)
+        self.assertEqual(["mean_x", "mean_y"], problem.outcome_names)
 
     def test_process_robust(self):
         pass
diff --git a/test/test_em_framework/test_parameter.py b/test/test_em_framework/test_parameter.py
index fe8dac540..c15cdbb9f 100644
--- a/test/test_em_framework/test_parameter.py
+++ b/test/test_em_framework/test_parameter.py
@@ -55,7 +55,7 @@ def test_comparison(self):
         name = "test"
         par1 = parameters.RealParameter(name, lower_bound, upper_bound, resolution)
 
-        name = "what?"
+        name = "what"
         par2 = parameters.RealParameter(name, lower_bound, upper_bound, resolution)
         self.assertFalse(par1 == par2)
 
diff --git a/test/test_em_framework/test_samplers.py b/test/test_em_framework/test_samplers.py
index a4d939fae..0d4522b53 100644
--- a/test/test_em_framework/test_samplers.py
+++ b/test/test_em_framework/test_samplers.py
@@ -23,9 +23,9 @@
 
 class SamplerTestCase(unittest.TestCase):
     uncertainties = [
-        RealParameter("1", 0, 10),
-        IntegerParameter("2", 0, 10),
-        CategoricalParameter("3", ["a", "b", "c"]),
+        RealParameter("a", 0, 10),
+        IntegerParameter("b", 0, 10),
+        CategoricalParameter("c", ["a", "b", "c"]),
     ]
 
     def _test_generate_designs(self, sampler):
@@ -37,9 +37,9 @@ def _test_generate_designs(self, sampler):
         for design in designs:
             actual_nr_designs += 1
 
-        self.assertIn("1", design, msg)
-        self.assertIn("2", design, msg)
-        self.assertIn("3", design, msg)
+        self.assertIn("a", design, msg)
+        self.assertIn("b", design, msg)
+        self.assertIn("c", design, msg)
         self.assertEqual(designs.n, actual_nr_designs, msg)
 
     def test_lhs_sampler(self):