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

Issue:
review and implement seaside restoration model
review Hesam's seaside restoration 
![image](https://user-images.githubusercontent.com/13950475/142280554-aa7ec4b8-34a5-4159-a2e1-40a219cbc7e7.png)



Discussion / hints:

Fix patch:
diff --git a/docs/source/modules.rst b/docs/source/modules.rst
index ee6b5ff92..98262622b 100644
--- a/docs/source/modules.rst
+++ b/docs/source/modules.rst
@@ -125,6 +125,31 @@ analyses/housingunitallocation
 ..  autoclass:: housingunitallocation.housingunitallocation.HousingUnitAllocation
     :members:
 
+analyses/indp
+==============================
+..  autoclass:: indp.indp.INDP
+    :members:
+..  autoclass:: indp.dislocationutils.DislocationUtil
+    :members:
+..  autoclass:: indp.indpcomponents.INDPComponents
+    :members:
+..  autoclass:: indp.indpresults.INDPResults
+    :members:
+..  autoclass:: indp.indputil.INDPUtil
+    :members:
+..  autoclass:: indp.infrastructurearc.InfrastructureArc
+    :members:
+..  autoclass:: indp.infrastructureinterdeparc.InfrastructureInterdepArc
+    :members:
+..  autoclass:: indp.infrastructurenetwork.InfrastructureNetwork
+    :members:
+..  autoclass:: indp.infrastructurenode.InfrastructureNode
+    :members:
+..  autoclass:: indp.infrastructurespace.InfrastructureSpace
+    :members:
+..  autoclass:: indp.infrastructureutil.InfrastructureUtil
+    :members:
+
 analyses/joplincge
 ==================
 ..  autoclass:: joplincge.joplincge.JoplinCGEModel
diff --git a/pyincore/analyses/epfrestoration/epfrestoration.py b/pyincore/analyses/epfrestoration/epfrestoration.py
index b1e750d72..9a9864131 100644
--- a/pyincore/analyses/epfrestoration/epfrestoration.py
+++ b/pyincore/analyses/epfrestoration/epfrestoration.py
@@ -55,7 +55,8 @@ def run(self):
         damage = self.get_input_dataset("damage").get_csv_reader()
         damage_result = AnalysisUtil.get_csv_table_rows(damage, ignore_first_row=False)
 
-        (inventory_restoration_map, pf_results, time_results, func_results) = self.electricpowerfacility_restoration(
+        (inventory_restoration_map, pf_results, time_results, func_results,
+         repair_times) = self.electricpowerfacility_restoration(
             inventory_list, damage_result, mapping_set, restoration_key, end_time, time_interval, pf_interval, 
             discretized_days)
 
@@ -66,6 +67,8 @@ def run(self):
         self.set_result_csv_data("time_results", pf_results, name="reptime_" + self.get_parameter("result_name"))
         self.set_result_csv_data("func_results", func_results,
                                  name=self.get_parameter("result_name") + "_discretized_restoration")
+        self.set_result_csv_data("repair_times", repair_times, name="full_reptime_" + self.get_parameter(
+            "result_name"))
 
         return True
 
@@ -112,6 +115,7 @@ def electricpowerfacility_restoration(self, inventory_list, damage_result, mappi
 
         time_results = []
         pf_results = []
+        repair_time = {}
 
         for mapping in mapping_set.mappings:
             # get restoration curves
@@ -142,6 +146,8 @@ def electricpowerfacility_restoration(self, inventory_list, damage_result, mappi
                     **new_dict
                 })
 
+            repair_time[restoration_curve_set.id] = restoration_curve_set.calculate_inverse_restoration_rates(time=0.99)
+
         # Compute discretized restoration
         func_result = []
         for dmg in damage_result:
@@ -162,7 +168,11 @@ def electricpowerfacility_restoration(self, inventory_list, damage_result, mappi
 
             func_result.append({"guid": guid, **result_dict})
 
-        return inventory_restoration_map, pf_results, time_results, func_result
+        repair_times = []
+        for inventory in inventory_restoration_map:
+            repair_times.append({"guid": inventory["guid"], **repair_time[inventory["restoration_id"]]})
+
+        return inventory_restoration_map, pf_results, time_results, func_result, repair_times
 
     def get_spec(self):
         return {
@@ -261,6 +271,13 @@ def get_spec(self):
                     'parent_type': '',
                     'description': 'A csv file recording discretized functionality over time',
                     'type': 'incore:epfDiscretizedRestorationFunc'
+                },
+                {
+                    'id': 'repair_times',
+                    'parent_type': '',
+                    'description': 'A csv file recording repair time at full functionality recovery for each guid '
+                                   'and limit state.',
+                    'type': 'incore:epfRepairTime'
                 }
 
             ]
diff --git a/pyincore/analyses/epfrestoration/epfrestorationutil.py b/pyincore/analyses/epfrestoration/epfrestorationutil.py
index f789a1195..6e02b9e85 100644
--- a/pyincore/analyses/epfrestoration/epfrestorationutil.py
+++ b/pyincore/analyses/epfrestoration/epfrestorationutil.py
@@ -2,25 +2,25 @@
 
 
 class EpfRestorationUtil:
-    def __init__(self, epfr: EpfRestoration):
+    def __init__(self, inventory_restoration_map, pf_results, time_results, time_interval, pf_interval, end_time):
 
-        # merge inventory_restoration_map with pf and time tables
-        inventory_restoration_map = epfr.get_output_dataset("inventory_restoration_map").get_dataframe_from_csv()
-        pf_results = epfr.get_output_dataset("pf_results").get_dataframe_from_csv()
-        time_results = epfr.get_output_dataset("time_results").get_dataframe_from_csv()
-        self.pf_results = inventory_restoration_map.merge(pf_results, on="restoration_id").set_index('guid')
-        self.time_results = inventory_restoration_map.merge(time_results, on="restoration_id").set_index('guid')
+        # merge inventory_restoration_map with pf and timetables
+        inventory_restoration_map_df = inventory_restoration_map.get_dataframe_from_csv()
+        pf_results_df = pf_results.get_dataframe_from_csv()
+        time_results_df = time_results.get_dataframe_from_csv()
+        self.pf_results_df = inventory_restoration_map_df.merge(pf_results_df, on="restoration_id").set_index('guid')
+        self.time_results_df = inventory_restoration_map_df.merge(time_results_df, on="restoration_id").set_index('guid')
 
-        self.time_interval = epfr.get_parameter("time_interval")
-        self.pf_interval = epfr.get_parameter("pf_interval")
-        self.end_time = epfr.get_parameter("end_time")
+        self.time_interval = time_interval
+        self.pf_interval = pf_interval
+        self.end_time = end_time
 
     def get_restoration_time(self, guid, damage_state="DS_0", pf=0.99):
         if pf > 1:
             raise ValueError("Percentage of functionality should not be larger than 1!")
 
         state = "time_" + damage_state.replace("DS", "PF")
-        df = self.pf_results.loc[guid].reset_index(drop=True)
+        df = self.pf_results_df.loc[guid].reset_index(drop=True)
         # round up and get the closest
         time = df.loc[(df["percentage_of_functionality"] >= pf) & (df["percentage_of_functionality"] <
                                                                    pf+self.pf_interval), state].values[0]
@@ -32,7 +32,7 @@ def get_percentage_func(self, guid, damage_state="DS_0", time=1):
             raise ValueError("restore time should not be larger than end time for restoration model!")
 
         state = damage_state.replace("DS", "PF")
-        df = self.time_results.loc[guid].reset_index(drop=True)
+        df = self.time_results_df.loc[guid].reset_index(drop=True)
         # round up and get the closest
         pf = df.loc[(df["time"] >= time) & df['time'] < time+self.time_interval, state].values[0]
 
diff --git a/pyincore/analyses/housingunitallocation/housingunitallocation.py b/pyincore/analyses/housingunitallocation/housingunitallocation.py
index 8041ec177..5e3222e5e 100644
--- a/pyincore/analyses/housingunitallocation/housingunitallocation.py
+++ b/pyincore/analyses/housingunitallocation/housingunitallocation.py
@@ -220,7 +220,7 @@ def merge_inventories(self, sorted_housing_unit: pd.DataFrame, sorted_infrastruc
                                   how='outer', left_on=["blockid", "randommergeorder"],
                                   right_on=["blockid", "randommergeorder"],
                                   sort=True, suffixes=("_x", "_y"),
-                                  copy=True, indicator=True, validate="1:1")
+                                  copy=True, indicator=True)
         huap_inventory = huap_inventory.rename(columns={"_merge": "aphumerge"})
         # check for duplicate columns from merge
         huap_inventory = self.compare_merges(sorted_housing_unit.columns,
diff --git a/pyincore/analyses/indp/__init__.py b/pyincore/analyses/indp/__init__.py
new file mode 100644
index 000000000..bb9c776d0
--- /dev/null
+++ b/pyincore/analyses/indp/__init__.py
@@ -0,0 +1,18 @@
+# Copyright (c) 2023 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.analyses.indp.indp import INDP
+from pyincore.analyses.indp.dislocationutils import DislocationUtil
+from pyincore.analyses.indp.indpcomponents import INDPComponents
+from pyincore.analyses.indp.indpresults import INDPResults
+from pyincore.analyses.indp.indputil import INDPUtil
+from pyincore.analyses.indp.infrastructurearc import InfrastructureArc
+from pyincore.analyses.indp.infrastructureinterdeparc import InfrastructureInterdepArc
+from pyincore.analyses.indp.infrastructurenetwork import InfrastructureNetwork
+from pyincore.analyses.indp.infrastructurenode import InfrastructureNode
+from pyincore.analyses.indp.infrastructurespace import InfrastructureSpace
+from pyincore.analyses.indp.infrastructureutil import InfrastructureUtil
diff --git a/pyincore/analyses/indp/dislocationutils.py b/pyincore/analyses/indp/dislocationutils.py
new file mode 100644
index 000000000..785bf464e
--- /dev/null
+++ b/pyincore/analyses/indp/dislocationutils.py
@@ -0,0 +1,158 @@
+"""
+These function compute the change of demand values over time based on population dislocation models for Seaside and
+dislocation time model for Lumberton.
+"""
+
+import math
+import os
+import pickle
+
+import numpy as np
+import pandas as pd
+
+
+class DislocationUtil:
+
+    @staticmethod
+    def create_dynamic_param(params, pop_dislocation, T=1, N=None):
+        """
+        This function computes the change of demand values over time based on population dislocation models.
+
+        Parameters
+        ----------
+        params : dict
+            Parameters that are needed to run the INDP optimization.
+        T : int, optional
+            Number of time steps to optimize over. T=1 shows an iINDP analysis, and T>1 shows a td-INDP. The default is
+            1.
+        N : :class:`~infrastructure.InfrastructureNetwork`
+            The object containing the network data.
+
+        Returns
+        -------
+         dynamic_params : dict
+             Dictionary of dynamic demand value for nodes
+        """
+        dynamic_param_dict = params['DYNAMIC_PARAMS']
+        return_type = dynamic_param_dict['RETURN']
+        dp_dict_col = ['time', 'node', 'current pop', 'total pop']
+        net_names = {'WATER': 1, 'GAS': 2, 'POWER': 3, 'TELECOME': 4}
+        dynamic_params = {}
+
+        output_file = dynamic_param_dict['TESTBED'] + '_pop_dislocation_demands_' + str(params['MAGNITUDE']) + 'yr.pkl'
+        if os.path.exists(output_file):
+            print("\nReading from file...")
+            with open(output_file, 'rb') as f:
+                dynamic_params = pickle.load(f)
+            return dynamic_params
+        for net in dynamic_param_dict['MAPPING'].keys():
+            nn = net_names[net]
+            mapping_data = dynamic_param_dict['MAPPING'][net]
+            dynamic_params[nn] = pd.DataFrame(columns=dp_dict_col)
+            if net == 'POWER':
+                for n, d in N.G.nodes(data=True):
+                    if n[1] != nn:
+                        continue
+                    guid = d['data']['inf_data'].guid
+                    # Find building in the service area of the node/substation
+                    try:
+                        serv_area = mapping_data[mapping_data['substations_guid'] == guid]
+                    except KeyError:
+                        serv_area = mapping_data[mapping_data['node_guid'] == guid]
+                    # compute dynamic_params
+                    num_dilocated = {t: 0 for t in range(T + 1)}
+                    total_pop_node = 0
+                    for _, bldg in serv_area.iterrows():
+                        try:
+                            pop_bldg_dict = pop_dislocation[pop_dislocation['guid'] == bldg['buildings_guid']]
+                        except KeyError:
+                            pop_bldg_dict = pop_dislocation[pop_dislocation['guid'] == bldg['bldg_guid']]
+                        for _, hh in pop_bldg_dict.iterrows():
+                            total_pop_node += hh['numprec'] if ~np.isnan(hh['numprec']) else 0
+                            if hh['dislocated']:
+                                # ..todo Lumebrton dislocation time model. Replace with that of Seaside when available
+                                return_time = DislocationUtil.disloc_time_mode(hh)
+                                for t in range(return_time):
+                                    if t <= T and return_type == 'step_function':
+                                        num_dilocated[t] += hh['numprec'] if ~np.isnan(hh['numprec']) else 0
+                                    elif t <= T and return_type == 'linear':
+                                        pass  # ..todo Add linear return here
+                    for t in range(T + 1):
+                        values = [t, n[0], total_pop_node - num_dilocated[t], total_pop_node]
+                        dynamic_params[n[1]] = dynamic_params[n[1]].append(dict(zip(dp_dict_col, values)),
+                                                                           ignore_index=True)
+            elif net == 'WATER':
+                node_pop = {}
+                for u, v, a in N.G.edges(data=True):
+                    if u[1] != nn or v[1] != nn:
+                        continue
+                    guid = a['data']['inf_data'].guid
+                    # Find building in the service area of the pipe
+                    serv_area = mapping_data[mapping_data['edge_guid'] == guid]
+                    start_node = u[0]
+                    if start_node not in node_pop.keys():
+                        node_pop[start_node] = {'total_pop_node': 0,
+                                                'num_dilocated': {t: 0 for t in range(T + 1)}}
+                    end_node = v[0]
+                    if end_node not in node_pop.keys():
+                        node_pop[end_node] = {'total_pop_node': 0,
+                                              'num_dilocated': {t: 0 for t in range(T + 1)}}
+                    # compute dynamic_params
+                    for _, bldg in serv_area.iterrows():
+                        pop_bldg_dict = pop_dislocation[pop_dislocation['guid'] == bldg['bldg_guid']]
+                        for _, hh in pop_bldg_dict.iterrows():
+                            # half of the arc's demand is assigned to each node
+                            # also, each arc is counted twice both as (u,v) and (v,u)
+                            node_pop[start_node]['total_pop_node'] += hh['numprec'] / 4 \
+                                if ~np.isnan(hh['numprec']) else 0
+                            node_pop[end_node]['total_pop_node'] += hh['numprec'] / 4 \
+                                if ~np.isnan(hh['numprec']) else 0
+                            if hh['dislocated']:
+                                # ..todo Lumebrton dislocation time model. Replace with that of Seaside when available
+                                return_time = DislocationUtil.disloc_time_mode(hh)
+                                for t in range(return_time):
+                                    if t <= T and return_type == 'step_function':
+                                        node_pop[start_node]['num_dilocated'][t] += hh['numprec'] / 4 \
+                                            if ~np.isnan(hh['numprec']) else 0
+                                        node_pop[end_node]['num_dilocated'][t] += hh['numprec'] / 4 \
+                                            if ~np.isnan(hh['numprec']) else 0
+                                    elif t <= T and return_type == 'linear':
+                                        pass  # ..todo Add linear return here
+                for n, val in node_pop.items():
+                    for t in range(T + 1):
+                        values = [t, n, val['total_pop_node'] - val['num_dilocated'][t],
+                                  val['total_pop_node']]
+                        dynamic_params[nn] = dynamic_params[nn].append(dict(zip(dp_dict_col, values)),
+                                                                       ignore_index=True)
+        with open(output_file, 'wb') as f:
+            pickle.dump(dynamic_params, f)
+        return dynamic_params
+
+    @staticmethod
+    def disloc_time_mode(household_data):
+        dt_params = {'DS0': 1.00, 'DS1': 2.33, 'DS2': 2.49, 'DS3': 3.62,
+                     'white': 0.78, 'black': 0.88, 'hispanic': 0.83,
+                     'income': -0.00, 'insurance': 1.06}
+        race_white = 1 if household_data['race'] == 1 else 0
+        race_balck = 1 if household_data['race'] == 2 else 0
+        hispan = household_data['hispan'] if ~np.isnan(household_data['hispan']) else 0
+        # ..todo verify that the explanatory variable correspond to columns in dt_params
+        # ..todo Replace random insurance assumption
+        linear_term = household_data['DS_0'] * dt_params['DS0'] + household_data['DS_1'] * dt_params['DS1'] + \
+            household_data['DS_2'] * dt_params['DS2'] + household_data['DS_3'] * dt_params[
+                          'DS3'] + race_white * dt_params['white'] + race_balck * dt_params['black'] + hispan * \
+            dt_params['hispanic'] + np.random.choice([0, 1], p=[.15, .85]) * dt_params['insurance']
+        # household_data['randincome']/1000*dt_params['income']+\#!!! income data
+        disloc_time = np.exp(linear_term)
+        return_time = math.ceil(disloc_time / 7)  # !!! assuming each time step is one week
+        return return_time
+
+    @staticmethod
+    def dynamic_parameters(N, original_N, t, dynamic_params):
+        for n, d in N.G.nodes(data=True):
+            data = dynamic_params[d['data']['inf_data'].net_id]
+            if d['data']['inf_data'].demand < 0:
+                current_pop = data.loc[(data['node'] == n[0]) & (data['time'] == t), 'current pop'].iloc[0]
+                total_pop = data.loc[(data['node'] == n[0]) & (data['time'] == t), 'total pop'].iloc[0]
+                original_demand = original_N.G.nodes[n]['data']['inf_data'].demand
+                d['data']['inf_data'].demand = original_demand * current_pop / total_pop
diff --git a/pyincore/analyses/indp/indp.py b/pyincore/analyses/indp/indp.py
new file mode 100644
index 000000000..a9dd43029
--- /dev/null
+++ b/pyincore/analyses/indp/indp.py
@@ -0,0 +1,981 @@
+# Copyright (c) 2023 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/
+
+
+import copy
+import os
+import sys
+import time
+
+import pandas as pd
+import pyomo.environ as pyo
+from pyomo.opt import SolverFactory, TerminationCondition
+from pyomo.util.infeasible import log_infeasible_constraints
+
+from pyincore import BaseAnalysis, NetworkDataset
+from pyincore.analyses.indp.dislocationutils import DislocationUtil
+from pyincore.analyses.indp.indpresults import INDPResults
+from pyincore.analyses.indp.indputil import INDPUtil
+from pyincore.analyses.indp.infrastructureutil import InfrastructureUtil
+
+from pyincore import globals as pyglobals
+
+
+class INDP(BaseAnalysis):
+    """
+        This class runs INDP or td-INDP for a given number of time steps and input parameters.This analysis takes a
+        decentralized approach to solve the Interdependent Network Design Problem (INDP), a family of
+        centralized Mixed-Integer Programming (MIP) models, which find the optimal restoration strategy of disrupted
+        networked systems subject to budget and operational constraints.
+
+        Contributors
+            | Science: Hesam Talebiyan
+            | Implementation: Chen Wang and NCSA IN-CORE Dev Team
+
+        Args:
+            incore_client (IncoreClient): Service authentication.
+        """
+
+    def __init__(self, incore_client):
+        super(INDP, self).__init__(incore_client)
+
+    def run(self):
+        # input parameters
+        network_type = self.get_parameter("network_type")
+        sample_range = self.get_parameter("sample_range")
+        MAGS = self.get_parameter("MAGS")
+        filter_sce = None
+        fail_sce_param = {
+            'TYPE': network_type,
+            'SAMPLE_RANGE': sample_range,
+            'MAGS': MAGS,
+            'FILTER_SCE': filter_sce,
+        }
+
+        RC = self.get_parameter('RC')
+        layers = self.get_parameter('layers')
+        method = self.get_parameter('method')
+
+        t_steps = self.get_parameter('t_steps')
+        if t_steps is None:
+            t_steps = 10
+
+        dynamic_params = None
+        dislocation_data_type = self.get_parameter("dislocation_data_type")
+        return_model = self.get_parameter("return_model")
+        testbed_name = self.get_parameter("testbed_name")
+
+        # more advanced; can hide from user for the first implementation
+        # this is for dynamically updating the demands
+        bldgs2elec_dataset = self.get_input_dataset("bldgs2elec")
+        bldgs2wter_dataset = self.get_input_dataset("bldgs2wter")
+
+        if dislocation_data_type is not None and return_model is not None and testbed_name is not \
+                None and bldgs2elec_dataset is not None and bldgs2wter_dataset is not None:
+            dynamic_params = {
+                "TYPE": dislocation_data_type,
+                "RETURN": return_model,
+                "TESTBED": testbed_name,
+                "MAPPING": {'POWER': bldgs2elec_dataset.get_dataframe_from_csv(low_memory=False),
+                            'WATER': bldgs2wter_dataset.get_dataframe_from_csv(low_memory=False)}
+            }
+
+        extra_commodity = self.get_parameter("extra_commodity")
+        time_resource = self.get_parameter("time_resource")
+        if time_resource is None:
+            time_resource = True
+
+        # save model or not; default to False
+        save_model = self.get_parameter("save_model")
+        if save_model is None:
+            save_model = False
+
+        action_result, cost_result, runtime_result = self.run_method(fail_sce_param, RC, layers, method=method,
+                                                                     t_steps=t_steps,
+                                                                     misc={'DYNAMIC_PARAMS': dynamic_params,
+                                                                           'EXTRA_COMMODITY': extra_commodity,
+                                                                           'TIME_RESOURCE': time_resource},
+                                                                     save_model=save_model)
+
+        self.set_result_csv_data("action", action_result, name="actions.csv")
+        self.set_result_csv_data("cost", cost_result, name="costs.csv")
+        self.set_result_csv_data("runtime", runtime_result, name="run_time.csv")
+
+    def run_method(self, fail_sce_param, v_r, layers, method, t_steps=10, misc=None, save_model=False):
+        """
+        This function runs restoration analysis based on INDP or td-INDP for different numbers of resources.
+
+        Args:
+            fail_sce_param (dict): information about damage scenarios.
+            v_r (list): number of resources, if this is a list of floats, each float is interpreted as a different
+            total number of resources, and INDP
+            is run given the total number of resources. If this is a list of lists of floats, each list is interpreted
+            as fixed upper bounds on the number of resources each layer can use (same for all time steps).
+            layers (list): List of layers.
+            method (str): Algorithm type.
+            t_steps (int): Number of time steps of the analysis.
+            misc (dict): A dictionary that contains miscellaneous data needed for the analysis
+            save_model (bool): Flag indicates if the model should be saved or not
+        Returns:
+
+        """
+
+        # input files
+        wf_repair_cost = self.get_input_dataset("wf_repair_cost").get_dataframe_from_csv(low_memory=False)
+        wf_repair_cost['budget'] = wf_repair_cost['budget'].str.split(',')
+        wf_repair_cost['repaircost'] = wf_repair_cost['repaircost'].str.split(',')
+        epf_repair_cost = self.get_input_dataset("epf_repair_cost").get_dataframe_from_csv(low_memory=False)
+        epf_repair_cost['budget'] = epf_repair_cost['budget'].str.split(',')
+        epf_repair_cost['repaircost'] = epf_repair_cost['repaircost'].str.split(',')
+
+        pipeline_restoration_time = self.get_input_dataset("pipeline_restoration_time").get_dataframe_from_csv(
+            low_memory=False)
+        pipeline_repair_cost = self.get_input_dataset("pipeline_repair_cost").get_dataframe_from_csv(low_memory=False)
+
+        powerline_supply_demand_info = self.get_input_dataset("powerline_supply_demand_info").get_dataframe_from_csv(
+            low_memory=False)
+        epf_supply_demand_info = self.get_input_dataset("epf_supply_demand_info").get_dataframe_from_csv(
+            low_memory=False)
+        power_network = NetworkDataset.from_dataset(self.get_input_dataset("power_network"))
+        power_arcs = power_network.links.get_dataframe_from_shapefile().merge(powerline_supply_demand_info, on="guid")
+        power_nodes = power_network.nodes.get_dataframe_from_shapefile().merge(epf_supply_demand_info, on="guid")
+
+        pipeline_supply_demand_info = self.get_input_dataset("pipeline_supply_demand_info").get_dataframe_from_csv(
+            low_memory=False)
+        wf_supply_demand_info = self.get_input_dataset("wf_supply_demand_info").get_dataframe_from_csv(low_memory=False)
+        water_network = NetworkDataset.from_dataset(self.get_input_dataset("water_network"))
+        water_arcs = water_network.links.get_dataframe_from_shapefile().merge(pipeline_supply_demand_info, on="guid")
+        water_nodes = water_network.nodes.get_dataframe_from_shapefile()
+        water_nodes = water_nodes.merge(wf_supply_demand_info, on="guid")
+
+        interdep = self.get_input_dataset("interdep").get_dataframe_from_csv(low_memory=False)
+
+        # get rid of distribution nodes
+        wf_failure_state_df = self.get_input_dataset("wf_failure_state").get_dataframe_from_csv(
+            low_memory=False).dropna()
+        wf_damage_state_df = self.get_input_dataset("wf_damage_state").get_dataframe_from_csv(low_memory=False).dropna()
+        pipeline_failure_state_df = self.get_input_dataset("pipeline_failure_state").get_dataframe_from_csv(
+            low_memory=False).dropna()
+        epf_failure_state_df = self.get_input_dataset("epf_failure_state").get_dataframe_from_csv(
+            low_memory=False).dropna()
+        epf_damage_state_df = self.get_input_dataset("epf_damage_state").get_dataframe_from_csv(
+            low_memory=False).dropna()
+
+        sample_range = self.get_parameter("sample_range")
+        initial_node = INDPUtil.generate_intial_node_failure_state(wf_failure_state_df, epf_failure_state_df,
+                                                                   water_nodes, power_nodes, sample_range)
+        initial_link = INDPUtil.generate_intial_link_failure_state(pipeline_failure_state_df,
+                                                                   water_arcs, power_arcs, sample_range)
+
+        pop_dislocation = self.get_input_dataset("pop_dislocation").get_dataframe_from_csv(low_memory=False)
+
+        wf_restoration_time = self.get_input_dataset("wf_restoration_time").get_dataframe_from_csv(low_memory=False)
+        wf_restoration_time = wf_restoration_time.merge(wf_damage_state_df, on="guid")
+
+        epf_restoration_time = self.get_input_dataset("epf_restoration_time").get_dataframe_from_csv(low_memory=False)
+        epf_restoration_time = epf_restoration_time.merge(epf_damage_state_df, on="guid")
+
+        # results
+        action_result = []
+        cost_result = []
+        runtime_result = []
+        for v_i, v in enumerate(v_r):
+            if method == 'INDP':
+                params = {"NUM_ITERATIONS": t_steps, "OUTPUT_DIR": 'indp_results', "V": v,
+                          "T": 1, 'L': layers, "ALGORITHM": "INDP"}
+            elif method == 'TDINDP':
+                params = {"NUM_ITERATIONS": t_steps, "OUTPUT_DIR": 'tdindp_results', "V": v,
+                          "T": t_steps, 'L': layers, "ALGORITHM": "INDP"}
+                if 'WINDOW_LENGTH' in misc.keys():
+                    params["WINDOW_LENGTH"] = misc['WINDOW_LENGTH']
+            else:
+                raise ValueError('Wrong method name: ' + method + '. We currently only support INDP and TDINDP as  '
+                                                                  'method name')
+
+            params['EXTRA_COMMODITY'] = misc['EXTRA_COMMODITY']
+            params['TIME_RESOURCE'] = misc['TIME_RESOURCE']
+            params['DYNAMIC_PARAMS'] = misc['DYNAMIC_PARAMS']
+            if misc['DYNAMIC_PARAMS']:
+                params['OUTPUT_DIR'] = 'dp_' + params['OUTPUT_DIR']
+
+            print('----Running for resources: ' + str(params['V']))
+            for m in fail_sce_param['MAGS']:
+                for i in fail_sce_param['SAMPLE_RANGE']:
+                    params["SIM_NUMBER"] = i
+                    params["MAGNITUDE"] = m
+
+                    print('---Running Magnitude ' + str(m) + ' sample ' + str(i) + '...')
+                    if params['TIME_RESOURCE']:
+                        print('Computing repair times...')
+
+                        wf_repair_cost_sample = wf_repair_cost.copy()
+                        wf_repair_cost_sample["budget"] = wf_repair_cost_sample['budget'].apply(lambda x: float(x[i]))
+                        wf_repair_cost_sample["repaircost"] = \
+                            wf_repair_cost_sample['repaircost'].apply(lambda x: float(x[i]))
+
+                        epf_repair_cost_sample = epf_repair_cost.copy()
+                        epf_repair_cost_sample["budget"] = epf_repair_cost_sample['budget'].apply(lambda x: float(x[i]))
+                        epf_repair_cost_sample["repaircost"] = \
+                            epf_repair_cost_sample['repaircost'].apply(lambda x: float(x[i]))
+
+                        # logic to read repair time
+                        wf_restoration_time_sample = pd.DataFrame()
+                        for index, row in wf_restoration_time.iterrows():
+                            failure_state = int(row["sample_damage_states"].split(",")[i].split("_")[1])  # DS_0,1,2,3,4
+                            if failure_state == 0:
+                                repairtime = 0
+                            else:
+                                repairtime = row["PF_" + str(failure_state - 1)]
+                            wf_restoration_time_sample = pd.concat([wf_restoration_time_sample,
+                                                                    pd.DataFrame([{"guid": row["guid"],
+                                                                                  "repairtime": repairtime}])],
+                                                                   ignore_index=True)
+                        epf_restoration_time_sample = pd.DataFrame()
+                        for index, row in epf_restoration_time.iterrows():
+                            failure_state = int(row["sample_damage_states"].split(",")[i].split("_")[1])  # DS_0,1,2,3,4
+                            if failure_state == 0:
+                                repairtime = 0
+                            else:
+                                repairtime = row["PF_" + str(failure_state - 1)]
+                            epf_restoration_time_sample = pd.concat([epf_restoration_time_sample,
+                                                                     pd.DataFrame([{"guid": row["guid"],
+                                                                                   "repairtime":repairtime}])],
+                                                                    ignore_index=True)
+
+                        water_nodes, water_arcs, power_nodes, power_arcs = \
+                            INDPUtil.time_resource_usage_curves(power_arcs, power_nodes, water_arcs, water_nodes,
+                                                                wf_restoration_time_sample, wf_repair_cost_sample,
+                                                                pipeline_restoration_time, pipeline_repair_cost,
+                                                                epf_restoration_time_sample, epf_repair_cost_sample)
+
+                    print("Initializing network...")
+                    params["N"] = INDPUtil.initialize_network(power_nodes, power_arcs, water_nodes, water_arcs,
+                                                              interdep, extra_commodity=params["EXTRA_COMMODITY"])
+
+                    if params['DYNAMIC_PARAMS']:
+                        print("Computing dynamic demand based on dislocation data...")
+                        dyn_dmnd = DislocationUtil.create_dynamic_param(params, pop_dislocation, N=params["N"],
+                                                                        T=params["NUM_ITERATIONS"])
+                        params['DYNAMIC_PARAMS']['DEMAND_DATA'] = dyn_dmnd
+
+                    if fail_sce_param['TYPE'] == 'from_csv':
+                        InfrastructureUtil.add_from_csv_failure_scenario(params["N"], sample=i,
+                                                                         initial_node=initial_node,
+                                                                         initial_link=initial_link)
+                    else:
+                        raise ValueError('Wrong failure scenario data type.')
+
+                    if params["ALGORITHM"] == "INDP":
+                        indp_results = self.run_indp(params, layers=params['L'], controlled_layers=params['L'],
+                                                     T=params["T"], save_model=save_model, print_cmd_line=False,
+                                                     co_location=False)
+                        for t in indp_results.results:
+                            actions = indp_results[t]['actions']
+                            costs = indp_results[t]['costs']
+                            runtimes = indp_results[t]['run_time']
+                            for a in actions:
+                                action_result.append({
+                                    "RC": str(v_i),
+                                    "layers": 'L' + str(len(layers)),
+                                    "magnitude": 'm' + str(m),
+                                    "sample_num": str(i),
+                                    "t": str(t),
+                                    "action": a,
+                                })
+
+                            runtime_result.append({
+                                "RC": str(v_i),
+                                "layers": 'L' + str(len(layers)),
+                                "magnitude": 'm' + str(m),
+                                "sample_num": str(i),
+                                "t": str(t),
+                                "runtime": runtimes,
+                            })
+
+                            cost_result.append({
+                                "RC": str(v_i),
+                                "layers": 'L' + str(len(layers)),
+                                "magnitude": 'm' + str(m),
+                                "sample_num": str(i),
+                                "t": str(t),
+                                "Space Prep": str(costs["Space Prep"]),
+                                "Arc": str(costs["Arc"]),
+                                "Node": str(costs["Node"]),
+                                "Over Supply": str(costs["Over Supply"]),
+                                "Under Supply": str(costs["Under Supply"]),
+                                "Flow": str(costs["Flow"]),
+                                "Total": str(costs["Total"]),
+                                "Under Supply Perc": str(costs["Under Supply Perc"]),
+                            })
+                    else:
+                        raise ValueError('Wrong algorithm type.')
+
+        return action_result, cost_result, runtime_result
+
+    def run_indp(self, params, layers=None, controlled_layers=None, functionality=None, T=1, save=True, suffix="",
+                 forced_actions=False, save_model=False, print_cmd_line=True, co_location=True):
+        """
+        This function runs iINDP (T=1) or td-INDP for a given number of time steps and input parameters.
+
+        Args:
+            params (dict): Parameters that are needed to run the INDP optimization.
+            layers (list): List of layers in the interdependent network. The default is 'None', which sets the list
+            to [1, 2, 3].
+            controlled_layers (list): List of layers that are included in the analysis. The default is 'None',
+            which  sets the list equal to layers.
+            functionality (dict): This dictionary is used to assign functionality values elements in the network
+            before  the analysis starts. The
+            default is 'None'.
+            T (int): Number of time steps to optimize over. T=1 shows an iINDP analysis, and T>1 shows a td-INDP.
+            The default is 1.
+            TODO save & suffice aare not exposed to outside should remove it
+            save (bool): If the results should be saved to file. The default is True.
+            suffix (str): The suffix that should be added to the output files when saved. The default is ''.
+
+            forced_actions (bool): If True, the optimizer is forced to repair at least one element in each time step.
+            The default is False.
+            TODO expose this parameter
+            save_model (bool): If the optimization model should be saved to file. The default is False.
+            TODO expose this parameter
+            print_cmd_line (bool): If full information about the analysis should be written to the console. The default
+            is True.
+            TODO expose this parameter
+            co_location (bool): If co-location and geographical interdependency should be considered in the analysis.
+            The default is True.
+
+        Returns:
+             indp_results (INDPResults): `~indputils.INDPResults` object containing the optimal restoration decisions.
+
+        """
+        # Initialize failure scenario.
+        global original_N
+        if functionality is None:
+            functionality = {}
+        if layers is None:
+            layers = [1, 2, 3]
+        if controlled_layers is None:
+            controlled_layers = layers
+
+        interdependent_net = params["N"]
+        if "NUM_ITERATIONS" not in params:
+            params["NUM_ITERATIONS"] = 1
+
+        out_dir_suffix_res = INDPUtil.get_resource_suffix(params)
+        indp_results = INDPResults(params["L"])
+        if T == 1:
+            print("--Running INDP (T=1) or iterative INDP.")
+            if print_cmd_line:
+                print("Num iters=", params["NUM_ITERATIONS"])
+
+            # Run INDP for 1 time step (original INDP).
+            output_dir = params["OUTPUT_DIR"] + '_L' + str(len(layers)) + '_m' + str(
+                params["MAGNITUDE"]) + "_v" + out_dir_suffix_res
+            # Initial calculations.
+            if params['DYNAMIC_PARAMS']:
+                original_N = copy.deepcopy(interdependent_net)  # !!! deepcopy
+                DislocationUtil.dynamic_parameters(interdependent_net, original_N, 0,
+                                                   params['DYNAMIC_PARAMS']['DEMAND_DATA'])
+            v_0 = {x: 0 for x in params["V"].keys()}
+            results = self.indp(interdependent_net, v_0, 1, layers, controlled_layers=controlled_layers,
+                                functionality=functionality, co_location=co_location)
+            indp_results = results[1]
+            if save_model:
+                INDPUtil.save_indp_model_to_file(results[0], output_dir + "/Model", 0)
+            for i in range(params["NUM_ITERATIONS"]):
+                print("-Time Step (iINDP)", i + 1, "/", params["NUM_ITERATIONS"])
+                if params['DYNAMIC_PARAMS']:
+                    DislocationUtil.dynamic_parameters(interdependent_net, original_N, i + 1,
+                                                       params['DYNAMIC_PARAMS']['DEMAND_DATA'])
+                results = self.indp(interdependent_net, params["V"], T, layers, controlled_layers=controlled_layers,
+                                    co_location=co_location, functionality=functionality)
+                indp_results.extend(results[1], t_offset=i + 1)
+                if save_model:
+                    INDPUtil.save_indp_model_to_file(results[0], output_dir + "/Model", i + 1)
+                # Modify network to account for recovery and calculate components.
+                INDPUtil.apply_recovery(interdependent_net, indp_results, i + 1)
+        else:
+            # td-INDP formulations. Includes "DELTA_T" parameter for sliding windows to increase efficiency.
+            # Edit 2/8/16: "Sliding window" now overlaps.
+            num_time_windows = 1
+            time_window_length = T
+            if "WINDOW_LENGTH" in params:
+                time_window_length = params["WINDOW_LENGTH"]
+                num_time_windows = T
+            output_dir = params["OUTPUT_DIR"] + '_L' + str(len(layers)) + "_m" + str(
+                params["MAGNITUDE"]) + "_v" + out_dir_suffix_res
+
+            print("Running td-INDP (T=" + str(T) + ", Window size=" + str(time_window_length) + ")")
+            # Initial percolation calculations.
+            v_0 = {x: 0 for x in params["V"].keys()}
+            results = self.indp(interdependent_net, v_0, 1, layers, controlled_layers=controlled_layers,
+                                functionality=functionality, co_location=co_location)
+            indp_results = results[1]
+            if save_model:
+                INDPUtil.save_indp_model_to_file(results[0], output_dir + "/Model", 0)
+            for n in range(num_time_windows):
+                print("-Time window (td-INDP)", n + 1, "/", num_time_windows)
+                functionality_t = {}
+                # Slide functionality matrix according to sliding time window.
+                if functionality:
+                    for t in functionality:
+                        if t in range(n, time_window_length + n + 1):
+                            functionality_t[t - n] = functionality[t]
+                    if len(functionality_t) < time_window_length + 1:
+                        diff = time_window_length + 1 - len(functionality_t)
+                        max_t = max(functionality_t.keys())
+                        for d in range(diff):
+                            functionality_t[max_t + d + 1] = functionality_t[max_t]
+                # Run td-INDP.
+                results = self.indp(interdependent_net, params["V"], time_window_length + 1, layers,
+                                    controlled_layers=controlled_layers, functionality=functionality_t,
+                                    co_location=co_location)
+                if save_model:
+                    INDPUtil.save_indp_model_to_file(results[0], output_dir + "/Model", n + 1)
+                if "WINDOW_LENGTH" in params:
+                    indp_results.extend(results[1], t_offset=n + 1, t_start=1, t_end=2)
+                    # Modify network for recovery actions and calculate components.
+                    INDPUtil.apply_recovery(interdependent_net, results[1], 1)
+                else:
+                    indp_results.extend(results[1], t_offset=0)
+                    for t in range(1, T):
+                        # Modify network to account for recovery actions.
+                        INDPUtil.apply_recovery(interdependent_net, indp_results, t)
+        # Save results of current simulation.
+        if save:
+            if not os.path.exists(output_dir + '/agents'):
+                os.makedirs(output_dir + '/agents')
+            indp_results.to_csv_layer(output_dir + '/agents', params["SIM_NUMBER"], suffix=suffix)
+
+        return indp_results
+
+    def indp(self, N, v_r, T=1, layers=None, controlled_layers=None, functionality=None, fixed_nodes=None,
+             print_cmd=True, co_location=True):
+        """
+        INDP optimization problem in Pyomo. It also solves td-INDP if T > 1.
+
+        Parameters
+        ----------
+        N : :class:`~infrastructure.InfrastructureNetwork`
+            An InfrastructureNetwork instance.
+        v_r : dict
+            Dictionary of the number of resources of different types in the analysis.
+            If the value is a scalar for a type, it shows the total number of resources of that type for all layers.
+            If the value is a list for a type, it shows the total number of resources of that type given to each layer.
+        T : int, optional
+            Number of time steps to optimize over. T=1 shows an iINDP analysis, and T>1 shows a td-INDP. The default is
+            1.
+        layers : list, optional
+            Layer IDs in N included in the optimization.
+        controlled_layers : list, optional
+            Layer IDs that can be recovered in this optimization. Used for decentralized optimization. The default is
+            None.
+        functionality : dict, optional
+            Dictionary of nodes to functionality values for non-controlled nodes.
+            Used for decentralized optimization. The default is None.
+        fixed_nodes : dict, optional
+            It fixes the functionality of given elements to a given value. The default is None.
+        print_cmd : bool, optional
+            If true, analysis information is written to the console. The default is True.
+        co_location : bool, optional
+            If false, exclude geographical interdependency from the optimization. The default is True.
+        Returns
+        -------
+        : list
+        A list of the form ``[m, results]`` for a successful optimization where m is the optimization model and
+            results is a :class:`~indputils.INDPResults` object generated using  :func:`collect_results`.
+            If :envvar:`solution_pool` is set to a number, the function returns ``[m, results,  sol_pool_results]``
+            where `sol_pool_results` is dictionary of solution that should be retrieved from the optimizer in
+            addition to the optimal one collected using :func:`collect_solution_pool`.
+
+        """
+        if functionality is None:
+            functionality = {}
+        if layers is None:
+            layers = [1, 2, 3]
+        if controlled_layers is None:
+            controlled_layers = layers
+
+        start_time = time.time()
+
+        m = pyo.ConcreteModel()
+        m.T = T
+        m.N = N
+        m.v_r = v_r
+        m.functionality = functionality
+
+        '''Sets and Dictionaries'''
+        g_prime_nodes = [n[0] for n in N.G.nodes(data=True) if n[1]['data']['inf_data'].net_id in layers]
+        g_prime = N.G.subgraph(g_prime_nodes)
+        # Nodes in controlled network.
+        m.n_hat_nodes = pyo.Set(
+            initialize=[n[0] for n in g_prime.nodes(data=True) if n[1]['data']['inf_data'].net_id in controlled_layers])
+        m.n_hat = g_prime.subgraph(m.n_hat_nodes.ordered_data())
+        # Damaged nodes in controlled network.
+        m.n_hat_prime_nodes = pyo.Set(
+            initialize=[n[0] for n in m.n_hat.nodes(data=True) if n[1]['data']['inf_data'].repaired == 0.0])
+        n_hat_prime = [n for n in m.n_hat.nodes(data=True) if n[1]['data']['inf_data'].repaired == 0.0]
+        # Arcs in controlled network.
+        m.a_hat = pyo.Set(
+            initialize=[(u, v) for u, v, a in g_prime.edges(data=True) if
+                        a['data']['inf_data'].layer in controlled_layers])
+        # Damaged arcs in whole network
+        m.a_prime = pyo.Set(
+            initialize=[(u, v) for u, v, a in g_prime.edges(data=True) if a['data']['inf_data'].functionality == 0.0])
+        a_prime = [(u, v, a) for u, v, a in g_prime.edges(data=True) if a['data']['inf_data'].functionality == 0.0]
+        # Damaged arcs in controlled network.
+        m.a_hat_prime = pyo.Set(initialize=[(u, v) for u, v, _ in a_prime if m.n_hat.has_node(u)
+                                            and m.n_hat.has_node(v)])
+        a_hat_prime = [(u, v, a) for u, v, a in a_prime if m.n_hat.has_node(u) and m.n_hat.has_node(v)]
+        # Sub-spaces
+        m.S = pyo.Set(initialize=N.S)
+
+        # Populate interdependencies. Add nodes to N' if they currently rely on a non-functional node.
+        m.interdep_nodes = {}
+        for u, v, a in g_prime.edges(data=True):
+            if not functionality:
+                if a['data']['inf_data'].is_interdep and g_prime.nodes[u]['data']['inf_data'].functionality == 0.0:
+                    # print "Dependency edge goes from:",u,"to",v
+                    if v not in m.interdep_nodes:
+                        m.interdep_nodes[v] = []
+                    m.interdep_nodes[v].append((u, a['data']['inf_data'].gamma))
+            else:
+                # Should populate m.n_hat with layers that are controlled. Then go through m.n_hat.edges(data=True)
+                # to find interdependencies.
+                for t in range(T):
+                    if t not in m.interdep_nodes:
+                        m.interdep_nodes[t] = {}
+                    if m.n_hat.has_node(v) and a['data']['inf_data'].is_interdep:
+                        if functionality[t][u] == 0.0:
+                            if v not in m.interdep_nodes[t]:
+                                m.interdep_nodes[t][v] = []
+                            m.interdep_nodes[t][v].append((u, a['data']['inf_data'].gamma))
+
+        '''Variables'''
+        m.time_step = pyo.Set(initialize=range(T))
+        # Add geographical space variables.
+        if co_location:
+            m.S_ids = pyo.Set(initialize=[s.id for s in m.S.value])
+            m.z = pyo.Var(m.S_ids, m.time_step, domain=pyo.Binary)
+        # Add functionality binary variables for each node in N'.
+        m.w = pyo.Var(m.n_hat_nodes, m.time_step, domain=pyo.Binary)
+        if T > 1:
+            m.w_tilde = pyo.Var(m.n_hat_nodes, m.time_step, domain=pyo.Binary)
+        # Add functionality binary variables for each arc in A'.
+        m.y = pyo.Var(m.a_hat_prime, m.time_step, domain=pyo.Binary)
+        if T > 1:
+            m.y_tilde = pyo.Var(m.a_hat_prime, m.time_step, domain=pyo.Binary)
+        # Add variables considering extra commodity in addition to the base one
+        node_com_idx = []
+        for n, d in m.n_hat.nodes(data=True):
+            node_com_idx.append((n, 'b'))
+            for key, val in d['data']['inf_data'].extra_com.items():
+                node_com_idx.append((n, key))
+        arc_com_idx = []
+        for u, v, a in m.n_hat.edges(data=True):
+            arc_com_idx.append((u, v, 'b'))
+            for key, val in a['data']['inf_data'].extra_com.items():
+                arc_com_idx.append((u, v, key))
+        # Add over/under-supply variables for each node.
+        m.delta_p = pyo.Var(node_com_idx, m.time_step, domain=pyo.NonNegativeReals)
+        m.delta_m = pyo.Var(node_com_idx, m.time_step, domain=pyo.NonNegativeReals)
+        # Add flow variables for each arc. (main commodity)
+        m.x = pyo.Var(arc_com_idx, m.time_step, domain=pyo.NonNegativeReals)
+
+        # Fix node values
+        if fixed_nodes:
+            for key, val in fixed_nodes.items():
+                if m.T == 1:
+                    m.w[key].fix(val)
+                else:
+                    m.w_tilde[key].fix(val)
+
+        '''Populate objective function'''
+        obj_func = 0
+        for t in range(m.T):
+            if co_location:
+                for s in N.S:
+                    obj_func += s.cost * m.z[s.id, t]
+            for u, v, a in a_hat_prime:
+                if T == 1:
+                    obj_func += (float(a['data']['inf_data'].reconstruction_cost) / 2.0) * m.y[u, v, t]
+                else:
+                    obj_func += (float(a['data']['inf_data'].reconstruction_cost) / 2.0) * m.y_tilde[u, v, t]
+            for n, d in n_hat_prime:
+                if T == 1:
+                    obj_func += d['data']['inf_data'].reconstruction_cost * m.w[n, t]
+                else:
+                    obj_func += d['data']['inf_data'].reconstruction_cost * m.w_tilde[n, t]
+            for n, d in m.n_hat.nodes(data=True):
+                obj_func += d['data']['inf_data'].oversupply_penalty * m.delta_p[n, 'b', t]
+                obj_func += d['data']['inf_data'].undersupply_penalty * m.delta_m[n, 'b', t]
+                for layer, val in d['data']['inf_data'].extra_com.items():
+                    obj_func += val['oversupply_penalty'] * m.delta_p[n, layer, t]
+                    obj_func += val['undersupply_penalty'] * m.delta_m[n, layer, t]
+            for u, v, a in m.n_hat.edges(data=True):
+                obj_func += a['data']['inf_data'].flow_cost * m.x[u, v, 'b', t]
+                for layer, val in a['data']['inf_data'].extra_com.items():
+                    obj_func += val['flow_cost'] * m.x[u, v, layer, t]
+        m.Obj = pyo.Objective(rule=obj_func, sense=pyo.minimize)
+
+        '''Constraints'''
+        # Time-dependent constraints.
+        if m.T > 1:
+            m.initial_state_node = pyo.Constraint(m.n_hat_prime_nodes,
+                                                  rule=INDPUtil.initial_state_node_rule,
+                                                  doc='Initialstate at node')
+            m.initial_state_arc = pyo.Constraint(m.a_hat_prime, rule=INDPUtil.initial_state_arc_rule,
+                                                 doc='Initial state at arc')
+        m.time_dependent_node = pyo.Constraint(m.n_hat_prime_nodes,
+                                               m.time_step,
+                                               rule=INDPUtil.time_dependent_node_rule,
+                                               doc='Time dependent recovery constraint at node')
+        m.time_dependent_arc = pyo.Constraint(m.a_hat_prime,
+                                              m.time_step,
+                                              rule=INDPUtil.time_dependent_arc_rule,
+                                              doc='Time dependent recovery constraint at arc')
+        # Enforce a_i,j to be fixed if a_j,i is fixed (and vice versa).
+        m.arc_equality = pyo.Constraint(m.a_hat_prime,
+                                        m.time_step,
+                                        rule=INDPUtil.arc_equality_rule,
+                                        doc='Arc reconstruction equality')
+
+        # Conservation of flow constraint. (2) in INDP paper.
+        m.flow_conserv_node = pyo.Constraint(m.delta_p_index_0,
+                                             m.time_step,
+                                             rule=INDPUtil.flow_conserv_node_rule,
+                                             doc='Flow conservation')
+        # Flow functionality constraints.
+        m.flow_in_functionality = pyo.Constraint(m.a_hat, m.time_step, rule=INDPUtil.flow_in_functionality_rule,
+                                                 doc='Flow In Functionality')
+        m.flow_out_functionality = pyo.Constraint(m.a_hat, m.time_step, rule=INDPUtil.flow_out_functionality_rule,
+                                                  doc='Flow Out Functionality')
+        m.flow_arc_functionality = pyo.Constraint(m.a_hat, m.time_step, rule=INDPUtil.flow_arc_functionality_rule,
+                                                  doc='Flow Arc Functionality')
+        # Resource availability constraints.
+        m.resource = pyo.Constraint(list(m.v_r.keys()), m.time_step, rule=INDPUtil.resource_rule,
+                                    doc='Resource availability')
+        # Interdependency constraints
+        m.interdependency = pyo.Constraint(m.n_hat_nodes, m.time_step, rule=INDPUtil.interdependency_rule,
+                                           doc='Interdependency')
+        # Geographic space constraints
+        if co_location:
+            m.node_geographic_space = pyo.Constraint(m.S_ids, m.n_hat_prime_nodes, m.time_step,
+                                                     rule=INDPUtil.node_geographic_space_rule,
+                                                     doc='Node Geographic space')
+            m.arc_geographic_space = pyo.Constraint(m.S_ids, m.a_hat_prime, m.time_step,
+                                                    rule=INDPUtil.arc_geographic_space_rule,
+                                                    doc='Arc Geographic space')
+
+        '''Solve'''
+        num_cont_vars = len([v for v in m.component_data_objects(pyo.Var) if v.domain == pyo.NonNegativeReals])
+        num_integer_vars = len([v for v in m.component_data_objects(pyo.Var) if v.domain == pyo.Binary])
+
+        solver_engine = self.get_parameter("solver_engine")
+        if solver_engine is None:
+            solver_engine = "scip"
+
+        solver_path = self.get_parameter("solver_path")
+        if solver_path is None:
+            solver_path = pyglobals.SCIP_PATH
+
+        solver_time_limit = self.get_parameter("solver_time_limit")
+
+        print(
+            "Solving... using %s solver (%d cont. vars, %d binary vars)" %
+            (solver_engine, num_cont_vars, num_integer_vars))
+
+        if solver_engine == "gurobi":
+            solver = SolverFactory(solver_engine, timelimit=solver_time_limit)
+        else:
+            solver = SolverFactory(solver_engine, timelimit=solver_time_limit, executable=solver_path)
+
+        solution = solver.solve(m)
+        run_time = time.time() - start_time
+
+        # Save results.
+        if solution.solver.termination_condition in [TerminationCondition.optimal, TerminationCondition.maxTimeLimit]:
+            if solution.solver.termination_condition == TerminationCondition.maxTimeLimit:
+                print('\nOptimizer time limit, gap = %1.3f\n' % solution.a.solution(0).gap)
+            results = INDPUtil.collect_results(m, controlled_layers, coloc=co_location)
+            results.add_run_time(t, run_time)
+            return [m, results]
+        else:
+            log_infeasible_constraints(m, log_expression=True, log_variables=True)
+            if solution.solver.termination_condition == TerminationCondition.infeasible:
+                print(solution.solver.termination_condition, ": SOLUTION NOT FOUND. (Check data and/or violated "
+                                                             "constraints in the infeasible_model.log).")
+            sys.exit()
+
+    def get_spec(self):
+        return {
+            'name': 'INDP',
+            'description': 'Interdependent Network Design Problem that models the restoration',
+            'input_parameters': [
+                {
+                    'id': 'network_type',
+                    'required': True,
+                    'description': 'type of the network, which is set to `from_csv` for Seaside networks. '
+                                   'e.g. from_csv, incore',
+                    'type': str
+                },
+                {
+                    'id': 'MAGS',
+                    'required': True,
+                    'description': 'sets the earthquake return period.',
+                    'type': list
+                },
+                {
+                    'id': 'sample_range',
+                    'required': True,
+                    'description': 'sets the range of sample scenarios to be analyzed',
+                    'type': range
+                },
+                {
+                    'id': 'dislocation_data_type',
+                    'required': False,
+                    'description': 'type of the dislocation data.',
+                    'type': str
+                },
+                {
+                    'id': 'return_model',
+                    'required': False,
+                    'description': 'type of the model for the return of the dislocated population. '
+                                   'Options: *step_function* and *linear*.',
+                    'type': str
+                },
+                {
+                    'id': 'testbed_name',
+                    'required': False,
+                    'description': 'sets the name of the testbed in analysis',
+                    'type': str
+                },
+                {
+                    'id': 'extra_commodity',
+                    'required': True,
+                    'description': 'multi-commodity parameters dict',
+                    'type': dict
+                },
+                {
+                    'id': 'RC',
+                    'required': True,
+                    'description': 'list of resource caps or the number of available resources in each step of the '
+                                   'analysis. Each item of the list is a dictionary whose items show the type of '
+                                   'resource and the available number of that type of resource. For example: '
+                                   '* If `network_type`=*from_csv*, you have two options:* if, for example, '
+                                   '`R_c`= [{"budget": 3}, {"budget": 6}], then the analysis is done for the cases '
+                                   'when there are 3 and 6 resources available of type "budget" '
+                                   '(total resource assignment).* if, for example, `R_c`= [{"budget": {1:1, 2:1}}, '
+                                   '{"budget": {1:1, 2:2}}, {"budget": {1:3, 2:3}}] and given there are 2 layers,'
+                                   ' then the analysis is done for the case where each layer gets 1 resource of '
+                                   'type "budget", AND the case where layer 1 gets 1 and layer 2 gets 2 resources of '
+                                   'type "budget", AND the case where each layer gets 3 resources of type '
+                                   '"budget" (Prescribed resource for each layer).',
+                    'type': list
+                },
+                {
+                    'id': 'layers',
+                    'required': True,
+                    'description': 'list of layers in the analysis',
+                    'type': list
+                },
+                {
+                    'id': 'method',
+                    'required': True,
+                    'description': 'There are two choices of method: 1. `INDP`: runs Interdependent Network '
+                                   'Restoration Problem (INDP). 2. `TDINDP`: runs time-dependent INDP (td-INDP).  In '
+                                   'both cases, if "TIME_RESOURCE" is True, then the repair time for each element '
+                                   'is considered in devising the restoration plans',
+                    'type': str,
+                },
+                {
+                    'id': 't_steps',
+                    'required': False,
+                    'description': 'Number of time steps of the analysis',
+                    'type': int
+                },
+                {
+                    'id': 'time_resource',
+                    'required': False,
+                    'description': 'if TIME_RESOURCE is True, then the repair time for each element is '
+                                   'considered in devising the restoration plans',
+                    'type': bool
+                },
+                {
+                    'id': 'save_model',
+                    'required': False,
+                    'description': 'If the optimization model should be saved to file. The default is False.',
+                    'type': bool
+                },
+                {
+                    'id': 'solver_engine',
+                    'required': False,
+                    'description': "Solver to use for optimization model. Such as gurobi/glpk/scip, default to scip.",
+                    'type': str
+                },
+                {
+                    'id': 'solver_path',
+                    'required': False,
+                    'description': "Solver to use for optimization model. Such as gurobi/glpk/scip, default to scip.",
+                    'type': str
+                },
+                {
+                    'id': 'solver_time_limit',
+                    'required': False,
+                    'description': "solver time limit in seconds",
+                    'type': int
+                },
+            ],
+            'input_datasets': [
+                {
+                    "id": "wf_repair_cost",
+                    "required": True,
+                    "description": "repair cost for each water facility",
+                    "type": ["incore:repairCost"]
+                },
+                {
+                    'id': 'wf_restoration_time',
+                    "required": True,
+                    'description': 'recording repair time at certain functionality recovery for each class '
+                                   'and limit state.',
+                    'type': ['incore:waterFacilityRepairTime']
+                },
+                {
+                    "id": "epf_repair_cost",
+                    "required": True,
+                    "description": "repair cost for each electric power facility",
+                    "type": ["incore:repairCost"]
+                },
+                {
+                    'id': 'epf_restoration_time',
+                    "required": True,
+                    'description': 'recording repair time at certain functionality recovery for each class '
+                                   'and limit state.',
+                    'type': ['incore:epfRepairTime']
+                },
+                {
+                    "id": "pipeline_repair_cost",
+                    "required": True,
+                    "description": "repair cost for each pipeline",
+                    "type": ["incore:pipelineRepairCost"]
+                },
+                {
+                    'id': 'pipeline_restoration_time',
+                    "required": True,
+                    'description': 'pipeline restoration times',
+                    'type': ['incore:pipelineRestorationVer1']
+                },
+                {
+                    'id': 'power_network',
+                    'required': True,
+                    'description': 'EPN Network Dataset',
+                    'type': ['incore:epnNetwork'],
+                },
+                {
+                    'id': 'water_network',
+                    'required': True,
+                    'description': 'Water Network Dataset',
+                    'type': ['incore:waterNetwork'],
+                },
+                {
+                    'id': 'powerline_supply_demand_info',
+                    'required': True,
+                    'description': 'Supply and demand information for powerlines',
+                    'type': ['incore:powerLineSupplyDemandInfo'],
+                },
+                {
+                    'id': 'epf_supply_demand_info',
+                    'required': True,
+                    'description': 'Supply and demand information for epfs',
+                    'type': ['incore:epfSupplyDemandInfo'],
+                },
+                {
+                    'id': 'wf_supply_demand_info',
+                    'required': True,
+                    'description': 'Supply and demand information for water facilities',
+                    'type': ['incore:waterFacilitySupplyDemandInfo'],
+                },
+                {
+                    'id': 'pipeline_supply_demand_info',
+                    'required': True,
+                    'description': 'Supply and demand information for water pipelines',
+                    'type': ['incore:pipelineSupplyDemandInfo'],
+                },
+                {
+                    "id": "interdep",
+                    "required": True,
+                    "description": "Interdepenency between water and electric power facilities",
+                    "type": ["incore:interdep"]
+                },
+                {
+                    "id": "wf_failure_state",
+                    "required": True,
+                    "description": "MCS failure state of water facilities",
+                    "type": ["incore:sampleFailureState"]
+                },
+                {
+                    "id": "wf_damage_state",
+                    "required": True,
+                    "description": "MCS damage state of water facilities",
+                    "type": ["incore:sampleDamageState"]
+                },
+                {
+                    "id": "pipeline_failure_state",
+                    "required": True,
+                    "description": "failure state of pipeline from pipeline functionality",
+                    "type": ["incore:sampleFailureState"]
+                },
+                {
+                    "id": "epf_failure_state",
+                    "required": True,
+                    "description": "MCS failure state of electric power facilities",
+                    "type": ["incore:sampleFailureState"]
+                },
+                {
+                    "id": "epf_damage_state",
+                    "required": True,
+                    "description": "MCS damage state of electric power facilities",
+                    "type": ["incore:sampleDamageState"]
+                },
+                {
+                    "id": "pop_dislocation",
+                    "required": True,
+                    "description": "Population dislocation output",
+                    "type": ["incore:popDislocation"]
+                },
+                {
+                    "id": "bldgs2elec",
+                    "required": False,
+                    "description": "relation between building and electric power facility",
+                    "type": ["incore:bldgs2elec"]
+                },
+                {
+                    "id": "bldgs2wter",
+                    "required": False,
+                    "description": "relation between building and water facility",
+                    "type": ["incore:bldgs2wter"]
+                },
+            ],
+            'output_datasets': [
+                {
+                    'id': 'action',
+                    'parent_type': '',
+                    'description': 'Restoration action plans',
+                    'type': 'incore:indpAction'
+                },
+                {
+                    'id': 'cost',
+                    'parent_type': '',
+                    'description': 'Restoration cost plans',
+                    'type': 'incore:indpCost'
+                },
+                {
+                    'id': 'runtime',
+                    'parent_type': '',
+                    'description': 'Restoration runtime plans',
+                    'type': 'incore:indpRuntime'
+                }
+            ]
+        }
diff --git a/pyincore/analyses/indp/indpcomponents.py b/pyincore/analyses/indp/indpcomponents.py
new file mode 100644
index 000000000..bc0b27f69
--- /dev/null
+++ b/pyincore/analyses/indp/indpcomponents.py
@@ -0,0 +1,129 @@
+# Copyright (c) 2023 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/
+
+import string
+
+
+class INDPComponents:
+    """
+    This class stores components of a network
+
+    Attributes
+    ----------
+    components : list
+        List of components of the network
+    num_components : int
+        Number of components of the network
+    gc_size : int
+        Size of the largest component of the network
+    """
+
+    def __init__(self):
+        self.components = []
+        self.num_components = 0
+        self.gc_size = 0
+
+    def add_component(self, members, excess_supply):
+        """
+        This function adds a components
+
+        Parameters
+        ----------
+        members : list
+            List of nodes in the component
+        excess_supply : float
+            Excess supply within the component
+
+        Returns
+        -------
+        None.
+
+        """
+        self.components.append((members, excess_supply))
+
+    def to_csv_string(self):
+        """
+        Convert the list of components to a string
+
+        Returns
+        -------
+            : str
+            List of components as a string
+
+        """
+        comp_strings = []
+        for c in self.components:
+            comp = c[0]
+            supp = c[1]
+            comp_string = "/".join(c[0])
+            comp_string += ":" + str(c[1])
+            comp_strings.append(comp_string)
+        return ",".join(comp_strings)
+
+    @classmethod
+    def calculate_components(clss, m, net, t=0, layers=None):
+        """
+        Find the components and the corresponding excess supply
+
+        Parameters
+        ----------
+        m : gurobi.Model
+            The object containing the solved optimization problem.
+        net : networkx.DiGraph
+            The networkx graph object that stores node, arc, and interdependency information
+        t : int
+            Time step. The default is zero.
+        layers : list
+            List of layers in the analysis
+
+        Returns
+        -------
+        indp_components : :class:`~INDPComponents`
+            The object containing the components
+
+        """
+        if layers is None:
+            layers = [1, 2, 3]
+        indp_components = INDPComponents()
+        components = net.get_clusters(layers[0])
+        indp_components.num_components = len(components)
+        indp_components.gc_size = net.gc_size(layers[0])
+        for c in components:
+            total_excess_supply = 0.0
+            members = []
+            for n in c:
+                members.append(str(n[0]) + "." + str(n[1]))
+                excess_supply = 0.0
+                excess_supply += m.getVarByName('delta+_' + str(n) + "," + str(t)).x
+                excess_supply += -m.getVarByName('delta-_' + str(n) + "," + str(t)).x
+                total_excess_supply += excess_supply
+            indp_components.add_component(members, total_excess_supply)
+        return indp_components
+
+    @classmethod
+    def from_csv_string(clss, csv_string):
+        """
+        This functions reads the components data from a string as a :class:`~INDPComponents` object
+
+        Parameters
+        ----------
+        csv_string : str
+            The string containing the component data
+
+        Returns
+        -------
+        indp_components : :class:`~INDPComponents`
+            The object containing the components
+
+        """
+        indp_components = INDPComponents()
+        comps = csv_string
+        for comp in comps:
+            data = string.split(comp, ":")
+            members = string.split(data[0], "/")
+            supp = str.strip(data[1])
+            indp_components.add_component(members, float(supp))
+        return indp_components
diff --git a/pyincore/analyses/indp/indpresults.py b/pyincore/analyses/indp/indpresults.py
new file mode 100644
index 000000000..43067d836
--- /dev/null
+++ b/pyincore/analyses/indp/indpresults.py
@@ -0,0 +1,356 @@
+# Copyright (c) 2023 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.analyses.indp.indpcomponents import INDPComponents
+import os
+
+
+class INDPResults:
+    """
+    This class saves INDP results including actions, costs, run time, and components
+
+    Attributes
+    ----------
+    results : dict
+        Dictionary containing INDP results including actions, costs, run time, and components
+    layers : list
+        List of layers in the analysis
+    results_layer : int
+        Dictionary containing INDP results for each layer including actions, costs, run time, and components
+    """
+    cost_types = ["Space Prep", "Arc", "Node", "Over Supply", "Under Supply", "Flow", "Total", "Under Supply Perc"]
+
+    def __init__(self, layers=None):
+        if layers is None:
+            layers = []
+        self.results = {}
+        self.layers = layers
+        self.results_layer = {layer: {} for layer in layers}
+
+    def __len__(self):
+        return len(self.results)
+
+    def __getitem__(self, index):
+        return self.results[index]
+
+    def extend(self, indp_result, t_offset=0, t_start=0, t_end=0):
+        """
+        This function extends the results to accommodate a new time step.
+
+        Parameters
+        ----------
+        indp_result : :class:`~INDPResults`
+            The current :class:`~INDPResults` object before extension
+        t_offset : int
+            The number of time steps that the current results should be shifted forward. The default is 0.
+        t_start : int
+            The starting time step. The default is 0.
+        t_end : int
+            The ending time step. The default is 0.
+
+        Returns
+        -------
+        None.
+
+        """
+        if t_end == 0:
+            t_end = len(indp_result)
+        for new_t, t in zip([x + t_offset for x in range(t_end - t_start)],
+                            [y + t_start for y in range(t_end - t_start)]):
+            self.results[new_t] = indp_result.results[t]
+        if self.layers:
+            if t_end == 0:
+                t_end = len(indp_result[self.layers[0]])
+            for layer in indp_result.results_layer.keys():
+                for new_t, t in zip([x + t_offset for x in range(t_end - t_start)],
+                                    [y + t_start for y in range(t_end - t_start)]):
+                    self.results_layer[layer][new_t] = indp_result.results_layer[layer][t]
+
+    def add_cost(self, t, cost_type, cost, cost_layer=None):
+        """
+        This function adds cost values to the results.
+
+        Parameters
+        ----------
+        t : int
+            The time steps to which the costs should be added
+        cost_type : str
+            The cost types that is added. The options are: "Space Prep", "Arc", "Node", "Over Supply", "Under Supply",
+            "Flow", "Total", "Under Supply Perc"
+        cost : float
+            The cost value that is added
+        cost_layer : dict
+            The cost value that is added for each layer. The default is None, which adds no value for layers.
+
+        Returns
+        -------
+        None.
+
+        """
+        if cost_layer is None:
+            cost_layer = {}
+        if t not in self.results:
+            self.results[t] = {
+                'costs': {"Space Prep": 0.0, "Arc": 0.0, "Node": 0.0, "Over Supply": 0.0, "Under Supply": 0.0,
+                          "Under Supply Perc": 0.0, "Flow": 0.0, "Total": 0.0}, 'actions': [], 'gc_size': 0,
+                'num_components': 0, 'components': INDPComponents(), 'run_time': 0.0}
+        self.results[t]['costs'][cost_type] = cost
+        if self.layers:
+            for layer in cost_layer.keys():
+                if t not in self.results_layer[layer]:
+                    self.results_layer[layer][t] = {
+                        'costs': {"Space Prep": 0.0, "Arc": 0.0, "Node": 0.0, "Over Supply": 0.0, "Under Supply": 0.0,
+                                  "Under Supply Perc": 0.0, "Flow": 0.0, "Total": 0.0}, 'actions': [], 'gc_size': 0,
+                        'num_components': 0, 'components': INDPComponents(), 'run_time': 0.0}
+                self.results_layer[layer][t]['costs'][cost_type] = cost_layer[layer]
+
+    def add_run_time(self, t, run_time, save_layer=True):
+        """
+        This function adds run time to the results.
+
+        Parameters
+        ----------
+        t : int
+            The time steps to which the run time should be added
+        run_time : float
+            The run time that is added
+        save_layer : bool
+            If the run time should be added for each layer. The default is True.
+
+        Returns
+        -------
+        None.
+
+        """
+        if t not in self.results:
+            self.results[t] = {
+                'costs': {"Space Prep": 0.0, "Arc": 0.0, "Node": 0.0, "Over Supply": 0.0, "Under Supply": 0.0,
+                          "Under Supply Perc": 0.0, "Flow": 0.0, "Total": 0.0}, 'actions': [], 'gc_size': 0,
+                'num_components': 0, 'components': INDPComponents(), 'run_time': 0.0}
+        self.results[t]['run_time'] = run_time
+        if self.layers and save_layer:
+            for layer in self.layers:
+                if t not in self.results_layer[layer]:
+                    self.results_layer[layer][t] = {
+                        'costs': {"Space Prep": 0.0, "Arc": 0.0, "Node": 0.0, "Over Supply": 0.0, "Under Supply": 0.0,
+                                  "Under Supply Perc": 0.0, "Flow": 0.0, "Total": 0.0}, 'actions': [], 'gc_size': 0,
+                        'num_components': 0, 'components': INDPComponents(), 'run_time': 0.0}
+                self.results_layer[layer][t]['run_time'] = run_time
+
+    def add_action(self, t, action, save_layer=True):
+        """
+        This function adds restoration actions to the results.
+
+        Parameters
+        ----------
+        t : int
+            The time steps to which the actions should be added
+        action : list
+            List of actions that are added
+        save_layer : bool
+            If the actions should be added for each layer. The default is True.
+
+        Returns
+        -------
+        None.
+
+        """
+        if t not in self.results:
+            self.results[t] = {
+                'costs': {"Space Prep": 0.0, "Arc": 0.0, "Node": 0.0, "Over Supply": 0.0, "Under Supply": 0.0,
+                          "Under Supply Perc": 0.0, "Flow": 0.0, "Total": 0.0}, 'actions': [], 'gc_size': 0,
+                'num_components': 0, 'components': INDPComponents(), 'run_time': 0.0}
+        self.results[t]['actions'].append(action)
+        if self.layers and save_layer:
+            action_layer = int(action[-1])
+            if t not in self.results_layer[action_layer]:
+                self.results_layer[action_layer][t] = {
+                    'costs': {"Space Prep": 0.0, "Arc": 0.0, "Node": 0.0, "Over Supply": 0.0, "Under Supply": 0.0,
+                              "Under Supply Perc": 0.0, "Flow": 0.0, "Total": 0.0}, 'actions': [], 'gc_size': 0,
+                    'num_components': 0, 'components': INDPComponents(), 'run_time': 0.0}
+            self.results_layer[action_layer][t]['actions'].append(action)
+
+    def add_gc_size(self, t, gc_size):
+        """
+        This function adds the giant component size to the results
+
+        Parameters
+        ----------
+        t : int
+            The time steps to which the giant component size should be added
+        gc_size : int
+            The giant component size that is added
+
+        Returns
+        -------
+        None.
+
+        """
+        if t not in self.results:
+            self.results[t] = {
+                'costs': {"Space Prep": 0.0, "Arc": 0.0, "Node": 0.0, "Over Supply": 0.0, "Under Supply": 0.0,
+                          "Under Supply Perc": 0.0, "Flow": 0.0, "Total": 0.0}, 'actions': [], 'gc_size': 0,
+                'num_components': 0, 'components': INDPComponents(), 'run_time': 0.0}
+        self.results[t]['gc_size'] = gc_size
+
+    def add_num_components(self, t, num_components):
+        """
+        This function adds the number of components to the results
+
+        Parameters
+        ----------
+        t : int
+            The time steps to which the number of components should be added
+        num_components : int
+            The number of components that is added
+
+        Returns
+        -------
+        None.
+
+        """
+        if t not in self.results:
+            self.results[t] = {
+                'costs': {"Space Prep": 0.0, "Arc": 0.0, "Node": 0.0, "Over Supply": 0.0, "Under Supply": 0.0,
+                          "Under Supply Perc": 0.0, "Flow": 0.0, "Total": 0.0}, 'actions': [], 'gc_size': 0,
+                'num_components': 0, 'components': INDPComponents(), 'run_time': 0.0}
+        self.results[t]['num_components'] = num_components
+
+    def add_components(self, t, components):
+        """
+        This function adds the components to the results
+
+        Parameters
+        ----------
+        t : int
+            The time steps to which the number of components should be added
+        components : list
+            The list of components that is added
+
+        Returns
+        -------
+        None.
+
+        """
+        if t not in self.results:
+            self.results[t] = {
+                'costs': {"Space Prep": 0.0, "Arc": 0.0, "Node": 0.0, "Over Supply": 0.0, "Under Supply": 0.0,
+                          "Under Supply Perc": 0.0, "Flow": 0.0, "Total": 0.0}, 'actions': [], 'gc_size': 0,
+                'num_components': 0, 'components': INDPComponents(), 'run_time': 0.0}
+        self.results[t]['components'] = components
+        self.add_num_components(t, components.num_components)
+        self.add_gc_size(t, components.gc_size)
+
+    def to_csv_layer(self, out_dir, sample_num=1, suffix=""):
+        """
+        This function writes the results to file for each layer. The file for each layer are distinguished by "_L" +
+        the layer number.
+
+        Parameters
+        ----------
+        out_dir : str
+            Output directory to which the results should be written.
+        sample_num : int
+            The sample number corresponding to the results, The default is 1.
+        suffix : str
+            The suffix of the file that is being written. The default is "".
+
+        Returns
+        -------
+        None.
+
+        """
+        for layer in self.layers:
+            action_file = out_dir + "/actions_" + str(sample_num) + "_L" + str(layer) + "_" + suffix + ".csv"
+            costs_file = out_dir + "/costs_" + str(sample_num) + "_L" + str(layer) + "_" + suffix + ".csv"
+            run_time_file = out_dir + "/run_time_" + str(sample_num) + "_L" + str(layer) + "_" + suffix + ".csv"
+            with open(action_file, 'w') as f:
+                f.write("t,action\n")
+                for t in self.results_layer[layer]:
+                    for a in self.results_layer[layer][t]['actions']:
+                        f.write(str(t) + "," + a + "\n")
+            with open(run_time_file, 'w') as f:
+                f.write("t,run_time\n")
+                for t in self.results_layer[layer]:
+                    f.write(str(t) + "," + str(self.results_layer[layer][t]['run_time']) + "\n")
+            with open(costs_file, 'w') as f:
+                f.write("t,Space Prep,Arc,Node,Over Supply,Under Supply,Flow,Total,Under Supply Perc\n")
+                for t in self.results_layer[layer]:
+                    costs = self.results_layer[layer][t]['costs']
+                    f.write(str(t) + "," + str(costs["Space Prep"]) + "," + str(costs["Arc"]) + "," + str(
+                        costs["Node"]) + "," + str(costs["Over Supply"]) + "," + str(costs["Under Supply"]) + "," + str(
+                        costs["Flow"]) + "," + str(costs["Total"]) + "," + str(costs["Under Supply Perc"]) + "\n")
+
+    @classmethod
+    def from_csv(clss, out_dir, sample_num=1, suffix=""):
+        """
+        This function reads the results from file.
+
+        Parameters
+        ----------
+        out_dir : str
+            Output directory from which the results should be read
+        sample_num : int
+            The sample number corresponding to the results, The default is 1.
+        suffix : str
+            The suffix of the file that is being read. The default is "".
+
+        Returns
+        -------
+        indp_result: :class:`~INDPResults`
+            The :class:`~INDPResults` object containing the read results
+
+        """
+        action_file = out_dir + "/actions_" + str(sample_num) + "_" + suffix + ".csv"
+        costs_file = out_dir + "/costs_" + str(sample_num) + "_" + suffix + ".csv"
+        perc_file = out_dir + "/percolation_" + str(sample_num) + "_" + suffix + ".csv"
+        comp_file = out_dir + "/components_" + str(sample_num) + "_" + suffix + ".csv"
+        run_time_file = out_dir + "/run_time_" + str(sample_num) + "_" + suffix + ".csv"
+        indp_result = INDPResults()
+        if os.path.isfile(action_file):  # ..todo: component-related results are not currently added to the results.
+            with open(action_file) as f:
+                lines = f.readlines()[1:]
+                for line in lines:
+                    data = line.strip().split(',')
+                    t = int(data[0])
+                    action = str.strip(data[1])
+                    indp_result.add_action(t, action)
+            with open(costs_file) as f:
+                lines = f.readlines()
+                cost_types = lines[0].strip().split(',')[1:]
+                for line in lines[1:]:
+                    data = line.strip().split(',')
+                    t = int(data[0])
+                    costs = data[1:]
+                    for ct in range(len(cost_types)):
+                        indp_result.add_cost(t, cost_types[ct], float(costs[ct]))
+            with open(run_time_file) as f:
+                lines = f.readlines()
+                for line in lines[1:]:
+                    data = line.strip().split(',')
+                    t = int(data[0])
+                    run_time = data[1]
+                    indp_result.add_run_time(t, run_time)
+            # with open(perc_file) as f:
+            #    lines=f.readlines()[1:]
+            #    for line in lines:
+            #        data=string.split(str.strip(line),",")
+            #        t=int(data[0])
+            #        indp_result.add_gc_size(t,int(data[1]))
+            #        indp_result.add_num_components(t,int(data[2]))
+            # with open(comp_file) as f:
+            #    lines=f.readlines()[1:]
+            #    for line in lines:
+            #        data=string.split(str.strip(line),",")
+            #        t=int(data[0])
+            #        comps=data[1:]
+            #        if comps[0]!='':
+            #            indp_result.add_components(t,INDPComponents.from_csv_string(comps))
+            #        else:
+            #            print("Caution: No component.")
+        else:
+            raise ValueError('File does not exist: ' + action_file)
+        return indp_result
diff --git a/pyincore/analyses/indp/indputil.py b/pyincore/analyses/indp/indputil.py
new file mode 100644
index 000000000..96a67b0af
--- /dev/null
+++ b/pyincore/analyses/indp/indputil.py
@@ -0,0 +1,643 @@
+# Copyright (c) 2023 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/
+
+import os
+
+import pyomo.environ as pyo
+from pyomo.core import value
+
+from pyincore.analyses.indp.infrastructureutil import InfrastructureUtil
+from pyincore.analyses.indp.indpresults import INDPResults
+import pandas as pd
+
+class INDPUtil:
+
+    @staticmethod
+    def get_resource_suffix(params):
+        """
+        This function generates the part of the suffix of result folders that pertains to resource cap(s).
+
+        Args:
+            params (dict): Parameters that are needed to run the INDP optimization.
+
+        Returns:
+            out_dir_suffix_res (str): The part of the suffix of result folders that pertains to resource cap(s).
+
+        """
+        out_dir_suffix_res = ''
+        for rc, val in params["V"].items():
+            if isinstance(val, int):
+                if rc != '':
+                    out_dir_suffix_res += rc[0] + str(val)
+                else:
+                    out_dir_suffix_res += str(val)
+            else:
+                out_dir_suffix_res += rc[0] + str(sum([lval for _, lval in val.items()])) + '_fixed_layer_Cap'
+        return out_dir_suffix_res
+
+    @staticmethod
+    def time_resource_usage_curves(power_arcs, power_nodes, water_arcs, water_nodes, wf_restoration_time_sample,
+                                   wf_repair_cost_sample, pipeline_restoration_time, pipeline_repair_cost,
+                                   epf_restoration_time_sample, epf_repair_cost_sample):
+        """
+        This module calculates the repair time for nodes and arcs for the current scenario based on their damage
+        state, and writes them to the input files of INDP. Currently, it is only compatible with NIST testbeds.
+
+        Args:
+            power_arcs (dataframe):
+            power_nodes (dataframe):
+            water_arcs (dataframe):
+            water_nodes (dataframe):
+            wf_restoration_time_sample (dataframe):
+            wf_repair_cost_sample (dataframe):
+            pipeline_restoration_time (dataframe):
+            pipeline_repair_cost (dataframe):
+            epf_restoration_time_sample (dataframe):
+            epf_repair_cost_sample (dataframe):
+
+        Returns:
+            water_nodes:
+            water_arcs:
+            power_nodes:
+            power_arcs:
+
+        """
+        _water_nodes = water_nodes.merge(wf_repair_cost_sample, on='guid', how='left').merge(
+            wf_restoration_time_sample, on='guid', how='left')
+        water_nodes['p_time'] = _water_nodes['repairtime']
+        water_nodes['p_time'].fillna(0, inplace=True)
+        water_nodes['p_budget'] = _water_nodes['budget']
+        water_nodes['p_budget'].fillna(0, inplace=True)
+        water_nodes['q'] = _water_nodes['repaircost']
+        water_nodes['q'].fillna(0, inplace=True)
+
+        _power_nodes = power_nodes.merge(epf_repair_cost_sample, on='guid', how='left').merge(
+            epf_restoration_time_sample, on='guid', how='left')
+        power_nodes['p_time'] = _power_nodes['repairtime']
+        power_nodes['p_time'].fillna(0, inplace=True)
+        power_nodes['p_budget'] = _power_nodes['budget']
+        power_nodes['p_budget'].fillna(0, inplace=True)
+        power_nodes['q'] = _power_nodes['repaircost']
+        power_nodes['q'].fillna(0, inplace=True)
+
+        _water_arcs = water_arcs.merge(pipeline_repair_cost, on='guid', how='left').merge(
+            pipeline_restoration_time, on='guid', how='left')
+        water_arcs['h_time'] = _water_arcs["repair_time"]
+        water_arcs['h_time'].fillna(0, inplace=True)
+        water_arcs['h_budget'] = _water_arcs['budget'].astype(float)
+        water_arcs['h_budget'].fillna(0, inplace=True)
+        water_arcs['f'] = _water_arcs['repaircost'].astype(float)
+        water_arcs['f'].fillna(0, inplace=True)
+
+        power_arcs['h_time'] = 0
+        power_arcs['h_budget'] = 0
+        power_arcs['f'] = 0
+
+        return water_nodes, water_arcs, power_nodes, power_arcs
+
+    @staticmethod
+    def initialize_network(power_nodes, power_arcs, water_nodes, water_arcs, interdep, cost_scale=1.0,
+                           extra_commodity=None):
+        """
+        This function initializes a :class:`~infrastructure.InfrastructureNetwork` object based on network data.
+
+        Args:
+            cost_scale (float): Scales the cost to improve efficiency. The default is 1.0:
+            extra_commodity (dict): Dictionary of commodities other than the default one for each layer of the
+            network. The default is 'None', which means that there is only one commodity per layer.
+
+        Returns:
+            interdep_net (class):`~infrastructure.InfrastructureNetwork` The object containing the network data.
+
+        """
+        interdep_net = InfrastructureUtil.load_infrastructure_array_format_extended(power_nodes,
+                                                                                    power_arcs,
+                                                                                    water_nodes,
+                                                                                    water_arcs,
+                                                                                    interdep,
+                                                                                    cost_scale=cost_scale,
+                                                                                    extra_commodity=extra_commodity)
+        return interdep_net
+
+    @staticmethod
+    def save_indp_model_to_file(model, out_model_dir, t, layer=0, suffix=''):
+        """
+        This function saves pyomo optimization model to file.
+
+        Parameters
+        ----------
+        model : Pyomo.Model
+            Pyomo optimization model
+        out_model_dir : str
+            Directory to which the models should be written
+        t : int
+            The time step corresponding to the model
+        layer : int
+            The layer number corresponding to the model. The default is 0, which means the model includes all layers in
+            the analysis
+        suffix : str
+            The suffix that should be added to files when saved. The default is ''.
+        Returns
+        -------
+        None.
+
+        """
+        if not os.path.exists(out_model_dir):
+            os.makedirs(out_model_dir)
+            # Write models to file
+        l_name = "/Model_t%d_l%d_%s.txt" % (t, layer, suffix)
+        file_id = open(out_model_dir + l_name, 'w')
+        model.pprint(ostream=file_id)
+        file_id.close()
+        # Write solution to file
+        s_name = "/Solution_t%d_l%d_%s.txt" % (t, layer, suffix)
+        file_id = open(out_model_dir + s_name, 'w')
+        for vv in model.component_data_objects(pyo.Var):
+            if vv.value:
+                file_id.write('%s %g\n' % (str(vv), vv.value))
+            else:
+                file_id.write('%s NONE\n' % (str(vv)))
+        file_id.write('Obj: %g' % value(model.Obj))
+        file_id.close()
+
+    @staticmethod
+    def apply_recovery(N, indp_results, t):
+        """
+        This function applies the restoration decisions (solution of INDP) to a Gurobi model by changing the state of
+        repaired elements to functional
+
+        Parameters
+        ----------
+        N : :class:`~infrastructure.InfrastructureNetwork`
+            The model of the interdependent network.
+        indp_results : INDPResults
+            A :class:`~indputils.INDPResults` object containing the optimal restoration decisions.
+        t : int
+            The time step to which the results should apply.
+
+        Returns
+        -------
+        None.
+
+        """
+        for action in indp_results[t]['actions']:
+            if "/" in action:
+                # Edge recovery action.
+                data = action.split("/")
+                src = tuple([int(x) for x in data[0].split(".")])
+                dst = tuple([int(x) for x in data[1].split(".")])
+                N.G[src][dst]['data']['inf_data'].functionality = 1.0
+            else:
+                # Node recovery action.
+                node = tuple([int(x) for x in action.split(".")])
+                # print "Applying recovery:",node
+                N.G.nodes[node]['data']['inf_data'].repaired = 1.0
+                N.G.nodes[node]['data']['inf_data'].functionality = 1.0
+
+    @staticmethod
+    def collect_results(model, controlled_layers, coloc=True):
+        """
+        This function computes the results (actions and costs) of the optimal results and writes them to a
+        :class:`~indputils.INDPResults` object.
+
+        Parameters
+        ----------
+        model : pyomo.model
+            The object containing the the solved optimization problem.
+        controlled_layers : list
+            Layer IDs that can be recovered in this optimization.
+        coloc : bool, optional
+            If false, exclude geographical interdependency from the results. The default is True.
+
+        Returns
+        -------
+        indp_results : INDPResults
+        A :class:`~indputils.INDPResults` object containing the optimal restoration decisions.
+
+        """
+        layers = controlled_layers
+        indp_results = INDPResults(layers)
+        # compute total demand of all layers and each layer
+        total_demand = 0.0
+        total_demand_layer = {layer: 0.0 for layer in layers}
+        for n, d in model.n_hat.nodes(data=True):
+            demand_value = d['data']['inf_data'].demand
+            if demand_value < 0:
+                total_demand += demand_value
+                total_demand_layer[n[1]] += demand_value
+        for t in range(model.T):
+            node_cost = 0.0
+            arc_cost = 0.0
+            flow_cost = 0.0
+            over_supp_cost = 0.0
+            under_supp_cost = 0.0
+            under_supp = 0.0
+            space_prep_cost = 0.0
+            node_cost_layer = {layer: 0.0 for layer in layers}
+            arc_cost_layer = {layer: 0.0 for layer in layers}
+            flow_cost_layer = {layer: 0.0 for layer in layers}
+            over_supp_cost_layer = {layer: 0.0 for layer in layers}
+            under_supp_cost_layer = {layer: 0.0 for layer in layers}
+            under_supp_layer = {layer: 0.0 for layer in layers}
+            space_prep_cost_layer = {layer: 0.0 for layer in layers}  # !!! populate this for each layer
+            # Record node recovery actions.
+            for n in model.n_hat_prime_nodes:
+                if model.T == 1:
+                    node_state = model.w[n, t].value
+                else:
+                    node_state = model.w_tilde[n, t].value
+                if node_state == 1:
+                    action = str(n[0]) + "." + str(n[1])
+                    indp_results.add_action(t, action)
+            # Record edge recovery actions.
+            for i, k, j, kb in model.a_hat_prime:
+                if model.T == 1:
+                    arc_state = model.y[i, k, j, kb, t].value
+                else:
+                    arc_state = model.y_tilde[i, k, j, kb, t].value
+                if arc_state == 1:
+                    action = str(i) + "." + str(k) + "/" + str(j) + "." + str(kb)
+                    indp_results.add_action(t, action)
+            # Calculate space preparation costs.
+            if coloc:
+                for s in model.S.value:
+                    if model.z[s.id, t].value:
+                        space_prep_cost += s.cost * model.z[s.id, t].value
+            indp_results.add_cost(t, "Space Prep", space_prep_cost, space_prep_cost_layer)
+            # Calculate arc preparation costs.
+            for i, k, j, kb in model.a_hat_prime:
+                a = model.n_hat[i, k][j, kb]['data']['inf_data']
+                if model.T == 1:
+                    arc_state = model.y[i, k, j, kb, t].value
+                else:
+                    arc_state = model.y_tilde[i, k, j, kb, t].value
+                arc_cost += (a.reconstruction_cost / 2.0) * arc_state
+                arc_cost_layer[k] += (a.reconstruction_cost / 2.0) * arc_state
+            indp_results.add_cost(t, "Arc", arc_cost, arc_cost_layer)
+            # Calculate node preparation costs.
+            for n in model.n_hat_prime_nodes:
+                d = model.n_hat.nodes[n]['data']['inf_data']
+                if model.T == 1:
+                    node_state = model.w[n, t].value
+                else:
+                    node_state = model.w_tilde[n, t].value
+                node_cost += d.reconstruction_cost * node_state
+                node_cost_layer[n[1]] += d.reconstruction_cost * node_state
+            indp_results.add_cost(t, "Node", node_cost, node_cost_layer)
+            # Calculate under/oversupply costs.
+            for n, d in model.n_hat.nodes(data=True):
+                over_supp_cost += d['data']['inf_data'].oversupply_penalty * model.delta_p[n, 'b', t].value
+                over_supp_cost_layer[n[1]] += d['data']['inf_data'].oversupply_penalty * model.delta_p[n, 'b', t].value
+                under_supp += model.delta_m[n, 'b', t].value
+                under_supp_layer[n[1]] += model.delta_m[n, 'b', t].value / total_demand_layer[n[1]]
+                under_supp_cost += d['data']['inf_data'].undersupply_penalty * model.delta_m[n, 'b', t].value
+                under_supp_cost_layer[n[1]] += d['data']['inf_data'].undersupply_penalty * model.delta_m[
+                    n, 'b', t].value
+            indp_results.add_cost(t, "Over Supply", over_supp_cost, over_supp_cost_layer)
+            indp_results.add_cost(t, "Under Supply", under_supp_cost, under_supp_cost_layer)
+            indp_results.add_cost(t, "Under Supply Perc", under_supp / total_demand, under_supp_layer)
+            # Calculate flow costs.
+            for i, k, j, kb in model.a_hat:
+                a = model.n_hat[i, k][j, kb]['data']['inf_data']
+                flow_cost += a.flow_cost * model.x[i, k, j, kb, 'b', t].value
+                flow_cost_layer[k] += a.flow_cost * model.x[i, k, j, kb, 'b', t].value
+            indp_results.add_cost(t, "Flow", flow_cost, flow_cost_layer)
+            # Calculate total costs.
+            total_lyr = {}
+            total_nd_lyr = {}
+            for layer in layers:
+                total_lyr[layer] = flow_cost_layer[layer] + arc_cost_layer[layer] + node_cost_layer[layer] + \
+                                   over_supp_cost_layer[layer] + under_supp_cost_layer[layer] + \
+                                   space_prep_cost_layer[layer]
+                total_nd_lyr[layer] = space_prep_cost_layer[layer] + arc_cost_layer[layer] + flow_cost \
+                                      + node_cost_layer[layer]
+            indp_results.add_cost(t, "Total", flow_cost + arc_cost + node_cost + over_supp_cost + under_supp_cost +
+                                  space_prep_cost, total_lyr)
+            indp_results.add_cost(t, "Total no disconnection", space_prep_cost + arc_cost + flow_cost + node_cost,
+                                  total_nd_lyr)
+        return indp_results
+
+    @staticmethod
+    def initial_state_node_rule(model, i, k):
+        return model.w[i, k, 0] == 0.0
+
+    @staticmethod
+    def initial_state_arc_rule(model, i, k, j, kb):
+        return model.y[i, k, j, kb, 0] == 0.0
+
+    @staticmethod
+    def time_dependent_node_rule(model, i, k, t):
+        if t > 0:
+            return sum(model.w_tilde[i, k, t_p] for t_p in range(1, t + 1) if t > 0) >= model.w[i, k, t]
+        else:
+            return pyo.Constraint.Skip
+
+    @staticmethod
+    def time_dependent_arc_rule(model, i, k, j, kb, t):
+        if t > 0:
+            return sum(model.y_tilde[i, k, j, kb, t_p] for t_p in range(1, t + 1) if t > 0) >= model.y[i, k, j, kb, t]
+        else:
+            return pyo.Constraint.Skip
+
+    @staticmethod
+    def arc_equality_rule(model, i, k, j, kb, t):
+        try:
+            if model.T == 1:
+                return model.y[i, k, j, kb, t] == model.y[j, kb, i, k, t]
+            else:
+                return model.y[i, k, j, kb, t] == model.y[j, kb, i, k, t], \
+                       model.y_tilde[i, k, j, kb, t] == model.y_tilde[j, kb, i, k, t]
+        except KeyError:
+            return pyo.Constraint.Skip
+
+    @staticmethod
+    def flow_conserv_node_rule(model, i, k, layer, t):
+        d = model.n_hat.nodes[(i, k)]
+        out_flow_constr = 0
+        in_flow_constr = 0
+        demand_constr = 0
+        for u, v, a in model.n_hat.out_edges((i, k), data=True):
+            if layer == 'b' or layer in a['data']['inf_data'].extra_com.keys():
+                out_flow_constr += model.x[u, v, layer, t]
+        for u, v, a in model.n_hat.in_edges((i, k), data=True):
+            if layer == 'b' or layer in a['data']['inf_data'].extra_com.keys():
+                in_flow_constr += model.x[u, v, layer, t]
+        if layer == 'b':
+            demand_constr += d['data']['inf_data'].demand - model.delta_p[i, k, layer, t] \
+                             + model.delta_m[i, k, layer, t]
+        else:
+            demand_constr += d['data']['inf_data'].extra_com[layer]['demand'] - model.delta_p[i, k, layer, t] + \
+                             model.delta_m[i, k, layer, t]
+        return out_flow_constr - in_flow_constr == demand_constr
+
+    @staticmethod
+    def flow_in_functionality_rule(model, i, k, j, kb, t):
+        if not model.functionality:
+            interdep_nodes_list = model.interdep_nodes.keys()  # Interdependent nodes with a damaged dependee node
+        else:
+            interdep_nodes_list = model.interdep_nodes[t].keys()  # Interdependent nodes with a damaged dependee node
+        a = model.n_hat[i, k][j, kb]['data']['inf_data']
+        lhs = model.x[i, k, j, kb, 'b', t]
+        for layer in a.extra_com.keys():
+            lhs += model.x[i, k, j, kb, layer, t]
+        if ((i, k) in model.n_hat_prime_nodes) | ((i, k) in interdep_nodes_list):
+            return lhs <= a.capacity * model.w[i, k, t]
+        else:
+            return lhs <= a.capacity * model.N.G.nodes[(i, k)]['data']['inf_data'].functionality
+
+    @staticmethod
+    def flow_out_functionality_rule(model, i, k, j, kb, t):
+        if not model.functionality:
+            interdep_nodes_list = model.interdep_nodes.keys()  # Interdependent nodes with a damaged dependee node
+        else:
+            interdep_nodes_list = model.interdep_nodes[t].keys()  # Interdependent nodes with a damaged dependee node
+        a = model.n_hat[i, k][j, kb]['data']['inf_data']
+        lhs = model.x[i, k, j, kb, 'b', t]
+        for layer in a.extra_com.keys():
+            lhs += model.x[i, k, j, kb, layer, t]
+        if ((j, kb) in model.n_hat_prime_nodes) | ((j, kb) in interdep_nodes_list):
+            return lhs <= a.capacity * model.w[j, kb, t]
+        else:
+            return lhs <= a.capacity * model.N.G.nodes[(j, kb)]['data']['inf_data'].functionality
+
+    @staticmethod
+    def flow_arc_functionality_rule(model, i, k, j, kb, t):
+        if not model.functionality:
+            interdep_nodes_list = model.interdep_nodes.keys()  # Interdependent nodes with a damaged dependee node
+        else:
+            interdep_nodes_list = model.interdep_nodes[t].keys()  # Interdependent nodes with a damaged dependee node
+        a = model.n_hat[i, k][j, kb]['data']['inf_data']
+        lhs = model.x[i, k, j, kb, 'b', t]
+        for layer in a.extra_com.keys():
+            lhs += model.x[i, k, j, kb, layer, t]
+        if (i, k, j, kb) in model.a_hat_prime:
+            return lhs <= a.capacity * model.y[i, k, j, kb, t]
+        else:
+            return lhs <= a.capacity * model.N.G[(i, k)][(j, kb)]['data']['inf_data'].functionality
+
+    @staticmethod
+    def resource_rule(model, rc, t):
+        resource_dict = model.v_r[rc]
+        is_sep_res = False
+        if isinstance(resource_dict, int):
+            total_resource = resource_dict
+        else:
+            is_sep_res = True
+            total_resource = sum([lval for _, lval in resource_dict.items()])
+            assert len(resource_dict.keys()) == len(model.layers), "The number of resource  values does not match the \
+                                                                    number of layers."
+        resource_left_constr = 0
+        if is_sep_res:
+            res_left_constr_sep = {key: 0 for key in resource_dict.keys()}
+        for i, k, j, kb in model.a_hat_prime:
+            a = model.n_hat[i, k][j, kb]['data']['inf_data']
+            idx_lyr = a.layer
+            res_use = 0.5 * a.resource_usage['h_' + rc]
+            if model.T == 1:
+                resource_left_constr += res_use * model.y[i, k, j, kb, t]
+                if is_sep_res:
+                    res_left_constr_sep[idx_lyr] += res_use * model.y[i, k, j, kb, t]
+            else:
+                resource_left_constr += res_use * model.y_tilde[i, k, j, kb, t]
+                if is_sep_res:
+                    res_left_constr_sep[idx_lyr] += res_use * model.y_tilde[i, k, j, kb, t]
+        for n in model.n_hat_prime_nodes:
+            idx_lyr = n[1]
+            d = model.n_hat.nodes[n]['data']['inf_data']
+            res_use = d.resource_usage['p_' + rc]
+            if model.T == 1:
+                resource_left_constr += res_use * model.w[n, t]
+                if is_sep_res:
+                    res_left_constr_sep[idx_lyr] += res_use * model.w[n, t]
+            else:
+                resource_left_constr += res_use * model.w_tilde[n, t]
+                if is_sep_res:
+                    res_left_constr_sep[idx_lyr] += res_use * model.w_tilde[n, t]
+        if not isinstance(resource_left_constr, int):
+            if not is_sep_res:
+                return resource_left_constr <= total_resource
+            else:
+                return resource_left_constr <= total_resource, [res_left_constr_sep[key] <= lval for key, lval in
+                                                                resource_dict]
+        else:
+            return pyo.Constraint.Skip
+
+    @staticmethod
+    def interdependency_rule(model, i, k, t):
+        if not model.functionality:
+            if (i, k) in model.interdep_nodes:
+                interdep_l_constr = 0
+                interdep_r_constr = 0
+                for interdep in model.interdep_nodes[(i, k)]:
+                    src = interdep[0]
+                    gamma = interdep[1]
+                    if not model.n_hat.has_node(src):
+                        interdep_l_constr += 0
+                    else:
+                        interdep_l_constr += model.w[src, t] * gamma
+                interdep_r_constr += model.w[i, k, t]
+                return interdep_l_constr >= interdep_r_constr
+        else:
+            if (i, k) in model.interdep_nodes[t]:
+                # print interdep_nodes[t]
+                interdep_l_constr = 0
+                interdep_r_constr = 0
+                for interdep in model.interdep_nodes[t][(i, k)]:
+                    src = interdep[0]
+                    gamma = interdep[1]
+                    if not model.n_hat.has_node(src):
+                        print("Forcing", str((i, k)), "to be 0 (dep. on", str(src), ")")
+                        interdep_l_constr += 0
+                    else:
+                        interdep_l_constr += model.w[src, t] * gamma
+                interdep_r_constr += model.w[i, k, t]
+                return interdep_l_constr >= interdep_r_constr
+        return pyo.Constraint.Skip
+
+    @staticmethod
+    def node_geographic_space_rule(model, s, i, k, t):
+        d = model.n_hat.nodes[(i, k)]['data']['inf_data']
+        if d.in_space(s):
+            if model.T == 1:
+                return model.w[(i, k), t] * d.in_space(s) <= model.z[s, t]
+            else:
+                return model.w_tilde[(i, k), t] * d.in_space(s) <= model.z[s, t]
+        return pyo.Constraint.Skip
+
+    @staticmethod
+    def arc_geographic_space_rule(model, s, i, k, j, kb, t):
+        a = model.n_hat[i, k][j, kb]['data']['inf_data']
+        if a.in_space(s):
+            if model.T == 1:
+                return model.y[i, k, j, kb, t] * a.in_space(s) <= model.z[s, t]
+            else:
+                return model.y_tilde[i, k, j, kb, t] * a.in_space(s) <= model.z[s, t]
+        return pyo.Constraint.Skip
+
+    @staticmethod
+    def collect_solution_pool(m, T, n_hat_prime, a_hat_prime):
+        """
+        This function collects the result (list of repaired nodes and arcs) for all feasible solutions in the
+        solution pool.
+
+        Parameters
+        ----------
+        m : gurobi.Model
+            The object containing the solved optimization problem.
+        T : int
+            Number of time steps in the optimization (T=1 for iINDP, and T>=1 for td-INDP).
+        n_hat_prime : list
+            List of damaged nodes in controlled networks.
+        a_hat_prime : list
+            List of damaged arcs in controlled networks.
+
+        Returns
+        -------
+        sol_pool_results : dict
+        A dictionary containing one dictionary per solution that contains list of repaired node and arcs in the
+        solution.
+
+        """
+        sol_pool_results = {}
+        current_sol_count = 0
+        for sol in range(m.SolCount):
+            m.setParam('SolutionNumber', sol)
+            # print(m.PoolObjVal)
+            sol_pool_results[sol] = {'nodes': [], 'arcs': []}
+            for t in range(T):
+                # Record node recovery actions.
+                for n, d in n_hat_prime:
+                    node_var = 'w_tilde_' + str(n) + "," + str(t)
+                    if T == 1:
+                        node_var = 'w_' + str(n) + "," + str(t)
+                    if round(m.getVarByName(node_var).xn) == 1:
+                        sol_pool_results[sol]['nodes'].append(n)
+                # Record edge recovery actions.
+                for u, v, a in a_hat_prime:
+                    arc_var = 'y_tilde_' + str(u) + "," + str(v) + "," + str(t)
+                    if T == 1:
+                        arc_var = 'y_' + str(u) + "," + str(v) + "," + str(t)
+                    if round(m.getVarByName(arc_var).x) == 1:
+                        sol_pool_results[sol]['arcs'].append((u, v))
+            if sol > 0 and sol_pool_results[sol] == sol_pool_results[current_sol_count]:
+                del sol_pool_results[sol]
+            elif sol > 0:
+                current_sol_count = sol
+        return sol_pool_results
+
+    @staticmethod
+    def generate_intial_node_failure_state(wf_failure_state_df, epf_failure_state_df, water_nodes,
+                                           power_nodes, sample_range):
+        # todo change those hard coding
+        network_name = {'Water': 1, 'Power': 3}
+        combined_node_failed_states = pd.DataFrame()
+        for node_fail_state, node_data, network_code in zip([wf_failure_state_df, epf_failure_state_df],
+                                                            [water_nodes, power_nodes],
+                                                            network_name.values()):
+            node_fail_state[[str(x) for x in sample_range]] = \
+                node_fail_state['failure'].str.split(',', expand=True).iloc[:, sample_range.start:sample_range.stop]
+            node_fail_state = node_fail_state.drop(columns=['failure'])
+            node_fail_state['name'] = 'nan'
+            for index, row in node_data.iterrows():
+                node_name = '(' + str(int(row['nodenwid'])) + ',' + str(network_code) + ')'
+                if not pd.isna(row['guid']):
+                    node_fail_state.loc[node_fail_state['guid'] == row['guid'], 'name'] = node_name
+                else:
+                    temp_dict = {**{'name': node_name, 'guid': 'nan'}, **{str(x): 1 for x in sample_range}}
+                    node_fail_state = pd.concat([node_fail_state, pd.DataFrame([temp_dict])], ignore_index=True)
+            combined_node_failed_states = pd.concat([combined_node_failed_states, node_fail_state])
+
+        return combined_node_failed_states.dropna(subset=["name"])
+
+    @staticmethod
+    def generate_intial_link_failure_state(pipeline_failure_state_df, water_arcs, power_arcs, sample_range):
+        # TODO remove hardcoded
+        network_name = {'Water': 1, 'Power': 3}
+        combined_arc_failed_states = pd.DataFrame()
+
+        powerline_failure_state_df = pd.DataFrame()
+        powerline_failure_state_df['name'] = 'nan'
+        for index, row in power_arcs.iterrows():
+            temp_dict = {str(x): 1 for x in sample_range}
+            temp_dict['guid'] = row['guid']
+            powerline_failure_state_df = pd.concat([powerline_failure_state_df, pd.DataFrame([temp_dict])],
+                                                   ignore_index=True)
+            s_node = int(row['fromnode'])
+            e_node = int(row['tonode'])
+            arc_name = '((' + str(s_node) + ',' + str(network_name["Power"]) + '),(' + str(e_node) + ',' + str(
+                network_name["Power"]) + '))'
+            powerline_failure_state_df.loc[powerline_failure_state_df['guid'] == row['guid'], 'name'] = arc_name
+
+        pipeline_failure_state_df[[str(x) for x in sample_range]] = \
+            pipeline_failure_state_df['failure'].str.split(',', expand=True).iloc[:,
+            sample_range.start:sample_range.stop]
+        pipeline_failure_state_df['name'] = 'nan'
+        pipeline_failure_state_df = pipeline_failure_state_df.drop(columns=['failure'])
+        for index, row in water_arcs.iterrows():
+            s_node = int(row['fromnode'])
+            e_node = int(row['tonode'])
+            arc_name = '((' + str(s_node) + ',' + str(network_name["Water"]) + '),(' + str(e_node) + ',' + str(
+                network_name["Water"]) + '))'
+            if not pd.isna(row['guid']):
+                pipeline_failure_state_df.loc[pipeline_failure_state_df['guid'] == row['guid'], 'name'] = arc_name
+            else:
+                temp_dict = {**{'name': arc_name, 'guid': 'nan'}, **{str(x): 1 for x in sample_range}}
+                pipeline_failure_state_df = pd.concat([pipeline_failure_state_df, pd.DataFrame([temp_dict])],
+                                                      ignore_index=True)
+
+        combined_arc_failed_states = pd.concat([combined_arc_failed_states, powerline_failure_state_df])
+        combined_arc_failed_states = pd.concat([combined_arc_failed_states, pipeline_failure_state_df])
+
+        return combined_arc_failed_states.dropna(subset=["name"])
+
+    @staticmethod
+    def generate_distribution_nodes(water_arc_df, water_nodes_df):
+        dist_nodes = list(set().union(water_arc_df['tonode'].unique(), water_arc_df['fromnode'].unique()))
+        for node in dist_nodes:
+            if node not in water_nodes_df["nodenwid"]:
+                water_nodes_df = pd.concat([water_nodes_df, pd.DataFrame([{'utilfcltyc': 'Distribution Node',
+                                                                          'nodenwid': node}])], ignore_index=True)
+        return water_nodes_df
diff --git a/pyincore/analyses/indp/infrastructurearc.py b/pyincore/analyses/indp/infrastructurearc.py
new file mode 100644
index 000000000..95cd6b00a
--- /dev/null
+++ b/pyincore/analyses/indp/infrastructurearc.py
@@ -0,0 +1,109 @@
+# Copyright (c) 2023 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/
+
+
+class InfrastructureArc(object):
+    """
+    This class models an arc in an infrastructure network
+
+    Attributes
+    ----------
+    source : int
+        Start (or head) node id
+    dest : int
+        End (or tail) node id
+    layer : int
+        The id of the layer of the network to which the arc belong
+    failure_probability : float
+        Failure probability of the arc
+    functionality : bool
+        Functionality state of the node
+    repaired : bool
+        If the arc is repaired or not
+    flow_cost : float
+        Unit cost of sending the main commodity through the arc
+    reconstruction_cost : float
+        Reconstruction cost of the arc
+    capacity : float
+        Maximum volume of the commodities that the arc can carry
+    space : int
+        The id of the geographical space where the arc is
+    resource_usage : dict
+        The dictionary that shows how many resource (of each resource type) is employed to repair the arc
+    extra_com : dict
+        The dictionary that shows flow_cost corresponding to commodities other than the main commodity
+    is_interdep : bool
+        If arc represent a normal arc (that carry commodity within a single layer) or physical interdependency between
+        nodes from different layers
+    """
+
+    def __init__(self, source, dest, layer, is_interdep=False):
+        self.source = source
+        self.dest = dest
+        self.layer = layer
+        self.failure_probability = 0.0
+        self.functionality = 1.0
+        self.repaired = 1.0
+        self.flow_cost = 0.0
+        self.reconstruction_cost = 0.0
+        self.resource_usage = {}
+        self.capacity = 0.0
+        self.space = 0
+        self.extra_com = {}
+        self.is_interdep = is_interdep
+
+    def set_extra_commodity(self, extra_commodity):
+        """
+        This function initialize the dictionary for the extra commodities
+
+        Parameters
+        ----------
+        extra_commodity : list
+            List of extra commodities
+
+        Returns
+        -------
+        None.
+
+        """
+        for ec in extra_commodity:
+            self.extra_com[ec] = {'flow_cost': 0}
+
+    def set_resource_usage(self, resource_names):
+        """
+        This function initialize the dictionary for resource usage per all resource types in the analysis
+
+        Parameters
+        ----------
+        resource_names : list
+            List of resource types
+
+        Returns
+        -------
+        None.
+
+        """
+        for rc in resource_names:
+            self.resource_usage[rc] = 0
+
+    def in_space(self, space_id):
+        """
+        This function checks if the arc is in a given space or not
+
+        Parameters
+        ----------
+        space_id :
+            The id of the space that is checked
+
+        Returns
+        -------
+            : bool
+            Returns 1 if the arc is in the space, and 0 otherwise.
+        """
+        if self.space == space_id:
+            return 1
+        else:
+            return 0
diff --git a/pyincore/analyses/indp/infrastructureinterdeparc.py b/pyincore/analyses/indp/infrastructureinterdeparc.py
new file mode 100644
index 000000000..37466d052
--- /dev/null
+++ b/pyincore/analyses/indp/infrastructureinterdeparc.py
@@ -0,0 +1,30 @@
+# Copyright (c) 2023 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.analyses.indp.infrastructurearc import InfrastructureArc
+
+
+class InfrastructureInterdepArc(InfrastructureArc):
+    """
+    This class models a physical interdependency between nodes from two different layers. This class inherits from
+    :class:`InfrastructureArc`, where `source` attributes corresponds to the dependee node, and `dest` corresponds to
+    the depender node. The depender node is non-functional if the corresponding dependee node is non-functional.
+
+    Attributes
+    ----------
+    source_layer : int
+        The id of the layer where the dependee node is
+    dest_layer : int
+        The id of the layer where the depender node is
+    gamma : float
+        The strength of the dependency, which is a number between 0 and 1.
+    """
+
+    def __init__(self, source, dest, source_layer, dest_layer, gamma):
+        super(InfrastructureInterdepArc, self).__init__(source, dest, source_layer, True)
+        self.source_layer = source_layer
+        self.dest_layer = dest_layer
+        self.gamma = gamma
diff --git a/pyincore/analyses/indp/infrastructurenetwork.py b/pyincore/analyses/indp/infrastructurenetwork.py
new file mode 100644
index 000000000..9ec1fc315
--- /dev/null
+++ b/pyincore/analyses/indp/infrastructurenetwork.py
@@ -0,0 +1,172 @@
+# Copyright (c) 2023 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/
+
+import networkx as nx
+
+
+class InfrastructureNetwork(object):
+    """
+    Stores information of the infrastructure network
+
+    Attributes
+    ----------
+    G : networkx.DiGraph
+        The networkx graph object that stores node, arc, and interdependency information
+    S : list
+        List of geographical spaces on which the network lays
+    id : int
+        Id of the network
+    """
+
+    def __init__(self, id):
+        self.G = nx.DiGraph()
+        self.S = []
+        self.id = id
+
+    def copy(self):
+        """
+        This function copies the current :class:`InfrastructureNetwork` object
+
+        Returns
+        -------
+        new_net: :class:`InfrastructureNetwork`
+            Copy of the current infrastructure network object
+        """
+        new_net = InfrastructureNetwork(self.id)
+        new_net.G = self.G.copy()
+        new_net.S = [s for s in self.S]
+        return new_net
+
+    def update_with_strategy(self, player_strategy):
+        """
+        This function modify the functionality of node and arc per a given strategy
+
+        Parameters
+        ----------
+        player_strategy : list
+            Given strategy, where the first list item shows the functionality of nodes, and the second one is for arcs
+
+        Returns
+        -------
+        None.
+
+        """
+        for q in player_strategy[0]:
+            node = q
+            strat = player_strategy[0][q]
+            self.G.node[q]['data']['inf_data'].repaired = round(strat['repair'])
+            self.G.node[q]['data']['inf_data'].functionality = round(strat['w'])
+        for q in player_strategy[1]:
+            src = q[0]
+            dst = q[1]
+            strat = player_strategy[1][q]
+            self.G[src][dst]['data']['inf_data'].repaired = round(strat['repair'])
+            self.G[src][dst]['data']['inf_data'].functionality = round(strat['y'])
+
+    def get_clusters(self, layer):
+        """
+        This function find the clusters in a layer of the network
+
+        Parameters
+        ----------
+        layer : int
+            The id of the desired layer
+
+        Returns
+        -------
+            : list
+            List of layer components
+
+        """
+        g_prime_nodes = [n[0] for n in self.G.nodes(data=True) if
+                         n[1]['data']['inf_data'].net_id == layer and n[1]['data']['inf_data'].functionality == 1.0]
+        g_prime = nx.DiGraph(self.G.subgraph(g_prime_nodes).copy())
+        g_prime.remove_edges_from(
+            [(u, v) for u, v, a in g_prime.edges(data=True) if a['data']['inf_data'].functionality == 0.0])
+        # print nx.connected_components(g_prime.to_undirected())
+        return list(nx.connected_components(g_prime.to_undirected()))
+
+    def gc_size(self, layer):
+        """
+        This function finds the size of the largest component in a layer of the network
+
+        Parameters
+        ----------
+        layer : int
+            The id of the desired layer
+
+        Returns
+        -------
+            : int
+            Size of the largest component in the layer
+        """
+        g_prime_nodes = [n[0] for n in self.G.nodes(data=True) if
+                         n[1]['data']['inf_data'].net_id == layer and n[1]['data']['inf_data'].functionality == 1.0]
+        g_prime = nx.Graph(self.G.subgraph(g_prime_nodes))
+        g_prime.remove_edges_from(
+            [(u, v) for u, v, a in g_prime.edges(data=True) if a['data']['inf_data'].functionality == 0.0])
+        cc = nx.connected_components(g_prime.to_undirected())
+        if cc:
+            # if len(list(cc)) == 1:
+            #    print "I'm here"
+            #    return len(list(cc)[0])
+            # cc_list=list(cc)
+            # print "Length",len(cc_list)
+            # if len(cc_list) == 1:
+            #    return len(cc_list[0])
+            return len(max(cc, key=len))
+        else:
+            return 0
+
+    def to_game_file(self, layers=None):
+        """
+        This function writes the multi-defender security games.
+
+        Parameters
+        ----------
+        layers : list
+            List of layers in the game.
+
+        Returns
+        -------
+        None.
+
+        """
+        if layers is None:
+            layers = [1, 3]
+        with open("../results/indp_gamefile.game") as f:
+            num_players = len(layers)
+            num_targets = len([n for n in self.G.nodes(data=True) if n[1]['data']['inf_data'].net_id in layers])
+            f.write(str(num_players) + "," + str(num_targets) + ",2\n")
+            for layer in range(len(layers)):
+                layer_nodes = [n for n in self.G.nodes(data=True) if n[1]['data']['inf_data'].net_id == layers[layer]]
+                for node in layer_nodes:
+                    def_values = [0.0] * len(layers)
+                    def_values[layer] = abs(node[1]['data']['inf_data'].demand)
+                    atk_values = sum(def_values)
+                    f.write(str(layer) + "\n")
+                    # f.write("0.0,1.0")
+                    for v in def_values:
+                        f.write(str(v) + "\n")
+                    f.write(str(atk_values) + "\n")
+
+    def to_csv(self, filename="infrastructure_adj.csv"):
+        """
+        This function writes the object to a csv file
+
+        Parameters
+        ----------
+        filename : str
+            Name of the file to which the network should be written
+
+        Returns
+        -------
+        None.
+
+        """
+        with open(filename, 'w') as f:
+            for u, v, a in self.G.edges(data=True):
+                f.write(str(u[0]) + "." + str(u[1]) + "," + str(v[0]) + "." + str(v[1]) + "\n")
diff --git a/pyincore/analyses/indp/infrastructurenode.py b/pyincore/analyses/indp/infrastructurenode.py
new file mode 100644
index 000000000..35278a043
--- /dev/null
+++ b/pyincore/analyses/indp/infrastructurenode.py
@@ -0,0 +1,127 @@
+# Copyright (c) 2023 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/
+
+
+class InfrastructureNode(object):
+    """
+    This class models a node in an infrastructure network
+
+    Attributes
+    ----------
+    id : int
+        Node id
+    net_id : int
+        The id of the layer of the network to which the node belong
+    local_id : int
+        Local node id
+    failure_probability : float
+        Failure probability of the node
+    functionality : bool
+        Functionality state of the node
+    repaired : bool
+        If the node is repaired or not
+    reconstruction_cost : float
+        Reconstruction cost of the node
+    oversupply_penalty : float
+        Penalty per supply unit  of the main commodity that is not used for the the node
+    undersupply_penalty : float
+        Penalty per demand unit  of the main commodity that is not satisfied for the the node
+    demand : float
+        Demand or supply value of the main commodity assigned to the node
+    space : int
+        The id of the geographical space where the node is
+    resource_usage : dict
+        The dictionary that shows how many resource (of each resource type) is employed to repair the node
+    extra_com : dict
+        The dictionary that shows demand, oversupply_penalty, and undersupply_penalty corresponding to commodities
+        other than the main commodity
+    """
+
+    def __init__(self, id, net_id, local_id=""):
+        self.id = id
+        self.net_id = net_id
+        if local_id == "":
+            self.local_id = id
+        self.local_id = local_id
+        self.failure_probability = 0.0
+        self.functionality = 1.0
+        self.repaired = 1.0
+        self.reconstruction_cost = 0.0
+        self.oversupply_penalty = 0.0
+        self.undersupply_penalty = 0.0
+        self.demand = 0.0
+        self.space = 0
+        self.resource_usage = {}
+        self.extra_com = {}
+
+    def set_failure_probability(self, failure_probability):
+        """
+        This function sets the failure probability of the node
+
+        Parameters
+        ----------
+        failure_probability : float
+            Assigned failure probability of node
+
+        Returns
+        -------
+        None.
+
+        """
+        self.failure_probability = failure_probability
+
+    def set_extra_commodity(self, extra_commodity):
+        """
+        This function initialize the dictionary for the extra commodities
+
+        Parameters
+        ----------
+        extra_commodity : list
+            List of extra commodities
+
+        Returns
+        -------
+        None.
+
+        """
+        for l in extra_commodity:
+            self.extra_com[l] = {'demand': 0, 'oversupply_penalty': 0, 'undersupply_penalty': 0}
+
+    def set_resource_usage(self, resource_names):
+        """
+        This function initialize the dictionary for resource usage per all resource types in the analysis
+
+        Parameters
+        ----------
+        resource_names : list
+            List of resource types
+
+        Returns
+        -------
+        None.
+
+        """
+        for rc in resource_names:
+            self.resource_usage[rc] = 0
+
+    def in_space(self, space_id):
+        """
+        This function checks if the node is in a given space or not
+
+        Parameters
+        ----------
+        space_id :
+            The id of the space that is checked
+
+        Returns
+        -------
+            : bool
+            Returns 1 if the node is in the space, and 0 otherwise.
+        """
+        if self.space == space_id:
+            return 1
+        else:
+            return 0
diff --git a/pyincore/analyses/indp/infrastructurespace.py b/pyincore/analyses/indp/infrastructurespace.py
new file mode 100644
index 000000000..5892a3630
--- /dev/null
+++ b/pyincore/analyses/indp/infrastructurespace.py
@@ -0,0 +1,22 @@
+# Copyright (c) 2023 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/
+
+
+class InfrastructureSpace(object):
+    """
+    This class models a geographical space.
+
+    Attributes
+    ----------
+    id : int
+        The id of the space
+    cost : float
+        The cost of preparing the space for a repair action
+    """
+
+    def __init__(self, id, cost):
+        self.id = id
+        self.cost = cost
diff --git a/pyincore/analyses/indp/infrastructureutil.py b/pyincore/analyses/indp/infrastructureutil.py
new file mode 100644
index 000000000..9375981b7
--- /dev/null
+++ b/pyincore/analyses/indp/infrastructureutil.py
@@ -0,0 +1,162 @@
+# Copyright (c) 2023 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.analyses.indp.infrastructurenode import InfrastructureNode
+from pyincore.analyses.indp.infrastructurearc import InfrastructureArc
+from pyincore.analyses.indp.infrastructureinterdeparc import InfrastructureInterdepArc
+from pyincore.analyses.indp.infrastructurenetwork import InfrastructureNetwork
+
+
+class InfrastructureUtil:
+
+    @staticmethod
+    def load_infrastructure_array_format_extended(power_nodes, power_arcs, water_nodes, water_arcs,
+                                                  interdep, cost_scale=1.0,
+                                                  extra_commodity=None):
+        """
+        This function reads the infrastructure network from file in the extended format
+
+        Args:
+            power_nodes:
+            power_arcs:
+            water_nodes:
+            water_arcs:
+            interdep:
+            cost_scale (float): The factor by which all cost values has to multiplied. The default is 1.0.
+            extra_commodity (dict): ( only for extended format of input data) List of extra-commodities in the
+            analysis. The default is None, which only considers a main commodity.
+
+        Returns:
+            G (class:`~infrastructure.InfrastructureNetwork`): The object containing the network data.
+
+        """
+        net_names = {'Water': 1, 'Gas': 2, 'Power': 3, 'Telecommunication': 4}  # !!!
+        G = InfrastructureNetwork("Test")
+        global_index = 0
+
+        for index, data in enumerate([power_nodes, water_nodes]):
+            if index == 0:
+                net = net_names["Power"]
+            else:
+                net = net_names["Water"]
+            for v in data.iterrows():
+                try:
+                    node_id = v[1]['ID']
+                except KeyError:
+                    node_id = v[1]['nodenwid']
+                n = InfrastructureNode(global_index, net, int(node_id))
+                G.G.add_node((n.local_id, n.net_id), data={'inf_data': n})
+                global_index += 1
+                node_main_data = G.G.nodes[(n.local_id, n.net_id)]['data']['inf_data']
+                node_main_data.reconstruction_cost = float(v[1]['q_ds_3']) * cost_scale
+                node_main_data.oversupply_penalty = float(v[1]['Mp']) * cost_scale
+                node_main_data.undersupply_penalty = float(v[1]['Mm']) * cost_scale
+                node_main_data.demand = float(v[1]['Demand'])
+                if 'guid' in v[1].index.values:
+                    node_main_data.guid = v[1]['guid']
+                resource_names = [x for x in list(v[1].index.values) if x[:2] == 'p_']
+                if len(resource_names) > 0:
+                    n.set_resource_usage(resource_names)
+                    for rc in resource_names:
+                        n.resource_usage[rc] = v[1][rc]
+                else:
+                    n.resource_usage['p_'] = 1
+                if extra_commodity:
+                    n.set_extra_commodity(extra_commodity[net])
+                    for layer in extra_commodity[net]:
+                        ext_com_data = G.G.nodes[(n.local_id, n.net_id)]['data']['inf_data'].extra_com[layer]
+                        ext_com_data['oversupply_penalty'] = float(v[1]['Mp_' + layer]) * cost_scale
+                        ext_com_data['undersupply_penalty'] = float(v[1]['Mm_' + layer]) * cost_scale
+                        ext_com_data['demand'] = float(v[1]['Demand_' + layer])
+
+        for index, data in enumerate([power_arcs, water_arcs]):
+            if index == 0:
+                net = net_names["Power"]
+            else:
+                net = net_names["Water"]
+            for v in data.iterrows():
+                try:
+                    start_id = v[1]['Start Node']
+                    end_id = v[1]['End Node']
+                except KeyError:
+                    start_id = v[1]['fromnode']
+                    end_id = v[1]['tonode']
+                for duplicate in range(2):
+                    if duplicate == 0:
+                        a = InfrastructureArc(int(start_id), int(end_id), net)
+                    elif duplicate == 1:
+                        a = InfrastructureArc(int(end_id), int(start_id), net)
+                    G.G.add_edge((a.source, a.layer), (a.dest, a.layer), data={'inf_data': a})
+                    arc_main_data = G.G[(a.source, a.layer)][(a.dest, a.layer)]['data']['inf_data']
+                    arc_main_data.flow_cost = float(v[1]['c']) * cost_scale
+                    arc_main_data.reconstruction_cost = float(v[1]['f']) * cost_scale
+                    arc_main_data.capacity = float(v[1]['u'])
+                    if 'guid' in v[1].index.values:
+                        arc_main_data.guid = v[1]['guid']
+                    resource_names = [x for x in list(v[1].index.values) if x[:2] == 'h_']
+                    if len(resource_names) > 0:
+                        a.set_resource_usage(resource_names)
+                        for rc in resource_names:
+                            a.resource_usage[rc] = v[1][rc]
+                    else:
+                        a.resource_usage['h_'] = 1
+                    if extra_commodity:
+                        a.set_extra_commodity(extra_commodity[net])
+                        for layer in extra_commodity[net]:
+                            ext_com_data = \
+                                G.G[(a.source, a.layer)][(a.dest, a.layer)]['data']['inf_data'].extra_com[layer]
+                            ext_com_data['flow_cost'] = float(v[1]['c_' + layer]) * cost_scale
+
+        for v in interdep.iterrows():
+            if v[1]['Type'] == 'Physical':
+                i = int(v[1]['Dependee Node'])
+                net_i = net_names[v[1]['Dependee Network']]
+                j = int(v[1]['Depender Node'])
+                net_j = net_names[v[1]['Depender Network']]
+                a = InfrastructureInterdepArc(i, j, net_i, net_j, gamma=1.0)
+                G.G.add_edge((a.source, a.source_layer), (a.dest, a.dest_layer), data={'inf_data': a})
+                if extra_commodity:
+                    a.set_extra_commodity(extra_commodity[net_i])
+                    a.set_extra_commodity(extra_commodity[net_j])
+
+        return G
+
+    @staticmethod
+    def add_from_csv_failure_scenario(G, sample, initial_node, initial_link):
+        """
+        This function reads initial damage data from file in the from_csv format, and applies it to the infrastructure
+        network. This format only considers one magnitude value (0), and there can be as many samples from that magnitude.
+
+        Parameters
+        ----------
+        G : :class:`~infrastructure.InfrastructureNetwork`
+             The object containing the network data.
+        sample : int
+            Sample number of the initial damage scenario,
+
+        Returns
+        -------
+        None.
+
+        """
+        for index, row in initial_node.iterrows():
+            raw_n = row["name"].split(",")
+            n = (int(raw_n[0].strip(' )(')), int(raw_n[1].strip(' )(')))
+            state = float(row[sample + 1])
+            G.G.nodes[n]['data']['inf_data'].functionality = state
+            G.G.nodes[n]['data']['inf_data'].repaired = state
+
+        for index, row in initial_link.iterrows():
+            raw_uv = row["name"].split(",")
+            u = (int(raw_uv[0].strip(' )(')), int(raw_uv[1].strip(' )(')))
+            v = (int(raw_uv[2].strip(' )(')), int(raw_uv[3].strip(' )(')))
+            state = float(row[sample + 1])
+            if state == 0.0:
+                G.G[u][v]['data']['inf_data'].functionality = state
+                G.G[u][v]['data']['inf_data'].repaired = state
+
+                G.G[v][u]['data']['inf_data'].functionality = state
+                G.G[v][u]['data']['inf_data'].repaired = state
diff --git a/pyincore/analyses/waterfacilityrestoration/waterfacilityrestoration.py b/pyincore/analyses/waterfacilityrestoration/waterfacilityrestoration.py
index 36cca9974..f864af626 100644
--- a/pyincore/analyses/waterfacilityrestoration/waterfacilityrestoration.py
+++ b/pyincore/analyses/waterfacilityrestoration/waterfacilityrestoration.py
@@ -53,12 +53,15 @@ def run(self):
         if discretized_days is None:
             discretized_days = [1, 3, 7, 30, 90]
 
-        damage = self.get_input_dataset("damage").get_csv_reader()
-        damage_result = AnalysisUtil.get_csv_table_rows(damage, ignore_first_row=False)
+        if self.get_input_dataset("damage") is not None:
+            damage = self.get_input_dataset("damage").get_csv_reader()
+            damage_result = AnalysisUtil.get_csv_table_rows(damage, ignore_first_row=False)
+        else:
+            damage_result = None
 
-        (inventory_restoration_map, pf_results, time_results, func_results) = self.waterfacility_restoration(
-            inventory_list, damage_result, mapping_set, restoration_key, end_time, time_interval, pf_interval,
-            discretized_days)
+        (inventory_restoration_map, pf_results, time_results, func_results, repair_times) = \
+            self.waterfacility_restoration(inventory_list, damage_result, mapping_set, restoration_key, end_time,
+                                           time_interval, pf_interval, discretized_days)
 
         self.set_result_csv_data("inventory_restoration_map", inventory_restoration_map,
                                  name="inventory_restoration_map_" + self.get_parameter("result_name"))
@@ -67,6 +70,8 @@ def run(self):
         self.set_result_csv_data("time_results", pf_results, name="reptime_" + self.get_parameter("result_name"))
         self.set_result_csv_data("func_results", func_results,
                                  name=self.get_parameter("result_name") + "_discretized_restoration")
+        self.set_result_csv_data("repair_times", repair_times, name="full_reptime_" + self.get_parameter(
+            "result_name"))
 
         return True
 
@@ -107,15 +112,18 @@ def waterfacility_restoration(self, inventory_list, damage_result, mapping_set,
 
             inventory_restoration_map.append({"guid": inventory['properties']['guid'],
                                               "restoration_id": restoration_set_id})
-            restoration_curve_set = restoration_sets[inventory["id"]]
+            if inventory["id"] in restoration_sets.keys():
+                restoration_curve_set = restoration_sets[inventory["id"]]
 
-            # For each facility, get the discretized restoration from the continuous curve
-            discretized_restoration = AnalysisUtil.get_discretized_restoration(restoration_curve_set, discretized_days)
-            inventory_class_map[inventory['properties']['guid']] = discretized_restoration
+                # For each facility, get the discretized restoration from the continuous curve
+                discretized_restoration = AnalysisUtil.get_discretized_restoration(restoration_curve_set, discretized_days)
+                inventory_class_map[inventory['properties']['guid']] = discretized_restoration
 
         time_results = []
         pf_results = []
 
+        repair_time = {}
+
         for mapping in mapping_set.mappings:
             # get restoration curves
             # if it's string:id; then need to fetch it from remote and cast to restorationcurveset object
@@ -144,28 +152,39 @@ def waterfacility_restoration(self, inventory_list, damage_result, mapping_set,
                     "percentage_of_functionality": p,
                     **new_dict
                 })
+
+            repair_time[restoration_curve_set.id] = \
+                restoration_curve_set.calculate_inverse_restoration_rates(time=0.99)
+
         # Compute discretized restoration
         func_result = []
 
-        for dmg in damage_result:
-            guid = dmg['guid']
+        if damage_result is not None:
+            for dmg in damage_result:
+                guid = dmg['guid']
+
+                # Dictionary of discretized restoration functionality
+                if guid in inventory_class_map.keys():
+                    rest_dict = inventory_class_map[guid]
 
-            # Dictionary of discretized restoration functionality
-            rest_dict = inventory_class_map[guid]
+                    ds_0, ds_1, ds_2, ds_3, ds_4 = dmg['DS_0'], dmg['DS_1'], dmg['DS_2'], dmg['DS_3'], dmg['DS_4']
+                    result_dict = {}
+                    for time in discretized_days:
+                        key = "day" + str(time)
+                        # Only compute if we have damage
+                        if ds_0:
+                            functionality = (rest_dict[key][0] * float(ds_0) + rest_dict[key][1] * float(ds_1) + rest_dict[
+                                key][2] * float(ds_2) + rest_dict[key][3] * float(ds_3) + rest_dict[key][4] * float(ds_4))
+                            result_dict.update({str(key): functionality})
 
-            ds_0, ds_1, ds_2, ds_3, ds_4 = dmg['DS_0'], dmg['DS_1'], dmg['DS_2'], dmg['DS_3'], dmg['DS_4']
-            result_dict = {}
-            for time in discretized_days:
-                key = "day" + str(time)
-                # Only compute if we have damage
-                if ds_0:
-                    functionality = (rest_dict[key][0] * float(ds_0) + rest_dict[key][1] * float(ds_1) + rest_dict[
-                        key][2] * float(ds_2) + rest_dict[key][3] * float(ds_3) + rest_dict[key][4] * float(ds_4))
-                    result_dict.update({str(key): functionality})
+                    func_result.append({"guid": guid, **result_dict})
 
-            func_result.append({"guid": guid, **result_dict})
+        repair_times = []
+        for inventory in inventory_restoration_map:
+            if inventory["restoration_id"] is not None:
+                repair_times.append({"guid": inventory["guid"], **repair_time[inventory["restoration_id"]]})
 
-        return inventory_restoration_map, pf_results, time_results, func_result
+        return inventory_restoration_map, pf_results, time_results, func_result, repair_times
 
     def get_spec(self):
         return {
@@ -226,7 +245,7 @@ def get_spec(self):
                 {
 
                     'id': 'damage',
-                    'required': True,
+                    'required': False,
                     'description': 'damage result that has damage intervals in it',
                     'type': ['ergo:waterFacilityDamageVer6']
                 }
@@ -258,6 +277,13 @@ def get_spec(self):
                     'parent_type': '',
                     'description': 'A csv file recording discretized functionality over time',
                     'type': 'incore:waterFacilityDiscretizedRestorationFunc'
+                },
+                {
+                    'id': 'repair_times',
+                    'parent_type': '',
+                    'description': 'A csv file recording repair time at full functionality recovery for each guid '
+                                   'and limit state.',
+                    'type': 'incore:waterFacilityRepairTime'
                 }
             ]
         }
diff --git a/pyincore/analyses/waterfacilityrestoration/waterfacilityrestorationutil.py b/pyincore/analyses/waterfacilityrestoration/waterfacilityrestorationutil.py
index d8b5290b8..345e2ec23 100644
--- a/pyincore/analyses/waterfacilityrestoration/waterfacilityrestorationutil.py
+++ b/pyincore/analyses/waterfacilityrestoration/waterfacilityrestorationutil.py
@@ -2,25 +2,26 @@
 
 
 class WaterFacilityRestorationUtil:
-    def __init__(self, wfr: WaterFacilityRestoration):
+    def __init__(self, inventory_restoration_map, pf_results, time_results, time_interval, pf_interval, end_time):
 
-        # merge inventory_restoration_map with pf and time tables
-        inventory_restoration_map = wfr.get_output_dataset("inventory_restoration_map").get_dataframe_from_csv()
-        pf_results = wfr.get_output_dataset("pf_results").get_dataframe_from_csv()
-        time_results = wfr.get_output_dataset("time_results").get_dataframe_from_csv()
-        self.pf_results = inventory_restoration_map.merge(pf_results, on="restoration_id").set_index('guid')
-        self.time_results = inventory_restoration_map.merge(time_results, on="restoration_id").set_index('guid')
+        # merge inventory_restoration_map with pf and timetables
+        inventory_restoration_map_df = inventory_restoration_map.get_dataframe_from_csv()
+        pf_results_df = pf_results.get_dataframe_from_csv()
+        time_results_df = time_results.get_dataframe_from_csv()
+        self.pf_results_df = inventory_restoration_map_df.merge(pf_results_df, on="restoration_id").set_index('guid')
+        self.time_results_df = inventory_restoration_map_df.merge(time_results_df, on="restoration_id").set_index(
+            'guid')
 
-        self.time_interval = wfr.get_parameter("time_interval")
-        self.pf_interval = wfr.get_parameter("pf_interval")
-        self.end_time = wfr.get_parameter("end_time")
+        self.time_interval = time_interval
+        self.pf_interval = pf_interval
+        self.end_time = end_time
 
     def get_restoration_time(self, guid, damage_state="DS_0", pf=0.99):
         if pf > 1:
             raise ValueError("Percentage of functionality should not be larger than 1!")
 
         state = "time_" + damage_state.replace("DS", "PF")
-        df = self.pf_results.loc[guid].reset_index(drop=True)
+        df = self.pf_results_df.loc[guid].reset_index(drop=True)
         # round up and get the closest
         time = df.loc[(df["percentage_of_functionality"] >= pf) & (df["percentage_of_functionality"] <
                                                                    pf+self.pf_interval), state].values[0]
@@ -32,7 +33,7 @@ def get_percentage_func(self, guid, damage_state="DS_0", time=1):
             raise ValueError("restore time should not be larger than end time for restoration model!")
 
         state = damage_state.replace("DS", "PF")
-        df = self.time_results.loc[guid].reset_index(drop=True)
+        df = self.time_results_df.loc[guid].reset_index(drop=True)
         # round up and get the closest
         pf = df.loc[(df["time"] >= time) & df['time'] < time+self.time_interval, state].values[0]
 
diff --git a/pyincore/dataset.py b/pyincore/dataset.py
index 1062e7b83..b6eed8ab7 100644
--- a/pyincore/dataset.py
+++ b/pyincore/dataset.py
@@ -256,6 +256,25 @@ def get_csv_reader(self):
 
         return self.readers["csv"]
 
+    def get_csv_reader_std(self):
+        """Utility method for reading different standard file formats: csv reader.
+
+         Returns:
+             obj: CSV reader.
+
+         """
+        if "csv" not in self.readers:
+            filename = self.local_file_path
+            if os.path.isdir(filename):
+                files = glob.glob(filename + "/*.csv")
+                if len(files) > 0:
+                    filename = files[0]
+
+            csvfile = open(filename, 'r')
+            return csv.reader(csvfile)
+
+        return self.readers["csv"]
+
     def get_file_path(self, type='csv'):
         """Utility method for reading different standard file formats: file path.
 
@@ -274,7 +293,7 @@ def get_file_path(self, type='csv'):
 
         return filename
 
-    def get_dataframe_from_csv(self, low_memory=True):
+    def get_dataframe_from_csv(self, low_memory=True, delimiter=None):
         """Utility method for reading different standard file formats: Pandas DataFrame from csv.
 
         Args:
@@ -288,7 +307,7 @@ def get_dataframe_from_csv(self, low_memory=True):
         filename = self.get_file_path('csv')
         df = pd.DataFrame()
         if os.path.isfile(filename):
-            df = pd.read_csv(filename, header="infer", low_memory=low_memory)
+            df = pd.read_csv(filename, header="infer", low_memory=low_memory, delimiter=delimiter)
         return df
 
     def get_dataframe_from_shapefile(self):
diff --git a/pyincore/globals.py b/pyincore/globals.py
index 19669eb37..432b6041f 100644
--- a/pyincore/globals.py
+++ b/pyincore/globals.py
@@ -40,5 +40,7 @@
 MAX_LOGIN_ATTEMPTS = 3
 
 IPOPT_PATH = shutil.which("ipopt")
+GLPK_PATH = shutil.which("glpsol")
+SCIP_PATH = shutil.which("scip")
 
 DAMAGE_PRECISION = 10
diff --git a/recipes/meta.yaml b/recipes/meta.yaml
index 0bf9a7d23..485eb33f5 100644
--- a/recipes/meta.yaml
+++ b/recipes/meta.yaml
@@ -40,6 +40,7 @@ requirements:
   run:
     - python>=3.6
     - ipopt>=3.11
+    - scip>=8.0.0
     - {{ pin_compatible('numpy') }}
     - fiona>=1.8.4
     - geopandas>=0.6.1
diff --git a/requirements.min b/requirements.min
index 44eaa227c..e65bfed30 100644
--- a/requirements.min
+++ b/requirements.min
@@ -1,6 +1,7 @@
 fiona>=1.8.4
 geopandas>=0.6.1
 ipopt>=3.11
+scip>=8.0.0
 matplotlib>=2.1.0
 networkx>=2.2
 numpy>=1.16.6,<2.0a0
diff --git a/requirements.sh b/requirements.sh
index 9e48abb09..02ecc7dad 100755
--- a/requirements.sh
+++ b/requirements.sh
@@ -5,7 +5,7 @@ set -o noglob
 #TESTPYPI="--extra-index-url https://test.pypi.org/simple/ --pre"
 
 MODULE="pyincore"
-CONDA_IMPORTS="ipopt>=3.11"
+CONDA_IMPORTS="ipopt>=3.11 scip>=8.0.0"
 
 # setup pip-compile to resolve dependencies, also is minimal version of python
 if [ ! -d /tmp/pyincore-requirements ]; then


Test patch:
diff --git a/tests/pyincore/analyses/epfrestoration/test_epfrestoration.py b/tests/pyincore/analyses/epfrestoration/test_epfrestoration.py
index ff61d9cf2..c412bca1e 100644
--- a/tests/pyincore/analyses/epfrestoration/test_epfrestoration.py
+++ b/tests/pyincore/analyses/epfrestoration/test_epfrestoration.py
@@ -60,7 +60,14 @@ def run_with_base_class():
     epf_rest.run_analysis()
 
     # test utility function
-    epf_rest_util = EpfRestorationUtil(epf_rest)
+    inventory_restoration_map = epf_rest.get_output_dataset("inventory_restoration_map")
+    pf_results = epf_rest.get_output_dataset("pf_results")
+    time_results = epf_rest.get_output_dataset("time_results")
+    time_interval = epf_rest.get_parameter("time_interval")
+    pf_interval = epf_rest.get_parameter("pf_interval")
+    end_time = epf_rest.get_parameter("end_time")
+    epf_rest_util = EpfRestorationUtil(inventory_restoration_map, pf_results, time_results, time_interval,
+                                       pf_interval, end_time)
     functionality = epf_rest_util.get_percentage_func(guid="60748fbd-67c3-4f8d-beb9-26685a53d3c5",
                                                       damage_state="DS_0", time=2.0)
     time = epf_rest_util.get_restoration_time(guid="60748fbd-67c3-4f8d-beb9-26685a53d3c5", damage_state="DS_1", pf=0.81)
diff --git a/tests/pyincore/analyses/indp/test_indp.py b/tests/pyincore/analyses/indp/test_indp.py
new file mode 100644
index 000000000..f186b5f9c
--- /dev/null
+++ b/tests/pyincore/analyses/indp/test_indp.py
@@ -0,0 +1,335 @@
+from pyincore import IncoreClient, Dataset, RestorationService, MappingSet, FragilityService, NetworkDataset
+from pyincore.analyses.buildingdamage import BuildingDamage
+from pyincore.analyses.epfrepaircost import EpfRepairCost
+from pyincore.analyses.epfrestoration import EpfRestoration
+from pyincore.analyses.housingunitallocation import HousingUnitAllocation
+from pyincore.analyses.indp import INDP
+import pyincore.globals as pyglobals
+from pyincore.analyses.pipelinerepaircost import PipelineRepairCost
+from pyincore.analyses.populationdislocation import PopulationDislocation
+from pyincore.analyses.waterfacilitydamage import WaterFacilityDamage
+from pyincore.analyses.epfdamage import EpfDamage
+from pyincore.analyses.pipelinedamagerepairrate import PipelineDamageRepairRate
+from pyincore.analyses.montecarlofailureprobability import MonteCarloFailureProbability
+from pyincore.analyses.pipelinefunctionality import PipelineFunctionality
+from pyincore.analyses.pipelinerestoration import PipelineRestoration
+from pyincore.analyses.waterfacilityrepaircost import WaterFacilityRepairCost
+from pyincore.analyses.waterfacilityrestoration import WaterFacilityRestoration
+from pyincore.dataservice import DataService
+
+
+def run_with_base_class():
+    dev_client = IncoreClient(pyglobals.INCORE_API_DEV_URL)
+    prod_client = IncoreClient()
+    dev_datasvc = DataService(dev_client)
+
+    hazard_type = "earthquake"
+    hazard_id = "5dfa3e36b9219c934b64c231"  # 1000 yr eq
+    hazard_id_dev = "5ba8f127ec2309043520906c"  # 1000 yr eq dev
+    num_cpu = 8
+    sim_number = 2
+    sample_range = range(0, sim_number)
+    result_name = "seaside_indp"
+
+    bldg_inv_id = "613ba5ef5d3b1d6461e8c415"  # prod
+    bldg_inv_id_dev = "64c7d058a62b20774f4107b5"  # dev
+    seed = 1111
+
+    power_network_dataset = Dataset.from_data_service("64ac73694e01de3af8fd8f2b", data_service=dev_datasvc)
+    power_network = NetworkDataset.from_dataset(power_network_dataset)
+    water_network_dataset = Dataset.from_data_service("64ad6abb4e01de3af8fe5201", data_service=dev_datasvc)
+    water_network = NetworkDataset.from_dataset(water_network_dataset)
+    water_facilities = water_network.nodes
+    epfs = power_network.nodes
+    pipeline = water_network.links
+
+    ###################################################
+    # water facility damage
+    ###################################################
+    wterfclty_dmg = WaterFacilityDamage(prod_client)
+    fragility_service = FragilityService(prod_client)
+    wterfclty_dmg.set_input_dataset("water_facilities", water_facilities)
+    mapping_id = "5d39e010b9219cc18bd0b0b6"  # 5 DS
+    mapping_set = MappingSet(fragility_service.get_mapping(mapping_id))
+    wterfclty_dmg.set_input_dataset("dfr3_mapping_set", mapping_set)
+    wterfclty_dmg.set_parameter("hazard_type", hazard_type)
+    wterfclty_dmg.set_parameter("fragility_key", "pga")
+    wterfclty_dmg.set_parameter("num_cpu", 4)
+    wterfclty_dmg.set_parameter("hazard_id", hazard_id)
+    wterfclty_dmg.set_parameter("result_name", result_name + "_wf_damage")
+    wterfclty_dmg.run_analysis()
+    wterfclty_dmg_result = wterfclty_dmg.get_output_dataset("result")
+
+    ###################################################
+    # water facility mcs
+    ###################################################
+    wterfclty_mc = MonteCarloFailureProbability(prod_client)
+    wterfclty_mc.set_input_dataset("damage", wterfclty_dmg_result)
+    wterfclty_mc.set_parameter("num_cpu", num_cpu)
+    wterfclty_mc.set_parameter("num_samples", sim_number)
+    wterfclty_mc.set_parameter("damage_interval_keys", ["DS_0", "DS_1", "DS_2", "DS_3", "DS_4"])
+    wterfclty_mc.set_parameter("failure_state_keys", ["DS_1", "DS_2", "DS_3", "DS_4"])
+    wterfclty_mc.set_parameter("result_name", result_name + "_wf")  # name of csv file with results
+    wterfclty_mc.run()
+    wterfclty_sample_failure_state = wterfclty_mc.get_output_dataset("sample_failure_state")
+    wterfclty_sample_damage_states = wterfclty_mc.get_output_dataset("sample_damage_states")
+
+    ###################################################
+    # water facility repair time
+    ###################################################
+    wterfclty_rest = WaterFacilityRestoration(prod_client)
+    restorationsvc = RestorationService(prod_client)
+    mapping_set = MappingSet(restorationsvc.get_mapping("61f075ee903e515036cee0a5"))
+    wterfclty_rest.set_input_dataset("water_facilities", water_facilities)
+    wterfclty_rest.set_input_dataset("dfr3_mapping_set", mapping_set)
+    wterfclty_rest.set_input_dataset('damage', wterfclty_dmg_result)
+    wterfclty_rest.set_parameter("result_name", result_name + "_wf_restoration")
+    wterfclty_rest.set_parameter("discretized_days", [1, 3, 7, 30, 90])
+    wterfclty_rest.set_parameter("restoration_key", "Restoration ID Code")
+    wterfclty_rest.set_parameter("end_time", 100.0)
+    wterfclty_rest.set_parameter("time_interval", 1.0)
+    wterfclty_rest.set_parameter("pf_interval", 0.01)
+    wterfclty_rest.run_analysis()
+    wf_restoration_time = wterfclty_rest.get_output_dataset("repair_times")
+
+    ###################################################
+    # water facility repair cost
+    ###################################################
+    wf_repair_cost = WaterFacilityRepairCost(prod_client)
+    wf_repair_cost.set_input_dataset("water_facilities", water_facilities)
+    wf_repair_cost.load_remote_input_dataset("replacement_cost", "64833bcdd3f39a26a0c8b147")
+    wf_repair_cost.set_input_dataset("sample_damage_states", wterfclty_sample_damage_states)
+    wf_repair_cost.load_remote_input_dataset("wf_dmg_ratios", "647e423d7ae18139d9758607")
+    wf_repair_cost.set_parameter("result_name", result_name + "_wf_repair_cost")
+    wf_repair_cost.set_parameter("num_cpu", 4)
+    wf_repair_cost.run_analysis()
+    wf_repair_cost_result = wf_repair_cost.get_output_dataset("result")
+
+    ###################################################
+    # epf damage
+    ###################################################
+    epf_dmg = EpfDamage(prod_client)
+    fragility_service = FragilityService(prod_client)
+    epf_dmg.set_input_dataset("epfs", epfs)
+    mapping_id = "64ac5f3ad2122d1f95f36356"  # 5 DS
+    mapping_set = MappingSet(fragility_service.get_mapping(mapping_id))
+    epf_dmg.set_input_dataset('dfr3_mapping_set', mapping_set)
+    epf_dmg.set_parameter("hazard_type", hazard_type)
+    epf_dmg.set_parameter("num_cpu", num_cpu)
+    epf_dmg.set_parameter('fragility_key', "pga")
+    epf_dmg.set_parameter("hazard_id", hazard_id)
+    epf_dmg.set_parameter("result_name", result_name + "_epf_dmg")
+    epf_dmg.run_analysis()
+    epf_dmg_result = epf_dmg.get_output_dataset("result")
+
+    ###################################################
+    # epf repair time
+    ###################################################
+    epf_rest = EpfRestoration(prod_client)
+    restorationsvc = RestorationService(prod_client)
+    mapping_set = MappingSet(restorationsvc.get_mapping("61f302e6e3a03e465500b3eb"))
+    epf_rest.set_input_dataset("epfs", epfs)
+    epf_rest.set_input_dataset("dfr3_mapping_set", mapping_set)
+    epf_rest.set_input_dataset('damage', epf_dmg_result)
+    epf_rest.set_parameter("result_name", result_name + "_epf_restoration")
+    epf_rest.set_parameter("discretized_days", [1, 3, 7, 30, 90])
+    epf_rest.set_parameter("restoration_key", "Restoration ID Code")
+    epf_rest.set_parameter("end_time", 100.0)
+    epf_rest.set_parameter("time_interval", 1.0)
+    epf_rest.set_parameter("pf_interval", 0.01)
+    epf_rest.run_analysis()
+    epf_restoration_time = epf_rest.get_output_dataset("repair_times")
+
+    ###################################################
+    # epf mcs
+    ###################################################
+    epf_mc = MonteCarloFailureProbability(prod_client)
+    epf_mc.set_input_dataset("damage", epf_dmg_result)
+    epf_mc.set_parameter("num_cpu", num_cpu)
+    epf_mc.set_parameter("num_samples", sim_number)
+    epf_mc.set_parameter("damage_interval_keys", ["DS_0", "DS_1", "DS_2", "DS_3", "DS_4"])
+    epf_mc.set_parameter("failure_state_keys", ["DS_1", "DS_2", "DS_3", "DS_4"])
+    epf_mc.set_parameter("result_name", result_name + "_epf")  # name of csv file with results
+    epf_mc.run()
+    epf_sample_failure_state = epf_mc.get_output_dataset("sample_failure_state")
+    epf_sample_damage_states = epf_mc.get_output_dataset("sample_damage_states")
+
+    ###################################################
+    # epf repair cost
+    ###################################################
+    epf_repair_cost = EpfRepairCost(prod_client)
+    epf_repair_cost.set_input_dataset("epfs", epfs)
+    epf_repair_cost.load_remote_input_dataset("replacement_cost", "647dff5b4dd25160127ca192")
+    epf_repair_cost.set_input_dataset("sample_damage_states", epf_sample_damage_states)
+    epf_repair_cost.load_remote_input_dataset("epf_dmg_ratios", "6483354b41181d20004efbd7")
+    epf_repair_cost.set_parameter("result_name", result_name + "_epf_repair_cost")
+    epf_repair_cost.set_parameter("num_cpu", 4)
+    epf_repair_cost.run_analysis()
+    epf_repair_cost_result = epf_repair_cost.get_output_dataset("result")
+
+    ###################################################
+    # pipeline repair rate damage
+    ###################################################
+    pipeline_dmg = PipelineDamageRepairRate(prod_client)
+    fragility_service = FragilityService(prod_client)
+    pipeline_dmg.set_input_dataset("pipeline", pipeline)
+    mapping_id = "5b47c227337d4a38464efea8"
+    mapping_set = MappingSet(fragility_service.get_mapping(mapping_id))
+    pipeline_dmg.set_input_dataset('dfr3_mapping_set', mapping_set)
+    pipeline_dmg.set_parameter("hazard_type", hazard_type)
+    pipeline_dmg.set_parameter("fragility_key", 'pgv')
+    pipeline_dmg.set_parameter("num_cpu", num_cpu)
+    pipeline_dmg.set_parameter("hazard_id", hazard_id)
+    pipeline_dmg.set_parameter("result_name", result_name + "_pipeline_dmg")
+    pipeline_dmg.run_analysis()
+    pipeline_dmg_result = pipeline_dmg.get_output_dataset("result")
+
+    ###################################################
+    # pipeline functionality
+    ###################################################
+    pipeline_func = PipelineFunctionality(prod_client)
+    pipeline_func.set_input_dataset("pipeline_repair_rate_damage", pipeline_dmg_result)
+    pipeline_func.set_parameter("result_name", result_name + "_pipeline")
+    pipeline_func.set_parameter("num_samples", sim_number)
+    pipeline_func.run_analysis()
+    pipeline_sample_failure_state = pipeline_func.get_output_dataset("sample_failure_state")
+
+    ###################################################
+    # pipeline repair time
+    ###################################################
+    pipeline_rest = PipelineRestoration(prod_client)
+    restorationsvc = RestorationService(prod_client)
+    mapping_set = MappingSet(restorationsvc.get_mapping("61f35f09903e515036cee106"))
+    pipeline_rest.set_input_dataset("pipeline", pipeline)
+    pipeline_rest.set_input_dataset("pipeline_damage", pipeline_dmg_result)
+    pipeline_rest.set_input_dataset("dfr3_mapping_set", mapping_set)
+    pipeline_rest.set_parameter("num_available_workers", 4)
+    pipeline_rest.set_parameter("result_name", result_name + "_pipeline_restoration")
+    pipeline_rest.run_analysis()
+    pipeline_restoration_time = pipeline_rest.get_output_dataset("pipeline_restoration")
+
+    ###################################################
+    # pipeline repair cost
+    ###################################################
+    pipeline_repair_cost = PipelineRepairCost(prod_client)
+    pipeline_repair_cost.set_input_dataset("pipeline", pipeline)
+    pipeline_repair_cost.load_remote_input_dataset("replacement_cost", "6480a2787ae18139d975e919")
+    pipeline_repair_cost.set_input_dataset("pipeline_dmg", pipeline_dmg_result)
+    pipeline_repair_cost.load_remote_input_dataset("pipeline_dmg_ratios", "6480a2d44dd25160127d2fcc")
+    pipeline_repair_cost.set_parameter("result_name", result_name + "_pipeline_repair_cost")
+    pipeline_repair_cost.set_parameter("num_cpu", 4)
+    pipeline_repair_cost.run_analysis()
+    pipeline_repair_cost_result = pipeline_repair_cost.get_output_dataset("result")
+
+    ###################################################
+    # building damage
+    ###################################################
+    bldg_dmg = BuildingDamage(dev_client)
+    fragility_service = FragilityService(dev_client)
+    bldg_dmg.load_remote_input_dataset("buildings", bldg_inv_id_dev)
+    mapping_id = "5e99c86d6129af000136defa"  # 4 DS dev
+    # mapping_id = "5d2789dbb9219c3c553c7977"  # 4 DS prod
+    mapping_set = MappingSet(fragility_service.get_mapping(mapping_id))
+    bldg_dmg.set_input_dataset('dfr3_mapping_set', mapping_set)
+    bldg_dmg.set_parameter("hazard_type", hazard_type)
+    bldg_dmg.set_parameter("num_cpu", 4)
+    bldg_dmg.set_parameter("hazard_id", hazard_id_dev)
+    bldg_dmg.set_parameter("result_name", result_name + "_bldg_dmg")
+    bldg_dmg.run_analysis()
+    building_dmg_result = bldg_dmg.get_output_dataset("ds_result")
+
+    ###################################################
+    # housing unit allocation
+    ###################################################
+    hua = HousingUnitAllocation(prod_client)
+    housing_unit_inv_id = "5d543087b9219c0689b98234"
+    address_point_inv_id = "5d542fefb9219c0689b981fb"
+    hua.load_remote_input_dataset("housing_unit_inventory", housing_unit_inv_id)
+    hua.load_remote_input_dataset("address_point_inventory", address_point_inv_id)
+    hua.load_remote_input_dataset("buildings", bldg_inv_id)
+    hua.set_parameter("result_name", result_name + "_hua")
+    hua.set_parameter("seed", seed)
+    hua.set_parameter("iterations", 1)
+    hua.run_analysis()
+    hua_result = hua.get_output_dataset("result")
+
+    ###################################################
+    # population dislocation
+    ###################################################
+    pop_dis = PopulationDislocation(prod_client)
+    pop_dis.set_input_dataset("building_dmg", building_dmg_result)
+    pop_dis.set_input_dataset("housing_unit_allocation", hua_result)
+    pop_dis.load_remote_input_dataset("block_group_data", "5d542bd8b9219c0689b90408")
+    pop_dis.load_remote_input_dataset("value_loss_param", "60354810e379f22e16560dbd")
+    pop_dis.set_parameter("result_name", result_name + "_popdislocation")
+    pop_dis.set_parameter("seed", seed)
+    pop_dis.run_analysis()
+    pop_dislocation_result = pop_dis.get_output_dataset("result")
+
+    ###################################################
+    # INDP
+    ###################################################
+    indp_analysis = INDP(dev_client)
+    indp_analysis.set_parameter("network_type", "from_csv")
+    indp_analysis.set_parameter("MAGS", [1000])
+    indp_analysis.set_parameter("sample_range", sample_range)
+    indp_analysis.set_parameter("dislocation_data_type", "incore")
+    indp_analysis.set_parameter("return_model", "step_function")
+    indp_analysis.set_parameter("testbed_name", "seaside")
+    indp_analysis.set_parameter("extra_commodity", {1: ["PW"], 3: []})
+    indp_analysis.set_parameter("RC", [{"budget": 240000, "time": 700}, {"budget": 300000, "time": 600}])
+    indp_analysis.set_parameter("layers", [1, 3])
+    indp_analysis.set_parameter("method", "INDP")
+    # indp_analysis.set_parameter("method", "TDINDP")
+    indp_analysis.set_parameter("t_steps", 10)
+    indp_analysis.set_parameter("time_resource", True)
+    # indp_analysis.set_parameter("save_model", False)
+    indp_analysis.set_parameter("save_model", True)
+
+    # scip
+    # indp_analysis.set_parameter("solver_engine", "scip") # recommended
+    # indp_analysis.set_parameter("solver_path", "/usr/local/bin/scip")
+
+    # cbc
+    # indp_analysis.set_parameter("solver_engine", "cbc")
+    # indp_analysis.set_parameter("solver_path", "/usr/local/bin/cbc")
+
+    # glpk
+    # indp_analysis.set_parameter("solver_engine", "glpk")
+    # indp_analysis.set_parameter("solver_path", "/usr/local/bin/glpsol")
+
+    # gurobi
+    # indp_analysis.set_parameter("solver_engine", "gurobi")
+
+    indp_analysis.set_parameter("solver_time_limit", 3600)  # if not set default to never timeout
+
+    indp_analysis.set_input_dataset("wf_restoration_time", wf_restoration_time)
+    indp_analysis.set_input_dataset("wf_repair_cost", wf_repair_cost_result)
+    indp_analysis.set_input_dataset("epf_restoration_time", epf_restoration_time)
+    indp_analysis.set_input_dataset("epf_repair_cost", epf_repair_cost_result)
+    indp_analysis.set_input_dataset("pipeline_restoration_time", pipeline_restoration_time)
+    indp_analysis.set_input_dataset("pipeline_repair_cost", pipeline_repair_cost_result)
+    indp_analysis.set_input_dataset("power_network", power_network_dataset)
+    indp_analysis.set_input_dataset("water_network", water_network_dataset)  # with distribution noes
+    indp_analysis.load_remote_input_dataset("powerline_supply_demand_info", "64ad8b434e01de3af8fea0ba")
+    indp_analysis.load_remote_input_dataset("epf_supply_demand_info", "64ad9ea54e01de3af8fea0f2")
+    indp_analysis.load_remote_input_dataset("wf_supply_demand_info", "64ad9e704e01de3af8fea0ec")
+    indp_analysis.load_remote_input_dataset("pipeline_supply_demand_info", "64ad9e274e01de3af8fea0e5")
+    indp_analysis.load_remote_input_dataset("interdep", "61c10104837ac508f9a178ef")
+    indp_analysis.set_input_dataset("wf_failure_state", wterfclty_sample_failure_state)
+    indp_analysis.set_input_dataset("wf_damage_state", wterfclty_sample_damage_states)
+    indp_analysis.set_input_dataset("pipeline_failure_state", pipeline_sample_failure_state)
+    indp_analysis.set_input_dataset("epf_failure_state", epf_sample_failure_state)
+    indp_analysis.set_input_dataset("epf_damage_state", epf_sample_damage_states)
+    indp_analysis.set_input_dataset("pop_dislocation", pop_dislocation_result)
+
+    # # optional inputs
+    # indp_analysis.load_remote_input_dataset("bldgs2elec", "61c10219837ac508f9a17904")
+    # indp_analysis.load_remote_input_dataset("bldgs2wter", "c102b0837ac508f9a1790a")
+
+    # Run Analysis
+    indp_analysis.run_analysis()
+
+
+if __name__ == "__main__":
+    run_with_base_class()
diff --git a/tests/pyincore/analyses/waterfacilityrestoration/test_waterfacilityrestoration.py b/tests/pyincore/analyses/waterfacilityrestoration/test_waterfacilityrestoration.py
index 3613fa528..6645dd46c 100644
--- a/tests/pyincore/analyses/waterfacilityrestoration/test_waterfacilityrestoration.py
+++ b/tests/pyincore/analyses/waterfacilityrestoration/test_waterfacilityrestoration.py
@@ -64,10 +64,17 @@ def run_with_base_class():
     wf_rest.run_analysis()
 
     # test utility function
-    wf_util = WaterFacilityRestorationUtil(wf_rest)
+    inventory_restoration_map = wf_rest.get_output_dataset("inventory_restoration_map")
+    pf_results = wf_rest.get_output_dataset("pf_results")
+    time_results = wf_rest.get_output_dataset("time_results")
+    time_interval = wf_rest.get_parameter("time_interval")
+    pf_interval = wf_rest.get_parameter("pf_interval")
+    end_time = wf_rest.get_parameter("end_time")
+    wf_util = WaterFacilityRestorationUtil(inventory_restoration_map, pf_results, time_results, time_interval,
+                                           pf_interval, end_time)
     functionality = wf_util.get_percentage_func(guid="e1bce78d-00a1-4605-95f3-3776ff907f73",
                                                 damage_state="DS_0", time=2.0)
-    time = wf_util.get_restoration_time(guid="e1bce78d-00a1-4605-95f3-3776ff907f73", damage_state="DS_1", pf=0.81)
+    time = wf_util.get_restoration_time(guid="e1bce78d-00a1-4605-95f3-3776ff907f73", damage_state="DS_1", pf=0.95)
     print(functionality, time)