Repository: IN-CORE/pyincore (language: python)

Issue:
Multiobjective retrofit optimization analysis
This code integrates the multiobjective retrofit optimization analysis by Yunjie and collaborators.

* Code was componentized into functions that implement the analysis and separate responsibilities
* Code only generates three output files insteads of >100, removing IO bounds from performance
* Code reproduces results of original notebook results

remove deprecated methods from pyIncore
remove deprecated methods (@deprecated) that are annotated in pyIncore.

 

remove the following scripts:

fragilitycurve.py

conditionalstandardfragilitycurve.py

customexpressionfragilitycurve.py

parametricfragilitycurve.py

periodbuildingfragilitycurve.py

periodstandardfragilitycurve.py

standardfragilitycurve.py

 

rename fragilitycurverefactored.py to fragilitycurve.py

 

rename fragilitycurveset.py's method calculate_limit_state_refactored_w_conversion to calculate_limit_state
remove deprecated methods from pyIncore
remove deprecated methods (@deprecated) that are annotated in pyIncore.

 

remove the following scripts:

fragilitycurve.py

conditionalstandardfragilitycurve.py

customexpressionfragilitycurve.py

parametricfragilitycurve.py

periodbuildingfragilitycurve.py

periodstandardfragilitycurve.py

standardfragilitycurve.py

 

rename fragilitycurverefactored.py to fragilitycurve.py

 

rename fragilitycurveset.py's method calculate_limit_state_refactored_w_conversion to calculate_limit_state
remove deprecated methods from pyIncore
remove deprecated methods (@deprecated) that are annotated in pyIncore.

 

remove the following scripts:

fragilitycurve.py

conditionalstandardfragilitycurve.py

customexpressionfragilitycurve.py

parametricfragilitycurve.py

periodbuildingfragilitycurve.py

periodstandardfragilitycurve.py

standardfragilitycurve.py

 

rename fragilitycurverefactored.py to fragilitycurve.py

 

rename fragilitycurveset.py's method calculate_limit_state_refactored_w_conversion to calculate_limit_state
pyincore util to convert population dislocation output for heatmap

 show the api response messages that the services return
When pyincore receives an error message from api, it just prints the error code (500 etc.) and ignores the message included. Below is an example of what pyincore shows to user

ERROR - client.py:return_http_response() - HTTPError: The server returned a 500 failure response code. You can find information about HTTP response status codes here: https://developer.mozilla.org/en-US/docs/Web/HTTP/Status
--
raise self._exception
requests.exceptions.HTTPError: 500 Server Error: Internal Server Error for url: https://incore-tst.ncsa.illinois.edu/hazard/api/earthquakes/5dfa4058b9219c934b64d495/values

Add github action to run unit tests
Here is the pytest line from the old build server:

pytest --cache-clear -s pyincore/test_client.py pyincore/test_dataservice.py pyincore/test_dfr3service.py pyincore/test_spaceservice.py pyincore/test_hazardservice.py --junitxml=test-reports/results.xml

A couple of things here. 

1. We need to probably setup a conda environment because requirements.txt has a reference to ipopt
2. We need to use the test user to connect to the services to run these tests
Add Joplin building functionality recovery time model
Implement only the original, empirical part
implement new get_building_period() method
now period and num_stories is a part of expression existing in the fragility curve;

rewrote get_building_period() and need to implement that in various place

<p>dmg_probability = selected_fragility_set.calculate_limit_state_refactored_w_conversion(<br />hval_dict, **building_args, period=xxx)</p>
<p>&nbsp;</p>
<p>in building damage this period =xxx seems to be not needed, cuz it's already a part of the equation of the fragility curves</p>
<p>make sure i test removing it</p>
<p>&nbsp;</p>
<p>using the&nbsp;5b47b34e337d4a36290754a0 as example&nbsp;</p>
<div class="table-wrap">


description | "default building period"
-- | --
fullName | null
expression | "(0.097) * math.pow(num_stories * (13.0), 0.624)"


</div>
<p>&nbsp;</p>
convert hua and pd outputs to json
Lisa's notebook does not run pyincore analysis. There is a csv file with retrofits Wangoutput_2020-05-28_sheet1.csv which is merged with IN-CORE_2ev3_SetupJoplin_FourInventories_2019-08-06_HUinventory.csv. There are two lines named Population dislocation which is compared to HUA dislocation. That might be the only code to be run in pyIncore.


Discussion / hints:

Fix patch:
diff --git a/CHANGELOG.md b/CHANGELOG.md
index f7df4847b..c88267c6c 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -5,6 +5,18 @@ All notable changes to this project will be documented in this file.
 The format is based on [Keep a Changelog](http://keepachangelog.com/)
 and this project adheres to [Semantic Versioning](http://semver.org/).
 
+## [1.2.0] - 2021-12-15
+
+### Added
+- Vacant household category to population dislocation output [#43](https://github.com/IN-CORE/pyincore/issues/43)
+- Input dataset with target functionality for Joplin empirical restoration [#56](https://github.com/IN-CORE/pyincore/issues/56)
+- Multi-objective retrofit optimization analysis [#19](https://github.com/IN-CORE/pyincore/issues/19)
+
+### Fixed
+- EPFDamage to remove deprecated LS/DS code and properly handle liquefaction [#32](https://github.com/IN-CORE/pyincore/issues/32)
+- Automatic build for pyincore documentation docker [#61](https://github.com/IN-CORE/pyincore/issues/61)
+- Fix get_building_period() method to work with new-format Fragility curves [#15](https://github.com/IN-CORE/pyincore/issues/15)
+
 ## [1.1.0] - 2021-10-27
 
 ### Added
diff --git a/Dockerfile b/Dockerfile
index 2069fe288..5900cb9b8 100644
--- a/Dockerfile
+++ b/Dockerfile
@@ -8,6 +8,7 @@ USER root
 # install packages
 WORKDIR /src
 COPY requirements.txt .
+ENV PATH "$MAMBA_ROOT_PREFIX/bin:$PATH"
 RUN micromamba install -y -n base -c conda-forge \
     sphinx sphinx_rtd_theme \
     -f requirements.txt
diff --git a/docs/docker.sh b/docs/docker.sh
index 77aaf38d9..c58a48922 100755
--- a/docs/docker.sh
+++ b/docs/docker.sh
@@ -13,4 +13,4 @@ BRANCH=${BRANCH:-"$(git rev-parse --abbrev-ref HEAD)"}
 cd ../
 
 # Build docker image
-$DEBUG docker build -t incore/doc/pyincore -f Dockerfile.docs .
+$DEBUG docker build -t incore/doc/pyincore -f Dockerfile .
diff --git a/docs/source/conf.py b/docs/source/conf.py
index ec04374f7..906f0ba5d 100644
--- a/docs/source/conf.py
+++ b/docs/source/conf.py
@@ -33,9 +33,9 @@
 author = ''
 
 # The short X.Y version
-version = '1.1'
+version = '1.2'
 # The full version, including alpha/beta/rc tags
-release = '1.1.0'
+release = '1.2.0'
 
 # -- General configuration ---------------------------------------------------
 
diff --git a/docs/source/modules.rst b/docs/source/modules.rst
index dc3f3e4b9..b9b1e2d7b 100644
--- a/docs/source/modules.rst
+++ b/docs/source/modules.rst
@@ -106,6 +106,11 @@ analyses/montecarlofailureprobability
 ..  autoclass:: montecarlofailureprobability.montecarlofailureprobability.MonteCarloFailureProbability
     :members:
 
+analyses/multiobjectiveretrofitoptimization
+=====================================
+..  autoclass:: multiobjectiveretrofitoptimization.multiobjectiveretrofitoptimization.MultiObjectiveRetrofitOptimization
+    :members:
+
 analyses/nonstructbuildingdamage
 ================================
 ..  autoclass:: nonstructbuildingdamage.nonstructbuildingdamage.NonStructBuildingDamage
diff --git a/pyincore/analyses/buildingdamage/buildingdamage.py b/pyincore/analyses/buildingdamage/buildingdamage.py
index 36db97f7b..99751bb35 100755
--- a/pyincore/analyses/buildingdamage/buildingdamage.py
+++ b/pyincore/analyses/buildingdamage/buildingdamage.py
@@ -168,11 +168,8 @@ def building_damage_analysis_bulk_input(self, buildings, retrofit_strategy, haza
             dmg_probability = dict()
             dmg_interval = dict()
             b_id = b["id"]
-            num_stories = b['properties']['no_stories']
             selected_fragility_set = fragility_sets[b_id]
 
-            building_period = selected_fragility_set.fragility_curves[0].get_building_period(num_stories)
-
             # TODO: Once all fragilities are migrated to new format, we can remove this condition
             if isinstance(selected_fragility_set.fragility_curves[0], FragilityCurve):
                 # Supports multiple demand types in same fragility
@@ -192,6 +189,9 @@ def building_damage_analysis_bulk_input(self, buildings, retrofit_strategy, haza
                 if not AnalysisUtil.do_hazard_values_have_errors(hazard_vals[i]["hazardValues"]):
                     building_args = selected_fragility_set.construct_expression_args_from_inventory(b)
 
+                    building_period = selected_fragility_set.fragility_curves[0].get_building_period(
+                        selected_fragility_set.fragility_curve_parameters, **building_args)
+
                     dmg_probability = selected_fragility_set.calculate_limit_state(
                         hval_dict, **building_args, period=building_period)
                     dmg_interval = selected_fragility_set.calculate_damage_interval(
diff --git a/pyincore/analyses/epfdamage/epfdamage.py b/pyincore/analyses/epfdamage/epfdamage.py
index 12ad9ef21..f1d8dd5bd 100644
--- a/pyincore/analyses/epfdamage/epfdamage.py
+++ b/pyincore/analyses/epfdamage/epfdamage.py
@@ -4,7 +4,6 @@
 # terms of the Mozilla Public License v2.0 which accompanies this distribution,
 # and is available at https://www.mozilla.org/en-US/MPL/2.0/
 
-import collections
 import concurrent.futures
 from itertools import repeat
 
@@ -51,11 +50,8 @@ def run(self):
         if self.get_parameter("use_hazard_uncertainty") is not None:
             use_hazard_uncertainty = self.get_parameter("use_hazard_uncertainty")
 
-        # Liquefaction
-        use_liquefaction = False
-        if self.get_parameter("use_liquefaction") is not None:
-            use_liquefaction = self.get_parameter("use_liquefaction")
-        liq_geology_dataset_id = self.get_parameter("liquefaction_geology_dataset_id")
+        if use_hazard_uncertainty:
+            raise ValueError("Uncertainty is not implemented yet.")
 
         user_defined_cpu = 1
 
@@ -75,10 +71,7 @@ def run(self):
         (ds_results, damage_results) = self.epf_damage_concurrent_future(self.epf_damage_analysis_bulk_input,
                                                                          num_workers,
                                                                          inventory_args, repeat(hazard_type),
-                                                                         repeat(hazard_dataset_id),
-                                                                         repeat(use_hazard_uncertainty),
-                                                                         repeat(use_liquefaction),
-                                                                         repeat(liq_geology_dataset_id))
+                                                                         repeat(hazard_dataset_id))
 
         self.set_result_csv_data("result", ds_results, name=self.get_parameter("result_name"))
         self.set_result_json_data("metadata", damage_results,
@@ -108,30 +101,47 @@ def epf_damage_concurrent_future(self, function_name, num_workers, *args):
 
         return output_ds, output_dmg
 
-    def epf_damage_analysis_bulk_input(self, epfs, hazard_type, hazard_dataset_id,
-                                       use_hazard_uncertainty, use_liquefaction, liq_geology_dataset_id):
+    def epf_damage_analysis_bulk_input(self, epfs, hazard_type, hazard_dataset_id):
         """Run analysis for multiple epfs.
 
         Args:
             epfs (list): Multiple epfs from input inventory set.
             hazard_type (str): A type of hazard exposure (earthquake, tsunami, tornado, or hurricane).
             hazard_dataset_id (str): An id of the hazard exposure.
-            use_hazard_uncertainty (bool):  Hazard uncertainty. True for using uncertainty when computing damage,
-                False otherwise.
-            use_liquefaction (bool): Liquefaction. True for using liquefaction information to modify the damage,
-                False otherwise.
-            liq_geology_dataset_id (str): geology_dataset_id (str): A dataset id for geology dataset for liquefaction.
 
         Returns:
             list: A list of ordered dictionaries with epf damage values and other data/metadata.
 
         """
+
+        use_liquefaction = False
+        liquefaction_available = False
+
         fragility_key = self.get_parameter("fragility_key")
 
         fragility_set = self.fragilitysvc.match_inventory(self.get_input_dataset("dfr3_mapping_set"), epfs,
                                                           fragility_key)
 
+        if hazard_type == "earthquake":
+            liquefaction_fragility_key = self.get_parameter("liquefaction_fragility_key")
+            if self.get_parameter("use_liquefaction") is True:
+                if liquefaction_fragility_key is None:
+                    liquefaction_fragility_key = self.DEFAULT_LIQ_FRAGILITY_KEY
+
+                use_liquefaction = self.get_parameter("use_liquefaction")
+
+                # Obtain the geology dataset
+                geology_dataset_id = self.get_parameter("liquefaction_geology_dataset_id")
+
+                if geology_dataset_id is not None:
+                    fragility_sets_liq = self.fragilitysvc.match_inventory( self.get_input_dataset("dfr3_mapping_set"),
+                                                                            epfs, liquefaction_fragility_key)
+
+                    if fragility_sets_liq is not None:
+                        liquefaction_available = True
+
         values_payload = []
+        values_payload_liq = []
         unmapped_epfs = []
         mapped_epfs = []
         for epf in epfs:
@@ -148,6 +158,17 @@ def epf_damage_analysis_bulk_input(self, epfs, hazard_type, hazard_dataset_id,
                 }
                 values_payload.append(value)
                 mapped_epfs.append(epf)
+
+                if liquefaction_available and epf["id"] in fragility_sets_liq:
+                    fragility_set_liq = fragility_sets_liq[epf["id"]]
+                    demands_liq = fragility_set_liq.demand_types
+                    units_liq = fragility_set_liq.demand_units
+                    value_liq = {
+                        "demands": demands_liq,
+                        "units": units_liq,
+                        "loc": loc
+                    }
+                    values_payload_liq.append(value_liq)
             else:
                 unmapped_epfs.append(epf)
 
@@ -163,6 +184,12 @@ def epf_damage_analysis_bulk_input(self, epfs, hazard_type, hazard_dataset_id,
         else:
             raise ValueError("Missing hazard type.")
 
+        liquefaction_resp = None
+        if liquefaction_available:
+            liquefaction_resp = self.hazardsvc.post_liquefaction_values(hazard_dataset_id,
+                                                                        geology_dataset_id,
+                                                                        values_payload_liq)
+
         ds_results = []
         damage_results = []
 
@@ -187,14 +214,39 @@ def epf_damage_analysis_bulk_input(self, epfs, hazard_type, hazard_dataset_id,
                 limit_states = selected_fragility_set.calculate_limit_state(hval_dict,
                                                                             inventory_type='electric_facility',
                                                                             **epf_args)
+
+                if liquefaction_resp is not None:
+                    fragility_set_liq = fragility_sets_liq[epf["id"]]
+
+                    if isinstance(fragility_set_liq.fragility_curves[0], FragilityCurve):
+                        liq_hazard_vals = AnalysisUtil.update_precision_of_lists(liquefaction_resp[i]["pgdValues"])
+                        liq_demand_types = liquefaction_resp[i]["demands"]
+                        liq_demand_units = liquefaction_resp[i]["units"]
+                        liquefaction_prob = liquefaction_resp[i]['liqProbability']
+
+                        hval_dict_liq = dict()
+
+                        for j, d in enumerate(fragility_set_liq.demand_types):
+                            hval_dict_liq[d] = liq_hazard_vals[j]
+
+                        facility_liq_args = fragility_set_liq.construct_expression_args_from_inventory(epf)
+                        pgd_limit_states = \
+                            fragility_set_liq.calculate_limit_state(
+                                hval_dict_liq, inventory_type="electric_facility",
+                                **facility_liq_args)
+                    else:
+                        raise ValueError("One of the fragilities is in deprecated format. "
+                                         "This should not happen If you are seeing this please report the issue.")
+
+                    limit_states = AnalysisUtil.adjust_limit_states_for_pgd(limit_states, pgd_limit_states)
+
+                dmg_interval = selected_fragility_set.calculate_damage_interval(limit_states,
+                                                                       hazard_type=hazard_type,
+                                                                       inventory_type='electric_facility')
             else:
                 raise ValueError("One of the fragilities is in deprecated format. This should not happen. If you are "
                                  "seeing this please report the issue.")
 
-            dmg_interval = selected_fragility_set.calculate_damage_interval(limit_states,
-                                                                            hazard_type=hazard_type,
-                                                                            inventory_type="electric_facility")
-
             ds_result["guid"] = epf["properties"]["guid"]
             ds_result.update(limit_states)
             ds_result.update(dmg_interval)
@@ -205,11 +257,21 @@ def epf_damage_analysis_bulk_input(self, epfs, hazard_type, hazard_dataset_id,
             damage_result["demandtypes"] = input_demand_types
             damage_result["demandunits"] = input_demand_units
             damage_result["hazardtype"] = hazard_type
-            damage_result["hazardval"] = hazard_val
-            damage_result['liq_fragility_id'] = None
-            damage_result['liqhaztype'] = None
-            damage_result['liqhazval'] = None
-            damage_result['liqprobability'] = None
+            damage_result["hazardvals"] = hazard_val
+
+            if hazard_type == "earthquake" and use_liquefaction is True:
+                if liquefaction_available:
+                    damage_result['liq_fragility_id'] = fragility_sets_liq[epf["id"]].id
+                    damage_result['liqdemandtypes'] = liq_demand_types
+                    damage_result['liqdemandunits'] = liq_demand_units
+                    damage_result['liqhazval'] = liq_hazard_vals
+                    damage_result['liqprobability'] = liquefaction_prob
+                else:
+                    damage_result['liq_fragility_id'] = None
+                    damage_result['liqdemandtypes'] = None
+                    damage_result['liqdemandunits'] = None
+                    damage_result['liqhazval'] = None
+                    damage_result['liqprobability'] = None
 
             ds_results.append(ds_result)
             damage_results.append(damage_result)
@@ -217,75 +279,6 @@ def epf_damage_analysis_bulk_input(self, epfs, hazard_type, hazard_dataset_id,
             i += 1
 
         #############################################################
-        # when there is liquefaction, limit state need to be modified
-        if hazard_type == 'earthquake' and use_liquefaction and liq_geology_dataset_id is not None:
-            liq_fragility_key = self.get_parameter("liquefaction_fragility_key")
-            if liq_fragility_key is None:
-                liq_fragility_key = self.DEFAULT_LIQ_FRAGILITY_KEY
-            liq_fragility_set = self.fragilitysvc.match_inventory(self.get_input_dataset("dfr3_mapping_set"),
-                                                                  epfs, liq_fragility_key)
-
-            liq_values_payload = []
-            for epf in epfs:
-                epf_id = epf["id"]
-                if epf_id in liq_fragility_set:
-                    location = GeoUtil.get_location(epf)
-                    liq_loc = str(location.y) + "," + str(location.x)
-                    liq_demands = liq_fragility_set[epf_id].demand_types
-                    liq_units = liq_fragility_set[epf_id].demand_units
-                    value = {
-                        "demands": liq_demands,
-                        "units": liq_units,
-                        "loc": liq_loc
-                    }
-                    liq_values_payload.append(value)
-
-            liquefaction_vals = self.hazardsvc.post_liquefaction_values(hazard_dataset_id,
-                                                                        geology_dataset_id=liq_geology_dataset_id,
-                                                                        payload=liq_values_payload)
-
-            i = 0
-            for liq_epf in epfs:
-                liq_epf_id = liq_epf["id"]
-                liq_input_demand_types = liquefaction_vals[i]["demands"][0]
-                liq_hazard_val = AnalysisUtil.update_precision(liquefaction_vals[i][liq_input_demand_types][0])
-                std_dev = 0.0
-                if use_hazard_uncertainty:
-                    raise ValueError("Uncertainty Not Implemented!")
-
-                liquefaction_prob = liquefaction_vals[i]['liqProbability']
-
-                selected_liq_fragility = liq_fragility_set[liq_epf_id]
-                pgd_limit_states = \
-                    selected_liq_fragility.calculate_limit_state_w_conversion(liq_hazard_val, std_dev=std_dev)
-
-                # match id and add liqhaztype, liqhazval, liqprobability field as well as rewrite limit
-                # states and dmg_interval
-                for ds_result in ds_results:
-                    if ds_result['guid'] == liq_epf["properties"]['guid']:
-                        # default EPF to be 4 ls and 5 ds
-                        if ['LS_0', 'LS_1', 'LS_2', 'LS_3'] in ds_result.keys():
-                            limit_states = {"LS_0": ds_result['LS_0'],
-                                            "LS_1": ds_result['LS_1'],
-                                            "LS_2": ds_result['LS_2'],
-                                            "LS_3": ds_result['LS_3']}
-                        else:
-                            raise ValueError("Do not support the limit state name")
-
-                        liq_limit_states = AnalysisUtil.adjust_limit_states_for_pgd(limit_states, pgd_limit_states)
-                        liq_dmg_interval = AnalysisUtil.calculate_damage_interval(liq_limit_states)
-
-                        ds_result.update(liq_limit_states)
-                        ds_result.update(liq_dmg_interval)
-
-                for damage_result in damage_results:
-                    if damage_result['guid'] == liq_epf["properties"]['guid']:
-                        damage_result['liq_fragility_id'] = liq_fragility_set.id
-                        damage_result['liqhaztype'] = liq_input_demand_types
-                        damage_result['liqhazval'] = liq_hazard_val
-                        damage_result['liqprobability'] = liquefaction_prob
-
-                i = i + 1
 
         # unmapped
         for epf in unmapped_epfs:
@@ -298,10 +291,12 @@ def epf_damage_analysis_bulk_input(self, epfs, hazard_type, hazard_dataset_id,
             damage_result['demandunits'] = None
             damage_result["hazardtype"] = None
             damage_result['hazardval'] = None
-            damage_result['liq_fragility_id'] = None
-            damage_result['liqhaztype'] = None
-            damage_result['liqhazval'] = None
-            damage_result['liqprobability'] = None
+            if hazard_type == "earthquake" and use_liquefaction is True:
+                damage_result['liq_fragility_id'] = None
+                damage_result['liqdemandtypes'] = None
+                damage_result['liqdemandunits'] = None
+                damage_result['liqhazval'] = None
+                damage_result['liqprobability'] = None
 
             ds_results.append(ds_result)
             damage_results.append(damage_result)
@@ -344,6 +339,12 @@ def get_spec(self):
                     'description': 'Fragility key to use in mapping dataset ()',
                     'type': str
                 },
+                {
+                    'id': 'liquefaction_fragility_key',
+                    'required': False,
+                    'description': 'Fragility key to use in liquefaction mapping dataset',
+                    'type': str
+                },
                 {
                     'id': 'use_liquefaction',
                     'required': False,
diff --git a/pyincore/analyses/joplinempiricalrestoration/joplinempiricalrestoration.py b/pyincore/analyses/joplinempiricalrestoration/joplinempiricalrestoration.py
index 63a1884f9..963fddaa0 100755
--- a/pyincore/analyses/joplinempiricalrestoration/joplinempiricalrestoration.py
+++ b/pyincore/analyses/joplinempiricalrestoration/joplinempiricalrestoration.py
@@ -42,6 +42,12 @@ def run(self):
         building_set = self.get_input_dataset("buildings").get_dataframe_from_shapefile()
         # Building damage dataset
         building_dmg = self.get_input_dataset("building_dmg").get_dataframe_from_csv(low_memory=False)
+        # Building functionality target level
+        building_target_fl = self.get_input_dataset("building_functionality_level")
+        if building_target_fl is not None:
+            building_target_fl = building_target_fl.get_dataframe_from_csv(low_memory=False)
+        else:
+            building_target_fl = None
 
         # merge and filter out archetypes > 5
         building_dmg_all = pd.merge(building_dmg, building_set, how="left", on="guid", copy=True, validate="1:1")
@@ -51,6 +57,15 @@ def run(self):
         building_func = building_func_5[["guid", "LS_0", "LS_1", "LS_2", "haz_expose"]].copy()
         building_func["targetFL"] = target_fl
 
+        if building_target_fl is not None:
+            building_func = pd.merge(building_func, building_target_fl,
+                                     how="left", on="guid", copy=True, validate="1:1")
+            # Replace NaN value from targetFL_y with targetFL_x value
+            building_func["targetFL"] = building_func["targetFL_y"].fillna(building_func["targetFL_x"])
+            # Drop merged columns
+            building_func = building_func.drop(["targetFL_x", "targetFL_y"], axis=1)
+            building_func = building_func.astype({"targetFL": "int64"})
+
         initial_func_level, restoration_days = self.get_restoration_days(seed_i, building_func)
         building_func["initialFL"] = initial_func_level
         building_func["restorDays"] = restoration_days
@@ -66,8 +81,9 @@ def get_restoration_days(self, seed_i, building_func):
 
         Args:
             seed_i (int): Seed for random number generator to ensure replication if run as part
-                    of a stochastic analysis, for example in connection with housing unit allocation analysis.
-            building_dmg (pd.DataFrame): Building damage dataset with guid, limit states, hazard exposure.
+                of a stochastic analysis, for example in connection with housing unit allocation analysis.
+            building_func (pd.DataFrame): Building damage dataset with guid, limit states, hazard exposure
+                and a target level column.
 
         Returns:
             np.array: Initial functionality level based on damage state
@@ -155,6 +171,13 @@ def get_spec(self):
                              "ergo:buildingInventory",
                              "ergo:nsBuildingInventoryDamage",
                              "ergo:nsBuildingInventoryDamageVer2"]
+                },
+                {
+                    "id": "building_functionality_level",
+                    "required": False,
+                    "description": "Functionality level per building. The target level defaults "
+                                   "to target_functionality_level parameter if building not in the dataset",
+                    "type": ["incore:buildingFuncTargetVer1"]
                 }
             ],
             "output_datasets": [
diff --git a/pyincore/analyses/multiobjectiveretrofitoptimization/__init__.py b/pyincore/analyses/multiobjectiveretrofitoptimization/__init__.py
new file mode 100644
index 000000000..aa4ac7f2d
--- /dev/null
+++ b/pyincore/analyses/multiobjectiveretrofitoptimization/__init__.py
@@ -0,0 +1,2 @@
+from pyincore.analyses.multiobjectiveretrofitoptimization.multiobjectiveretrofitoptimization import \
+    MultiObjectiveRetrofitOptimization
diff --git a/pyincore/analyses/multiobjectiveretrofitoptimization/multiobjectiveretrofitoptimization.py b/pyincore/analyses/multiobjectiveretrofitoptimization/multiobjectiveretrofitoptimization.py
new file mode 100644
index 000000000..0fd04657e
--- /dev/null
+++ b/pyincore/analyses/multiobjectiveretrofitoptimization/multiobjectiveretrofitoptimization.py
@@ -0,0 +1,1440 @@
+# Copyright (c) 2021 University of Illinois and others. All rights reserved.
+#
+# This program and the accompanying materials are made available under the
+# terms of the Mozilla Public License v2.0 which accompanies this distribution,
+# and is available at https://www.mozilla.org/en-US/MPL/2.0/
+
+from pyincore import BaseAnalysis
+import pandas as pd
+import numpy as np
+import time
+from pandas import DataFrame
+from pyomo.environ import ConcreteModel, Set, Var, Param, Objective, Constraint
+from pyomo.environ import quicksum, minimize, maximize, NonNegativeReals, Any
+from pyomo.environ import sum_product
+import pyomo.environ as pyo
+from pyomo.opt import SolverStatus, TerminationCondition
+from pyomo.util.infeasible import log_infeasible_constraints
+
+
+class MultiObjectiveRetrofitOptimization(BaseAnalysis):
+    """
+        This analysis computes a series of linear programming models for single- and multi-objective
+        optimization related to the effect of extreme weather on a community in terms of three objective functions.
+        The three objectives used in this program are to minimize economic loss, minimize population dislocation,
+        and maximize building functionality. The computation proceeds by iteratively solving constrained linear
+        models using epsilon steps.
+
+        The output of the computation a collection of optimal resource allocations.
+
+        Contributors
+            | Science: Charles Nicholson, Yunjie Wen
+            | Implementation: Dale Cochran , Tarun Adluri , Jorge Duarte, Santiago Núñez-Corrales, Diego Calderon
+                              and NCSA IN-CORE Dev Team
+
+        Related publications
+
+        Args:
+
+            incore_client (IncoreClient): Service authentication.
+
+        """
+
+    # Column descriptors
+    __Q_col = 'Q_t_hat'
+    __Q_rs_col = 'Q_t_hat_rs'
+    __SC_col = 'Sc'
+    __SC_rs_col = 'Sc_rs'
+
+    __budget_default = 0.2
+
+    def __init__(self, incore_client):
+        super(MultiObjectiveRetrofitOptimization, self).__init__(incore_client)
+
+    def run(self):
+        """Execute the multiobjective retrofit optimization analysis using parameters and input data."""
+        # Read parameters
+        model_solver = self.get_parameter('model_solver')
+        num_epsilon_steps = self.get_parameter('num_epsilon_steps')
+
+        budget_available = self.__budget_default
+        if self.get_parameter('max_budget') != 'default':
+            budget_available = self.get_parameter('budget_available')
+
+        inactive_submodels = []
+
+        in_subm = self.get_parameter('inactive_submodels')
+
+        if in_subm is not None:
+            inactive_submodels = in_subm
+
+        # Perform code scaling
+        scaling_factor = 1.0
+        if self.get_parameter('scale_data'):
+            scaling_factor = self.get_parameter('scaling_factor')
+
+        building_related_data = self.get_input_dataset('building_related_data').get_dataframe_from_csv()
+        strategy_costs = self.get_input_dataset('strategy_costs_data').get_dataframe_from_csv()
+
+        # Convert Z columns to text in both datasets
+        building_related_data['Z'] = building_related_data['Z'].astype(str)
+        strategy_costs['Z'] = strategy_costs['Z'].astype(str)
+
+        self.multiobjective_retrofit_optimization_model(model_solver, num_epsilon_steps, budget_available,
+                                                        scaling_factor, inactive_submodels, building_related_data,
+                                                        strategy_costs)
+
+    def multiobjective_retrofit_optimization_model(self, model_solver, num_epsilon_steps, budget_available,
+                                                   scaling_factor, inactive_submodels, building_related_data,
+                                                   strategy_costs):
+        """Performs the computation of the model.
+
+        Args:
+
+            model_solver (str): model solver to use for analysis
+            num_epsilon_steps (int): number of epsilon values for the multistep optimization algorithm
+            budget_available (float): budget constraint of the optimization analysis
+            scaling_factor (float): scaling factor for Q and Sc matrices
+            inactive_submodels (list): submodels to avoid during the computation
+            building_related_data (pd.DataFrame): building repairs after a disaster event
+            strategy_costs (pd.DataFrame): strategy cost data per building
+
+        """
+        # Setup the stream that will collect data
+        self.ostr = open('pyincore.txt', 'w')
+
+        # Setup the model
+        model, sum_sc = self.configure_model(budget_available, scaling_factor, building_related_data,
+                                             strategy_costs)
+        self.configure_model_objectives(model)
+        print("With constraints model")
+        self.configure_model_retrofit_costs(model) # Suspicious
+
+        # Choose the solver setting
+        if model_solver == "gurobi" or model_solver is None:
+            model_solver_setting = pyo.SolverFactory('gurobi', solver_io="python")
+        else:
+            model_solver_setting = pyo.SolverFactory(model_solver)
+
+        # Solve each model individually
+        print("With individual model")
+        obj_list = self.solve_individual_models(model, model_solver_setting, sum_sc)
+        print("Epsilon values")
+        self.configure_min_max_epsilon_values(model, obj_list, num_epsilon_steps)
+
+        print("Epsilon model")
+        xresults_df, yresults_df = self.solve_epsilon_models(model, model_solver_setting, inactive_submodels)
+
+        df_list = self.compute_optimal_results(inactive_submodels, xresults_df, yresults_df)
+
+        self.set_result_csv_data("optimal_solution_dv_x", df_list[0], name="optimal_solution_dv_x",
+                                 source="dataframe")
+        self.set_result_csv_data("optimal_solution_dv_y", df_list[1], name="optimal_solution_dv_y",
+                                 source="dataframe")
+        return True
+
+    def configure_model(self, budget_available, scaling_factor, building_related_data, strategy_costs):
+        """ Configure the base model to perform the multiobjective optimization.
+
+        Args:
+
+            budget_available (float): available budget
+            scaling_factor (float): value to scale monetary input data
+            building_related_data (DataFrame): table containing building functionality data
+            strategy_costs (DataFrame): table containing retrofit strategy costs data
+
+        Returns:
+            ConcreteModel: a base, parameterized cost/functionality model
+
+        """
+        # Rescale data
+        if scaling_factor != 1.0:
+            building_related_data[self.__Q_col] = \
+                building_related_data[self.__Q_col].map(lambda a: a / scaling_factor)
+            strategy_costs[self.__SC_col] = strategy_costs[self.__SC_col].map(lambda a: a / scaling_factor)
+
+        # Setup pyomo
+        model = ConcreteModel()
+
+        model.Z = Set(initialize=building_related_data.Z.unique())
+        model.S = Set(initialize=building_related_data.S.unique())
+        model.K = Set(initialize=building_related_data.K.unique())
+        model.K_prime = Set(initialize=strategy_costs["K'"].unique())
+
+        zsk = []
+        for y in range(len(building_related_data)):
+            i = building_related_data.loc[y, 'Z']  # Identify the i ∈ Z value.
+            j = building_related_data.loc[y, 'S']  # Identify the j ∈ S value.
+            k = building_related_data.loc[y, 'K']  # Identify the k ∈ K value.
+            zsk.append((i, j, k))  # Add the combination to the list.
+        zsk = sorted(set(zsk), key=zsk.index)  # Convert the list to an ordered set for Pyomo.
+        model.ZSK = Set(initialize=zsk)  # Define and initialize the ZSK set in Pyomo.
+
+        zs = []
+        for y in range(len(building_related_data)):
+            i = building_related_data.loc[y, 'Z']  # Identify the i ∈ Z value.
+            j = building_related_data.loc[y, 'S']  # Identify the j ∈ S value.
+            zs.append((i, j))  # Add the combination to the list.
+        zs = sorted(set(zs), key=zs.index)  # Convert the list to an ordered set for Pyomo.
+        model.ZS = Set(initialize=zs)  # Define and initialize the ZS set in Pyomo.
+
+        kk_prime = []
+        for y in range(len(strategy_costs)):
+            k = strategy_costs.loc[y, 'K']  # Identify the k ∈ K value.
+            k_prime = strategy_costs.loc[y, "K'"]  # Identify the k ∈ K value.
+            kk_prime.append((k, k_prime))  # Add the combination to the list.
+        kk_prime = sorted(set(kk_prime), key=kk_prime.index)  # Convert the list to an ordered set for Pyomo.
+        model.KK_prime = Set(initialize=kk_prime)  # Define and initialize the KK_prime set in Pyomo.
+
+        k_primek = []
+        for y in range(len(strategy_costs)):
+            k = strategy_costs.loc[y, 'K']  # Identify the k ∈ K value.
+            k_prime = strategy_costs.loc[y, "K'"]  # Identify the k ∈ K value.
+            if k_prime <= k:
+                k_primek.append((k_prime, k))  # Add the combination to the list.
+        k_primek = sorted(set(k_primek), key=k_primek.index)  # Convert the list to an ordered set for Pyomo.
+        model.K_primeK = Set(initialize=k_primek)  # Define and initialize the K_primeK set in Pyomo.
+
+        # Define the set of all ZSKK' combinations:
+        zskk_prime = []
+        for y in range(len(strategy_costs)):
+            i = strategy_costs.loc[y, 'Z']  # Identify the i ∈ Z value.
+            j = strategy_costs.loc[y, 'S']  # Identify the j ∈ S value.
+            k = strategy_costs.loc[y, 'K']  # Identify the k ∈ K value.
+            k_prime = strategy_costs.loc[y, "K'"]  # Identify the k ∈ K value.
+            zskk_prime.append((i, j, k, k_prime))  # Add the combination to the list.
+        zskk_prime = sorted(set(zskk_prime), key=zskk_prime.index)  # Convert the list to an ordered set for Pyomo.
+        model.ZSKK_prime = Set(initialize=zskk_prime)  # Define and initialize the ZSKK_prime set in Pyomo.
+        model.zskk_prime = zskk_prime # Use the redundant index for later reference
+
+        ####################################################################################################
+        # DEFINE VARIABLES AND PARAMETERS:
+        ####################################################################################################
+        # Declare the decision variable x_ijk (total # buildings in zone i of structure type
+        # j at code level k after retrofitting):
+        model.x_ijk = Var(model.ZSK, within=NonNegativeReals)
+
+        # Declare the decision variable y_ijkk_prime
+        # (total # buildings in zone i of structure type j retrofitted from code level k to code level k_prime):
+        model.y_ijkk_prime = Var(model.ZSKK_prime, within=NonNegativeReals)
+
+        # Declare economic loss cost parameter l_ijk:
+        model.l_ijk = Param(model.ZSK, within=NonNegativeReals, mutable=True)
+        for y in range(len(building_related_data)):
+            i = building_related_data.loc[y, 'Z']
+            j = building_related_data.loc[y, 'S']
+            k = building_related_data.loc[y, 'K']
+            model.l_ijk[i, j, k] = building_related_data.loc[y, 'l']
+
+        # Declare dislocation parameter d_ijk:
+        model.d_ijk = Param(model.ZSK, within=NonNegativeReals, mutable=True)
+        for y in range(len(building_related_data)):
+            i = building_related_data.loc[y, 'Z']
+            j = building_related_data.loc[y, 'S']
+            k = building_related_data.loc[y, 'K']
+            model.d_ijk[i, j, k] = building_related_data.loc[y, 'd_ijk']
+
+        # Declare the number of buildings parameter b_ijk:
+        model.b_ijk = Param(model.ZSK, within=NonNegativeReals, mutable=True)
+        for y in range(len(building_related_data)):
+            i = building_related_data.loc[y, 'Z']
+            j = building_related_data.loc[y, 'S']
+            k = building_related_data.loc[y, 'K']
+            model.b_ijk[i, j, k] = building_related_data.loc[y, 'b']
+
+        # Declare the building functionality parameter Q_t_hat:
+        model.Q_t_hat = Param(model.ZSK, within=NonNegativeReals, mutable=True)
+        for y in range(len(building_related_data)):
+            i = building_related_data.loc[y, 'Z']
+            j = building_related_data.loc[y, 'S']
+            k = building_related_data.loc[y, 'K']
+            model.Q_t_hat[i, j, k] = building_related_data.loc[y, 'Q_t_hat']
+
+        # Declare the retrofit cost parameter Sc_ijkk':
+        model.Sc_ijkk_prime = Param(model.ZSKK_prime, within=NonNegativeReals, mutable=True)
+        for y in range(len(strategy_costs)):
+            i = strategy_costs.loc[y, 'Z']
+            j = strategy_costs.loc[y, 'S']
+            k = strategy_costs.loc[y, 'K']
+            k_prime = strategy_costs.loc[y, "K'"]
+            model.Sc_ijkk_prime[i, j, k, k_prime] = strategy_costs.loc[y, 'Sc']
+
+        ####################################################################################################
+        # DECLARE THE TOTAL MAX BUDGET AND TOTAL AVAILABLE BUDGET:
+        ####################################################################################################
+        # Define the total max budget based on user's input:
+        model.B = Param(mutable=True, within=NonNegativeReals)
+
+        if budget_available == self.__budget_default:
+            sumSc = quicksum(
+                pyo.value(model.Sc_ijkk_prime[i, j, k, 3]) * pyo.value(model.b_ijk[i, j, k]) for i, j, k in model.ZSK)
+        else:
+            sumSc = budget_available
+        # Define the total available budget based on user's input:
+        model.B = sumSc * budget_available
+
+        return model, sumSc
+
+    def configure_model_objectives(self, model):
+        """ Configure the model by adding objectives
+
+        Args:
+            model (ConcreteModel): a base cost/functionality model
+
+        Returns:
+            ConcreteModel: a model extended with objective functions
+
+        """
+        model.objective_1 = Objective(rule=self.obj_economic, sense=minimize)
+        model.econ_loss = Param(mutable=True, within=NonNegativeReals)  # ,default=10000000000)
+
+        model.objective_2 = Objective(rule=self.obj_dislocation, sense=minimize)
+        model.dislocation = Param(mutable=True, within=NonNegativeReals)  # ,default=30000)
+
+        model.objective_3 = Objective(rule=self.obj_functionality, sense=maximize)
+        model.functionality = Param(mutable=True, within=NonNegativeReals)  # ,default=1)
+
+    def configure_model_retrofit_costs(self, model):
+        model.retrofit_budget_constraint = Constraint(rule=self.retrofit_cost_rule)
+        model.number_buildings_ij_constraint = Constraint(model.ZS, rule=self.number_buildings_ij_rule)
+        model.a = Param(mutable=True)
+        model.c = Param(mutable=True)
+        model.building_level_constraint = Constraint(model.ZSK, rule=self.building_level_rule)
+
+    def solve_individual_models(self, model, model_solver_setting, sum_sc):
+        print("Max Budget: $", sum_sc)
+        print("Available Budget: $", pyo.value(model.B))
+        print("")
+
+        rlist_obj_1 = self.solve_model_1(model, model_solver_setting)
+        rlist_obj_2 = self.solve_model_2(model, model_solver_setting)
+        rlist_obj_3 = self.solve_model_3(model, model_solver_setting)
+
+        values_list = [sum_sc, pyo.value(model.B)] + rlist_obj_1 + rlist_obj_2 + rlist_obj_3
+
+        no_epsilon_constr_init_results = pd.DataFrame(data={
+            "Label": ["Max Budget: $", "Budget (20% of max)", "Economic loss min epsilon (optimal value)",
+                      "Dislocation when optimizing Economic Loss", "Functionality when optimizing Economic Loss",
+                      "Dislocation min epsilon (optimal value)", "Economic Loss when optimizing Dislocation",
+                      "Functionality when optimizing Dislocation", "Functionality max epsilon (optimal value)",
+                      "Economic Loss when optimizing Functionality", "Dislocation when optimizing Functionality"]
+            , "Value": values_list})
+        filename = 'no_epsilon_constr_init_results_' + str(time.strftime("%m-%d-%Y")) + '.csv'
+        no_epsilon_constr_init_results.to_csv(filename)
+
+        return rlist_obj_1 + rlist_obj_2 + rlist_obj_3
+
+    def solve_model_1(self, model, model_solver_setting):
+        starttime = time.time()
+        print("Initial solve for objective function 1 starting.")
+        # Activate objective function 1 (minimize economic loss) and deactivate others:
+        model.objective_1.activate()
+        model.objective_2.deactivate()
+        model.objective_3.deactivate()
+
+        # Solve the model:
+        results = model_solver_setting.solve(model)
+
+        # Save the results if the solver returns an optimal solution:
+        if (results.solver.status == SolverStatus.ok) and (
+                results.solver.termination_condition == TerminationCondition.optimal):
+            self.extract_optimization_results(model)
+            obj_1_min_epsilon = pyo.value(model.objective_1)  # Save the optimal economic loss value.
+            obj_2_value_1 = pyo.value(model.dislocation)  # Save the dislocation value when optimizing economic loss.
+            obj_3_value_1 = pyo.value(
+                model.functionality)  # Save the functionality value when optimizing economic loss.
+            print("Initial solve for objective function 1 complete.")
+            print("Economic Loss: ", pyo.value(model.econ_loss),
+                  "Dislocation: ", pyo.value(model.dislocation),
+                  "Functionality: ", pyo.value(model.functionality))
+
+            print("Economic loss min epsilon (optimal value):", obj_1_min_epsilon)
+            print("Dislocation when optimizing Economic Loss:", obj_2_value_1)
+            print("Functionality when optimizing Economic Loss:", obj_3_value_1)
+            print("")
+            results_list = [obj_1_min_epsilon, obj_2_value_1, obj_3_value_1]
+        else:
+            print("Not Optimal")
+            results_list = [0, 0, 0]
+
+        endtime = time.time()
+        elapsedtime = endtime - starttime
+        print("Elapsed time for initial obj 1 solve: ", elapsedtime)
+        print("")
+
+        return results_list
+
+    def solve_model_2(self, model, model_solver_setting):
+        starttime = time.time()
+        print("Initial solve for objective function 2 starting.")
+        # Activate objective function 2 (minimize dislocation) and deactivate others:
+        model.objective_1.deactivate()
+        model.objective_2.activate()
+        model.objective_3.deactivate()
+
+        # Solve the model:
+        results = model_solver_setting.solve(model)
+
+        # Save the results if the solver returns an optimal solution:
+        if (results.solver.status == SolverStatus.ok) and (
+                results.solver.termination_condition == TerminationCondition.optimal):
+            self.extract_optimization_results(model)
+            obj_2_min_epsilon = pyo.value(model.objective_2)  # Save the optimal dislocation value.
+            obj_1_value_2 = pyo.value(model.econ_loss)  # Save the economic loss value when optimizing dislocation.
+            obj_3_value_2 = pyo.value(model.functionality)  # Save the functionality value when optimizing dislocation.
+            print("Initial solve for objective function 2 complete.")
+            print("Economic Loss: ", pyo.value(model.econ_loss),
+                  "Dislocation: ", pyo.value(model.dislocation),
+                  "Functionality: ", pyo.value(model.functionality))
+
+            print("Dislocation min epsilon (optimal value):", obj_2_min_epsilon)
+            print("Economic Loss when optimizing Dislocation:", obj_1_value_2)
+            print("Functionality when optimizing Dislocation:", obj_3_value_2)
+            print("")
+            results_list = [obj_2_min_epsilon, obj_1_value_2, obj_3_value_2]
+        else:
+            print("Not Optimal")
+            results_list = [0, 0, 0]
+
+        endtime = time.time()
+        elapsedtime = endtime - starttime
+        print("Elapsed time for initial obj 2 solve: ", elapsedtime)
+        print("")
+
+        return results_list
+
+    def solve_model_3(self, model, model_solver_setting):
+        starttime = time.time()
+        print("Initial solve for objective function 3 starting.")
+        # Activate objective function 3 (maximize functionality) and deactivate others:
+        model.objective_1.deactivate()
+        model.objective_2.deactivate()
+        model.objective_3.activate()
+
+        # Solve the model:
+        results = model_solver_setting.solve(model)
+
+        # Save the results if the solver returns an optimal solution:
+        if (results.solver.status == SolverStatus.ok) and (
+                results.solver.termination_condition == TerminationCondition.optimal):
+            self.extract_optimization_results(model)
+            obj_3_max_epsilon = pyo.value(model.objective_3)  # Save the optimal functionality value.
+            obj_1_value_3 = pyo.value(model.econ_loss)  # Save the economic loss value when optimizing functionality.
+            obj_2_value_3 = pyo.value(model.dislocation)  # Save the dislocation value when optimizing functionality.
+            print("Initial solve for objective function 3 complete.")
+            print("Economic Loss: ", pyo.value(model.econ_loss),
+                  "Dislocation: ", pyo.value(model.dislocation),
+                  "Functionality: ", pyo.value(model.functionality))
+
+            print("Functionality max epsilon (optimal value):", obj_3_max_epsilon)
+            print("Economic Loss when optimizing Functionality:", obj_1_value_3)
+            print("Dislocation when optimizing Functionality:", obj_2_value_3)
+            results_list = [obj_3_max_epsilon, obj_1_value_3, obj_2_value_3]
+        else:
+            print("Not Optimal")
+            results_list = [0, 0, 0]
+
+        endtime = time.time()
+        elapsedtime = endtime - starttime
+        print("Elapsed time for initial obj 3 solve: ", elapsedtime)
+
+        return results_list
+
+    def configure_min_max_epsilon_values(self, model, objs, num_epsilon_steps):
+        # Define positions in the incoming list
+        __obj_1_value_2_pos = 4
+        __obj_1_value_3_pos = 7
+        __obj_1_min_epsilon_pos = 0
+        __obj_2_value_1_pos = 1
+        __obj_2_value_3_pos = 8
+        __obj_2_min_epsilon_pos = 3
+        __obj_3_value_1_pos = 2
+        __obj_3_value_2_pos = 5
+        __obj_3_max_epsilon_pos = 6
+        model.econ_loss_max = Param(within=NonNegativeReals, initialize=max(objs[__obj_1_value_2_pos],
+                                                                            objs[__obj_1_value_3_pos]))
+        model.econ_loss_min = Param(within=NonNegativeReals, initialize=objs[__obj_1_min_epsilon_pos])
+        model.dislocation_max = Param(within=NonNegativeReals, initialize=max(objs[__obj_2_value_1_pos],
+                                                                              objs[__obj_2_value_3_pos]))
+        model.dislocation_min = Param(within=NonNegativeReals, initialize=objs[__obj_2_min_epsilon_pos])
+        model.functionality_max = Param(within=NonNegativeReals, initialize=objs[__obj_3_max_epsilon_pos])
+        model.functionality_min = Param(within=NonNegativeReals, initialize=min(objs[__obj_3_value_1_pos],
+                                                                                objs[__obj_3_value_2_pos]))
+
+        model.econ_loss_step = Param(within=NonNegativeReals,
+                                     initialize=(pyo.value(model.econ_loss_max) - pyo.value(model.econ_loss_min)) * (
+                                                 1 / (num_epsilon_steps - 1)))
+        model.dislocation_step = Param(within=NonNegativeReals, initialize=(pyo.value(
+            model.dislocation_max) - pyo.value(model.dislocation_min)) * (1 / (num_epsilon_steps - 1)))
+        model.functionality_step = Param(within=NonNegativeReals, initialize=(pyo.value(
+            model.functionality_max) - pyo.value(model.functionality_min)) * (1 / (num_epsilon_steps - 1)))
+
+    def solve_epsilon_models(self, model, model_solver_setting, inactive_submodels):
+        xresults_df = pd.DataFrame()
+        yresults_df = pd.DataFrame()
+
+        if 1 not in inactive_submodels:
+           self.solve_epsilon_model_1(model, model_solver_setting)
+
+        if 2 not in inactive_submodels:
+            self.solve_epsilon_model_2(model, model_solver_setting)
+
+        if 3 not in inactive_submodels:
+            self.solve_epsilon_model_3(model, model_solver_setting)
+
+        if 4 not in inactive_submodels:
+            self.solve_epsilon_model_4(model, model_solver_setting)
+
+        if 5 not in inactive_submodels:
+            self.solve_epsilon_model_5(model, model_solver_setting)
+
+        if 6 not in inactive_submodels:
+            self.solve_epsilon_model_6(model, model_solver_setting)
+
+        if 7 not in inactive_submodels:
+            xresults_df, yresults_df = self.solve_epsilon_model_7(model, model_solver_setting, xresults_df,
+                                                                      yresults_df)
+
+        if 8 not in inactive_submodels:
+            xresults_df, yresults_df = self.solve_epsilon_model_8(model, model_solver_setting, xresults_df,
+                                                                      yresults_df)
+
+        if 9 not in inactive_submodels:
+            xresults_df, yresults_df = self.solve_epsilon_model_9(model, model_solver_setting, xresults_df,
+                                                                      yresults_df)
+
+        return xresults_df, yresults_df
+
+    def solve_epsilon_model_1(self, model, model_solver_setting):
+        starttime = time.time()
+        print("****OPTIMIZING ECONOMIC LOSS SUBJECT TO POPULATION DISLOCATION EPSILON CONSTRAINTS****")
+        # Activate objective function 1 and deactivate others:
+        model.objective_1.activate()
+        model.objective_2.deactivate()
+        model.objective_3.deactivate()
+
+        # Add objective function 2 as epsilon constraint of objective function 1:
+        # Dataframe to store optimization results:
+        obj_1_2_epsilon_results = pd.DataFrame(
+            columns=['Economic Loss(Million Dollars)', 'Dislocation Value', 'Functionality Value'])
+
+        # Parameter for objective 2 (dislocation) epsilon value:
+        model.obj_2_e = Param(mutable=True, within=NonNegativeReals)
+        # For each dislocation epsilon value (starting at min) set as constraint and solve model:
+        counter = 0
+        for e in np.arange(pyo.value(model.dislocation_min), pyo.value(model.dislocation_max) + 0.000001,
+                           pyo.value(model.dislocation_step)):
+            counter += 1
+            print("Step ", counter, ": ", e)
+
+            model.obj_2_e = e  # Set the model parameter to the epsilon value.
+            model.add_component("objective_2_constraint",
+                                Constraint(expr=sum_product(model.d_ijk, model.x_ijk) <= pyo.value(
+                                    model.obj_2_e)))  # Add the epsilon constraint.
+            # Solve the model:
+            results = model_solver_setting.solve(model)
+
+            # Save the results if the solver returns an optimal solution:
+            if (results.solver.status == SolverStatus.ok) and (
+                    results.solver.termination_condition == TerminationCondition.optimal):
+                self.extract_optimization_results(model)
+                # Add objective (economic loss), dislocation, and functionality values to results dataframe:
+                obj_1_2_epsilon_results.loc[counter - 1, 'Economic Loss(Million Dollars)'] = pyo.value(
+                    model.econ_loss)  # Save the optimal economic loss value.
+                obj_1_2_epsilon_results.loc[counter - 1, 'Dislocation Value'] = pyo.value(
+                    model.dislocation)  # Save the resulting dislocation value.
+                obj_1_2_epsilon_results.loc[counter - 1, 'Functionality Value'] = pyo.value(
+                    model.functionality)  # Save the resulting functionality value.
+            else:
+                print(results.solver.termination_condition)
+                log_infeasible_constraints(model)
+
+            # Remove the given epsilon constraint from the model now that optimization is complete:
+            model.del_component(model.objective_2_constraint)
+
+        # Remove the dislocation epsilon parameter from the model:
+        model.del_component(model.obj_2_e)
+
+        # Display and save the results of optimizing economic loss subject to dislocation epsilon constraints:
+        print(obj_1_2_epsilon_results)
+        filename = 'obj_1_2_epsilon_results_' + str(time.strftime("%m-%d-%Y")) + '.csv'
+        obj_1_2_epsilon_results.to_csv(filename)
+        endtime = time.time()
+        elapsedtime = endtime - starttime
+        print("Elapsed time: ", elapsedtime)
+
+    def solve_epsilon_model_2(self, model, model_solver_setting):
+        starttime = time.time()
+        print("****OPTIMIZING ECONOMIC LOSS SUBJECT TO BUILDING FUNCTIONALITY EPSILON CONSTRAINTS****")
+        # Activate objective function 1 and deactivate others:
+        model.objective_1.activate()
+        model.objective_2.deactivate()
+        model.objective_3.deactivate()
+
+        # Add objective function 3 as epsilon constraint of objective function 1:
+        # Dataframe to store optimization results:
+        obj_1_3_epsilon_results = pd.DataFrame(
+            columns=['Economic Loss(Million Dollars)', 'Dislocation Value', 'Functionality Value'])
+
+        # Parameter for objective 3 (functionality) epsilon value:
+        model.obj_3_e = Param(mutable=True, within=NonNegativeReals)
+        counter = 0
+        # For each functionality epsilon value (starting at min) set as constraint and solve model:
+        # Adding 0.0000000000001 to the maximum value allows np.arange() to include the maximum functionality.
+        for e in np.arange(pyo.value(model.functionality_min),
+                           pyo.value(model.functionality_max) + 0.000000000000000001,
+                           pyo.value(model.functionality_step)):
+            counter += 1
+            print("Step ", counter, ": ", e)
+
+            model.obj_3_e = e  # Set the model parameter to the epsilon value.
+            model.add_component("objective_3_constraint",
+                                Constraint(expr=sum_product(model.Q_t_hat, model.x_ijk) >= pyo.value(
+                                    model.obj_3_e)))  # Add the epsilon constraint.
+            # Solve the model:
+            results = model_solver_setting.solve(model)
+
+            # Save the results if the solver returns an optimal solution:
+            if (results.solver.status == SolverStatus.ok) and (
+                    results.solver.termination_condition == TerminationCondition.optimal):
+                self.extract_optimization_results(model)
+                # Add objective (economic loss), dislocation, and functionality values to results dataframe:
+                obj_1_3_epsilon_results.loc[counter - 1, 'Economic Loss(Million Dollars)'] = pyo.value(
+                    model.econ_loss)  # Save the optimal economic loss value.
+                obj_1_3_epsilon_results.loc[counter - 1, 'Dislocation Value'] = pyo.value(
+                    model.dislocation)  # Save the resulting dislocation value.
+                obj_1_3_epsilon_results.loc[counter - 1, 'Functionality Value'] = pyo.value(
+                    model.functionality)  # Save the resulting functionality value.
+            else:
+                print(results.solver.termination_condition)
+                log_infeasible_constraints(model)
+
+            # Remove the given epsilon constraint from the model now that optimization is complete:
+            model.del_component(model.objective_3_constraint)
+
+        # Remove the functionality epsilon parameter from the model:
+        model.del_component(model.obj_3_e)
+
+        # Display and save the results of optimizing economic loss subject to functionality epsilon constraints:
+        print(obj_1_3_epsilon_results)
+        filename = 'obj_1_3_epsilon_results_' + str(time.strftime("%m-%d-%Y")) + '.csv'
+        obj_1_3_epsilon_results.to_csv(filename)
+        endtime = time.time()
+        elapsedtime = endtime - starttime
+        print("Elapsed time: ", elapsedtime)
+
+    def solve_epsilon_model_3(self, model, model_solver_setting):
+        starttime = time.time()
+        print("****OPTIMIZING POPULATION DISLOCATION SUBJECT TO ECONOMIC LOSS EPSILON CONSTRAINTS****")
+        # Activate objective function 2 and deactivate others:
+        model.objective_1.deactivate()
+        model.objective_2.activate()
+        model.objective_3.deactivate()
+
+        # Add objective function 1 as epsilon constraint of objective function 2:
+        # Dataframe to store optimization results:
+        obj_2_1_epsilon_results = pd.DataFrame(
+            columns=['Economic Loss(Million Dollars)', 'Dislocation Value', 'Functionality Value'])
+
+        # Parameter for objective 1 (economic loss) epsilon value:
+        model.obj_1_e = Param(mutable=True, within=NonNegativeReals)
+        # For each economic loss epsilon value (starting at min) set as constraint and solve model:
+        # Adding 1 to the maximum value allows np.arange() to include the maximum economic loss.
+        counter = 0
+        for e in np.arange(pyo.value(model.econ_loss_min), pyo.value(model.econ_loss_max) + 0.1,
+                           pyo.value(model.econ_loss_step)):
+            counter += 1
+            print("Step ", counter, ": ", e)
+
+            model.obj_1_e = e  # Set the model parameter to the epsilon value.
+            model.add_component("objective_1_constraint",
+                                Constraint(expr=sum_product(model.l_ijk, model.x_ijk) <= pyo.value(
+                                    model.obj_1_e)))  # Set epsilon constraint.
+            # Solve the model:
+            results = model_solver_setting.solve(model)
+
+            # Save the results if the solver returns an optimal solution:
+            if (results.solver.status == SolverStatus.ok) and (
+                    results.solver.termination_condition == TerminationCondition.optimal):
+                self.extract_optimization_results(model)
+                # Add objective (dislocation), economic loss, and functionality values to results dataframe:
+                obj_2_1_epsilon_results.loc[counter - 1, 'Dislocation Value'] = pyo.value(
+                    model.dislocation)  # Save the optimal dislocation value.
+                obj_2_1_epsilon_results.loc[counter - 1, 'Economic Loss(Million Dollars)'] = pyo.value(
+                    model.econ_loss)  # Save the resulting economic loss value.
+                obj_2_1_epsilon_results.loc[counter - 1, 'Functionality Value'] = pyo.value(
+                    model.functionality)  # Save the resulting functionality value.
+            else:
+                print(results.solver.termination_condition)
+                log_infeasible_constraints(model)
+
+            # Remove the given epsilon constraint from the model now that optimization is complete:
+            model.del_component(model.objective_1_constraint)
+
+        # Remove the functionality epsilon parameter from the model:
+        model.del_component(model.obj_1_e)
+
+        # Display the results of optimizing dislocation subject to economic loss epsilon constraints:
+        print(obj_2_1_epsilon_results)
+        filename = 'obj_2_1_epsilon_results_' + str(time.strftime("%m-%d-%Y")) + '.csv'
+        obj_2_1_epsilon_results.to_csv(filename)
+        endtime = time.time()
+        elapsedtime = endtime - starttime
+        print("Elapsed time: ", elapsedtime)
+
+    def solve_epsilon_model_4(self, model, model_solver_setting):
+        starttime = time.time()
+        print("****OPTIMIZING POPULATION DISLOCATION SUBJECT TO BUILDING FUNCTIONALITY EPSILON CONSTRAINTS****")
+        # Activate objective function 2 and deactivate others:
+        model.objective_1.deactivate()
+        model.objective_2.activate()
+        model.objective_3.deactivate()
+
+        # Add objective function 3 as epsilon constraint of objective function 2:
+        # Dataframe to store optimization results:
+        obj_2_3_epsilon_results = pd.DataFrame(
+            columns=['Economic Loss(Million Dollars)', 'Dislocation Value', 'Functionality Value'])
+
+        # Parameter for objective 3 (functionality) epsilon value:
+        model.obj_3_e = Param(mutable=True, within=NonNegativeReals)
+        # For each functionality epsilon value (starting at min) set as constraint and solve model:
+        # Adding 0.0000000000001 to the maximum value allows np.arange() to include the maximum functionality.
+        counter = 0
+        for e in np.arange(pyo.value(model.functionality_min), pyo.value(model.functionality_max) + 0.0000000000001,
+                           pyo.value(model.functionality_step)):
+            counter += 1
+            print("Step ", counter, ": ", e)
+
+            model.obj_3_e = e  # Set the model parameter to the epsilon value.
+            model.add_component("objective_3_constraint",
+                                Constraint(expr=sum_product(model.Q_t_hat, model.x_ijk) >= (
+                                    pyo.value(model.obj_3_e))))  # Set the epsilon constraint.
+            # Solve the model:
+            results = model_solver_setting.solve(model)
+
+            # Save the results if the solver returns an optimal solution:
+            if (results.solver.status == SolverStatus.ok) and (
+                    results.solver.termination_condition == TerminationCondition.optimal):
+                self.extract_optimization_results(model)
+                # Add objective (dislocation), economic loss, and functionality values to results dataframe:
+                obj_2_3_epsilon_results.loc[counter - 1, 'Dislocation Value'] = pyo.value(
+                    model.dislocation)  # Save the optimal dislocation value.
+                obj_2_3_epsilon_results.loc[counter - 1, 'Economic Loss(Million Dollars)'] = pyo.value(
+                    model.econ_loss)  # Save the resulting economic loss value.
+                obj_2_3_epsilon_results.loc[counter - 1, 'Functionality Value'] = pyo.value(
+                    model.functionality)  # Save the resulting functionality value.
+            else:
+                print(results.solver.termination_condition)
+                log_infeasible_constraints(model)
+
+            # Remove the given epsilon constraint from the model now that optimization is complete:
+            model.del_component(model.objective_3_constraint)
+
+        # Remove the functionality epsilon parameter from the model:
+        model.del_component(model.obj_3_e)
+
+        # Display the results of optimizing dislocation subject to functionality epsilon constraints:
+        print(obj_2_3_epsilon_results)
+        filename = 'obj_2_3_epsilon_results_' + str(time.strftime("%m-%d-%Y")) + '.csv'
+        obj_2_3_epsilon_results.to_csv(filename)
+        endtime = time.time()
+        elapsedtime = endtime - starttime
+        print("Elapsed time: ", elapsedtime)
+
+    def solve_epsilon_model_5(self, model, model_solver_setting):
+        starttime = time.time()
+        print("****OPTIMIZING BUILDING FUNCTIONALITY SUBJECT TO ECONOMIC LOSS EPSILON CONSTRAINTS****")
+        # Activate objective function 3 and deactivate others:
+        model.objective_1.deactivate()
+        model.objective_2.deactivate()
+        model.objective_3.activate()
+
+        # Add objective function 1 as epsilon constraint of objective function 3:
+        # Dataframe to store optimization results:
+        obj_3_1_epsilon_results = pd.DataFrame(
+            columns=['Economic Loss(Million Dollars)', 'Dislocation Value', 'Functionality Value'])
+
+        # Parameter for objective 1 (economic loss) epsilon value:
+        model.obj_1_e = Param(mutable=True, within=NonNegativeReals)
+        # For each economic loss epsilon value (starting at min) set as constraint and solve model:
+        # Adding 1 to the maximum value allows np.arange() to include the maximum economic loss.
+        counter = 0
+        for e in np.arange(pyo.value(model.econ_loss_min), pyo.value(model.econ_loss_max) + 0.000001,
+                           pyo.value(model.econ_loss_step)):
+            counter += 1
+            print("Step ", counter, ": ", e)
+
+            model.obj_1_e = e
+            model.add_component("objective_1_constraint",
+                                Constraint(expr=sum_product(model.l_ijk, model.x_ijk) <= pyo.value(
+                                    model.obj_1_e)))  # Set the epsilon constraint.
+            # Solve the model:
+            results = model_solver_setting.solve(model)
+
+            # Save the results if the solver returns an optimal solution:
+            if (results.solver.status == SolverStatus.ok) and (
+                    results.solver.termination_condition == TerminationCondition.optimal):
+                self.extract_optimization_results(model)
+                # Add objective (functionality), economic loss, and dislocation values to results dataframe:
+                obj_3_1_epsilon_results.loc[counter - 1, 'Functionality Value'] = pyo.value(
+                    model.functionality)  # Save the optimal functionality value.
+                obj_3_1_epsilon_results.loc[counter - 1, 'Economic Loss(Million Dollars)'] = pyo.value(
+                    model.econ_loss)  # Save the resulting economic loss value.
+                obj_3_1_epsilon_results.loc[counter - 1, 'Dislocation Value'] = pyo.value(
+                    model.dislocation)  # Save the resulting dislocation value.
+            else:
+                print(results.solver.termination_condition)
+                log_infeasible_constraints(model)
+
+            # Remove the given epsilon constraint from the model now that optimization is complete:
+            model.del_component(model.objective_1_constraint)
+
+        # Remove the economic loss epsilon parameter from the model:
+        model.del_component(model.obj_1_e)
+
+        # Display the results of optimizing functionality subject to economic loss epsilon constraints:
+        print(obj_3_1_epsilon_results)
+        filename = 'obj_3_1_epsilon_results_' + str(time.strftime("%m-%d-%Y")) + '.csv'
+        obj_3_1_epsilon_results.to_csv(filename)
+        endtime = time.time()
+        elapsedtime = endtime - starttime
+        print("Elapsed time: ", elapsedtime)
+
+    def solve_epsilon_model_6(self, model, model_solver_setting):
+        starttime = time.time()
+        print("****OPTIMIZING BUILDING FUNCTIONALITY SUBJECT TO POPULATION DISLOCATION EPSILON CONSTRAINTS****")
+        # Activate objective function 3 (functionality) and deactivate others:
+        model.objective_1.deactivate()
+        model.objective_2.deactivate()
+        model.objective_3.activate()
+
+        # Add objective function 2 as epsilon constraint of objective function 3:
+        # Dataframe to store optimization results:
+        obj_3_2_epsilon_results = pd.DataFrame(
+            columns=['Economic Loss(Million Dollars)', 'Dislocation Value', 'Functionality Value'])
+
+        # Parameter for objective 2 (dislocation) epsilon value:
+        model.obj_2_e = Param(mutable=True, within=NonNegativeReals)
+        # For each dislocation epsilon value (starting at min) set as constraint and solve model:
+        # Adding 1 to the maximum value allows np.arange() to include the maximum dislocation.
+        counter = 0
+        for e in np.arange(pyo.value(model.dislocation_min), pyo.value(model.dislocation_max) + 0.000001,
+                           pyo.value(model.dislocation_step)):
+            counter += 1
+            print("Step ", counter, ": ", e)
+
+            model.obj_2_e = e
+            model.add_component("objective_2_constraint",
+                                Constraint(expr=sum_product(model.d_ijk, model.x_ijk) <= pyo.value(
+                                    model.obj_2_e)))  # Set the epsilon constraint.
+            # Solve the model:
+            results = model_solver_setting.solve(model)
+
+            # Save the results if the solver returns an optimal solution:
+            if (results.solver.status == SolverStatus.ok) and (
+                    results.solver.termination_condition == TerminationCondition.optimal):
+                self.extract_optimization_results(model)
+                # Add objective (functionality), economic loss, and dislocation values to results dataframe:
+                obj_3_2_epsilon_results.loc[counter - 1, 'Functionality Value'] = pyo.value(
+                    model.functionality)  # Save the optimal functionality value.
+                obj_3_2_epsilon_results.loc[counter - 1, 'Economic Loss(Million Dollars)'] = pyo.value(
+                    model.econ_loss)  # Save the resulting economic loss value.
+                obj_3_2_epsilon_results.loc[counter - 1, 'Dislocation Value'] = pyo.value(
+                    model.dislocation)  # Save the resulting dislocation value.
+            else:
+                print(results.solver.termination_condition)
+                log_infeasible_constraints(model)
+
+            # Remove the given epsilon constraint from the model now that optimization is complete:
+            model.del_component(model.objective_2_constraint)
+
+        # Remove the economic loss epsilon parameter from the model:
+        model.del_component(model.obj_2_e)
+
+        # Display the results of optimizing functionality subject to dislocation epsilon constraints:
+        print(obj_3_2_epsilon_results)
+        filename = 'obj_3_2_epsilon_results_' + str(time.strftime("%m-%d-%Y")) + '.csv'
+        obj_3_2_epsilon_results.to_csv(filename)
+        endtime = time.time()
+        elapsedtime = endtime - starttime
+        print("Elapsed time: ", elapsedtime)
+
+    def solve_epsilon_model_7(self, model, model_solver_setting, xresults_df, yresults_df):
+        starttime = time.time()
+        print(
+            "****OPTIMIZING ECONOMIC LOSS SUBJECT TO POPULATION DISLOCATION "
+            "AND BUILDING FUNCTIONALITY EPSILON CONSTRAINTS****")
+        model.objective_1.activate()
+        model.objective_2.deactivate()
+        model.objective_3.deactivate()
+
+        # Add objective functions 2 and 3 as epsilon constraints of objective function 1:
+        # Dataframe to store optimization results:
+        obj_1_23_epsilon_results = pd.DataFrame(
+            columns=['Economic Loss(Million Dollars)', 'Dislocation Value', 'Functionality Value'])
+
+        # Parameter for objective 2 (dislocation) and objective 3 (functionality) epsilon values:
+        model.obj_2_e = Param(mutable=True, within=NonNegativeReals)
+        model.obj_3_e = Param(mutable=True, within=NonNegativeReals)
+
+        # Generate datasets for this epsilon analysis
+        epsilon7_xresult_df = pd.DataFrame()
+        epsilon7_yresult_df = pd.DataFrame()
+
+        counter = 0
+        for e in np.arange(pyo.value(model.dislocation_min), pyo.value(model.dislocation_max) + 0.000001,
+                           pyo.value(model.dislocation_step)):
+            model.obj_2_e = e
+            model.add_component("objective_2_constraint",
+                                Constraint(expr=sum_product(model.d_ijk, model.x_ijk) <= pyo.value(model.obj_2_e)))
+            for e2 in np.arange(pyo.value(model.functionality_min),
+                                pyo.value(model.functionality_max) + 0.0000000000001,
+                                pyo.value(model.functionality_step)):
+                counter += 1
+                model.obj_3_e = e2
+                print("Step ", counter, " e: ", e, "  e2:", e2)
+                #             obj_1_23_epsilon_results.loc[counter-1,'Dislocation Epsilon']=e
+                #             obj_1_23_epsilon_results.loc[counter-1,'Functionality Epsilon']=e2
+
+                model.add_component("objective_3_constraint",
+                                    Constraint(
+                                        expr=sum_product(model.Q_t_hat, model.x_ijk) >= pyo.value(model.obj_3_e)))
+                # Solve the model:
+                results = model_solver_setting.solve(model)
+
+                # Save the results if the solver returns an optimal solution:
+                if (results.solver.status == SolverStatus.ok) and (
+                        results.solver.termination_condition == TerminationCondition.optimal):
+                    self.extract_optimization_results(model)
+                    budget_used = quicksum(pyo.value(model.y_ijkk_prime[i, j, k, k_prime]) * pyo.value(
+                        model.Sc_ijkk_prime[i, j, k, k_prime]) for (i, j, k, k_prime) in model.ZSKK_prime)
+                    percent_budget_used = (budget_used / pyo.value(
+                        model.B)) * 100  # Record percentage of available budget used.
+                    # Add objective (economic loss), dislocation, and functionality values to results dataframe:
+                    # TODO: optimize this search
+                    obj_1_23_epsilon_results.loc[counter - 1, 'Functionality Value'] = pyo.value(model.functionality)
+                    obj_1_23_epsilon_results.loc[counter - 1, 'Economic Loss(Million Dollars)'] = pyo.value(
+                        model.econ_loss)
+                    obj_1_23_epsilon_results.loc[counter - 1, 'Dislocation Value'] = pyo.value(model.dislocation)
+
+                    # Extract results per each variable and convert to dataframe
+                    x_data: dict = model.x_ijk.extract_values()
+                    newx_df = self.assemble_dataframe_from_solution("x_ijk", x_data, counter)
+
+                    y_data: dict = model.y_ijkk_prime.extract_values()
+                    newy_df = self.assemble_dataframe_from_solution("y_ijkk_prime", y_data, counter)
+
+                    # Append to local analysis result
+                    epsilon7_xresult_df = epsilon7_xresult_df.append(newx_df)
+                    epsilon7_yresult_df = epsilon7_yresult_df.append(newy_df)
+                else:
+                    print(results.solver.termination_condition)
+                    log_infeasible_constraints(model)
+
+                # Remove the given epsilon constraint from the model now that optimization is complete:
+                model.del_component(model.objective_3_constraint)
+
+            # Remove the given epsilon constraint from the model now that optimization is complete:
+            model.del_component(model.objective_2_constraint)
+
+        # Remove the economic loss epsilon parameter from the model:
+        model.del_component(model.obj_2_e)
+        model.del_component(model.obj_3_e)
+
+        print(obj_1_23_epsilon_results)
+
+        # Drop rows with infeasible results:
+        initial_length = len(obj_1_23_epsilon_results)
+        obj_1_23_data = obj_1_23_epsilon_results.dropna(axis=0, how='any')
+        print("Infeasible rows dropped: ", initial_length - len(obj_1_23_data), " rows.")
+
+        # Save results:
+        filename = '7_epsilon_results.csv'
+        obj_1_23_data.index += 1
+        obj_1_23_data.to_csv(filename)
+
+        endtime = time.time()
+        elapsedtime = endtime - starttime
+        print("Elapsed time: ", elapsedtime)
+
+        epsilon7_xresult_df['Epsilon'] = 7
+        epsilon7_yresult_df['Epsilon'] = 7
+
+        return xresults_df.append(epsilon7_xresult_df), yresults_df.append(epsilon7_yresult_df)
+
+    def solve_epsilon_model_8(self, model, model_solver_setting, xresults_df, yresults_df):
+        starttime = time.time()
+        print(
+            "****OPTIMIZING POPULATION DISLOCATION SUBJECT TO "
+            "ECONOMIC LOSS AND BUILDING FUNCTIONALITY EPSILON CONSTRAINTS****")
+        model.objective_1.deactivate()
+        model.objective_2.activate()
+        model.objective_3.deactivate()
+
+        # Add objective functions 1 and 3 as epsilon constraints of objective function 2:
+        # Dataframe to store optimization results:
+        obj_2_13_epsilon_results = pd.DataFrame(
+            columns=['Economic Loss(Million Dollars)', 'Dislocation Value', 'Functionality Value'])
+
+        # Parameter for objective 1 (economic loss) and objective 3 (functionality) epsilon values:
+        model.obj_1_e = Param(mutable=True, within=NonNegativeReals)
+        model.obj_3_e = Param(mutable=True, within=NonNegativeReals)
+
+        # Generate datasets for this epsilon analysis
+        epsilon8_xresult_df = pd.DataFrame()
+        epsilon8_yresult_df = pd.DataFrame()
+
+        counter = 0
+        for e in np.arange(pyo.value(model.econ_loss_min), pyo.value(model.econ_loss_max) + 0.000001,
+                           pyo.value(model.econ_loss_step)):
+            model.obj_1_e = e
+            model.add_component("objective_1_constraint",
+                                Constraint(expr=sum_product(model.l_ijk, model.x_ijk) <= pyo.value(model.obj_1_e)))
+            for e2 in np.arange(pyo.value(model.functionality_min), pyo.value(model.functionality_max) +
+                                0.0000000000001, pyo.value(model.functionality_step)):
+                counter += 1
+                model.obj_3_e = e2
+                print("Step ", counter, " e: ", e, "  e2:", e2)
+                #             obj_2_13_epsilon_results.loc[counter-1,'Economic Loss Epsilon']=e
+                #             obj_2_13_epsilon_results.loc[counter-1,'Functionality Epsilon']=e2
+                model.add_component("objective_3_constraint",
+                                    Constraint(
+                                        expr=sum_product(model.Q_t_hat, model.x_ijk) >= pyo.value(model.obj_3_e)))
+                # Solve the model:
+                results = model_solver_setting.solve(model)
+
+                # Save the results if the solver returns an optimal solution:
+                if (results.solver.status == SolverStatus.ok) and (
+                        results.solver.termination_condition == TerminationCondition.optimal):
+                    self.extract_optimization_results(model)
+                    budget_used = quicksum(pyo.value(model.y_ijkk_prime[i, j, k, k_prime]) * pyo.value(
+                        model.Sc_ijkk_prime[i, j, k, k_prime]) for (i, j, k, k_prime) in model.ZSKK_prime)
+                    percent_budget_used = (budget_used / pyo.value(
+                        model.B)) * 100  # Record percentage of available budget used.
+                    # Add objective (dislocation), economic loss, and functionality values to results dataframe:
+                    # TODO: optimize this search
+                    obj_2_13_epsilon_results.loc[counter - 1, 'Functionality Value'] = pyo.value(model.functionality)
+                    obj_2_13_epsilon_results.loc[counter - 1, 'Economic Loss(Million Dollars)'] = pyo.value(
+                        model.econ_loss)
+                    obj_2_13_epsilon_results.loc[counter - 1, 'Dislocation Value'] = pyo.value(model.dislocation)
+
+                    # Extract results per each variable and convert to dataframe
+                    x_data: dict = model.x_ijk.extract_values()
+                    newx_df = self.assemble_dataframe_from_solution("x_ijk", x_data, counter)
+
+                    y_data: dict = model.y_ijkk_prime.extract_values()
+                    newy_df = self.assemble_dataframe_from_solution("y_ijkk_prime", y_data, counter)
+
+                    # Append to local analysis result
+                    epsilon8_xresult_df = epsilon8_xresult_df.append(newx_df)
+                    epsilon8_yresult_df = epsilon8_yresult_df.append(newy_df)
+                else:
+                    print(results.solver.termination_condition)
+                    log_infeasible_constraints(model)
+
+                # Remove the given epsilon constraint from the model now that optimization is complete:
+                model.del_component(model.objective_3_constraint)
+
+            # Remove the given epsilon constraint from the model now that optimization is complete:
+            model.del_component(model.objective_1_constraint)
+
+        # Remove the economic loss epsilon parameter from the model:
+        model.del_component(model.obj_1_e)
+        model.del_component(model.obj_3_e)
+
+        print(obj_2_13_epsilon_results)
+
+        # Drop rows with infeasible results:
+        initial_length = len(obj_2_13_epsilon_results)
+        obj_2_13_data = obj_2_13_epsilon_results.dropna(axis=0, how='any')
+        print("Infeasible rows dropped: ", initial_length - len(obj_2_13_data), " rows.")
+
+        # Save results:
+        filename = '8_epsilon_results.csv'
+        obj_2_13_data.index += 1
+        obj_2_13_data.to_csv(filename)
+
+        endtime = time.time()
+        elapsedtime = endtime - starttime
+        print("Elapsed time: ", elapsedtime)
+
+        epsilon8_xresult_df['Epsilon'] = 8
+        epsilon8_yresult_df['Epsilon'] = 8
+
+        return xresults_df.append(epsilon8_xresult_df), yresults_df.append(epsilon8_yresult_df)
+
+    def solve_epsilon_model_9(self, model, model_solver_setting, xresults_df, yresults_df):
+        starttime = time.time()
+        print(
+            "****OPTIMIZING BUILDING FUNCTIONALITY SUBJECT TO "
+            "ECONOMIC LOSS AND POPULATION DISLOCATION EPSILON CONSTRAINTS****")
+        model.objective_1.deactivate()
+        model.objective_2.deactivate()
+        model.objective_3.activate()
+
+        # Add objective functions 1 and 2 as epsilon constraints of objective function 3:
+        # Dataframe to store optimization results:
+        obj_3_12_epsilon_results = pd.DataFrame(
+            columns=['Economic Loss(Million Dollars)', 'Dislocation Value', 'Functionality Value'])
+
+        # Parameter for objective 2 (dislocation) and objective 1 (economic loss) epsilon values:
+        model.obj_1_e = Param(mutable=True, within=NonNegativeReals)
+        model.obj_2_e = Param(mutable=True, within=NonNegativeReals)
+
+        # Generate datasets for this epsilon analysis
+        epsilon9_xresult_df = pd.DataFrame()
+        epsilon9_yresult_df = pd.DataFrame()
+
+        counter = 0
+        for e in np.arange(pyo.value(model.econ_loss_min), pyo.value(model.econ_loss_max),
+                           pyo.value(model.econ_loss_step)):
+            model.obj_1_e = e
+            model.add_component("objective_1_constraint",
+                                Constraint(expr=sum_product(model.l_ijk, model.x_ijk) <= pyo.value(model.obj_1_e)))
+            for e2 in np.arange(pyo.value(model.dislocation_min), pyo.value(model.dislocation_max),
+                                pyo.value(model.dislocation_step)):
+                counter += 1
+                model.obj_2_e = e2
+                print("Step ", counter, " e: ", e, "  e2:", e2)
+                #             obj_3_12_epsilon_results.loc[counter-1,'Economic Loss Epsilon']=e
+                #             obj_3_12_epsilon_results.loc[counter-1,'Dislocation Epsilon']=e2
+                model.add_component("objective_2_constraint",
+                                    Constraint(
+                                        expr=sum_product(model.d_ijk, model.x_ijk) <= (pyo.value(model.obj_2_e))))
+                # Solve the model:
+                results = model_solver_setting.solve(model)
+
+                # Save the results if the solver returns an optimal solution:
+                if (results.solver.status == SolverStatus.ok) and (
+                        results.solver.termination_condition == TerminationCondition.optimal):
+                    self.extract_optimization_results(model)
+                    budget_used = quicksum(pyo.value(model.y_ijkk_prime[i, j, k, k_prime]) * pyo.value(
+                        model.Sc_ijkk_prime[i, j, k, k_prime]) for (i, j, k, k_prime) in model.ZSKK_prime)
+                    percent_budget_used = (budget_used / pyo.value(
+                        model.B)) * 100  # Record percentage of available budget used.
+                    # Add objective (functionality), economic loss, and dislocation values to results dataframe:
+                    # TODO: optimize this search
+                    obj_3_12_epsilon_results.loc[counter - 1, 'Functionality Value'] = pyo.value(model.functionality)
+                    obj_3_12_epsilon_results.loc[counter - 1, 'Economic Loss(Million Dollars)'] = pyo.value(
+                        model.econ_loss)
+                    obj_3_12_epsilon_results.loc[counter - 1, 'Dislocation Value'] = pyo.value(model.dislocation)
+
+                    # Extract results per each variable and convert to dataframe
+                    x_data: dict = model.x_ijk.extract_values()
+                    newx_df = self.assemble_dataframe_from_solution("x_ijk", x_data, counter)
+
+                    y_data: dict = model.y_ijkk_prime.extract_values()
+                    newy_df = self.assemble_dataframe_from_solution("y_ijkk_prime", y_data, counter)
+
+                    # Append to local analysis result
+                    epsilon9_xresult_df = epsilon9_xresult_df.append(newx_df)
+                    epsilon9_yresult_df = epsilon9_yresult_df.append(newy_df)
+                else:
+                    print(results.solver.termination_condition)
+                    log_infeasible_constraints(model)
+
+                # Remove the given epsilon constraint from the model now that optimization is complete:
+                model.del_component(model.objective_2_constraint)
+
+            # Remove the given epsilon constraint from the model now that optimization is complete:
+            model.del_component(model.objective_1_constraint)
+
+        # Remove the economic loss epsilon parameter from the model:
+        model.del_component(model.obj_1_e)
+        model.del_component(model.obj_2_e)
+
+        print(obj_3_12_epsilon_results)
+
+        # Drop rows with infeasible results:
+        initial_length = len(obj_3_12_epsilon_results)
+        obj_3_12_data = obj_3_12_epsilon_results.dropna(axis=0, how='any')
+        print("Infeasible rows dropped: ", initial_length - len(obj_3_12_data), " rows.")
+
+        # Save results:
+        filename = '9_epsilon_results.csv'
+        obj_3_12_data.index += 1
+        obj_3_12_data.to_csv(filename)
+
+        endtime = time.time()
+        elapsedtime = endtime - starttime
+        print("Elapsed time: ", elapsedtime)
+
+        epsilon9_xresult_df['Epsilon'] = 9
+        epsilon9_yresult_df['Epsilon'] = 9
+
+        return xresults_df.append(epsilon9_xresult_df), yresults_df.append(epsilon9_yresult_df)
+
+    def compute_optimal_results(self, inactive_submodels, xresults_df, yresults_df):
+        # Fixed for the moment, will be expanded to the full model set in later iterations
+        epsilon_models = [7, 8, 9]
+
+        xresults_list: list = []
+        yresults_list: list = []
+
+        for k in epsilon_models:
+            if k not in inactive_submodels:
+                results = pd.read_csv(str(k) + '_epsilon_results.csv', usecols=['Economic Loss(Million Dollars)',
+                                                                             'Dislocation Value',
+                                                                             'Functionality Value'],
+                                      low_memory=False)
+                list_loss = results['Economic Loss(Million Dollars)'].values.tolist()
+                list_dislocation = results['Dislocation Value'].values.tolist()
+                list_func = results['Functionality Value'].values.tolist()
+                zipped_list = self.optimal_points(list_loss, list_dislocation, list_func)
+                optimal = pd.DataFrame(zipped_list, columns=['Iteration', 'Economic Loss(Million Dollars)',
+                                                             'Dislocation Value', 'Functionality Value'])
+
+                # Select only results corresponding to current epsilon step
+                epsilon_xresults_df = xresults_df[xresults_df['Epsilon'] == k]
+                epsilon_yresults_df = yresults_df[yresults_df['Epsilon'] == k]
+
+                # Filter results depending on optimal epsilon model
+                opt_xresults_df = epsilon_xresults_df[epsilon_xresults_df['Iteration'].isin(optimal['Iteration'])]
+                opt_yresults_df = epsilon_yresults_df[epsilon_yresults_df['Iteration'].isin(optimal['Iteration'])]
+
+                # Append for later construction of final dataset
+                xresults_list.append(opt_xresults_df)
+                yresults_list.append(opt_yresults_df)
+
+        # Construct the final dataframe per variable
+        return [pd.concat(xresults_list), pd.concat(yresults_list)]
+
+    # Objective functions
+    @staticmethod
+    def obj_economic(model):
+        # return(sum_product(model.l_ijk,model.x_ijk))
+        return quicksum(model.l_ijk[i, j, k] * model.x_ijk[i, j, k] for (i, j, k) in model.ZSK)
+
+    @staticmethod
+    def obj_dislocation(model):
+        # return(sum_product(model.d_ijk,model.x_ijk))
+        return quicksum(model.d_ijk[i, j, k] * model.x_ijk[i, j, k] for (i, j, k) in model.ZSK)
+
+    @staticmethod
+    def obj_functionality(model):
+        # return(sum_product(model.Q_t_hat,model.x_ijk))
+        return quicksum(model.Q_t_hat[i, j, k] * model.x_ijk[i, j, k] for (i, j, k) in model.ZSK)
+
+    @staticmethod
+    def retrofit_cost_rule(model):
+        return (None,
+                quicksum(
+                    model.Sc_ijkk_prime[i, j, k, k_prime] * model.y_ijkk_prime[i, j, k, k_prime] for
+                    (i, j, k, k_prime)
+                    in model.ZSKK_prime),
+                pyo.value(model.B))
+
+    @staticmethod
+    def number_buildings_ij_rule(model, i, j):
+        return (quicksum(pyo.value(model.b_ijk[i, j, k]) for k in model.K),
+                quicksum(model.x_ijk[i, j, k] for k in model.K),
+                quicksum(pyo.value(model.b_ijk[i, j, k]) for k in model.K))
+
+    @staticmethod
+    def building_level_rule(model, i, j, k):
+        model.a = quicksum(model.y_ijkk_prime[i, j, k_prime, k] for k_prime in model.K_prime if
+                           (i, j, k_prime, k) in model.zskk_prime)
+        model.c = quicksum(model.y_ijkk_prime[i, j, k, k_prime] for k_prime in model.K_prime if
+                           (i, j, k, k_prime) in model.zskk_prime)
+        return (pyo.value(model.b_ijk[i, j, k]),
+                model.x_ijk[i, j, k] + quicksum(model.y_ijkk_prime[i, j, k, k_prime] for k_prime in model.K_prime if
+                                                (i, j, k, k_prime) in model.zskk_prime) -
+                quicksum(model.y_ijkk_prime[i, j, k_prime, k] for k_prime in model.K_prime if
+                         (i, j, k_prime, k) in model.zskk_prime),
+                pyo.value(model.b_ijk[i, j, k]))
+
+    @staticmethod
+    def extract_optimization_results(model):
+        model.econ_loss = quicksum(
+            pyo.value(model.l_ijk[i, j, k]) * pyo.value(model.x_ijk[i, j, k]) for (i, j, k) in model.ZSK)
+        model.dislocation = quicksum(
+            pyo.value(model.d_ijk[i, j, k]) * pyo.value(model.x_ijk[i, j, k]) for (i, j, k) in model.ZSK)
+        model.functionality = quicksum(
+            pyo.value(model.Q_t_hat[i, j, k]) * pyo.value(model.x_ijk[i, j, k]) for (i, j, k) in model.ZSK)
+
+    @staticmethod
+    def assemble_dataframe_from_solution(variable, sol_dict, iteration):
+        x_index_lists = list(map(list, zip(*sol_dict.keys())))
+        x_dict = {
+            "Z": x_index_lists[0],
+            "S": x_index_lists[1],
+            "K": x_index_lists[2],
+        }
+
+        if variable == "y_ijkk_prime":
+            x_dict["K'"] = x_index_lists[3]
+
+        x_dict[variable] = sol_dict.values()
+
+        df = pd.DataFrame(x_dict)
+        df['Iteration'] = iteration
+
+        return df
+
+    @staticmethod
+    def optimal_points(list_loss, list_dislocation, list_func):
+        l_com = []
+        for i in range(0, len(list_loss)):
+            l_com.append(tuple([list_loss[i], list_dislocation[i], list_func[i]]))
+        # remove the repeated the points
+        set_com = set(l_com)
+
+        # change the set to list
+        list_temp = []
+        for x in set_com:
+            list_temp.append(x)
+
+        final_f = []
+        f_temp = []
+        for i in range(len(list_temp)):
+            for j in range(i + 1, len(list_temp)):
+                if list_temp[i][0] == list_temp[j][0] and list_temp[i][1] == list_temp[j][1]:
+                    if list_temp[i][2] > list_temp[j][2]:
+                        f_temp.append(list_temp[j])
+
+                    elif list_temp[j][2] > list_temp[i][2]:
+                        f_temp.append(list_temp[i])
+        final_f = list(set(list_temp) - set(f_temp))
+
+        # find the optimal value for dislocation if loss and functionality are same as the last loop
+        final_d = []
+        d_temp = []
+        for i in range(len(final_f)):
+            for j in range(i + 1, len(final_f)):
+                if final_f[i][0] == final_f[j][0] and final_f[i][2] == final_f[j][2]:
+                    if final_f[i][1] < final_f[j][1]:
+                        d_temp.append(final_f[j])
+                    elif final_f[i][1] > final_f[j][1]:
+                        d_temp.append(final_f[i])
+
+        final_d = list(set(final_f) - set(d_temp))
+
+        # find the optimal value for loss if dislocation and functionality are same as the last loop
+        final_l = []
+        l_temp = []
+        for i in range(len(final_d)):
+            for j in range(i + 1, len(final_d)):
+                if final_d[i][1] == final_d[j][1] and final_d[i][2] == final_d[j][2]:
+                    if final_d[i][0] < final_d[j][0]:
+                        l_temp.append(final_d[j])
+                    elif final_d[i][0] > final_d[j][0]:
+                        l_temp.append(final_d[i])
+        final_l = list(set(final_d) - set(l_temp))
+
+        loss_optimal = []
+        dislocation_optimal = []
+        func_optimal = []
+        for x in final_l:
+            loss_optimal.append(x[0])
+            dislocation_optimal.append(x[1])
+            func_optimal.append(x[2])
+        comb1 = []
+        for i in range(len(list_loss)):
+            comb1.append([list_loss[i], list_dislocation[i], list_func[i]])
+
+        comb2 = []
+        for i in range(len(loss_optimal)):
+            comb2.append([loss_optimal[i], dislocation_optimal[i],
+                          func_optimal[i]])
+
+        optimal_index = {}
+        for i in range(len(comb2)):
+            if comb2[i] not in comb1:
+                pass
+            else:
+                optimal_index[comb1.index(comb2[i])] = comb2[i]
+        optimal_index_list = list(optimal_index.keys())
+        zipped_list = list(zip(optimal_index_list, loss_optimal, dislocation_optimal, func_optimal))
+        return zipped_list
+
+    def get_spec(self):
+        """Get specifications of the multiobjective retrofit optimization model.
+
+        Returns:
+            obj: A JSON object of specifications of the multiobjective retrofit optimization model.
+
+        """
+        return {
+            'name': 'multiobjective-retrofit-optimization',
+            'description': 'Multiobjective retrofit optimization model',
+            'input_parameters': [
+                {
+                    'id': 'result_name',
+                    'required': False,
+                    'description': 'Result CSV dataset name',
+                    'type': str
+                },
+                {
+                    'id': 'model_solver',
+                    'required': False,
+                    'description': 'Choice of the model solver to use. Gurobi is the default solver.',
+                    'type': str
+                },
+                {
+                    'id': 'num_epsilon_steps',
+                    'required': True,
+                    'description': 'Number of epsilon values to evaluate',
+                    'type': int
+                },
+                {
+                    'id': 'max_budget',
+                    'required': True,
+                    'description': 'Selection of maximum possible budget',
+                    'type': str
+                },
+                {
+                    'id': 'budget_available',
+                    'required': False,
+                    'description': 'Custom budget value',
+                    'type': float
+                },
+                {
+                    'id': 'inactive_submodels',
+                    'required': False,
+                    'description': 'Identifier of submodels to inactivate during analysis',
+                    'type': [int]
+                },
+                {
+                    'id': 'scale_data',
+                    'required': True,
+                    'description': 'Choice for scaling data',
+                    'type': bool
+                },
+                {
+                    'id': 'scaling_factor',
+                    'required': False,
+                    'description': 'Custom scaling factor',
+                    'type': float
+                },
+            ],
+            'input_datasets': [
+                {
+                    'id': 'building_related_data',
+                    'required': True,
+                    'description': 'A csv file with building related data required to evaluate retrofit strategies',
+                    'type': 'incore:multiobjectiveBuildingRelatedData'
+                },
+                {
+                    'id': 'strategy_costs_data',
+                    'required': True,
+                    'description': 'A csv file with strategy cost data'
+                                   'per building',
+                    'type': 'incore:multiobjectiveStrategyCosts'
+                },
+            ],
+            'output_datasets': [
+                {
+                    'id': 'optimal_solution_dv_x',
+                    'parent_type': '',
+                    'description': 'Optimal solution for decision variable x',
+                    'type': 'incore:multiobjectiveOptimalSolutionX'
+                },
+                {
+                    'id': 'optimal_solution_dv_y',
+                    'parent_type': '',
+                    'description': 'Optimal solution for decision variable y with initial and final retrofitted '
+                                   'strategies',
+                    'type': 'incore:multiobjectiveOptimalSolutionY'
+                }
+            ]
+        }
diff --git a/pyincore/globals.py b/pyincore/globals.py
index dfa90e3d4..1153e2aa3 100644
--- a/pyincore/globals.py
+++ b/pyincore/globals.py
@@ -10,7 +10,7 @@
 import os
 import shutil
 
-PACKAGE_VERSION = "1.1.0"
+PACKAGE_VERSION = "1.2.0"
 
 INCORE_API_PROD_URL = "https://incore.ncsa.illinois.edu"
 INCORE_API_DEV_URL = "https://incore-dev.ncsa.illinois.edu"
diff --git a/pyincore/models/fragilitycurve.py b/pyincore/models/fragilitycurve.py
index 3045d399d..c71988e4e 100644
--- a/pyincore/models/fragilitycurve.py
+++ b/pyincore/models/fragilitycurve.py
@@ -101,31 +101,33 @@ def calculate_limit_state_probability(self, hazard_values: dict, fragility_curve
 
         return probability
 
-    def get_building_period(self, num_stories):
-        """Get building period from the fragility curve.
+    def get_building_period(self, fragility_curve_parameters, **kwargs):
+        """
+                Get building period from the fragility curve.
+
         Args:
-            num_stories (int): Number of building stories.
+            fragility_curve_parameters (dict): Fragility curve parameters.
+            **kwargs: Keyword arguments.
+
         Returns:
             float: Building period.
+
         """
+
         period = 0.0
-        return period
+        num_stories = 1.0
+        for parameter in fragility_curve_parameters:
+            # if default exists, use default
+            if parameter["name"] == "num_stories" and "expression" in parameter and parameter["expression"] is not None:
+                num_stories = evaluateexpression.evaluate(parameter["expression"])
 
-    # def get_building_period(self, fragility_curve_parameters, **kwargs):
-    #     period = 0.0
-    #     num_stories = 1.0
-    #     for parameter in fragility_curve_parameters:
-    #         # if default exists, use default
-    #         if parameter["name"] == "num_stories" and "expression" in parameter and parameter["expression"] is not None:
-    #             num_stories = evaluateexpression.evaluate(parameter["expression"])
-    #
-    #         # if exist in building inventory
-    #         for kwargs_key, kwargs_value in kwargs.items():
-    #             if kwargs_key.lower() == "num_stories" and kwargs_value is not None and kwargs_value > 0:
-    #                 num_stories = kwargs_value
-    #
-    #         # calculate period
-    #         if parameter["name"] == "period" and "expression" in parameter and parameter["expression"] is not None:
-    #             period = evaluateexpression.evaluate(parameter["expression"], {"num_stories": num_stories})
-    #
-    #     return period
+            # if exist in building inventory
+            for kwargs_key, kwargs_value in kwargs.items():
+                if kwargs_key.lower() == "num_stories" and kwargs_value is not None and kwargs_value > 0:
+                    num_stories = kwargs_value
+
+            # calculate period
+            if parameter["name"] == "period" and "expression" in parameter and parameter["expression"] is not None:
+                period = evaluateexpression.evaluate(parameter["expression"], {"num_stories": num_stories})
+
+        return period
diff --git a/pyincore/utils/analysisutil.py b/pyincore/utils/analysisutil.py
index d6c60d2bd..53154a061 100644
--- a/pyincore/utils/analysisutil.py
+++ b/pyincore/utils/analysisutil.py
@@ -11,7 +11,6 @@
 from typing import List, Dict
 from collections import Counter
 
-from deprecated.sphinx import deprecated
 from decimal import getcontext, Decimal
 
 from pyincore import DataService
@@ -69,27 +68,6 @@ def float_list_to_decimal(num_list: list):
     def float_dict_to_decimal(num_dict: dict):
         return {key: Decimal(str(num_dict[key])) for key in num_dict}
 
-    @staticmethod
-    @deprecated(version="0.9.0", reason="Use calculate_damage_interval in fragilitycurveset class instead.")
-    def calculate_damage_interval(damage):
-        output = collections.OrderedDict()
-
-        # convert damage to decimal dictionary
-        damage = AnalysisUtil.float_dict_to_decimal(damage)
-
-        if len(damage) == 4:
-            output['ds-none'] = 1 - damage['ls-slight']
-            output['ds-slight'] = damage['ls-slight'] - damage['ls-moderat']
-            output['ds-moderat'] = damage['ls-moderat'] - damage['ls-extensi']
-            output['ds-extensi'] = damage['ls-extensi'] - damage['ls-complet']
-            output['ds-complet'] = damage['ls-complet']
-        elif len(damage) == 3:
-            output['insignific'] = 1 - damage['immocc']
-            output['moderate'] = damage['immocc'] - damage['lifesfty']
-            output['heavy'] = damage['lifesfty'] - damage['collprev']
-            output['complete'] = damage['collprev']
-
-        return output
 
     @staticmethod
     def calculate_mean_damage(dmg_ratio_tbl, dmg_intervals,
@@ -251,12 +229,10 @@ def adjust_limit_states_for_pgd(limit_states, pgd_limit_states):
             adj_limit_states = collections.OrderedDict()
 
             for key, value in limit_states.items():
-                adj_limit_states[key] = limit_states[key] + pgd_limit_states[
-                    key] - \
-                                        (limit_states[key] * pgd_limit_states[
-                                            key])
+                adj_limit_states[key] = limit_states[key] + pgd_limit_states[key] - \
+                                        (limit_states[key] * pgd_limit_states[key])
 
-            return adj_limit_states
+            return AnalysisUtil.update_precision_of_dicts(adj_limit_states)
 
         except KeyError as e:
             print('Mismatched keys encountered in the limit states')
@@ -401,7 +377,10 @@ def get_hazard_demand_type(building, fragility_set, hazard_type):
         if hazard_type.lower() == "earthquake":
             num_stories = building[PROPERTIES][BLDG_STORIES]
             # Get building period from the fragility if possible
-            building_period = fragility_set.fragility_curves[0].get_building_period(num_stories)
+
+            building_args = fragility_set.construct_expression_args_from_inventory(building)
+            building_period = fragility_set.fragility_curves[0].get_building_period(
+                fragility_set.fragility_curve_parameters, **building_args)
 
             if fragility_hazard_type.endswith('sa') and fragility_hazard_type != 'sa':
                 # This fixes a bug where demand type is in a format similar to 1.0 Sec Sa
@@ -463,7 +442,9 @@ def get_hazard_demand_types(building, fragility_set, hazard_type):
                 demand_type = demand_type.lower()
                 num_stories = building[PROPERTIES][BLDG_STORIES]
                 # Get building period from the fragility if possible
-                building_period = fragility_set.fragility_curves[0].get_building_period(num_stories)
+                building_args = fragility_set.construct_expression_args_from_inventory(building)
+                building_period = fragility_set.fragility_curves[0].get_building_period(
+                    fragility_set.fragility_curve_parameters, **building_args)
 
                 # TODO: There might be a bug here as this is not handling SD
                 if demand_type.endswith('sa'):
diff --git a/pyincore/utils/popdisloutputprocess.py b/pyincore/utils/popdisloutputprocess.py
index 3d32d92e2..aeb9ea622 100644
--- a/pyincore/utils/popdisloutputprocess.py
+++ b/pyincore/utils/popdisloutputprocess.py
@@ -17,35 +17,62 @@ class PopDislOutputProcess:
         pop_disl_result (obj): IN-CORE dataset for Joplin Population Dislocation (PD) results.
         pop_disl_result_path (obj): A fallback for the case that Joplin PD object is not provided.
             For example a user wants to directly pass in csv files, a path to PD results.
-        filter_on (bool): A flag to filter all data, default True counts only Joplin buildings
+        filter_name (str): A string to filter data by name, default empty. Example: filter_name="Joplin" for Joplin
+            inventory, other is Duquesne etc. Name must be valid.
+        filter_guid (bool): A flag to filter all data, default True counts only Joplin buildings.
+        vacant_disl (bool): A flag to include vacant (Vacant for tenure) dislocation
 
     """
     HUPD_CATEGORIES = ["household_characteristics",
                        "household_dislocated",
                        "total_households",
-                       "%_household_dislocated",
+                       "percent_household_dislocated",
                        "population_dislocated",
                        "total_population",
-                       "%_population_dislocated"
+                       "percent_population_dislocated"
                        ]
 
-    def __init__(self, pop_disl_result, pop_disl_result_path=None, filter_on=True):
+    def __init__(self, pop_disl_result, pop_disl_result_path=None,
+                 filter_name=None, filter_guid=True, vacant_disl=True):
         if pop_disl_result_path:
             pd_result = pd.read_csv(pop_disl_result_path, low_memory=False)
         else:
             pd_result = pop_disl_result.get_dataframe_from_csv(low_memory=False)
         pd_result["geometry"] = pd_result["geometry"].apply(wkt.loads)
 
-        pd_result_shp = None
         # keep only inventory with guid; filter for Joplin since only Joplin inventory has guids
-        if filter_on:
-            pd_result = pd_result[(pd_result["guid"].notnull()) &
-                                  (pd_result["numprec"].notnull())]
-            # only keep guid and dislocated
-            pd_result_shp = pd_result[(pd_result["dislocated"]) &
-                                      (pd_result["guid"].notnull()) &
-                                      (pd_result["numprec"].notnull())]
-        self.pop_disl_result = pd_result
+        if filter_guid:
+            if filter_name:
+                pd_result_flag = pd_result[(pd_result["guid"].notnull()) &
+                                           (pd_result["numprec"].notnull()) &
+                                           (pd_result["plcname10"] == filter_name)]
+                # only keep guid, place and dislocated
+                pd_result_shp = pd_result[(pd_result["dislocated"]) &
+                                          (pd_result["guid"].notnull()) &
+                                          (pd_result["numprec"].notnull()) &
+                                          (pd_result["plcname10"] == filter_name)]
+            else:
+                pd_result_flag = pd_result[(pd_result["guid"].notnull()) &
+                                           (pd_result["numprec"].notnull())]
+                # only keep guid and dislocated
+                pd_result_shp = pd_result[(pd_result["dislocated"]) &
+                                          (pd_result["guid"].notnull()) &
+                                          (pd_result["numprec"].notnull())]
+        else:
+            if filter_name:
+                pd_result_flag = pd_result[(pd_result["numprec"].notnull()) &
+                                           (pd_result["plcname10"] == filter_name)]
+                # only keep guid, place and dislocated
+                pd_result_shp = pd_result[(pd_result["dislocated"]) &
+                                          (pd_result["numprec"].notnull()) &
+                                          (pd_result["plcname10"] == filter_name)]
+            else:
+                pd_result_flag = pd_result[(pd_result["numprec"].notnull())]
+                # only keep guid and dislocated
+                pd_result_shp = pd_result[(pd_result["dislocated"]) &
+                                          (pd_result["numprec"].notnull())]
+        self.vacant_disl = vacant_disl
+        self.pop_disl_result = pd_result_flag
         self.pop_disl_result_shp = pd_result_shp
 
     def get_heatmap_shp(self, filename="pop-disl-numprec.shp"):
@@ -93,7 +120,8 @@ def pd_by_race(self, filename_json=None):
         # The numbering follows the Community description notebook
         # 0 - Vacant HU No Race Ethnicity Data, 1 - Not Hispanic/White, 2 - Not Hispanic/Black
         # 3 - Not Hispanic/Other race, 4 - Hispanic, 5 - No Race or Ethnicity Data
-        race_categories = ["Not Hispanic/White",
+        race_categories = ["Vacant HU No Race or Ethnicity Data",
+                           "Not Hispanic/White",
                            "Not Hispanic/Black",
                            "Not Hispanic/Other Race",
                            "Hispanic",
@@ -109,54 +137,57 @@ def pd_by_race(self, filename_json=None):
         huapd.loc[(huapd["hispan"] == 1), "hua_re"] = "4"
         huapd.loc[(huapd["gqtype"] >= 1), "hua_re"] = "5"
         hua_vals = huapd["hua_re"].value_counts()
+
         hua_tot = []
-        for i in range(len(race_categories)):
+        for i in range(len(race_categories) - 1):
             try:
                 hua_tot.append(int(hua_vals[str(i)]))
             except Exception:
                 hua_tot.append(0)
-        hua_tot.append(int(sum(hua_tot[1:])))
+        hua_tot.append(int(sum(hua_tot)))
 
         pop_tot = []
-        for i in range(len(race_categories)):
+        for i in range(len(race_categories) - 1):
             pop_tot.append(int(huapd["numprec"].where(huapd["hua_re"] == str(i)).sum()))
-        pop_tot.append(int(sum(pop_tot[1:])))
+        pop_tot.append(int(sum(pop_tot)))
 
         # Dislocated by race and ethnicity
-        huapd["hud_re"] = "0"
+        huapd["hud_re"] = ""
+        huapd.loc[huapd["dislocated"], "hud_re"] = "0"
         huapd.loc[(huapd["race"] == 1) & (huapd["hispan"] == 0) & huapd["dislocated"], "hud_re"] = "1"
         huapd.loc[(huapd["race"] == 2) & (huapd["hispan"] == 0) & huapd["dislocated"], "hud_re"] = "2"
         huapd.loc[(huapd["race"].isin([3, 4, 5, 6, 7])) & (huapd["hispan"] == 0) & huapd["dislocated"], "hud_re"] = "3"
         huapd.loc[(huapd["hispan"] == 1) & huapd["dislocated"], "hud_re"] = "4"
         huapd.loc[(huapd["gqtype"] >= 1) & huapd["dislocated"], "hud_re"] = "5"
         hud_vals = huapd["hud_re"].value_counts()
+
         hua_disl = []
-        for i in range(len(race_categories)):
+        for i in range(len(race_categories) - 1):
             try:
                 hua_disl.append(int(hud_vals[str(i)]))
             except Exception:
                 hua_disl.append(0)
-        hua_disl.append(int(sum(hua_disl[1:])))
+        hua_disl.append(int(sum(hua_disl)))
 
         pd_disl = []
-        for i in range(len(race_categories)):
+        for i in range(len(race_categories) - 1):
             pd_disl.append(int(huapd["numprec"].where(huapd["hud_re"] == str(i)).sum()))
-        pd_disl.append(int(sum(pd_disl[1:])))
+        pd_disl.append(int(sum(pd_disl)))
 
         pd_by_race_json = []
         for i in range(len(race_categories)):
             huapd_race = {}
             huapd_race[self.HUPD_CATEGORIES[0]] = race_categories[i]
-            huapd_race[self.HUPD_CATEGORIES[1]] = hua_disl[i + 1]
-            huapd_race[self.HUPD_CATEGORIES[2]] = hua_tot[i + 1]
-            if hua_tot[i + 1]:
-                huapd_race[self.HUPD_CATEGORIES[3]] = 100 * (hua_disl[i + 1] / hua_tot[i + 1])
+            huapd_race[self.HUPD_CATEGORIES[1]] = hua_disl[i]
+            huapd_race[self.HUPD_CATEGORIES[2]] = hua_tot[i]
+            if hua_tot[i]:
+                huapd_race[self.HUPD_CATEGORIES[3]] = 100 * (hua_disl[i] / hua_tot[i])
             else:
                 huapd_race[self.HUPD_CATEGORIES[3]] = None
-            huapd_race[self.HUPD_CATEGORIES[4]] = pd_disl[i + 1]
-            huapd_race[self.HUPD_CATEGORIES[5]] = pop_tot[i + 1]
-            if pop_tot[i + 1]:
-                huapd_race[self.HUPD_CATEGORIES[6]] = 100 * (pd_disl[i + 1] / pop_tot[i + 1])
+            huapd_race[self.HUPD_CATEGORIES[4]] = pd_disl[i]
+            huapd_race[self.HUPD_CATEGORIES[5]] = pop_tot[i]
+            if pop_tot[i]:
+                huapd_race[self.HUPD_CATEGORIES[6]] = 100 * (pd_disl[i] / pop_tot[i])
             else:
                 huapd_race[self.HUPD_CATEGORIES[6]] = None
             pd_by_race_json.append(huapd_race)
@@ -319,7 +350,8 @@ def pd_by_tenure(self, filename_json=None):
         pop_tot.append(int(sum(pop_tot[1:])))
 
         # Dislocated by tenure
-        huapd["hud_tnr"] = "0"
+        huapd["hud_tnr"] = ""
+        huapd.loc[huapd["dislocated"], "hud_tnr"] = "0"
         huapd.loc[(huapd["ownershp"] == 1.0) & huapd["dislocated"], "hud_tnr"] = "1"
         huapd.loc[(huapd["ownershp"] == 2.0) & huapd["dislocated"], "hud_tnr"] = "2"
         huapd.loc[(huapd["gqtype"] == 3) & huapd["dislocated"], "hud_tnr"] = "3"
@@ -334,6 +366,11 @@ def pd_by_tenure(self, filename_json=None):
                 hua_disl.append(int(hud_vals[str(i)]))
             except Exception:
                 hua_disl.append(0)
+        # If vacant_disl is False the Vacant places do not dislocate (set to 0).
+        for i in range(len(tenure_categories)):
+            if not self.vacant_disl and "Vacant" in tenure_categories[i]:
+                hua_disl[i + 1] = 0
+
         hua_disl.append(int(sum(hua_disl[1:])))
 
         pd_disl = []
@@ -414,7 +451,8 @@ def pd_by_housing(self, filename_json=None):
         pop_tot.append(int(sum(pop_tot[1:])))
 
         # Dislocated by household
-        huapd["hud_house"] = "0"
+        huapd["hud_house"] = ""
+        huapd.loc[huapd["dislocated"], "hud_house"] = "0"
         huapd.loc[(huapd["huestimate"] == 1.0) & huapd["dislocated"], "hud_house"] = "1"
         huapd.loc[(huapd["huestimate"] > 1.0) & huapd["dislocated"], "hud_house"] = "2"
         hud_vals = huapd["hud_house"].value_counts()
@@ -478,40 +516,35 @@ def pd_total(self, filename_json=None):
         hud_vals = hud["dislocated"].value_counts()
         hua_disl = [int(hud_vals[False]), int(hud_vals[True])]
 
-        pd_disl = []
-        pd_disl.append(int(hud["numprec"].where(hud["dislocated"] == 0).sum()))
-        pd_disl.append(int(hud["numprec"].where(hud["dislocated"] == 1).sum()))
+        pd_disl = [int(hud["numprec"].where(hud["dislocated"] == 0).sum()),
+                   int(hud["numprec"].where(hud["dislocated"] == 1).sum())]
 
         hua_tot = sum(hua_disl)
         pop_tot = sum(pd_disl)
 
+        no_hua_tot = {"households": None, "percent_of_households": None}
         hua_disl_tot = {}
+        hua_disl_tot["dislocated"] = no_hua_tot
+        hua_disl_tot["not_dislocated"] = no_hua_tot
+        hua_disl_tot["total"] = no_hua_tot
         if hua_tot:
             hua_disl_tot["dislocated"] = {"households": hua_disl[1],
-                                          "%_of_households": 100 * (hua_disl[1]/hua_tot)}
+                                          "percent_of_households": 100 * (hua_disl[1]/hua_tot)}
             hua_disl_tot["not_dislocated"] = {"households": hua_tot - hua_disl[1],
-                                              "%_of_households": 100 * ((hua_tot - hua_disl[1])/hua_tot)}
-            hua_disl_tot["total"] = {"households": hua_tot, "%_of_households": 100}
-        else:
-            hua_disl_tot["dislocated"] = {"households": None,
-                                          "%_of_households": None}
-            hua_disl_tot["not_dislocated"] = {"households": None,
-                                              "%_of_households": None}
-            hua_disl_tot["total"] = {"households": None, "%_of_households": None}
+                                              "percent_of_households": 100 * ((hua_tot - hua_disl[1])/hua_tot)}
+            hua_disl_tot["total"] = {"households": hua_tot, "percent_of_households": 100}
 
+        no_pop_tot = {"population": None, "percent_of_population": None}
         pop_disl_tot = {}
+        hua_disl_tot["dislocated"] = no_pop_tot
+        hua_disl_tot["not_dislocated"] = no_pop_tot
+        hua_disl_tot["total"] = no_pop_tot
         if pop_tot:
             pop_disl_tot["dislocated"] = {"population": pd_disl[1],
-                                          "%_of_population": 100 * (pd_disl[1]/pop_tot)}
+                                          "percent_of_population": 100 * (pd_disl[1]/pop_tot)}
             pop_disl_tot["not_dislocated"] = {"population": pop_tot - pd_disl[1],
-                                              "%_of_population": 100 * ((pop_tot - pd_disl[1])/pop_tot)}
-            pop_disl_tot["total"] = {"population": pop_tot, "%_of_population": 100}
-        else:
-            pop_disl_tot["dislocated"] = {"population": None,
-                                          "%_of_population": None}
-            pop_disl_tot["not_dislocated"] = {"population": None,
-                                              "%_of_population": None}
-            pop_disl_tot["total"] = {"population": None, "%_of_population": None}
+                                              "percent_of_population": 100 * ((pop_tot - pd_disl[1])/pop_tot)}
+            pop_disl_tot["total"] = {"population": pop_tot, "percent_of_population": 100}
 
         pd_total_json = {"household_dislocation_in_total": hua_disl_tot,
                          "population_dislocation_in_total": pop_disl_tot}
diff --git a/recipes/meta.yaml b/recipes/meta.yaml
index ac07cf5ef..0fedce7c0 100644
--- a/recipes/meta.yaml
+++ b/recipes/meta.yaml
@@ -1,12 +1,12 @@
 {% set name = "pyincore" %}
-{% set version = "1.0.0" %}
+{% set version = "1.2.0" %}
  
 package:
   name: {{ name|lower }}
   version: {{ version }}
  
 source:
-  path: /path/to/pyincore
+  path: /path/to/pyincore/recipes
   # PR: https://opensource.ncsa.illinois.edu/bitbucket/projects/INCORE1/repos/pyincore/pull-requests/12/overview
  
 build:
diff --git a/scripts/release-packages.yml b/scripts/release-packages.yml
index 519ca72cc..13631a4eb 100644
--- a/scripts/release-packages.yml
+++ b/scripts/release-packages.yml
@@ -1,7 +1,7 @@
 # List of pyincore analyses to either include or exclude for release
 # If provided, 'excludes' takes preference over 'includeOnly'
 
-version: '1.1.0'
+version: '1.2.0'
 
 includeOnly:
   - buildingdamage
diff --git a/setup.py b/setup.py
index 408405c0e..84407fdf3 100644
--- a/setup.py
+++ b/setup.py
@@ -8,7 +8,7 @@
 
 setup(
     name='pyincore',
-    version='1.1.0',
+    version='1.2.0',
     packages=find_packages(where=".", exclude=["*.tests", "*.tests.*", "tests.*", "tests"]),
     include_package_data=True,
     package_data={


Test patch:
diff --git a/tests/pyincore/analyses/epfdamage/test_epfdamage.py b/tests/pyincore/analyses/epfdamage/test_epfdamage.py
index dd9e286df..02d27fcf6 100644
--- a/tests/pyincore/analyses/epfdamage/test_epfdamage.py
+++ b/tests/pyincore/analyses/epfdamage/test_epfdamage.py
@@ -6,83 +6,62 @@
 def run_with_base_class():
     client = IncoreClient(pyglobals.INCORE_API_DEV_URL)
 
+    # ## Memphis EPF Damage with Earthquake ##
     hazard_type_eq = "earthquake"
+    hazard_id_eq = "5b902cb273c3371e1236b36b"
+    epf_dataset_id = "6189c103d5b02930aa3efc35"
+    # mapping_id = "61980b11e32da63f4b9d86f4"
+    mapping_id = "61980cf5e32da63f4b9d86f5"  # PGA and PGD
 
-    hazard_id_eq = "5eebcbb08f80fe3899ad6039"
+    use_liquefaction = True
+    use_hazard_uncertainty = False
+    liquefaction_geology_dataset_id = "5a284f53c7d30d13bc08249c"
 
-    epf_dataset_id = "5eebcaa17a00803abc85ec11"
+    # Run epf damage
+    epf_dmg_eq_memphis = EpfDamage(client)
+    epf_dmg_eq_memphis.load_remote_input_dataset("epfs", epf_dataset_id)
 
+    # Load fragility mapping
+    fragility_service = FragilityService(client)
+    mapping_set = MappingSet(fragility_service.get_mapping(mapping_id))
+    epf_dmg_eq_memphis.set_input_dataset('dfr3_mapping_set', mapping_set)
+
+    epf_dmg_eq_memphis.set_parameter("result_name", "memphis_eq_epf_dmg_result")
+    epf_dmg_eq_memphis.set_parameter("hazard_type", hazard_type_eq)
+    epf_dmg_eq_memphis.set_parameter("hazard_id", hazard_id_eq)
+    epf_dmg_eq_memphis.set_parameter("use_liquefaction", use_liquefaction)
+    epf_dmg_eq_memphis.set_parameter("use_hazard_uncertainty", use_hazard_uncertainty)
+    epf_dmg_eq_memphis.set_parameter("liquefaction_geology_dataset_id", liquefaction_geology_dataset_id)
+    epf_dmg_eq_memphis.set_parameter("num_cpu", 1)
+    # Run Analysis
+    epf_dmg_eq_memphis.run_analysis()
+
+    # ##  Seaside EPF Damage with Earthquake ##
+    hazard_type_eq = "earthquake"
+    hazard_id_eq = "5eebcbb08f80fe3899ad6039"
+    epf_dataset_id = "5eebcaa17a00803abc85ec11"
     # Earthquake mapping
-    mapping_id = "5eebcc13e7226233ce4ef0d7"
+    mapping_id = "5eebcc13e7226233ce4ef0d7"  # only PGA
 
     # Run epf damage
     epf_dmg_eq = EpfDamage(client)
-
     epf_dmg_eq.load_remote_input_dataset("epfs", epf_dataset_id)
 
     # Load fragility mapping
     fragility_service = FragilityService(client)
     mapping_set = MappingSet(fragility_service.get_mapping(mapping_id))
     epf_dmg_eq.set_input_dataset('dfr3_mapping_set', mapping_set)
-
-    epf_dmg_eq.set_parameter("result_name", "earthquake_epf_dmg_result")
+    
+    epf_dmg_eq.set_parameter("result_name", "seaside_eq_epf_dmg_result")
     epf_dmg_eq.set_parameter("hazard_type", hazard_type_eq)
     epf_dmg_eq.set_parameter("hazard_id", hazard_id_eq)
     epf_dmg_eq.set_parameter("num_cpu", 1)
-
     # Run Analysis
     epf_dmg_eq.run_analysis()
 
+    # ## Seaside Tsunami EPF Damage ##
     hazard_type_tsu = "tsunami"
     hazard_id_tsu = "5bc9eaf7f7b08533c7e610e1"
-
-    epf_dataset_id = "5eebcaa17a00803abc85ec11"
-
-    # Tsunami mapping
-    mapping_id = "5eebce11e7226233ce4ef305"
-
-    # Run epf damage
-    epf_dmg_tsu = EpfDamage(client)
-    epf_dmg_tsu.load_remote_input_dataset("epfs", epf_dataset_id)
-
-    # Load fragility mapping
-    fragility_service = FragilityService(client)
-    mapping_set = MappingSet(fragility_service.get_mapping(mapping_id))
-    epf_dmg_tsu.set_input_dataset('dfr3_mapping_set', mapping_set)
-
-    epf_dmg_tsu.set_parameter("fragility_key", "Non-Retrofit inundationDepth Fragility ID Code")
-    epf_dmg_tsu.set_parameter("result_name", "tsunami_epf_dmg_result")
-    epf_dmg_tsu.set_parameter("hazard_type", hazard_type_tsu)
-    epf_dmg_tsu.set_parameter("hazard_id", hazard_id_tsu)
-    epf_dmg_tsu.set_parameter("num_cpu", 1)
-
-    # Run Analysis
-    epf_dmg_tsu.run_analysis()
-
-    ##################################
-    # test refactored
-
-    # Earthquake mapping
-    mapping_id = "605bf99b618178207f64754b"
-
-    # Run epf damage
-    epf_dmg_eq_re = EpfDamage(client)
-
-    epf_dmg_eq_re.load_remote_input_dataset("epfs", epf_dataset_id)
-
-    # Load fragility mapping
-    fragility_service = FragilityService(client)
-    mapping_set = MappingSet(fragility_service.get_mapping(mapping_id))
-    epf_dmg_eq_re.set_input_dataset('dfr3_mapping_set', mapping_set)
-
-    epf_dmg_eq_re.set_parameter("result_name", "refactored_earthquake_epf_dmg_result")
-    epf_dmg_eq_re.set_parameter("hazard_type", hazard_type_eq)
-    epf_dmg_eq_re.set_parameter("hazard_id", hazard_id_eq)
-    epf_dmg_eq_re.set_parameter("num_cpu", 1)
-
-    # Run Analysis
-    epf_dmg_eq_re.run_analysis()
-
     epf_dataset_id = "5eebcaa17a00803abc85ec11"
     # Tsunami mapping
     mapping_id = "605bf9e1618178207f6475c6"
@@ -90,6 +69,8 @@ def run_with_base_class():
     # Run epf damage
     epf_dmg_tsu_re = EpfDamage(client)
     epf_dmg_tsu_re.load_remote_input_dataset("epfs", epf_dataset_id)
+    epf_dmg_eq = EpfDamage(client)
+    epf_dmg_eq.load_remote_input_dataset("epfs", epf_dataset_id)
 
     # Load fragility mapping
     fragility_service = FragilityService(client)
@@ -97,11 +78,12 @@ def run_with_base_class():
     epf_dmg_tsu_re.set_input_dataset('dfr3_mapping_set', mapping_set)
 
     epf_dmg_tsu_re.set_parameter("fragility_key", "Non-Retrofit inundationDepth Fragility ID Code")
-    epf_dmg_tsu_re.set_parameter("result_name", "refactored_tsunami_epf_dmg_result")
+    epf_dmg_tsu_re.set_parameter("result_name", "seaside_tsunami_epf_dmg_result")
     epf_dmg_tsu_re.set_parameter("hazard_type", hazard_type_tsu)
     epf_dmg_tsu_re.set_parameter("hazard_id", hazard_id_tsu)
     epf_dmg_tsu_re.set_parameter("num_cpu", 1)
 
+    epf_dmg_tsu_re.set_input_dataset('dfr3_mapping_set', mapping_set)
     # Run Analysis
     epf_dmg_tsu_re.run_analysis()
 
diff --git a/tests/pyincore/analyses/joplinempiricalrestoration/test_joplinempiricalrestoration.py b/tests/pyincore/analyses/joplinempiricalrestoration/test_joplinempiricalrestoration.py
index c85d8db50..af13a182c 100644
--- a/tests/pyincore/analyses/joplinempiricalrestoration/test_joplinempiricalrestoration.py
+++ b/tests/pyincore/analyses/joplinempiricalrestoration/test_joplinempiricalrestoration.py
@@ -40,10 +40,13 @@ def run_with_base_class():
     building_dmg_result.get_dataframe_from_csv()
 
     restoration = JoplinEmpiricalRestoration(client)
+    # Building final target level dataset defining custom final restoration values
+    building_fl_id = "61a68979d5b02930aa43ae39"
 
     restoration.load_remote_input_dataset("buildings", bldg_dataset_id)
     # restoration.load_remote_input_dataset("building_dmg", building_dmg_result)
     restoration.set_input_dataset("building_dmg", building_dmg_result)
+    restoration.load_remote_input_dataset("building_functionality_level", building_fl_id)
 
     result_name = "Joplin_empirical_restoration_result"
     restoration.set_parameter("result_name", result_name)
diff --git a/tests/pyincore/analyses/multiobjectiveretrofitoptimization/test_multiobjectiveretrofitoptimization.py b/tests/pyincore/analyses/multiobjectiveretrofitoptimization/test_multiobjectiveretrofitoptimization.py
new file mode 100644
index 000000000..19cd42c30
--- /dev/null
+++ b/tests/pyincore/analyses/multiobjectiveretrofitoptimization/test_multiobjectiveretrofitoptimization.py
@@ -0,0 +1,25 @@
+from pyincore import IncoreClient
+from pyincore.analyses.multiobjectiveretrofitoptimization import MultiObjectiveRetrofitOptimization
+import pyincore.globals as pyglobals
+
+
+def run_base_analysis():
+    client = IncoreClient(pyglobals.INCORE_API_PROD_URL)
+    retrofit_optimization = MultiObjectiveRetrofitOptimization(client)
+
+    building_related_data = "6193ef5b6bee8c1fac5c915e"
+    strategy_costs_data = "6193efa69340a2170d51f495"
+
+    retrofit_optimization.set_parameter("model_solver", "ipopt")
+    retrofit_optimization.set_parameter("num_epsilon_steps", 10)
+    retrofit_optimization.set_parameter("max_budget", "default")
+    retrofit_optimization.set_parameter("scale_data", False)
+
+    retrofit_optimization.load_remote_input_dataset("building_related_data", building_related_data)
+    retrofit_optimization.load_remote_input_dataset("strategy_costs_data", strategy_costs_data)
+
+    retrofit_optimization.run_analysis()
+
+
+if __name__ == '__main__':
+    run_base_analysis()
diff --git a/tests/pyincore/utils/test_csvoutputjson.py b/tests/pyincore/utils/test_csvoutputjson.py
new file mode 100644
index 000000000..790c3bedb
--- /dev/null
+++ b/tests/pyincore/utils/test_csvoutputjson.py
@@ -0,0 +1,33 @@
+# This program and the accompanying materials are made available under the
+# terms of the Mozilla Public License v2.0 which accompanies this distribution,
+# and is available at https://www.mozilla.org/en-US/MPL/2.0/
+
+from pyincore.utils.cgeoutputprocess import CGEOutputProcess
+import os
+
+PYINCOREPATH = "path-to-pyincore"
+TESTAPATH = "pyincore/tests/pyincore/analyses/"
+
+def run_convert_cge_json():
+    # run the JoplinCGE analysis first to get results, csv files
+    cge_json = CGEOutputProcess()
+    filepath = os.path.join(PYINCOREPATH, TESTAPATH, "joplincge")
+
+    cge_json.get_cge_household_count(None,
+                                     os.path.join(filepath, "joplin-pop-disl-results.csv"),
+                                     "cge_total_household_count.json")
+    cge_json.get_cge_gross_income(None,
+                                  os.path.join(filepath, "gross-income.csv"),
+                                  "cge_total_household_income.json")
+    cge_json.get_cge_employment(None, None,
+                                os.path.join(filepath, "pre-disaster-factor-demand.csv"),
+                                os.path.join(filepath, "post-disaster-factor-demand.csv"),
+                                "cge_employment.json")
+    cge_json.get_cge_domestic_supply(None,
+                                     os.path.join(filepath, "domestic-supply.csv"),
+                                     "cge_domestic_supply.json")
+    return True
+
+
+if __name__ == '__main__':
+    run_convert_cge_json()
diff --git a/tests/pyincore/utils/test_popdisloutputprocess.py b/tests/pyincore/utils/test_popdisloutputprocess.py
index 16784c30b..0e5367f56 100644
--- a/tests/pyincore/utils/test_popdisloutputprocess.py
+++ b/tests/pyincore/utils/test_popdisloutputprocess.py
@@ -40,7 +40,7 @@ def run_convert_pd_json_chained(client):
     # incore-dev
     building_dmg = "602d96e4b1db9c28aeeebdce"  # dev Joplin
     # building_dmg = "602d975db1db9c28aeeebe35" # 15 guids test - dev Joplin
-    housing_unit_alloc = "61563545483ecb19e4304c2a"  # dev Joplin 2ev3
+    housing_unit_alloc = "6177089781ce0838fbacd6a8"  # dev Joplin 2ev3
     bg_data = "5df7cb0b425e0b00092c9464"  # Joplin 2ev2
     value_loss = "602d508fb1db9c28aeedb2a5"
 
@@ -70,6 +70,10 @@ def run_convert_pd_json_chained(client):
     pd_result = run_convert_pd_json_chained(client)
 
     pd_process = PopDislOutputProcess(pd_result)
+    # pd_process = PopDislOutputProcess(pd_result,
+    #                                   filter_name="Joplin",
+    #                                   filter_guid=False,
+    #                                   vacant_disl=False)
 
     pd_process.pd_by_race("PD_by_race.json")
     pd_process.pd_by_income("PD_by_income.json")
@@ -86,7 +90,11 @@ def run_convert_pd_json_chained(client):
     # testpath = "/Users/<user>/<path_to_pyincore>/pyincore/tests/pyincore/utils"
     if testpath:
         pd_process = PopDislOutputProcess(None, os.path.join(testpath, "joplin-pop-disl-results.csv"))
-
+        # pd_process = PopDislOutputProcess(None,
+        #                                   os.path.join(testpath, "joplin-pop-disl-results.csv"),
+        #                                   filter_name="Joplin",
+        #                                   filter_guid=True,
+        #                                   vacant_disl=True)
         pd_process.pd_by_race("PD_by_race.json")
         pd_process.pd_by_income("PD_by_income.json")
         pd_process.pd_by_tenure("PD_by_tenure.json")