>>> # module imports needed
>>> import numpy as np
>>> import matplotlib.pyplot as plt
>>> import matplotlib.patches as patches
>>>
>>> def plot_feasibility(solutions, solver, samples=200, test_set_size=10):
... # seed the random number generator for deterministic sampling
... rng = np.random.default_rng(123456)
...
... # nominal uncertain parameter realizations
... nom_vals = np.array([10, 1635])
...
... # sample points from box uncertainty set of specified test size
... point_samples = np.empty((samples, 2))
... point_samples[0] = nom_vals
... point_samples[1:] = rng.uniform(
... low=nom_vals * (1 - test_set_size / 100),
... high=nom_vals * (1 + test_set_size / 100),
... size=(samples - 1, 2),
... )
...
... costs = np.empty((len(solutions), point_samples.shape[0]), dtype=float)
... mdl = build_model()
... for sol_idx, (size, sol) in enumerate(solutions.items()):
... # fix the first-stage variables
... mdl.Vhat.fix(sol[0])
... mdl.A.fix(sol[1])
...
... for pt_idx, pt in enumerate(point_samples):
... # update parameter realization to sampled point
... mdl.kR.set_value(pt[0])
... mdl.U.set_value(pt[1])
...
... # update the values of the operational variables
... initialize_model(mdl, Vhat=sol[0], A=sol[1])
...
... # try solving the model to inspect for feasibility
... res = solver.solve(mdl, load_solutions=False)
... if pyo.check_optimal_termination(res):
... mdl.solutions.load_from(res)
... costs[sol_idx, pt_idx] = pyo.value(mdl.obj)
... else:
... costs[sol_idx, pt_idx] = np.nan
...
... # now generate the plot(s)
... fig, axs = plt.subplots(
... figsize=(0.5 * (len(solutions) - 1) + 5 * len(solutions), 4),
... ncols=len(solutions),
... squeeze=False,
... sharey=True,
... layout="constrained",
... )
... for sol_idx, (size, ax) in enumerate(zip(solutions, axs[0])):
... # track realizations for which solution feasible
... is_feas = ~np.isnan(costs[sol_idx])
...
... # realizations under which design is feasible
... # (colored by objective)
... plot = ax.scatter(
... point_samples[is_feas][:, 0],
... point_samples[is_feas][:, 1],
... c=costs[sol_idx, is_feas],
... vmin=np.nanmin(costs),
... vmax=np.nanmax(costs),
... cmap="plasma_r",
... marker="o",
... )
... # realizations under which design is infeasible
... ax.scatter(
... point_samples[~is_feas][:, 0],
... point_samples[~is_feas][:, 1],
... color="none",
... edgecolors="black",
... label="infeasible",
... marker="^",
... )
... if size != 0:
... # boundary of the box uncertainty set mapped to the design
... rect = patches.Rectangle(
... nom_vals * (1 - size / 100),
... *tuple(nom_vals * 2 * size / 100),
... facecolor="none",
... edgecolor="black",
... linestyle="dashed",
... label=f"{size}% box set",
... )
... ax.add_patch(rect)
...
... ax.legend(bbox_to_anchor=(1, -0.15), loc="upper right")
... ax.set_xlabel(r"$k_\mathrm{R}$ (per hr)")
... ax.set_ylabel("$U$ (kJ/sqm-h-K)")
...
... is_in_set = np.logical_and(
... np.all(nom_vals * (1 - size / 100) <= point_samples, axis=1),
... np.all(point_samples <= nom_vals * (1 + size / 100), axis=1),
... )
... feas_in_set = np.logical_and(is_feas, is_in_set)
...
... # add plot title summarizing statistics of the results
... ax.set_title(
... f"Solution for {size}% box set\n"
... "Avg ± SD objective "
... f"{costs[sol_idx, is_feas].mean():.2f} ± {costs[sol_idx, is_feas].std():.2f}\n"
... f"Feas. for {feas_in_set.sum()}/{is_in_set.sum()} samples in set\n"
... f"Feas. for {is_feas.sum()}/{len(point_samples)} samples overall"
... )
...
... cbar = fig.colorbar(plot, ax=axs.ravel().tolist(), pad=0.03)
... cbar.ax.set_ylabel("Objective ($/yr)")
...
... plt.show()
...