{
"cells": [
{
"cell_type": "markdown",
"id": "4fd4792f",
"metadata": {},
"source": [
"# Substation Clustering for Multi-Root Wind Farms\n",
"\n",
"This notebook demonstrates how `OptiWindNet` uses **substation clustering** (`clusterize()`) to partition multi-substation wind farms into single-root sub-fields.\n",
"\n",
"## Key Motivation & Concepts\n",
"\n",
"1. **Enabling Meta-Heuristics (HGS & LKH-3) for Multi-Root Farms**:\n",
" CVRP algorithms like Hybrid Genetic Search (`hgs_cvrp`) and Lin-Kernighan-Helsgaun (`lkh3`) are single-depot solvers. They cannot solve multi-substation instances as a single global problem. `OptiWindNet` automatically uses `clusterize()` to decompose multi-root farms into single-depot sub-problems before passing them to HGS or LKH-3.\n",
"\n",
"2. **Avoiding Cross-Root Rings in Ringed Topologies**:\n",
" When optimizing a multi-substation site for a **ringed topology** with `MILPRouter(model_options=ModelOptions(topology='ringed'))`, the MILP formulation can form **cross-root rings** (rings that start at Substation A and return to Substation B). By running `clusterize()` first and optimizing each cluster independently, every ring is guaranteed to return to the *same* substation root (single-root rings)."
]
},
{
"cell_type": "markdown",
"id": "3308b673",
"metadata": {},
"source": [
"## Load Modules & Data"
]
},
{
"cell_type": "code",
"execution_count": 1,
"id": "464c6146",
"metadata": {},
"outputs": [],
"source": [
"from optiwindnet.importer import load_repository\n",
"from optiwindnet.mesh import make_planar_embedding\n",
"from optiwindnet.interarraylib import S_from_G, as_normalized, G_from_S, rings_from_S\n",
"from optiwindnet.clustering import clusterize\n",
"from optiwindnet.baselines.hgs import hgs_cvrp\n",
"from optiwindnet.api import WindFarmNetwork, MILPRouter, ModelOptions\n",
"from optiwindnet.pathfinding import PathFinder\n",
"from optiwindnet.svg import svgplot\n",
"import networkx as nx\n",
"import numpy as np\n",
"\n",
"%config InlineBackend.figure_formats = ['svg']"
]
},
{
"cell_type": "markdown",
"id": "cc60fccc",
"metadata": {},
"source": [
"## 1. Graphically Visualizing Substation Clusters\n",
"\n",
"Let's load two multi-substation wind farm layouts (**Borssele** with 2 substations and **Hornsea One** with 3 substations) and inspect how `clusterize()` partitions their layout meshes into distinct sub-fields. `svgplot()` colors turbines by their `subtree` node attribute, so assigning each substation cluster a distinct `subtree` value makes the partition visible."
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "61b10b19",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Borssele cluster sizes: [96, 77]\n"
]
},
{
"data": {
"image/svg+xml": [
""
],
"text/plain": [
""
]
},
"execution_count": 2,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"def color_clusters(L, clusters):\n",
" subtree = {\n",
" node: cluster_id for cluster_id, nodes in enumerate(clusters) for node in nodes\n",
" }\n",
" nx.set_node_attributes(L, subtree, name='subtree')\n",
"\n",
"\n",
"locations = load_repository()\n",
"\n",
"# Borssele (2 substations)\n",
"L_borssele = locations.borssele\n",
"P_borssele, A_borssele = make_planar_embedding(L_borssele)\n",
"clusters_borssele = clusterize(as_normalized(A_borssele), capacity=7)\n",
"\n",
"color_clusters(L_borssele, clusters_borssele)\n",
"\n",
"print(f'Borssele cluster sizes: {[len(nodes) for nodes in clusters_borssele]}')\n",
"svgplot(L_borssele)"
]
},
{
"cell_type": "code",
"execution_count": 3,
"id": "421df73b",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Hornsea One cluster sizes: [56, 49, 69]\n"
]
},
{
"data": {
"image/svg+xml": [
""
],
"text/plain": [
""
]
},
"execution_count": 3,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Hornsea One (3 substations)\n",
"L_hornsea = locations.hornsea\n",
"P_hornsea, A_hornsea = make_planar_embedding(L_hornsea)\n",
"clusters_hornsea = clusterize(as_normalized(A_hornsea), capacity=7)\n",
"\n",
"color_clusters(L_hornsea, clusters_hornsea)\n",
"\n",
"print(f'Hornsea One cluster sizes: {[len(nodes) for nodes in clusters_hornsea]}')\n",
"svgplot(L_hornsea)"
]
},
{
"cell_type": "markdown",
"id": "d9272371",
"metadata": {},
"source": [
"## 2. Multi-Root Meta-Heuristic Routing (`HGSRouter`)\n",
"\n",
"Because HGS solves single-depot CVRPs, `HGSRouter` and `hgs_cvrp()` use `clusterize()` under the hood to solve multi-root instances smoothly."
]
},
{
"cell_type": "code",
"execution_count": 4,
"id": "082138d5",
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
""
],
"text/plain": [
""
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# HGSRouter automatically decomposes multi-root layouts via clusterize()\n",
"S_hgs = hgs_cvrp(as_normalized(A_borssele), capacity=7, time_limit=1.0)\n",
"G_hgs = G_from_S(S_hgs, A_borssele)\n",
"H_hgs = PathFinder(G_hgs, P_borssele, A_borssele).create_detours()\n",
"svgplot(H_hgs)"
]
},
{
"cell_type": "markdown",
"id": "33c2876c",
"metadata": {},
"source": [
"## 3. Preventing Cross-Root Rings in Ringed Topologies\n",
"\n",
"A ringed MILP is free to use different roots at the two ends of a ring. On **Neart na Gaoithe**, the smallest packaged multi-substation site, an unrestricted capacity-4 model produces such a bridge. Both cases below use Gurobi with a 20-second limit per MILP solve. Ring endpoints are checked after routing by reconstructing the topology from `G`, because `PathFinder` may re-hook a tentative feeder.\n",
"\n",
"### Without `clusterize()`\n",
"\n",
"The unrestricted routed network joins both substations in one connected component. The SVG therefore shows a ring bridging the two green substation squares."
]
},
{
"cell_type": "code",
"execution_count": 5,
"id": "2d751569",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Set parameter WLSAccessID\n",
"Set parameter WLSSecret\n",
"Set parameter LicenseID to value 937681\n",
"Set parameter MIPFocus to value 1\n",
"Academic license 937681 - for non-commercial use only - registered to ma___@dtu.dk\n",
"Roots in each routed component: [[-2, -1]]\n",
"Cross-root rings: [((-2, -1), [44, 3, 37, 10, 48, 24, 23, 43])]\n"
]
},
{
"data": {
"image/svg+xml": [
""
],
"text/plain": [
""
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"ring_site = WindFarmNetwork(cables=4, L=locations.neart)\n",
"\n",
"\n",
"def make_ringed_milp():\n",
" return MILPRouter(\n",
" solver_name='gurobi',\n",
" time_limit=20,\n",
" mip_gap=0.0,\n",
" model_options=ModelOptions(topology='ringed'),\n",
" )\n",
"\n",
"\n",
"def roots_by_component(S):\n",
" roots = set(range(-S.graph['R'], 0))\n",
" return [\n",
" sorted(roots.intersection(component))\n",
" for component in nx.connected_components(S)\n",
" if roots.intersection(component)\n",
" ]\n",
"\n",
"\n",
"ring_site.optimize(router=make_ringed_milp())\n",
"S_unclustered_routed = S_from_G(ring_site.G)\n",
"rings_unclustered = rings_from_S(S_unclustered_routed)\n",
"bridges_unclustered = [\n",
" (roots, terminals) for roots, terminals in rings_unclustered if roots[0] != roots[1]\n",
"]\n",
"\n",
"print('Roots in each routed component:', roots_by_component(S_unclustered_routed))\n",
"print('Cross-root rings:', bridges_unclustered)\n",
"svgplot(ring_site.G)"
]
},
{
"cell_type": "markdown",
"id": "f9ad63db",
"metadata": {},
"source": [
"### With `clusterize()`\n",
"\n",
"`clusterize()` assigns each turbine to one substation. Each cluster is then built, optimized, and routed as a fresh single-root site. The already-routed cluster graphs are remapped and merged only for visualization, so neither the MILP nor `PathFinder` can reach the other root. This also avoids modifying the available-links graph or its `diagonals` metadata."
]
},
{
"cell_type": "code",
"execution_count": 6,
"id": "471162ef",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Set parameter WLSAccessID\n",
"Set parameter WLSSecret\n",
"Set parameter LicenseID to value 937681\n",
"Set parameter MIPFocus to value 1\n",
"Academic license 937681 - for non-commercial use only - registered to ma___@dtu.dk\n",
"Set parameter WLSAccessID\n",
"Set parameter WLSSecret\n",
"Set parameter LicenseID to value 937681\n",
"Set parameter MIPFocus to value 1\n",
"Academic license 937681 - for non-commercial use only - registered to ma___@dtu.dk\n",
"Cluster sizes: [23, 31]\n",
"Roots in each routed component: [[-2], [-1]]\n",
"Cross-root rings: []\n"
]
},
{
"data": {
"image/svg+xml": [
""
],
"text/plain": [
""
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"def make_cluster_site(site, cluster_id, cluster):\n",
" L = site.L\n",
" R = L.graph['R']\n",
" VertexC = L.graph['VertexC']\n",
" border = L.graph.get('border', np.array([], dtype=int))\n",
" obstacles = L.graph.get('obstacles', [])\n",
" return WindFarmNetwork(\n",
" cables=site.cables,\n",
" turbinesC=VertexC[sorted(cluster)],\n",
" substationsC=VertexC[[-R + cluster_id]],\n",
" borderC=VertexC[border],\n",
" obstacleC_=[VertexC[obstacle] for obstacle in obstacles],\n",
" name=f'{L.graph[\"name\"]} — cluster {cluster_id + 1}',\n",
" )\n",
"\n",
"\n",
"def node_coord(G, node):\n",
" fnT = G.graph.get('fnT')\n",
" coord_id = fnT[node] if fnT is not None else node\n",
" return G.graph['VertexC'][coord_id]\n",
"\n",
"\n",
"def merge_cluster_routes(site, solved):\n",
" \"\"\"Merge already-routed single-root graphs without rerunning PathFinder.\"\"\"\n",
" L = site.L\n",
" R, T, B = (L.graph[key] for key in 'RTB')\n",
" contours = []\n",
" detours = []\n",
" for terminals, cluster_id, G_local in solved:\n",
" T_local, B_local = (G_local.graph[key] for key in 'TB')\n",
" C_local = G_local.graph.get('C', 0)\n",
" contours.extend(\n",
" (cluster_id, node, node_coord(G_local, node))\n",
" for node in range(T_local + B_local, T_local + B_local + C_local)\n",
" )\n",
" detours.extend(\n",
" (cluster_id, node, node_coord(G_local, node))\n",
" for node in range(\n",
" T_local + B_local + C_local,\n",
" T_local + B_local + C_local + G_local.graph.get('D', 0),\n",
" )\n",
" )\n",
"\n",
" aux_coord = [coord for *_, coord in contours + detours]\n",
" VertexC = L.graph['VertexC']\n",
" merged_VertexC = np.vstack((VertexC[: T + B], *aux_coord, VertexC[-R:]))\n",
" G = nx.Graph()\n",
" G.graph.update(\n",
" L.graph\n",
" | {\n",
" 'VertexC': merged_VertexC,\n",
" 'capacity': site.cables_capacity,\n",
" 'topology': 'ringed',\n",
" 'C': len(contours),\n",
" 'D': len(detours),\n",
" 'cables': site.cables,\n",
" 'creator': 'clusterize + MILP.gurobi',\n",
" 'has_loads': True,\n",
" }\n",
" )\n",
" aux_map = {\n",
" (cluster_id, node): T + B + i\n",
" for i, (cluster_id, node, _) in enumerate(contours + detours)\n",
" }\n",
"\n",
" subtree_offset = 0\n",
" for terminals, cluster_id, G_local in solved:\n",
" mapping = {local: global_ for local, global_ in enumerate(terminals)}\n",
" mapping[-1] = -R + cluster_id\n",
" mapping.update(\n",
" {\n",
" node: global_\n",
" for (owner, node), global_ in aux_map.items()\n",
" if owner == cluster_id\n",
" }\n",
" )\n",
" local_subtrees = {\n",
" data['subtree']\n",
" for node, data in G_local.nodes.items()\n",
" if node >= 0 and 'subtree' in data\n",
" }\n",
" for node, data in G_local.nodes.items():\n",
" data = data.copy()\n",
" if 'subtree' in data:\n",
" data['subtree'] += subtree_offset\n",
" G.add_node(mapping[node], **data)\n",
" G.add_edges_from(\n",
" (mapping[u], mapping[v], data.copy())\n",
" for u, v, data in G_local.edges(data=True)\n",
" )\n",
" subtree_offset += len(local_subtrees)\n",
"\n",
" G.graph['max_load'] = max(\n",
" (data['load'] for _, _, data in G.edges(data=True)), default=0\n",
" )\n",
" return G\n",
"\n",
"\n",
"clusters_ring = clusterize(\n",
" as_normalized(ring_site.A), capacity=ring_site.cables_capacity\n",
")\n",
"solved_clusters = []\n",
"for cluster_id, cluster in enumerate(clusters_ring):\n",
" cluster_site = make_cluster_site(ring_site, cluster_id, cluster)\n",
" cluster_site.optimize(router=make_ringed_milp())\n",
" solved_clusters.append((sorted(cluster), cluster_id, cluster_site.G))\n",
"\n",
"G_clustered = merge_cluster_routes(ring_site, solved_clusters)\n",
"S_clustered_routed = S_from_G(G_clustered)\n",
"rings_clustered = rings_from_S(S_clustered_routed)\n",
"bridges_clustered = [\n",
" (roots, terminals) for roots, terminals in rings_clustered if roots[0] != roots[1]\n",
"]\n",
"\n",
"print('Cluster sizes:', [len(cluster) for cluster in clusters_ring])\n",
"print('Roots in each routed component:', roots_by_component(S_clustered_routed))\n",
"print('Cross-root rings:', bridges_clustered)\n",
"svgplot(G_clustered)"
]
},
{
"cell_type": "markdown",
"id": "90ba28f5",
"metadata": {},
"source": [
"## Summary\n",
"\n",
"* **Meta-heuristics**: `HGSRouter` and `LKH` rely on `clusterize()` to handle multi-substation farms.\n",
"* **Cluster visualization**: assigning the cluster index to each turbine's `subtree` attribute lets `svgplot()` display the partition directly.\n",
"* **Ringed MILP**: Neart's unrestricted capacity-4 solution bridges roots, while independently optimizing and routing each cluster guarantees single-root rings."
]
}
],
"metadata": {
"language_info": {
"name": "python"
}
},
"nbformat": 4,
"nbformat_minor": 5
}