In the previous post, we defined robustness using the satisficing metric where (1) reliability should be at least 98%, (2) restriction frequency should be not more than 10% and (3) worst-case cost of drought mitigation action should not be more than 10% of annual net volumetric income. To calculate the robustness of these set of actions (portfolios) against the effects of challenging states of the world (SOWs) on the initial set of actions, we once again re-simulated them to discover how they fail.
In this penultimate post, we will perform simple sensitivity analysis across the average performance of all sixty-nine portfolios of actions to understand which uncertainties control the performance of each utility (Raleigh, Durham and Cary) and the regions across all uncertain SOWs.
Calculating average performance across 100 DU SOWs
First, we create a new folder to hold the output of the next few post-processing steps. Navigate to the WaterPaths/ folder and create a folder called post_processing. Now, let’s calculate the average performance of each of the sixty-nine initial portfolios across the 100 DU SOWs that we previously simulated them over. The code for this can be found in the post_processing_code folder under gather_objs.py file and should look like this:
# -*- coding: utf-8 -*-
"""
Created on Mon April 26 2022 11:12
@author: Lillian Bei Jia Lau
Organizes output objectives by mean across RDMs and across solutions
"""
import numpy as np
obj_names = ['REL_C', 'RF_C', 'INF_NPC_C', 'PFC_C', 'WCC_C', \
'REL_D', 'RF_D', 'INF_NPC_D', 'PFC_D', 'WCC_D', \
'REL_R', 'RF_R', 'INF_NPC_R', 'PFC_R', 'WCC_R', \
'REL_reg', 'RF_reg', 'INF_NPC_reg', 'PFC_reg', 'WCC_reg']
'''
Performs regional minimax
'''
def minimax(N_SOLNS, objs):
for i in range(N_SOLNS):
for j in range(5):
if j == 0:
objs[i,15] = np.min([objs[i,0],objs[i,5], objs[i,10]])
else:
objs[i, (j+15)] = np.max([objs[i,j],objs[i,j+5], objs[i,j+10]])
return objs
'''
Calculates the mean performance acorss all SOWs
'''
def mean_performance_across_rdms(objs_by_rdm_dir, N_RDMS, N_SOLNS):
objs_matrix = np.zeros((N_SOLNS,20,N_RDMS), dtype='float')
objs_means = np.zeros((N_SOLNS,20), dtype='float')
for i in range(N_RDMS):
filepathname = objs_by_rdm_dir + str(i) + '_sols0_to_' + str(N_SOLNS) + '.csv'
objs_file = np.loadtxt(filepathname, delimiter=",")
objs_matrix[:,:15,i] = objs_file
objs_file_wRegional = minimax(N_SOLNS, objs_matrix[:,:,i])
objs_matrix[:,:,i] = objs_file_wRegional
array_has_nan = np.isnan(np.mean(objs_matrix[:,3,i]))
if(array_has_nan == True):
print('NaN found at RDM ', str(i))
for n in range(N_SOLNS):
for n_objs in range(20):
objs_means[n,n_objs] = np.mean(objs_matrix[n,n_objs,:])
return objs_means
'''
Calculates the mean performance acorss all SOWs
'''
def mean_performance_across_solns(objs_by_rdm_dir, N_RDMS, N_SOLNS):
objs_matrix = np.zeros((N_SOLNS,20,N_RDMS), dtype='float')
objs_means = np.zeros((N_RDMS,20), dtype='float')
for i in range(N_RDMS):
filepathname = objs_by_rdm_dir + str(i) + '_sols0_to_' + str(N_SOLNS) + '.csv'
objs_file = np.loadtxt(filepathname, delimiter=",")
objs_matrix[:,:15,i] = objs_file
objs_file_wRegional = minimax(N_SOLNS, objs_matrix[:,:,i])
objs_matrix[:,:,i] = objs_file_wRegional
array_has_nan = np.isnan(np.mean(objs_matrix[:,3,i]))
if(array_has_nan == True):
print('NaN found at RDM ', str(i))
for n in range(N_RDMS):
for n_objs in range(20):
objs_means[n,n_objs] = np.mean(objs_matrix[:,n_objs,n])
return objs_means
# change number of solutions available depending on the number of solutions
# that you identified
N_SOLNS = 69
N_RDMS = 100
# change the filepaths
objs_by_rdm_dir = '/yourFilePath/WaterPaths/output/Objectives_RDM'
objs_og_dir = '/yourFilePath/WaterPaths/'
fileoutpath = '/yourFilePath/WaterPaths/post_processing/'
fileoutpath_byRDMs = fileoutpath + 'meanObjs_acrossRDM.csv'
fileoutpath_bySoln = fileoutpath + 'meanObjs_acrossSoln.csv'
# should have shape (N_SOLNS, 20)
objs_byRDM = mean_performance_across_rdms(objs_by_rdm_dir, N_RDMS, N_SOLNS)
# should have shape (N_RDMS, 20)
objs_bySoln = mean_performance_across_solns(objs_by_rdm_dir, N_RDMS, N_SOLNS)
np.savetxt(fileoutpath_byRDMs, objs_byRDM, delimiter=",")
np.savetxt(fileoutpath_bySoln, objs_bySoln, delimiter=",")
This will output two .csv files: meanObjs_acrossRDM.csv contains the average performance of each of the sixty-nine objectives evaluated over 100 DU SOWs, while meanObjs_acrossSoln.csv contains the average performance of all solutions within one SOW. Take some time to understand this difference, as this will be important when performing sensitivity analysis and scenario discovery.
Calculate the robustness of each portfolio to deep uncertainty
Now, let’s identify how each of these solutions perform under a set of more challenging SOWs. Within post_processing_code/, identify the file called calc_robustness_across_rdms.py. It should look like this:
# -*- coding: utf-8 -*-
"""
Created on Mon April 26 2022 11:12
@author: Lillian Bei Jia Lau
Calculates the fraction of RDMs over which each perturbed version of the solution meets all four satisficing criteria
"""
import numpy as np
import pandas as pd
obj_names = ['REL_C', 'RF_C', 'INF_NPC_C', 'PFC_C', 'WCC_C', \
'REL_D', 'RF_D', 'INF_NPC_D', 'PFC_D', 'WCC_D', \
'REL_R', 'RF_R', 'INF_NPC_R', 'PFC_R', 'WCC_R', \
'REL_reg', 'RF_reg', 'INF_NPC_reg', 'PFC_reg', 'WCC_reg']
utilities = ['Cary', 'Durham', 'Raleigh', 'Regional']
'''
Performs regional minimax
'''
def minimax(N_SOLNS, objs):
for i in range(N_SOLNS):
for j in range(5):
if j == 0:
objs[i,15] = np.min([objs[i,0],objs[i,5], objs[i,10]])
else:
objs[i, (j+15)] = np.max([objs[i,j],objs[i,j+5], objs[i,j+10]])
return objs
'''
For each rdm, identify if the perturbed solution version x satisfies the satisficing criteria
'''
def satisficing(df_objs):
for i in range(4):
df_objs['satisficing_C'] = (df_objs['REL_C'] >= 0.98).astype(int) *\
(df_objs['WCC_C'] <= 0.10).astype(int) *\
(df_objs['RF_C'] <= 0.10).astype(int)
df_objs['satisficing_D'] = (df_objs['REL_D'] >= 0.98).astype(int) *\
(df_objs['WCC_D'] <= 0.10).astype(int) *\
(df_objs['RF_D'] <= 0.10).astype(int)
df_objs['satisficing_R'] = (df_objs['REL_R'] >= 0.98).astype(int) *\
(df_objs['WCC_R'] <= 0.10).astype(int) *\
(df_objs['RF_R'] <= 0.10).astype(int)
df_objs['satisficing_reg'] = np.max(df_objs.iloc[:, 20:23])
return df_objs
def calc_robustness(objs_by_rdm_dir, N_RDMS, N_SOLNS):
# matrix structure: (N_SOLNS, N_OBJS, N_RDMS)
objs_matrix = np.zeros((N_SOLNS,20,N_RDMS), dtype='float')
satisficing_matrix = np.zeros((N_SOLNS,4,N_RDMS), dtype='float')
solution_robustness = np.zeros((N_SOLNS,4), dtype='float')
for i in range(N_RDMS):
# get one perturbed instance's behavior over all RDMs
filepathname = objs_by_rdm_dir + str(i) + '_sols0_to_' + str(N_SOLNS) + '.csv'
objs_file = np.loadtxt(filepathname, delimiter=",")
objs_matrix[:,:15,i] = objs_file
objs_file_wRegional = minimax(N_SOLNS, objs_matrix[:,:,i])
objs_matrix[:,:,i] = objs_file_wRegional
# NaN check
array_has_nan = np.isnan(np.mean(objs_matrix[:,3,i]))
if(array_has_nan == True):
print('NaN found at RDM ', str(i))
# for the perturbed instances
for r in range(N_RDMS):
df_curr_rdm = pd.DataFrame(objs_matrix[:, :, r], columns = obj_names)
df_curr_rdm = satisficing(df_curr_rdm)
satisficing_matrix[:N_SOLNS,:,r] = df_curr_rdm.iloc[:,20:24]
for n in range(N_SOLNS):
solution_robustness[n,0] = np.sum(satisficing_matrix[n,0,:])/N_RDMS
solution_robustness[n,1] = np.sum(satisficing_matrix[n,1,:])/N_RDMS
solution_robustness[n,2] = np.sum(satisficing_matrix[n,2,:])/N_RDMS
solution_robustness[:,3] = np.min(solution_robustness[:,:3], axis=1)
return solution_robustness
'''
Change number of solutions available depending on the number of solutions
that you identified and the number of SOWs that you are evaluating them over.
'''
N_RDMS = 100
N_SOLNS = 69
objs_by_rdm_dir = '/scratch/lbl59/blog/WaterPaths/output/Objectives_RDM'
fileoutpath_robustness = '/scratch/lbl59/blog/WaterPaths/post_processing/' + \
'robustness_' + str(N_RDMS) + '_SOWs.csv'
robustness = calc_robustness(objs_by_rdm_dir, N_RDMS, N_SOLNS)
np.savetxt(fileoutpath_robustness, robustness, delimiter=",")
When you run this script from the terminal, you should have a .csv file called ‘robustness_100_SOWs.csv‘ appear in your post_processing/ folder. Now, let’s get onto some sensitivity analysis!
Delta moment-independent sensitivity analysis
The Delta moment-independent (DMI) method (Borgonovo, 2007) is sensitivity analysis method that compares the entire probability distribution of both input and output parameters to estimate the sensitivity of the output to a specific input parameter. It is one of many global sensitivity analysis methods, which in itself is one of two main categories of sensitivity analysis that enables the assessment of the degree to which uncertainty in model inputs map to the degree of uncertainty in model output. For purposes of this test case, the DMI is preferable as it does not rely on any one statistical moment (variance, mean and skew) to describe the dependence of model output to its input parameters. It is also time-sensitive, reflecting the current state of knowledge within the system, which philosophically pairs well with our use of the ROF triggers. More information on alternative global sensitivity methods can be found here.
Within the context of our test case, we will be using the DMI method to identify uncertainties in our decision variables that most strongly influence our performance over the 100 DU SOWs. To perform DMI sensitivity analysis, first navigate to the post_processing/ folder. Within the folder, create two sub-folders called delta_output_DV/ and delta_output_DUF/. This is where all your DMI output will be stored. Next, locate the delta_sensitivity.py file within the post_processing_code/ folder. The code should look similar to the script provided below:
import sys
from SALib.analyze import delta
from SALib.util import read_param_file
import numpy as np
import pandas as pd
'''
Finds the upper and lower bounds of input parameters
'''
def find_bounds(input_file):
bounds = np.zeros((input_file.shape[1],2), dtype=float)
for i in range(input_file.shape[1]):
bounds[i,0] = min(input_file[:,i])
bounds[i,1] = max(input_file[:,i])
return bounds
'''
Performs delta moment-independent sensitivity analysis
Source: https://github.com/SALib/SALib/tree/main/examples/delta
'''
def delta_sensitivity(dec_vars, measured_outcomes, names, mo_names, bounds, rdm, mode):
X = dec_vars
Y = measured_outcomes
problem = {
'num_vars': int(dec_vars.shape[1]),
'names': names,
'bounds': bounds
}
for i in range(measured_outcomes.shape[1]):
mo_label = mo_names[i]
if i == 2 and mode == 'objs':
break
else:
filename = '../post_processing/delta_output_' + rdm + '/S1_' + mo_label + '.csv'
S1 = delta.analyze(problem, X, Y[mo_label].values, num_resamples=10, conf_level=0.95, print_to_console=False)
numpy_S1 = np.array(S1["S1"])
fileout = pd.DataFrame([names, numpy_S1], index = None, columns = None)
fileout.to_csv(filename, sep=",")
'''
0 - Name all file headers and compSol to be analyzed
'''
obj_names = ['REL_C', 'RF_C', 'INF_NPC_C', 'PFC_C', 'WCC_C', \
'REL_D', 'RF_D', 'INF_NPC_D', 'PFC_D', 'WCC_D', \
'REL_R', 'RF_R', 'INF_NPC_R', 'PFC_R', 'WCC_R', \
'REL_reg', 'RF_reg', 'INF_NPC_reg', 'PFC_reg', 'WCC_reg']
dv_names = ['RT_C', 'RT_D', 'RT_R', 'TT_D', 'TT_R', 'LMA_C', 'LMA_D', 'LMA_R',\
'RC_C', 'RC_D', 'RC_R', 'IT_C', 'IT_D', 'IT_R', 'IP_C', 'IP_D', \
'IP_R', 'INF_C', 'INF_D', 'INF_R']
rdm_headers_dmp = ['Cary restr. eff', 'Durham restr. eff', 'Raleigh restr. eff']
rdm_headers_utilities = ['Demand growth\nmultiplier', 'Bond term\nmultiplier', \
'Bond interest\nrate multiplier', 'Infrastructure interest\nrate multiplier']
rdm_headers_ws = ['Streamflow amp', 'SCR PT', 'SCR CT', 'NRR PT', 'NRR CT', 'CR Low PT', 'CR Low CT',\
'CR High PT', 'CR High CT', 'WR1 PT', 'WR1 CT', 'WR2 PT', 'WR2 CT',\
'DR PT', 'DR CT', 'FR PT', 'FR CT']
duf_names = ['Cary restr. eff', 'Durham restr. eff', 'Raleigh restr. eff', 'Demand growth\nmultiplier',\
'Bond term\nmultiplier', 'Bond interest\nrate multiplier', 'Infrastructure interest\nrate multiplier',\
'Streamflow amp\nmultiplier', 'SCR PT\nmultiplier', 'SCR CT\nmultiplier', 'NRR PT\nmultiplier',\
'NRR CT\nmultiplier', 'CR Low PT\nmultiplier', 'CR Low CT\nmultiplier', 'CR High PT\nmultiplier',\
'CR High CT\nmultiplier', 'WR1 PT\nmultiplier', 'WR1 CT\nmultiplier', 'WR2 PT\nmultiplier',\
'WR2 CT\nmultiplier', 'DR PT\nmultiplier', 'DR CT\nmultiplier', 'FR PT\nmultiplier', 'FR CT\nmultiplier',\
'DR PT\nmultiplier', 'DR CT\nmultiplier', 'FR PT\nmultiplier', 'FR CT\nmultiplier']
utilities = ['Cary', 'Durham', 'Raleigh', 'Regional']
N_RDMS = 100
N_SOLNS = 69
'''
1 - Load DU factor files and DV files
'''
# change to your own filepath
rdm_factors_directory = '/yourFilePath/WaterPaths/TestFiles/'
rdm_dmp_filename = rdm_factors_directory + 'rdm_dmp_test_problem_reeval.csv'
rdm_utilities_filename = rdm_factors_directory + 'rdm_utilities_test_problem_reeval.csv'
rdm_watersources_filename = rdm_factors_directory + 'rdm_water_sources_test_problem_reeval.csv'
rdm_dmp = pd.read_csv(rdm_dmp_filename, sep=",", names=rdm_headers_dmp)
rdm_utilities = pd.read_csv(rdm_utilities_filename, sep=",", names=rdm_headers_utilities)
rdm_ws_all = np.loadtxt(rdm_watersources_filename, delimiter=",")
rdm_ws = pd.DataFrame(rdm_ws_all[:,:17], columns=rdm_headers_ws)
dufs = pd.concat([rdm_dmp, rdm_utilities, rdm_ws], axis=1, ignore_index=True)
duf_numpy = dufs.to_numpy()
# change to your own filepath
dv_directory = '/yourFilePath/WaterPaths/'
dv_filename = dv_directory + 'NC_dvs_all_noheader.csv'
dvs = np.loadtxt(dv_filename, delimiter=",")
'''
2 - Get bounds for DU factors and DVs
'''
duf_bounds = find_bounds(duf_numpy)
dv_bounds = find_bounds(dvs)
'''
3 - Load robustness file and objectives file
'''
# change to your own filepath
main_dir = '/yourFilePath/WaterPaths/post_processing/'
robustness_filename = main_dir + 'robustness_100_SOWs.csv'
robustness_arr = np.loadtxt(robustness_filename, delimiter=",")
robustness_df = pd.DataFrame(robustness_arr, columns=utilities)
objs_mean_rdm_filename = main_dir + 'meanObjs_acrossRDM.csv'
objs_mean_rdm_arr = np.loadtxt(objs_mean_rdm_filename, delimiter=",")
objs_mean_rdm_df = pd.DataFrame(objs_mean_rdm_arr, columns=obj_names)
objs_mean_soln_filename = main_dir + 'meanObjs_acrossSoln.csv'
objs_mean_soln_arr = np.loadtxt(objs_mean_soln_filename, delimiter=",")
objs_mean_soln_df = pd.DataFrame(objs_mean_soln_arr, columns=obj_names)
# to change depending on whether DV or DUF is being analyzed
dec_vars = dvs
measured_outcomes = objs_mean_rdm_df
names = dv_names
mo_names = obj_names
bounds = dv_bounds
rdm = 'DV'
mode = 'objs'
###
delta_sensitivity(dec_vars, measured_outcomes, names, mo_names, bounds, rdm, mode)
The code above identifies the sensitivity of the average values of all sixty-nine performance objective sets over all 100 deeply-uncertain SOWs to the decision variables. This is why we use the meanObjs_acrossRDM.csv file – this file contains sixty-nine mean values of the performance objectives, where each set of performance objectives inversely maps to their corresponding portfolio of actions.
To identify the sensitivity of performance objectives to the DU factors, change lines 121 to 127 to the following:
# to change depending on whether DV or DUF is being analyzed
dec_vars = duf_numpy[:100,:]
measured_outcomes = objs_mean_soln_df
names = duf_names
mo_names = obj_names
bounds = duf_bounds[:100,:]
rdm = 'DUF'
mode = 'objs'
###
The code above identifies the sensitivity of the average values of all twenty performance objectives over each of the sixty-nine different portfolios to the set of deeply uncertain hydroclimatic and demand scenarios. This is why we use the meanObjs_acrossSoln.csv file – this file contains one hundred mean values of the performance objectives, where each set of performance objectives inversely maps to their corresponding DU SOW.
Great job so far! Now let’s visualize the sensitivity of our output to our input parameters using heatmaps. Before being able to visualize each utility’s performance sensitivity, we must first organize the sensitivity indices of the decision variables into a file containing the indices over all objectives for each utility. The gather_delta.py script does this. Simply change the value of mode on line 11 to ‘DUF‘ to gather the indices for the DU factors.
"""
Created on Tue April 26 2022 16:12
@author: Lillian Bei Jia Lau
Gathers the delta sensitivity indices into files per utility
"""
import numpy as np
import pandas as pd
mode = 'DV'
main_dir = '/yourFilePath/WaterPaths/post_processing/delta_output_' + mode + '/'
utilities = ['_C', '_D', '_R', '_reg']
objs = ['REL', 'RF', 'INF_NPC', 'PFC', 'WCC']
utilities_full = ['Cary', 'Durham', 'Raleigh', 'Regional']
dv_names = ['RT_C', 'RT_D', 'RT_R', 'TT_D', 'TT_R', 'LMA_C', 'LMA_D', 'LMA_R',\
'RC_C', 'RC_D', 'RC_R', 'IT_C', 'IT_D', 'IT_R', 'IP_C', 'IP_D', \
'IP_R', 'INF_C', 'INF_D', 'INF_R']
duf_names = ['Cary restr. eff', 'Durham restr. eff', 'Raleigh restr. eff', 'Demand growth\nmultiplier',\
'Bond term\nmultiplier', 'Bond interest\nrate multiplier', 'Infrastructure interest\nrate multiplier',\
'Streamflow amp\nmultiplier', 'SCR PT\nmultiplier', 'SCR CT\nmultiplier', 'NRR PT\nmultiplier',\
'NRR CT\nmultiplier', 'CR Low PT\nmultiplier', 'CR Low CT\nmultiplier', 'CR High PT\nmultiplier',\
'CR High CT\nmultiplier', 'WR1 PT\nmultiplier', 'WR1 CT\nmultiplier', 'WR2 PT\nmultiplier',\
'WR2 CT\nmultiplier', 'DR PT\nmultiplier', 'DR CT\nmultiplier', 'FR PT\nmultiplier', 'FR CT\nmultiplier',\
'DR PT\nmultiplier', 'DR CT\nmultiplier', 'FR PT\nmultiplier', 'FR CT\nmultiplier']
s1_dv_cary = np.zeros((len(objs), len(dv_names)), dtype=float)
s1_dv_durham = np.zeros((len(objs), len(dv_names)), dtype=float)
s1_dv_raleigh = np.zeros((len(objs), len(dv_names)), dtype=float)
s1_dv_regional = np.zeros((len(objs), len(dv_names)), dtype=float)
s1_dv_dict = {
'_C': s1_dv_cary,
'_D': s1_dv_durham,
'_R': s1_dv_raleigh,
'_reg': s1_dv_regional
}
s1_duf_cary = np.zeros((len(objs), len(duf_names)), dtype=float)
s1_duf_durham = np.zeros((len(objs), len(duf_names)), dtype=float)
s1_duf_raleigh = np.zeros((len(objs), len(duf_names)), dtype=float)
s1_duf_regional = np.zeros((len(objs), len(duf_names)), dtype=float)
s1_duf_dict = {
'_C': s1_duf_cary,
'_D': s1_duf_durham,
'_R': s1_duf_raleigh,
'_reg': s1_duf_regional
}
for i in range(len(utilities)):
s1_util = []
hdrs = []
if mode == 'DV':
s1_util = s1_dv_dict[utilities[i]]
hdrs = dv_names
elif mode == 'DUF':
s1_util = s1_duf_dict[utilities[i]]
hdrs = duf_names
for j in range(len(objs)):
curr_file = main_dir + 'S1_' + objs[j] + utilities[i] + '.csv'
s1_util[j, :] = pd.read_csv(curr_file, sep=',', skiprows=2, header=None).iloc[0,1:]
s1_util_df = pd.DataFrame(s1_util, columns=hdrs)
out_filepath = main_dir + utilities_full[i] + '.csv'
s1_util_df.to_csv(out_filepath, sep=',', index=False)
Now, let’s plot our heatmaps! The code to do so can be found in sensitivity_heatmap.py, and should look similar to the code provided below:
import numpy as np
import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt
from mpl_toolkits.axes_grid1 import AxesGrid
sns.set_theme()
# change depending on compromise solution and whether its sensitivity to DUF or DVs
mode = 'DUF'
rot = 90 # if DV use 0; if DUF use 45
main_dir = '/YourFilePath/WaterPaths/post_processing/delta_output_' + mode + '/'
c_filepath = main_dir + 'Cary.csv'
d_filepath = main_dir + 'Durham.csv'
r_filepath = main_dir + 'Raleigh.csv'
reg_filepath = main_dir + 'Regional.csv'
cary = pd.read_csv(c_filepath, index_col=False, header=0)
durham = pd.read_csv(d_filepath, index_col=False, header=0)
raleigh = pd.read_csv(r_filepath, index_col=False, header=0)
regional = pd.read_csv(reg_filepath, index_col=False, header=0)
savefig_dir = '/YourFilePath/WaterPaths/post_processing/'
savefig_name = savefig_dir + 'dmi_heatmap_' + mode + '.svg'
grid_kws = {"height_ratios": (0.20, 0.20, 0.20, 0.20, .02), "hspace": 0.5}
f, (ax1, ax2, ax3, ax4, cbar_ax) = plt.subplots(5, figsize=(15, 20), gridspec_kw=grid_kws)
plt.subplots_adjust(top = 0.95, bottom = 0.05,
hspace = 0, wspace = 0.05)
y_objs=['REL', 'RF', 'INPC', 'PFC', 'WCC']
x_dvs=['$RT_{C}$', '$RT_{D}$', '$RT_{R}$', '$TT_{D}$', '$TT_{R}$', '$LM_{C}$', '$LM_{D}$', '$LM_{R}$',\
'$RC_{C}$', '$RC_{D}$', '$RC_{R}$', '$IT_{C}$', '$IT_{D}$', '$IT_{R}$', '$IP_{C}$', \
'$IP_{D}$', '$IP_{R}$','$INF_{C}$', '$INF_{D}$', '$INF_{R}$']
x_dufs = ['Cary\nrestr. eff', 'Durham\nrestr. eff', 'Raleigh\nrestr. eff', 'Dem. growth\nmultiplier',\
'Bond term\nmultiplier', 'Bond interest\nrate multiplier', 'Infra. interest\nrate multiplier',\
'Streamflow amp\nmultiplier', 'SCR PT\nmultiplier', 'SCR CT\nmultiplier', 'NRR PT\nmultiplier',\
'NRR CT\nmultiplier', 'CR Low PT\nmultiplier', 'CR Low CT\nmultiplier', 'CR High PT\nmultiplier',\
'CR High CT\nmultiplier', 'WR1 PT\nmultiplier', 'WR1 CT\nmultiplier', 'WR2 PT\nmultiplier',\
'WR2 CT\nmultiplier', 'DR PT\nmultiplier', 'DR CT\nmultiplier', 'FR PT\nmultiplier', 'FR CT\nmultiplier',\
'DR PT\nmultiplier', 'DR CT\nmultiplier', 'FR PT\nmultiplier', 'FR CT\nmultiplier']
x_labs = []
if mode == 'DV':
x_labs = x_dvs
elif mode == 'DUF':
x_labs = x_dufs
plt.rc('xtick', labelsize=1)
plt.rc('ytick', labelsize=3)
plt.rc('axes', labelsize=5)
plt.rc('axes', titlesize=14)
#ax1 = fig.add_subplot(411)
ax1.set_title("Cary")
sns.heatmap(cary, linewidths=.05, cmap="YlOrBr", xticklabels=[],
yticklabels=y_objs, ax=ax1, cbar=False)
ax1.set_yticklabels(y_objs, rotation=0)
#ax2 = fig.add_subplot(412)
ax2.set_title("Durham")
sns.heatmap(durham, linewidths=.05, cmap="YlOrBr", xticklabels=[],
yticklabels=y_objs, ax=ax2, cbar=False)
ax2.set_yticklabels(y_objs, rotation=0)
#ax3 = fig.add_subplot(413)
ax3.set_title("Raleigh")
sns.heatmap(raleigh, linewidths=.05, cmap="YlOrBr", xticklabels=[],
yticklabels=y_objs, ax=ax3, cbar=False)
ax3.set_yticklabels(y_objs, rotation=0)
#ax4 = fig.add_subplot(414, fontsize=10)
ax4.set_title("Regional")
ax4 = sns.heatmap(regional, linewidths=.05, cmap="YlOrBr", xticklabels=x_labs,
yticklabels=y_objs, ax=ax4, cbar=True, cbar_ax=cbar_ax,
cbar_kws={'orientation': 'horizontal'}) # change depending on whether analyzing DUF or DV
ax4.set_xticklabels(x_labs, rotation=rot, fontsize=10)
ax4.set_yticklabels(y_objs, rotation=0)
plt.savefig(savefig_name)
Running this for the sensitivity to decision variables and DU factors will generate the following images:
In the figure above, the color of each box represents the sensitivity of a performance objective (y-axis) to a specific decision variable (x-axis). It is interesting to note that the restriction triggers (RT) of all utilities strongly influence each of their individual and regional reliability and restriction frequency. This indicates the potential for regional conflict, as possible errors is operating one utility’s restriction trigger may adversely affect other utilities’ reliabilities and ability to maintain full control over their own use of water-use restrictions. Furthermore, Raleigh’s performance is sensitive to more decision variables than its remaining two counterparts, with it’s worst-case cost (WCC) being affected most by Cary’s infrastructure investments. This observation highlights the importance of careful cooperation between a region’s member utilities to ensure that all partners play their part in maintaining both their own and their counterparts’ performance.
In this next figure, we observe that uncertainty in demand growth is the only DU factor that significantly drives changes in individual and regional performance. This finding can thus help utilities to focus on demand management programs, or formulate operations and management policies that enable them to more quickly adapt to changes in consumer and industrial demand growth.
Overall, in this post, we have performed a simple sensitivity analysis to identify uncertainties in the decision variables and DU factors that control regional and individual performance. All the code for processing the output data can be found in this GitHub repository here. In the next post, we will end the MORDM blogpost series by performing scenario discovery to map regions of success and failure as defined by our robustness metrics.
References
Borgonovo, E. (2007). A new uncertainty importance measure. Reliability Engineering &Amp; System Safety, 92(6), 771-784. doi: 10.1016/j.ress.2006.04.015
Herman, J. D., Reed, P. M., Zeff, H. B., & Characklis, G. W. (2015). How should robustness be defined for water systems planning under change? Journal of Water Resources Planning and Management, 141(10), 04015012. doi:10.1061/(asce)wr.1943-5452.0000509
Reed, P.M., Hadjimichael, A., Malek, K., Karimi, T., Vernon, C.R., Srikrishnan, V., Gupta, R.S., Gold, D.F., Lee, B., Keller, K., Thurber, T.B, & Rice, J.S. (2022). Addressing Uncertainty in Multisector Dynamics Research [Book]. Zenodo. https://doi.org/10.5281/zenodo.6110623
Water Programming: A Collaborative Research Blog 