Exploring the Near-Optimal Feasible Space¶
Here, we use modelling to generate alternatives (MGA) to explore the near-opimal feasible space of a simple one-node model in the style of model.energy (the same as in the previous MGA example).
While we explored the use of n.optimize.optimize_mga() in the previous example, PyPSA includes a few additional functions making it convenient to work with near-optimal spaces in a geometric way. While the optimize_mga function solves a network with an alternative objective function (given by the weights argument), n.optimize_mga_in_direction is a simple alternative useful for exploring trade-offs between multiple alternative objectives. It takes two key arguments: a dictionary of dimensions, specifying multiple alternative objectives in the same format as the previously mentioned weights, and a direction representing a vector in the coordinate space defined by dimensions.
In this example, we explore the trade-offs between wind and solar expansion in a simple renewables-based system. Therefore, we define our dimensions (or alternative objectives) as total installed wind capacity and total installed solar capacity in MW, respectively. These dimensions can also be seen as defining a projection of near-optimal space of the model down to two dimensions.
:::{note} See also this research article for a background on dimension reduction and approximation of near-optimal spaces. :::
import matplotlib.pyplot as plt
from matplotlib.collections import LineCollection
from scipy.spatial import ConvexHull, convex_hull_plot_2d
import pypsa
Before any near-optimal analysis can start, the model has to be solved to optimality first. This is in order to find the minimum system cost, which is later used to define a system cost bound for near-optimality (i.e. "slack", usually in percentage of minimum system cost).
n = pypsa.examples.model_energy()
n.snapshots = n.snapshots[::3] # Crude reduction of time resolution for speed
n.generators.at["load shedding", "marginal_cost"] = 20000
n.optimize(solver_name="gurobi")
Individual dimensions are specified just like the weights argument for the n.optimize.optimize_mga function: a nested dictionary of components, variables and finally a series or dictionary giving weights to individual components.
Here, we define two dimensions, and give them the user-defined names "wind" and "solar" (the keys in the outer dictionary). These names are arbitrary and don't have to match names of underlaying components.
dimensions = {
"wind": {"Generator": {"p_nom": {"wind": 1}}},
"solar": {"Generator": {"p_nom": {"solar": 1}}},
}
Suppose we are interested in minimizing wind capacity alone, and also interested in minimizing the total capacity of wind and solar. These correspond to the directions (-1, 0) and (-1, -1) in wind-solar coordinate space, respectively. Let's see what is possible within a 5% total system cost slack.
_, _, min_wind = n.optimize.optimize_mga_in_direction(
dimensions=dimensions,
direction={
"wind": -1,
"solar": 0,
}, # Use coordinate names matching those given in `dimensions`
slack=0.05,
)
_, _, min_wind_solar = n.optimize.optimize_mga_in_direction(
dimensions=dimensions,
direction={"wind": -1, "solar": -1},
slack=0.05,
)
display(min_wind)
display(min_wind_solar)
Results are returned as dictionaries represting coordinates in the wind-solar space. The dictionaries use the same user-given coordinate names that we passed in the dimensions argument. We see that minimizing wind alone leads to a large installation of solar, while minimizing both jointly leads to a more balanced solution.
In order to get a full picture of the near-optimal space (projected to the wind and solar dimensions), it is useful to optimize in many directions and consider the convex hull of the resulting points.
Of course, you can use optimize_mga_in_direction manually to do this. For convenience, PyPSA also provides the optimize_mga_in_multiple_directions function. This function takes a list of directions (or a pandas DataFrame with rows for directions), and uses built-in Python multiprocessing functionality to optimize in these directions in parallel. In the event of any individual optimisations failing, only the successful ones are returned.
While the user is free to select directions themselves, PyPSA also includes some simple functions (for instance, pypsa.optimization.mga.generate_directions_randomber of random directions in the format expected by optimize_mga_in_multiple_directions`.
directions = pypsa.optimization.mga.generate_directions_random(
dimensions.keys(), 10, seed=0
)
dirs, pts = n.optimize.optimize_mga_in_multiple_directions(
dimensions=dimensions,
directions=directions,
slack=0.05,
solver_name="gurobi",
# Solve in up to 8 directions in parallel.
max_parallel=8,
)
The results are returned in the form of two pandas DataFrames; one for directions and one for resulting points / coordinates. Using scipy, we can compute and plot the convex hull of the points we generated.
display(pts.head() / 1e3)
fig, ax = plt.subplots(figsize=(8, 6))
convex_hull_plot_2d(ConvexHull(pts / 1e3), ax)
ax.axis("equal")
ax.set_xlabel("Wind capacity (GW)")
ax.set_ylabel("Solar capacity (GW)")
ax.set_title("Near-optimal space in wind-solar coordinate space")
plt.show()
While random directions generally work well, they are not always perfectly spaced out. In two dimensions, we can easily generate evenly spaced directions (points on a circle); implemented in PyPSA in pypsa.optimization.mga.generate_directions_evenly_spaced. In higher dimensions, generating evenly spaced directions is not trivial, but Halton sequences (pypsa.optimization.mga.generate_directions_halton) do a better job than randomly generated directions.
n_points = 20
directions_random = pypsa.optimization.mga.generate_directions_random(
dimensions.keys(), n_points, seed=0
)
directions_even = pypsa.optimization.mga.generate_directions_evenly_spaced(
dimensions.keys(), n_points
)
directions_halton = pypsa.optimization.mga.generate_directions_halton(
dimensions.keys(), n_points, seed=0
)
fig, axs = plt.subplots(1, 3, figsize=(12, 3), sharey=True)
directions_random.plot(kind="scatter", x="wind", y="solar", ax=axs[0])
directions_even.plot(kind="scatter", x="wind", y="solar", ax=axs[1])
directions_halton.plot(kind="scatter", x="wind", y="solar", ax=axs[2])
for ax in axs:
ax.set_xlim(-1.1, 1.1)
ax.set_ylim(-1.1, 1.1)
ax.axis("equal")
plt.show()
fig, axs = plt.subplots(1, 3, figsize=(12, 3), sharey=True)
for ax, directions in zip(axs, [directions_random, directions_even, directions_halton]):
_, pts = n.optimize.optimize_mga_in_multiple_directions(
dimensions=dimensions,
directions=directions,
slack=0.05,
solver_name="gurobi",
max_parallel=8,
)
convex_hull_plot_2d(ConvexHull(pts / 1e3), ax)
ax.axis("equal")
ax.set_xlabel("Wind capacity (GW)")
ax.set_ylabel("Solar capacity (GW)")
plt.show()
While more difficult to visualize, more than two dimensions can be used to analyse near-optimal spaces. Suppose we are interested in the mutual trade-offs between wind, solar and energy storage. Storage, as a sum of battery and hydrogen storage, is easily defined as a third dimension. In this case, since we are comparing different kinds of components with heterogenous units, it makes sense to scale our dimensions by capital cost in order to get comparable investment numbers.
dimensions = {
"wind": {"Generator": {"p_nom": {"wind": n.generators.at["wind", "capital_cost"]}}},
"solar": {
"Generator": {"p_nom": {"solar": n.generators.at["solar", "capital_cost"]}}
},
"storage": {
"StorageUnit": {
"p_nom": {
"battery storage": n.storage_units.at["battery storage", "capital_cost"]
}
},
"Store": {
"e_nom": {
"hydrogen storage": n.stores.at["hydrogen storage", "capital_cost"]
}
},
"Link": {"p_nom": n.links.capital_cost},
},
}
directions = pypsa.optimization.mga.generate_directions_halton(
dimensions.keys(), 50, seed=0
)
dirs, pts = n.optimize.optimize_mga_in_multiple_directions(
dimensions=dimensions,
directions=directions,
slack=0.05,
solver_name="gurobi",
max_parallel=8,
)
fig, ax = plt.subplots(figsize=(8, 6))
ax.scatter(
pts["wind"] / 1e6, pts["solar"] / 1e6, c=pts["storage"] / 1e6, cmap="viridis"
)
hull = ConvexHull(pts[["wind", "solar"]] / 1e6)
line_segments = [hull.points[simplex] for simplex in hull.simplices]
ax.add_collection(LineCollection(line_segments, colors="k", linestyle="solid"))
ax.axis("equal")
ax.set_xlabel("Wind investment (mEUR)")
ax.set_ylabel("Solar investment (mEUR)")
cbar = plt.colorbar(ax.collections[0], ax=ax)
cbar.set_label("Storage investment (mEUR)")
In the above example, we see that low investment in renewables necessitates higher investment in storage, while high investment in renewables excludes (via the total system cost bound) a high investment in storage.