In my last blog post, I showed how to run an MORDM experiment using Rhodium. This process included the multi-objective optimization to an assumed state of the world (SOW) as well as the re-evaluation of the Pareto-approximate solutions on alternative SOWs, before using sensitivity and classification tools such as PRIM and CART for the scenario discovery analysis. However, there may be cases where you have to run the optimization and re-evaluation outside of Rhodium, for instance if your model is in another programming language than Python. There are two ways you can do this while still using Rhodium for the scenario discovery. The first option is to run the model through the Executioner. Another option is to run the model separately and import the output into the same format as is generated by Rhodium for post-analysis. I will explain the second method here for the fish game described in my last post.

The first step is to read the decision variables and objectives from the optimization into 2D arrays. Then the uncertainties, levers and responses can be defined as before, except they no longer have to be associated with an object of the class ‘Model‘.

# read in output of optimization variables = np.loadtxt('FishGame/FishGame.resultfile',usecols=[0,1]) objectives = np.loadtxt('FishGame/FishGame.resultfile',usecols=[2,3]) # make maximization objectives positive maxIndices = [0] objectives[:,maxIndices] = -objectives[:,maxIndices] # define X of XLRM framework uncertainties = [UniformUncertainty("a", 1.5, 4.0), UniformUncertainty("b0", 0.25, 0.67)] # define L of XLRM framework levers = [RealLever("vars", 0.0, 1.0, length=2)] # define R of XLRM framework responses = [Response("NPVharvest", Response.MAXIMIZE), Response("std_z", Response.MINIMIZE)]

Note: If you are interested in using Rhodium’s plotting tools to visualize the results of the optimization, you can still make the uncertainties, levers and responses attributes of a model object. However, you will have to create a model function to instantiate the model. This is sloppy, but you can fake this by just creating a function that takes in the decision variables and model parameters, and returns the objective values, but doesn’t actually perform any calculations.

def fishGame(vars, a = 1.75, # rate of prey growth b0 = 0.6, # initial rate of predator growth F = 0, # rate of change of radiative forcing per unit time S = 0.5): # climate sensitivity) NPVharvest = None std_z = None return (NPVharvest, std_z) model = Model(fishGame) # define all parameters to the model that we will be studying model.parameters = [Parameter("vars"), Parameter("a"), Parameter("b0"), Parameter("F"), Parameter("S")]

If using Rhodium for the optimization, this function would actually perform the desired calculation and Platypus could be used for the optimization. Since we have already performed the optimization, we just need to reformat the output of the optimization into that used by Rhodium for the RDM analysis. This can be done by mimicking the output structure that would be returned by the function ‘optimize‘.

# find number of solutions nsolns = np.shape(objectives)[0] # properly format output of optimization output = DataSet() for i in range(nsolns): env = OrderedDict() offset = 0 for lever in levers: if lever.length == 1: env[lever.name] = list(variables[i,:]) else: env[lever.name] = list(variables[i,offset:offset+lever.length]) offset += lever.length for j, response in enumerate(responses): env[response.name] = objectives[i,j] output.append(env) # write output to file with open("FishGame/FishGame_data.txt","w") as f: json.dump(output, f)

Next we need to read in the uncertain parameters that were sampled for the re-evaluation and format the results of the re-evaluation into the same format as would be output by calling ‘evaluate‘ within Rhodium. Below is an example with the first solution (soln_index=0).

# read in LH samples of uncertain parameters and determine # of samples LHsamples = np.loadtxt('FishGame/LHsamples.txt') nsamples = np.shape(LHsamples)[0] # load policy from optimization soln_index = 0 policy = output[soln_index] # load its objective values from re-evaluation and make maximization objectives positive objectives = np.loadtxt('FishGame/MORDMreeval/FishGame_Soln' + str(soln_index+1) + '.obj') objectives[:,maxIndices] = -objectives[:,maxIndices] # convert re-evaluation output to proper format results = DataSet() for j in range(nsamples): env = OrderedDict() offset = 0 for k, uncertainty in enumerate(uncertainties): env[uncertainty.name] = LHsamples[j,k] for k, response in enumerate(responses): env[response.name] = objectives[j,k] for lever in levers: if lever.length == 1: env[lever.name] = list(variables[soln_index,:]) else: env[lever.name] = list(variables[soln_index,offset:offset+lever.length]) offset += lever.length results.append(env) # write results to file with open("FishGame/FishGame_Soln" + str(soln_index+1) + "_reeval.txt","w") as f: json.dump(results, f)

Finally, you have to define the metrics.

# calculate M of XLRM framework metric = ["Profitable" if v["NPVharvest"] >= 3.0 else "Unprofitable" for v in results]

Then you can run PRIM and CART. This requires defining the names, or ‘keys’, of the uncertain parameters. If you created a fake model object, you can pass ‘include=model.uncertainties.keys()’ to the functions Prim() and Cart(). If not, you have to create your own list of ‘keys’ as I do below.

keys = [] for i in range(len(uncertainties)): keys.append(uncertainties[i].name) # run PRIM and CART on metrics p = Prim(results, metric, include=keys, coi="Profitable") box = p.find_box() box.show_details() plt.show() c = Cart(results, metrics[j], include=keys) c.print_tree(coi="Profitable") c.show_tree() plt.show()

The above code creates the following two figures.

If you had run the analysis using Sobol samples, you could use the SALib wrapper to calculate sensitivity indices and make bar charts or radial convergence plots of the results. (Note: My previous post did not show how to make these plots, but has since been updated. Check it out here.)

import seaborn as sns from SALib.analyze import sobol from SALib.util import read_param_file # Read the parameter range file and Sobol samples problem = read_param_file('FishGame/uncertain_params.txt') param_values = np.loadtxt('FishGame/SobolSamples.txt') # Load the first solution Y = np.loadtxt('FishGame/SobolReeval/FishGame_Soln' + (soln_index+1) + '.obj') # Evaluate sensitivity to the first objective, NPVharvest obj_index = 0 Si = sobol.analyze(problem, Y[:,obj_index], calc_second_order=True, conf_level=0.95, print_to_console=False) pretty_result = get_pretty_result(Si) sns.set() fig1 = pretty_result.plot() fig2 = pretty_result.plot_sobol(threshold=0.01,groups={"Prey Growth Parameters" : ["a"], "Predator Growth Parameters" : ["b0"]}) def get_pretty_result(result): pretty_result = SAResult(result["names"] if "names" in result else problem["names"]) if "S1" in result: pretty_result["S1"] = {k : float(v) for k, v in zip(problem["names"], result["S1"])} if "S1_conf" in result: pretty_result["S1_conf"] = {k : float(v) for k, v in zip(problem["names"], result["S1_conf"])} if "ST" in result: pretty_result["ST"] = {k : float(v) for k, v in zip(problem["names"], result["ST"])} if "ST_conf" in result: pretty_result["ST_conf"] = {k : float(v) for k, v in zip(problem["names"], result["ST_conf"])} if "S2" in result: pretty_result["S2"] = _S2_to_dict(result["S2"], problem) if "S2_conf" in result: pretty_result["S2_conf"] = _S2_to_dict(result["S2_conf"], problem) if "delta" in result: pretty_result["delta"] = {k : float(v) for k, v in zip(problem["names"], result["delta"])} if "delta_conf" in result: pretty_result["delta_conf"] = {k : float(v) for k, v in zip(problem["names"], result["delta_conf"])} if "vi" in result: pretty_result["vi"] = {k : float(v) for k, v in zip(problem["names"], result["vi"])} if "vi_std" in result: pretty_result["vi_std"] = {k : float(v) for k, v in zip(problem["names"], result["vi_std"])} if "dgsm" in result: pretty_result["dgsm"] = {k : float(v) for k, v in zip(problem["names"], result["dgsm"])} if "dgsm_conf" in result: pretty_result["dgsm_conf"] = {k : float(v) for k, v in zip(problem["names"], result["dgsm_conf"])} if "mu" in result: pretty_result["mu"] = {k : float(v) for k, v in zip(result["names"], result["mu"])} if "mu_star" in result: pretty_result["mu_star"] = {k : float(v) for k, v in zip(result["names"], result["mu_star"])} if "mu_star_conf" in result: pretty_result["mu_star_conf"] = {k : float(v) for k, v in zip(result["names"], result["mu_star_conf"])} if "sigma" in result: pretty_result["sigma"] = {k : float(v) for k, v in zip(result["names"], result["sigma"])} return pretty_result def _S2_to_dict(matrix, problem): result = {} names = list(problem["names"]) for i in range(problem["num_vars"]): for j in range(i+1, problem["num_vars"]): if names[i] not in result: result[names[i]] = {} if names[j] not in result: result[names[j]] = {} result[names[i]][names[j]] = result[names[j]][names[i]] = float(matrix[i][j]) return result

So don’t feel like you need to run your optimization and re-evaluation in Python in order to use Rhodium!

Pingback: Water Programming Blog Guide (Part I) – Water Programming: A Collaborative Research Blog