Global Optimization of Offshore Wind Farm Collection Systems

TL;DR

This paper presents a MILP-based global optimization method for offshore wind farm cable layouts, incorporating electrical losses and using an iterative approach, demonstrated on real-world large-scale cases with promising results.

Contribution

Introduces a novel MILP model with pre-processing and iterative solving for optimizing offshore wind farm cable layouts considering electrical losses.

Findings

Feasible solutions found in seconds.

Near-optimal solutions obtained within hours.

Method outperforms existing state-of-the-art approaches.

Abstract

A mathematical program for global optimization of the cable layout of Offshore Wind Farms (OWFs) is presented. The model consists on a Mixed Integer Linear Program (MILP). Modern branch-and-cut solvers are able to solve large-scale instances, defined by more than hundred Wind Turbines (WTs), and a reasonable number of Offshore Substations (OSSs). In addition to the MILP model to optimize total cable length or initial investment, a pre-processing strategy is proposed in order to incorporate total electrical power losses into the objective function. High fidelity models are adapted to calculate cables current capacities, spatial currents. The MILP model is embedded in an iterative algorithmic framework, consisting in solving a sequence of problems with increasing size of the search space. The search space is defined as a set of underlying candidate arcs. The applicability of the method is illustrated through 10 case studies of real-world large-scale wind farms. Results show that: (i) feasible points can quickly be obtained in seconds, (ii) points near the global optimum with an imposed maximum tolerance, are calculable in reasonable computational time in the order of hours, and (iii) the proposed method compares favorably against a state-of-the art method available in literature.

Tables

Table 1. TABLE I: Main inputs parameters for sensitivity analysis

OWF	Obj.	$r$ [$\%$]	$c_e$ [€ /$\mathrm{MWh}$]	$m$	$P_n$ [$\mathrm{MW}$]	$n_w$	$n_o$	$\varphi$	$\eta$	$V_n$ [$\mathrm{kV}$]	$𝑻$	$𝑼$	${𝑪}_{𝒄}+{𝑪}_{𝒑}$ [M€ /$\mathrm{km}$]
WDS	IP	$5$	$40$	$30$	$3.6$	$108$	$1$	$10$	$1$	$33$	$\{1,2,3\}$	$\{7,10,13\}$	$\{0.36,0.58,0.90\}$
TH	I	-	-	-	$3$	$100$	$1$	$10$	$1$	$33$	$\{4,5\}$	$\{7,10\}$	$\{0.44,0.62\}$

Table 2. TABLE II: In-detail results summary

	$k_f$/		${\upsilon }_f$/		CT1		CT2		Obj	${ϵ}_{k_f}$/${ϵ}_{k_o}$	${ϵ}_{k_g}$
	$k_o$		${\upsilon }_o$		Time	Unit	Time	Unit	[M€ ]	[$\%$]	[$\%$]
W	$k_f\leftarrow$	$1$	$v_f\leftarrow$	$5$	$7.89$	s	$36.93$	s	$72.84$	$0$	$38.62$
D	$k_o\leftarrow$	$1$	$v_o\leftarrow$	$15$	$0.08$	h	$1.01$	h	$45.70$	$0.46$	$2.17$
S	$k_o\leftarrow$	$2$	$v_o\leftarrow$	$20$	$0.11$	h	$0.36$	h	$44.93$	$0.47$	$0.50$
	$k_o\leftarrow$	$3$	$v_o\leftarrow$	$25$	$0.14$	h	$0.16$	h	$44.93$	$0.50$	$0.50$
T	$k_f\leftarrow$	$1$	$v_f\leftarrow$	$5$	$6.65$	s	$10.37$	s	$40$	$0$	$33.60$
H	$k_o\leftarrow$	$1$	$v_o\leftarrow$	$15$	$15.80$	s	$1.02$	h	$26.64$	$0.30$	$0.30$
	$k_o\leftarrow$	$2$	$v_o\leftarrow$	$20$	$30.22$	s	$0.75$	h	$26.64$	$0.30$	$0.30$

Table 3. TABLE III: Main inputs parameters for benchmarking

Ins.	OWF	Obj.	$r$ [$\%$]	$c_e$ [€ /$\mathrm{MWh}$]	$m$	$P_n$ [$\mathrm{MW}$]	$n_w$	$n_o$	$\varphi$	$\eta$	$V_n$ [$\mathrm{kV}$]	$𝑻$	$𝑼$	${𝑪}_{𝒄}+{𝑪}_{𝒑}$ [M€ /$\mathrm{km}$]
1	HR1	I	-	-	-	$2$	$80$	$1$	$10$	$1$	$33$	$\{6,7,8\}$	$\{7,11,13\}$	$\{0.37,0.39,0.43\}$
2	HR1	I	-	-	-	$2$	$80$	$1$	$10$	$1$	$33$	$\{9,10\}$	$\{7,12\}$	$\{0.44,0.45\}$
3	HR1	I	-	-	-	$2$	$80$	$1$	$10$	$1$	$33$	$\{4,5\}$	$\{10,14\}$	$\{0.44,0.62\}$
4	O	I	-	-	-	$5$	$30$	$1$	$4$	$1$	$33$	$\{11,12\}$	$\{5,10\}$	$\{0.41,0.61\}$
5	O	I	-	-	-	$5$	$30$	$1$	$4$	$1$	$33$	$\{13,14\}$	$\{4,9\}$	$\{0.38,0.63\}$
6	DT	I	-	-	-	$3.6$	$80$	$1$	$10$	$1$	$33$	$\{6,7,8\}$	$\{4,6,8\}$	$\{0.37,0.39,0.43\}$
7	DT	I	-	-	-	$3.6$	$80$	$1$	$10$	$1$	$33$	$\{4,5\}$	$\{6,8\}$	$\{0.44,0.62\}$
8	TH	I	-	-	-	$3$	$100$	$1$	$10$	$1$	$33$	$\{13,14\}$	$\{7,15\}$	$\{0.38,0.63\}$
9	TH	I	-	-	-	$3$	$100$	$1$	$10$	$1$	$33$	$\{4,5\}$	$\{7,10\}$	$\{0.44,0.62\}$
\cdashline1-15 10	LA	IP	$5$	$40$	$30$	$3.6$	$175$	$2$	$10$	$1$	$33$	$\{1,2,3\}$	$\{7,10,13\}$	$\{0.36,0.58,0.90\}$

Table 4. TABLE IV: Results for benchmarking

Ins.	OWF	Obj.[9] [M€ ]	Obj. [M€ ]	Diff.[M€ ]	${ϵ}_{k_g}$[9] [%]	${ϵ}_{k_g}$ [%]	Diff.[%]	Time[9] [min]	Time [min]	Diff.[min]
1	HR1	$19.44$	$19.44$	$0$	$0.01$	$0.01$	$0$	$30$	$1.57$	$28.43$
2	HR1	$22.61$	$22.61$	$0$	$0.01$	$0.01$	$0$	$30$	$1$	$29$
3	HR1	$23.48$	$23.48$	$0$	$0.17$	$0.01$	$0.16$	$1440$	$7.26$	$1432.74$
4	O	$8.05$	$8.05$	$0$	$0$	$0.01$	$-0.01$	$0.42$	$0.46$	$-0.04$
5	O	$8.36$	$8.36$	$0$	$0$	$0.01$	$-0.01$	$1.95$	$0.45$	$1.50$
6	DT	$38.98$	$38.98$	$0$	$4.85$	$0.01$	$4.84$	$30$	$1.88$	$28.12$
7	DT	$50.38$	$49.83$	$0.55$	$7.36$	$0.01$	$7.35$	$10$	$2.65$	$7.35$
8	TH	$22.34$	$22.31$	$0.03$	$3.37$	$0.7$	$2.67$	$5$	$627.60$	$-622.60$
9	TH	$26.64$	$26.64$	$0$	$2.51$	$0.3$	$2.21$	$1440$	$107.02$	$1332.98$
\cdashline1-11 10	LA	-	$68.97$	-	-	$0.5$	-	-	$748$	-

Equations

$$u_t=\lfloor \frac{S_{r_t}=\sqrt{3}\cdot V_n\cdot I_t}{P_n}\rfloor \forall t\in \mathbf{T}$$

$$\mathbf{U}=\left\{1,\cdots ,u_{|\mathbf{T}|}\right\}$$

$${\vec{I}}_{ij,t}^k=\frac{k\cdot P_n}{\sqrt{3}\cdot V_n}\cdot \cosh ({\vec{\gamma }}_t\cdot d_{ij})-\frac{V_n}{{\overrightarrow{Z_c}}_t}\cdot \sinh ({\vec{\gamma }}_t\cdot d_{ij})$$

$$S_{ij,t}^k=\sqrt{3}\cdot V_n\cdot |{\vec{I}}_{ij,t}^k|$$

$$\sum _{i\in \mathbf{N}}\sum _{k=1}^{f(i)}k\cdot y_{ij}^k\cdot p^{\omega }-\sum _{i\in {\mathbf{N}}_{\mathbf{w}}}\sum _{k=1}^{f(i)}k\cdot y_{ji}^k\cdot p^{\omega }=p^{\omega }$$

$$\forall j\in {\mathbf{N}}_{\mathbf{w}}\wedge \omega \in {\mathbf{\Omega }}^{\mu }\wedge \mu \in \mathbf{M}$$

$${\vec{I}}_{ij,t}^{\omega ,k}=\frac{f_{ij}^{\omega ,k}}{\sqrt{3}\cdot V_n}\cdot \cosh ({\vec{\gamma }}_t\cdot d_{ij})-\frac{V_n}{{\overrightarrow{Z_c}}_t}\cdot \sinh ({\vec{\gamma }}_t\cdot d_{ij})$$

$$l_{ij,t}^{\mu ,k}\approx 3\cdot (1+{\lambda }_1+{\lambda }_2)\cdot \sum _{\omega \in {\mathbf{\Omega }}^{\mu }}\omega \cdot R_t\cdot d_{ij}\cdot |{\vec{I}}_{ij,t}^{\omega ,k}{|}^2$$

$$+3\cdot |{\mathbf{\Omega }}^{\mu }|\cdot \omega \cdot W_{d_t}\cdot d_{ij}$$

$$c_{ij}^k=\text{min}$$

$$\sum _{t\in \mathbf{T}}{x}_{ij,t}=1$$

$$x_{ij,t}\cdot (S_{ij,t}^k-S_{r_t})\le 0\forall t\in \mathbf{T}$$

$$x_{ij,t}\in \left\{0,1\right\}\forall t\in \mathbf{T}$$

$$c_{c_t}=a_{p_t}+b_{p_t}\cdot e^{{\left(\frac{c_{p_t}{S}_{r_t}}{10^8}\right)}^2}$$

$$\min \sum _{i\in \mathbf{N}}\sum _{j\in {\mathbf{N}}_{\mathbf{w}}}\sum _{k=1}^{f(i)}{c}_{ij}^k\cdot y_{ij}^k$$

$$\sum _{i\in {\mathbf{N}}_{\mathbf{o}}}{\sigma }_i=n_w$$

$$\sum _{j\in {\mathbf{N}}_{\mathbf{w}}}\sum _{k=1}^{f(j)}k\cdot y_{ij}^k={\sigma }_i\forall i\in {\mathbf{N}}_{\mathbf{o}}$$

$$\sum _{j\in {\mathbf{N}}_{\mathbf{w}}}\sum _{k=1}^{f(j)}{y}_{ij}^k\le \varphi \forall i\in {\mathbf{N}}_{\mathbf{o}}$$

$$\sum _{i\in \mathbf{N}}\sum _{k=1}^{f(i)}{y}_{ij}^k=1\forall j\in {\mathbf{N}}_{\mathbf{w}}$$

$$\sum _{i\in \mathbf{N}}\sum _{k=1}^{f(i)}k\cdot y_{ij}^k-\sum _{i\in {\mathbf{N}}_{\mathbf{w}}}\sum _{k=1}^{f(i)}k\cdot y_{ji}^k=1\forall j\in {\mathbf{N}}_{\mathbf{w}}$$

$$x_{ij}+x_{ji}+x_{uv}+x_{vu}\le 1\forall \left\{(i,j),(u,v)\right\}\in \mathbf{\chi }$$

$$\sum _{k=1}^{f(i)}{y}_{ij}^k-x_{ij}\le 0\forall (i,j)\in {\mathbf{A}}_{\mathbf{r}}$$

$$-\sum _{i\in \mathbf{N}}\sum _{k=v+1}^{f(i)}\lfloor \frac{k-1}{v}\rfloor \cdot y_{ij}^k+\sum _{i\in {\mathbf{N}}_{\mathbf{w}}}\sum _{k=v}^{f(i)}{y}_{ji}^k\le 0$$

$$\forall v=\left\{2,\cdots ,U-1\right\}\wedge j\in {\mathbf{N}}_{\mathbf{w}}$$

$$x_{ij}\in \left\{0,1\right\}{y}_{ij}^k\in \left\{0,1\right\}$$

$$\forall (i,j)\in {\mathbf{A}}_{\mathbf{r}}\wedge k\in \left\{1,\cdots ,f(i)\right\}$$

$$0\le {\sigma }_i\le \eta \cdot \lceil \frac{n_w}{n_o}\rceil {\sigma }_i\in {\mathbb{Z}}_{+}\forall i\in {\mathbf{N}}_{\mathbf{o}}$$

Full text

Nomenclature

Acronyms

OWF(s)

Offshore Wind Farm(s).

WT(s)

Wind Turbines(s).

OSS(s)

Offshore Substation(s).

BIP

Binary Integer Programming.

MILP

Mixed Integer Linear Programming.

MIQP

Mixed Integer Quadratic Programming.

MINLP

Mixed Integer Nonlinear Programming.

SCETM

Single-Core Equivalent Thermal Model.

NP

Non-Deterministic Polynomial.

L

Length.

LP

Length plus total Power losses.

I

Initial investment.

IP

Investment plus total Power losses.

WDS

West of Duddon Sands.

TH

Thanet.

LA

London Array.

HR1

Horns Rev 1.

O

Ormonde.

DT

DanTysk.

Parameters (non-sets)

$n_{w}$

Number of wind turbines.

$n_{o}$

Number of offshore substations.

$m$

Total number of years.

$a_{ij}$

Arc connecting point $i$ to $j$ $(i,j)$.

$d_{ij}$

Euclidean norm for arc $(i,j)$.

$u_{t}$

Capacity of cable $t$ in number of wind turbines.

$c_{c_{t}}$

Metric capital cost of cable $t$.

$c_{p_{t}}$

Metric installation cost of cable $t$.

$I_{t}$

Steady-state rated current of cable $t$.

$V_{n}$

Nominal line-to-line voltage of the system.

$P_{n}$

Nominal power of the wind turbines.

$S_{r_{t}}$

Nominal power of cable $t$.

$\vec{\gamma}_{t}$

Propagation constant of cable $t$.

${\vec{Z_{c}}}_{t}$

Characteristic impedance of cable $t$.

$\vec{z}_{t}$

Metric series impedance of cable $t$.

$\vec{y}_{t}$

Metric admittance of cable $t$.

$\vec{I}_{ij,t}^{\,k}$

Nominal phasor current of arc $(i,j)$ using cable

$t$, when $k$ wind turbines are connected.

$S_{ij,t}^{k}$

Nominal power of arc $(i,j)$ using cable $t$, when

$k$ wind turbines are connected.

$p^{\omega}$

Power produced by a wind turbine at $\omega$ hour-slot.

$f^{\omega,k}_{ij}$

Power flow in arc $(i,j)$ at $\omega$ hour-slot,

when $k$ wind turbines are connected.

$\vec{I}_{ij,t}^{\,\omega,k}$

Phasor current through arc $(i,j)$ using cable $t$, at

$\omega$ hour-slot, when $k$ wind turbines are connected.

$\lambda_{1}$

Screen losses factor.

$\lambda_{2}$

Armouring losses factor.

$W_{d_{t}}$

Metric dielectric loss of cable $t$.

$R_{t}$

Metric electrical resistance of cable $t$.

$l^{\mu,k}_{ij,t}$

Annual total power losses through arc $(i,j)$ using

cable $t$, at year $\mu$, when $k$ wind turbines are co-

nnected.

$U$

Capacity of the biggest cable available given as

maximum number of supportable wind turbines.

$r$

Discount rate.

$c_{ij}^{k}$

Metric cost of arc $(i,j)$, when $k$ wind turbines

are connected.

$\phi$

Maximum number of feeders per offshore substa-

tion.

$\eta$

Loading symmetry coefficient for offshore subs-

tations.

$\upsilon$

Number of wind turbines arcs set to a wind turbi-

ne.

$\upsilon_{f}$

Number of wind turbines arcs set to a wind turbi-

ne for the feasibility problem.

$\upsilon_{f_{min}}$

Algorithm parameters for the feasibility problem.

$\upsilon_{f_{\delta}},\upsilon_{f_{max}}$

$\upsilon_{o}$

Number of wind turbines arcs set to a wind turbi-

ne for the global optimization problem.

$\upsilon_{o_{min}}$

Algorithm parameters for the global opt. problem.

$\upsilon_{o_{\delta}},\upsilon_{o_{max}}$

$\epsilon$

Required relative optimality gap.

Parameters (sets)

$\bm{N_{o}}$

Set of offshore substations.

$\bm{N_{w}}$

Set of wind turbines.

$\bm{N}$

Set of offshore substations and wind turbines.

$G$

Weighted directed graph.

$\bm{A}$

Set of available arcs.

$\bm{D}$

Set of arcs weights.

$\bm{T}$

Set of available cables.

$\bm{U}$

Set of cables capacities in wind turbines number.

$\bm{C_{c}}$

Set of cables capital expenditures costs.

$\bm{C_{p}}$

Set of cables installation costs.

$\bm{\Omega}^{\mu}$

Set of hours-slot for a year $\mu$.

$\bm{M}$

Set of operational years.

$G_{r}$

First reduced graph.

$\bm{A_{r}}$

Set of first reduced arcs.

$\bm{\chi}$

Set of crossing pairs arcs.

$\bm{\Upsilon}_{i}$

Set of wind turbines connected to $i$.

$G^{{}^{\prime}}_{r}$

Second reduced graph.

$\bm{A^{{}^{\prime}}_{r}}$

Set of second reduced arcs.

Variables

$x_{ij}$

Binary variable to activate arc $(i,j)$.

$x_{ij,t}$

Binary variable to select optimum cable type $t$

for arc $(i,j)$.

$y^{k}_{ij}$

Binary variable to activate arc $(i,j)$, when $k$ wind

turbines are connected.

$\sigma_{i}$

Integer variable of number of wind turbines co-

nnected to offshore substation $i$.

Optimization output (non-sets)

$\epsilon_{k_{f}}$

Calculated gap at iteration $k_{f}$.

$\epsilon_{k_{o}}$

Calculated gap at iteration $k_{o}$.

$\epsilon_{k_{g}}$

Recalculated gap at iteration $k_{g}$.

Optimization output (sets)

$\bm{I}$

Set of first feasible solution.

$\bm{O}_{k_{o}}$

Set of feasible solution at iteration $k_{o}$.

$\bm{\Gamma}_{k_{o}}$

Set of candidate arcs at iteration $k_{o}$.

$\bm{Z}_{k_{o}}$

Set of active variables $x_{ij}=1$ of the problem

defined in the iteration $k_{o}$.

Subscripts

$i$

Element in the set $\bm{N}$.

$j$

Element in the set $\bm{N}$.

$ij$

Arc $(i,j)$ with tail at $i$ and head at $j$.

$k_{f}$

Iteration for the feasibility problem.

$k_{o}$

Iteration for the global optimization problem.

$k_{g}$

Iteration for the general problem.

$t$

Cable type in the set $\bm{T}$.

Superscripts

$k$

Number of turbines connected in $(i,j)$.

$\omega$

Hour-slot $\in$ $\bm{\Omega}^{\mu}$.

$\mu$

Year $\in$ $\bm{M}$.

I Introduction

Offshore wind energy represents a backbone technology towards the transition to power systems fully based on renewable energy.

After the invention, and experimentation age during the 1990s, the commercialization and development period is fundamentally focused on turning this technology into not only an environmentally sustainable paradigm, but also financially competitive compared to other classic, and emergent types of energy generation.

The share of Offshore Wind Farms (OWFs) has increased almost five times in the last seven years [1], reaching a globally installed power of nearly $19\text{\,}\mathrm{G}\mathrm{W}$. OWF projects are capital intensive, having large values of operating leverage where the required electrical infrastructure costs can raise up to 15% compared to the total system costs [2].

The electrical collection system is the set of electrical infrastructure components (AC submarine cables, switchgears, transformers, protection, and control units, etc.). This is required to interconnect the Wind Turbines (WTs) with each other and the Offshore Substation (OSS), guaranteeing an effective, reliable, and efficient collection of energy to the export infrastructure.

Between 2018 and 2028 more than $19,000\text{\,}\mathrm{km}$ of cables for collection systems are prognosed to be installed in UK only, with an estimated worth of £5.36bn [3]. Economies of scale pushes the development of large-scale OWFs, having more than 80 WTs while increasing their rated power.

The collection system design and optimization problem has been studied with accentuated focus in the last ten years [4, 5]. Finding the global optimum of this problem is generally NP-hard [6]. Four big clusters of methods for tackling this problem can be established: heuristics, metaheuristics, global optimization with mathematical formulations, and hybrids, such as matheuristics.

Global optimization encompasses a large set of different alternatives to model the cable layout problem, like Binary Integer Programming (BIP) [7], Mixed Integer Linear Programming (MILP) [8, 9, 10, 11, 12], MILP with decomposition techniques for stochastic programming [13, 14], Mixed Integer Quadratic Programming (MIQP) [15, 16], and Mixed Integer Non-Linear Programming (MINLP) [17, 18].

Important advances on mathematical modelling are provided in [7], but disregarding fundamental practical considerations for OWFs, such as cables crossings, multiple OSSs, maximum number of feeders per OSS, wind power stochasticity, among others. The application focuses on global optimization for medium-scale OWFs (with the largest problem instance of $57$ units). Similarly, [11] and [12] are tackling small-scale projects (with less than or equal to $30$ WTs), without proposing any strategies for scaling the applications to large-scale problems and, as in [7], do not include practical considerations in their modelling.

Larger instances are designed to optimality in [8], however, following a Planar Open Vehicle Routing Problem approach (no branching), restricting the cables set size, and ignoring power losses. Large-scale OWFs are tackled in [9] and [10], combining a MILP flow-based model with up to four heuristics, considering power losses, and other practical applications.

The works [14, 13], and [16] provide remarkable advances on stochastic optimization for problems in this context. Different stochastic scenarios are supported, accounting for wind power variability, and cables failure. Distinctive theoretical strategies to accelerate convergence are applied and compared. Nonetheless, some simplifications are incorporated, such as iterative processes to calculate total electrical power losses, combined with experts analysis to estimate their impact on the cable layout. Case studies are limited to small-scale applications.

By means of explicit formulation of electrical losses economic costs in the model, a MINLP program is proposed in [18], where, additionally, clustering algorithms are used prior the execution of the program into a commercial solver. Finally, losses can be also computed in the accurate quadratic form as in [15], but at the expense of a decrease in computational efficiency.

Each of these mathematical formulations impose certain limitations about the physics modelling options. For instance, using flow-based MILP makes it more difficult to include the quadratic active power losses explicitly into the objective function. The commonly used power flow equations solved with e.g. the Newton-Raphson method cannot be considered in MILP or MIQP formulations.

The side effects of more flexible modelling formulations are the compromise of solver functionality and performance. The proper balance between solution method and complexity on modelling, represents one of the main challenges for the OWF developer, and trade-offs have to be adopted within certain assumptions. As a generalization, linear-based formulations are computationally more efficient than quadratic or non-linear.

Since the trend in OWFs is to deploy large-scale projects, focus is directed into this aspect. To the best of the authors’ knowledge, only [8], [9], [18], and [19] have tackled collection systems for large-scale OWFs with $80$ - $100$ WTs using global optimization. Only one cable type and no-branching at WTs nodes is considered by Bauer et. al. [8], where heuristics have also been proposed. A matheuristic framework is developed by Fischetti et. al. [9], including the total electrical power losses in the objective function using pre-processing strategies, and flow formulation. A MINLP flow formulation is used in [18], having sets of linear constraints but with explicit linearized losses inclusion in the objective; a fixed neighboring search is implemented, and pre-clustering of WTs to OSS is required, limiting artificially the search space, leading to sub-optimal points. No losses have been included in the optimization in [19], and for large-scale problems, computational experiments point out that feasible solutions may not be found. Additionally, only a single cable is supported.

The main contributions of this work are: (i) development, testing, and application of a mathematical model to quickly find feasible points for large-scale OWF instances (with more than $100$ WTs). This is possible by proposing a MILP model with reduced constraints and variables than previously used flow formulations for collection systems in OWFs. (ii) combination of an algorithmic framework with mathematical formulation for obtaining global optimum solution points (or near to it) in reasonable computational time. Additionally, improvements on the complexity and fidelity of modelling aspects for integration with a linear formulation, such as total power losses calculated using time series, and capacitive currents, have been considered. The method does not required pre-clustering of WTs to OSSs, as this is tackled intrinsically in the model.

In Section II, the modelling aspects are explained in details; followed by Section III, where the optimization model is formulated including objective function and constraints definition. The whole framework description is presented in Section IV, where the full mathematical formulation is compacted into a single main iterative algorithm. Computational experiments are performed in Section V, and the work is finalized with the conclusions in Section VI.

II Modelling aspects

II-A Problem graph representation

The aim of the optimization is to design the cable layout of the collection system for OWFs, i.e., to interconnect the $n_{w}$ WTs to the available OSSs, $n_{o}$.

In graph theory, a OWF is represented by a forest with roots at the OSSs. A few OWF developers opt for systems with loops. The inclusion of loops may increase reliability, but raises other technical challenges as cables sizing, and generally these systems are designed following heuristic and experts rules. Radial layouts are the common practice in the industry, since array cables failures are rather rare as they are buried in the seabed [20]. Therefore, the desired result is $n_{o}$ spanning-trees, which minimize the required objective function, while satisfying the operational and topological constraints.

Let the OSSs define the set $\bm{N_{o}}=\left\{1,\cdots,n_{o}\right\}$. Likewise, for the WTs, let $\bm{N_{w}}=\left\{n_{o}+1,\cdots,n_{o}+n_{w}\right\}$. In this way, each one of the OSSs and WTs (modelled as points in the space) have associated a unique identifier $i$, such as $i\in\bm{N}=\bm{N_{o}}\cup\bm{N_{w}}$. The Euclidean norm between the positions of the points $i$ and $j$, is defined as $d_{ij}$. The aforementioned inputs are condensed as a weighted directed graph $G(\bm{N},\bm{A},\bm{D})$, where $\bm{N}$ represents the vertex set, $\bm{A}$ the set of available arcs arranged as a pair-set, and $\bm{D}$ the set of associated weights for each element $a_{ij}\in\bm{A}$, where $i\in\bm{N}$ $\wedge$ $j\in\bm{N}$. For instance, for $a_{ij}=(i,j)$, $d=d_{ij}$, where $d\in\bm{D}$. In general, $G(\bm{N},\bm{A},\bm{D})$ is a complete directed graph.

Additionally, a predefined list of available cables types is required to interconnect the WTs towards the OSSs. Let the set of cables be $\bm{T}$ and let the capacity of a cable $t\in\bm{T}$ be $u_{t}$ measured in terms of number of supportable WTs connected downstream. Hence, let $\bm{U}$ be the set of capacities sorted as in $\bm{T}$ (see the definition in Section II-B).

Furthermore, each cable type $t$ has a cost per unit of length, $c_{c_{t}}$, in such a way that $u_{t}$ and $c_{c_{t}}$ describe a positive correlation, following an exponential regression model. The set of metric capital expenditures is defined as $\bm{C_{c}}$. Similarly the set of metric installation costs is defined by $\bm{C_{p}}$.

After defining the graph representation of the problem, the underlying variables associated to the calculation of the desired output, are unequivocally established. Let $x_{ij}$ represent a binary variable that is one if the arc between the vertex $i$ and $j$ is selected in the solution, and zero otherwise. Likewise, the binary variable $y^{k}_{ij}$ models the $k$ number of WTs connected downstream from $j$, including the wind turbine at node $j$ (under the condition that $x_{ij}=1$). Finally, the integer variable $\sigma_{i}$ represents the number of WTs connected to the OSS $i$.

II-B Cable capacity

The current capacity $I_{t}$ of a cable $t$ is calculated using the model given in [21]. This method represents the most common industrial practice, and it is based on a Single-Core Equivalent Thermal Model (SCETM) as generalized in [22] for single-core and three-core cables.

[TABLE]

[TABLE]

The set of cable capacities in terms of number of supportable WTs is defined in (1) and (2), where $P_{n}$ represents the nominal power of an individual WT, and $V_{n}$ the line-to-line nominal voltage of the system.

The cables capacity must be assessed in terms of calculated maximum conductor temperature, and the maximum allowable value by design (usually considered $90\text{\,}\mathrm{\SIUnitSymbolCelsius}$). In this manuscript, the conservative approach described in [21] is followed. This model neglects the slow thermal time constant of submarine cables, and the high variability of offshore wind power, since it assumes nominal steady state conditions. On the other hand, this practical assumption brings along as a result a robust design, while avoiding power curtailment at any instant.

II-C Arcs nominal power

The power flow model of a transmission line solving the resultant differential equations considering parameters distributed uniformly throughout the length of the cable [23], is implemented with the expressions:

[TABLE]

[TABLE]

The model (3) provides the highest complexity on physics modelling possible to apply given the flexibility of the proposed mathematical formulation (see Section III). This provides more accuracy for securely determining the cable type in an active arc $x_{ij}$. The characteristic impedance is calculated as ${\vec{Z_{c}}}_{t}=\sqrt{\vec{z}_{t}}/\vec{y}_{t}$, and the propagation constant, $\vec{\gamma}_{t}=\sqrt{\vec{z}_{t}}\cdot\vec{y}_{t}$. The series impedance is represented by $\vec{z}_{t}$, and the admittance by $\vec{y}_{t}$. The maximum power flowing through the arc $(i,j)$, given $k$ turbines connected downstream ($y^{k}_{ij}$), is calculated as per (4), accounting for the worst-case scenario, as the current in the arc is strictly increasing with length, and the value $\vec{I}_{ij,t}^{\,k}$ is calculated at the extreme of it when using cable $t$.

II-D Power flow and total power losses

A transportation model, accurate enough for radial systems, is implemented following [4]. As explained in Section II-B, the set $\bm{U}$ assumes nominal steady state conditions, describing a practical conservative criterion.

[TABLE]

The transportation model is generalized considering the temporal dimension (5), where $\bm{M}=\left\{1,\cdots,m\right\}$, is the set of operational years with upper limit in the project lifetime, $\bm{\Omega}^{\mu}$ is the set of hours-slot for a year $\mu$, and $\omega$ a specific hour-slot in $\bm{\Omega}^{\mu}$. Let $p^{w}$ be the power in $\mathrm{M}\mathrm{W}$ produced by one WT in that hour-slot, simulated with [24].

In this way, let define the auxiliary variable $f^{\omega,k}_{ij}=k\cdot y^{k}_{ij}\cdot p^{\omega}$ as the power flow ($\mathrm{M}\mathrm{W}$) in arc $(i,j)$ when $k$ WTs are connected downstream (including the one in $j$), in time instant $\omega$. This ignores wake losses, a reasonable assumption since micrositing optimization techniques are considered already applied [25].

[TABLE]

[TABLE]

Including the capacitive currents, (6) expresses the current at the end of the arc (with respect to $i$), with magnitude $|\vec{I}_{ij,t}^{\,\omega,k}|$. The annual total power losses $l^{\mu,k}_{ij,t}$ are calculated with (7). The factor $(1+\lambda_{1}+\lambda_{2})$ accounts for the screen and armouring losses and $W_{d_{t}}$ is the dielectric loss per unit length for the insulation surrounding the conductor in $\mathrm{W}\mathrm{/}\mathrm{m}$[21], while the constant 3 is for the three-phase system. This value must be scaled in $\mathrm{M}\mathrm{W}\mathrm{h}$.

III Optimization Model

The proposed optimization model described in this Section is able to cope with an arbitrary number of WTs, $n_{w}$, and similarly any reasonable number of OSSs, $n_{o}$.

The underlying mathematical formulation is inspired by the formulations and analysis proposed in [7] and [26] with additional constraints stemming from the nature of the problem, and with the objective to improve its tractability.

In Section II-A the binary variables $x_{ij}$ and $y^{k}_{ij}$, and the integer variable $\sigma_{i}$ are defined. They refer to an active arc, the number $k$ of WTs connect to that active arc, and the number of WTs connected to OSS $i$, respectively.

To increase the computational efficiency, the number of variables is reduced as follows. The capacity of the biggest cable is calculated as $U=\max{\bm{U}}$, therefore the possible maximum value of $k$ for $i\in\bm{N_{o}}$ is equal to $f(i)=U$, while for $i\in\bm{N_{w}}$ is $f(i)=U-1$. This acknowledges that the biggest cable available could be only used at maximum capacity when is connected from a OSS.

Analogously, the set of variables $x_{ij}$, where $i\in\bm{N_{w}}$, and $j\in\bm{N_{o}}$ are intrinsically discarded, considering the nature of the power flow, i.e., the OSSs collects the energy from the WTs and not the other way around. Lastly, since the export system is outside the scope of this article, all the arcs between OSSs are disregarded, i.e., $x_{ij}=0$ $\forall i\in\bm{N_{o}}$ $\wedge$ $j\in\bm{N_{o}}$.

The graph $G(\bm{N},\bm{A},\bm{D})$ is reduced to $G_{r}(\bm{N},\bm{A_{r}},\bm{D_{r}})$ after this stage. While including the OSSs in the optimization could offer some flexibility, their locations are typically decided already in the development process. Furthermore, their locations are strongly driven by the distance to the onshore connection points, hence is deemed as plausible assumption to consider fixed location of the OSSs.

III-A Cost coefficients

Note that the previously defined decision variables $x_{ij}$ and $y^{k}_{ij}$, do not include any information related to the cable type selected in a given arc.

This is because the cable type selection process is handled in a pre-processing stage, given that all the required data is present, and the task is totally independent to any other part of the desired tree(s) [7].

A pre-processing strategy allows integrating more complex power flow and total electrical power losses models, increasing the accuracy without compromising the computational efficiency. Along these lines, the distributed model explained in Section II-C is harmonized with the main optimization problem.

At the same time, the aforementioned point allows for a power flow estimation in a conservative fashion, i.e., overestimating the incoming power flow by neglecting the total power losses downstream.

In the pre-processing state, for the case of $y^{k}_{ij}$, the length of the arc is known ($d_{ij}$), and the number of WTs connected by it is also defined (by $k$). No more inputs are required for this task.

Hence, for each $y^{k}_{ij}$, the sub-problem defined by (8) to (11) is solved beforehand and independently by enumeration.

[TABLE]

In the objective function (8), the first term ($(c_{c_{t}}+c_{p_{t}})\cdot d_{ij}$) is for capital expenditure ($c_{c_{t}}$) and installation costs ($c_{p_{t}}$) of cable $t$.

The values for $c_{c_{t}}\in\bm{C_{c}}$ are obtained from the exponential regression function, given as

[TABLE]

where $a_{p_{t}}$, $b_{p_{t}}$, and $c_{p_{t}}$ are coefficients dependent on the nominal voltage of cable type $t\in{\bm{T}}$, $S_{r_{t}}$ is the rated power of $t$ in $\mathrm{V}\mathrm{A}$ (also depending on the rated line to line voltage level, $V_{n}$, see the definition in Section II-B). The cables capital expenses cost function (12) is extracted from [27], being based on a comprehensive semi-empirical cost survey, applying data fitting techniques.

The second term in the summation part of (8) accounts for the discounted cash flow of the economic losses caused by the energy dissipation in the cables, the parameters $m$, $l^{\mu,k}_{ij,t}$ (see Section II-D), ${c}_{e}$, and $r$, represents the project lifetime (years), total power losses at year $\mu$ for cable $t$ ($\mathrm{M}\mathrm{W}\mathrm{h}$), cost of energy (€ /$\mathrm{M}\mathrm{W}\mathrm{h}$), and discount rate (p.u.) respectively.

The objective function can be simplified by zeroing any of its terms, such as, if only the total length is minimized (L), the first term in (8) has to be replaced uniquely by $d_{ij}$, while the other term in the same equation is dropped. Similarly, if only the total initial investment (I) is targeted, only the first is kept.

Therefore, the set of single objectives available in the model are: length (L), length plus total power losses (LP)—monetizing lengths by assuming the same capital and installation costs for all cables types, initial investment (I), and initial investment plus total power losses (IP).

Likewise, (9) ensures that exactly one cable type is selected, while (10) guarantees that the capacity of cable $t$ is not exceeded; $S_{ij,t}^{k}$ is the power through arc $(i,j)$ when $k$ turbines are connected in $j$ using cable $t$ defined in Section II-C. The binary variable for selecting a cable type $t$ for arc $(i,j)$ is defined in (11).

The sub-problem from (8) to (11) seeks to find the cable type $t$ to be used for the arc $(i,j)$, which minimizes the objective (8).

III-B Objective function

After solving the multiple sub-problems related to cable selection and cost evaluation (maximum $U\cdot|\bm{N}|^{2}$ problems) from (8), a cost value $c_{ij}^{k}$ is associated to each $y_{ij}^{k}$ variable. The linear objective function of the main mathematical model is then

[TABLE]

III-C Constraints

In order to present the solution connecting all WTs between each other and to the OSSs, the following constraint is added

[TABLE]

Constraint (14) models the full OWF to be divided into multiple disconnected trees (forest) with $\sigma_{i}$ being the number of WTs associated to a OSS $i$. Hence, the total amount of WTs ($n_{w}$) are integrated into the electrical system.

To guarantee full connectivity in OSS $i$, the next constraint is added

[TABLE]

Note that (14) and (15) are combined in the case of only one OSS.

To limit the maximum number of feeders per OSS ($\phi$), it is used:

[TABLE]

To simultaneously ensure a tree topology, ensure that only one cable type used per arc, and to define the head-tail convention, the next expression is included into the model

[TABLE]

The flow conservation, which also avoids disconnected solutions, is considered by means of one linear equality per wind turbine

[TABLE]

The set $\bm{\chi}$ stores pairs of arcs $\left\{(i,j),(u,v)\right\}$, which are crossing each other. Excluding crossing arcs in the solution is ensured by the simultaneous application of the following linear inequalities

[TABLE]

[TABLE]

The no-crossing cables restriction is a practical requirement in order to avoid hot-spots, and potential single-points of failure caused by overlapping cables [8]. Constraint (19) exhaustively lists all combinations of crossings arcs, including also the corresponding inverse elements. The constraints in (20) ensure that no active arcs are crossing or overlapping between each other. These constraints thus link the variables $y^{k}_{ij}$ and $x_{ij}$.

Cables crossings are detected based on a procedure of slopes evaluation. Two arcs are crossing if the crossing point is inside of the lines, but not if this point is located at the extremes of the lines or beyond in the lines’ projections.

[TABLE]

Constraint (21) represents a set of valid inequalities, initially proposed in [7], to tighten the mathematical model. Given an active arc $y^{k}_{ij}$, the maximum number of active arcs rooted in $j$ and connecting $v$ WTs, is expressed by $\left\lfloor\frac{k-1}{v}\right\rfloor$, hence the constraint restricts the maximum number feasible arcs, reducing the search space without excluding valid solutions to the problem.

[TABLE]

[TABLE]

Constraints (22) and (23) define the nature of the formulation by the variables definition, a MILP.

Note that variables $\sigma_{i}$ are limited in their upper bounds to avoid uneven loading of OSSs (in case $\eta=1$, otherwise $1<\eta\leq n_{o}$). Equally rated OSSs bring benefits like design standardization, and decreasing of the dependency upon a single transformation unit for transporting the generated power.

To summarize, the complete formulation of the main MILP model consists of the objective function (13) and the constraints defined in (14) - (23).

The base formulation presented so far has a maximum number of binary variables equal to $|\bm{N}|^{2}+U\cdot|\bm{N}|^{2}$, integer variables number equal to $|\bm{N_{o}}|$ (linear in function of $n_{o}$), and constraints (excluding the crossing constraints and valid inequalities) of $1+2\cdot|\bm{N_{o}}|+2\cdot|\bm{N_{w}}|$. Flow formulations, such as the one proposed in [9], have more variables ($2\cdot|\bm{N}|^{2}+U\cdot|\bm{N}|^{2}$) and constraints ($|\bm{N}|^{2}+2\cdot|\bm{N_{w}}|+|\bm{N_{o}}|$); integer and binary variables are quadratic in function of the problem size. This fact along with the addition of valid inequalities may explain why the model from (13) to (23) is often more efficient to solve.

Further simplifications to the model are presented in Section IV-A.

IV Optimization framework

IV-A Candidate arcs

The reduced graph $G_{r}$ is obtained after performing the described considerations in the introduction of Section III.

However, given the NP-Hard nature of this problem, which is similar to a Capacitated Minimum Spanning Tree with additional constraints [28, 11], more reductions are required.

The limitations for successfully finding feasible points and high quality solutions, using solely mathematical models and commercial solvers, is demonstrated in [9].

For large-scale OWFs (with more than $100$ WTs) the computing time for robust global optimization solvers generally becomes notoriously long. Likewise, in general, solution times become unpredictable, while very large memory requirements are demanded to build the branch-and-cut tree. Besides, the constraints generation must be done with special care (the full set of crossing constraints has a combinatorial nature) to increase computational efficiency.

To make the formulation more flexible and implementable, a further operation to the graph $G_{r}$ is proposed. The function $f(i,G_{r},\upsilon)$ calculates the set $\bm{\Upsilon}_{i}$, defined as the $\upsilon$-closest WTs to $i$. In other words, it is intuitively considered that a WT will be connected to one of the WTs in its vicinity. Therefore, by systematically applying $f(i,G_{r},\upsilon)$ to each $i\in\bm{N_{w}}$, the reduced graph $G^{{}^{\prime}}_{r}$ is found. The set $\bm{A^{{}^{\prime}}_{r}}$ contains the candidate arcs to the solution of the problem.

With this strategy, the maximum number of variables is reduced to $|\bm{N_{o}}|+(U+1)\cdot|\bm{N_{o}}|\cdot|\bm{N_{w}}|+U\cdot\upsilon\cdot|\bm{N_{w}}|$. Additionally, the number of crossing constraints decreases dramatically as well.

Overall, the arcs set transformation follows $\bm{A}\rightarrow\bm{A_{r}}\rightarrow\bm{A^{{}^{\prime}}_{r}}$.

All $(ij)$ indexed variables and constraints, presented in the Section III, must be adapted to this reduction process.

IV-B The Algorithm

The main algorithm defining the full framework with the mathematical model is presented in Algorithm 1.

From line 1 to 12 the task is to efficiently solve a feasibility problem. The idea is to subsequently increase $\upsilon$ from an initial value $\upsilon=\upsilon_{f}=\upsilon_{f_{min}}$ to a maximum value $\upsilon=\upsilon_{f}=\upsilon_{f_{max}}$, with steps $\upsilon_{f_{\delta}}$, until a feasible point is found. If this is achieved in iteration $k_{f}$, the first task is terminated with a feasible point $\bm{I}$. Conversely, if the model is infeasible, the candidate arcs set is augmented with $\upsilon_{f_{\delta}}$ units, and the process is taken to the iteration $k_{f}+1$, where a new trial is attempted.

In order to formulate the MILP model, the cost coefficients calculation from (8) to (11) is omitted by setting them equal to zero, and the black-box MILP solver terminates when the first feasible point is found.

At this point, the Algorithm 1 requires as parameters $\upsilon_{f_{min}}$, $\upsilon_{f_{\delta}}$, and $\upsilon_{f_{max}}$. The greater $\upsilon_{f_{min}}$ and $\upsilon_{f_{\delta}}$ the less efficient the feasibility problem becomes, however, increasing the odds to defining a feasible instance of the problem promptly.

Likewise, from line 13 to 35, the global optimization task is performed. The target is to obtain a feasible point with a given relative gap $\epsilon$, expressed as the relative difference of the best feasible point ($\tau$) minus the best achievable value objective ($\kappa$), with respect to $\tau$. These values are indexed by iteration number. Similarly to the feasibility task, the iterative process increases the candidate arcs set from $\upsilon=\upsilon_{o}=\upsilon_{o_{min}}$ to $\upsilon=\upsilon_{o}=\upsilon_{o_{max}}$, with steps $\upsilon_{o_{\delta}}$.

The termination criterion is when the set $\bm{Z}_{k_{o}}$ of active variables $x_{ij}=1$ of the problem defined in the iteration $k_{o}$, is a subset of the arcs set $\bm{A^{{}^{\prime}}_{r}}$ defined in the previous iteration $k_{o}-1$ ($\bm{\Gamma}_{k_{o}-1}$).

In this way, it is inferred that it is not longer necessary to increase $\upsilon_{o}$, as the optimum variables have been already provided in the previous iteration.

To guarantee along the process a monotonously decreasing value of the objective function, in iteration $k_{o}$, the mathematical model is warm-started with the feasible solution found in $k_{o}-1$ ($\bm{O}_{k_{o}-1}$). This strategy may help in shortening the convergence time for the sub-instance $k_{o}$.

Conceptually, Algorithm 1 intends to determine a reduced search space, where the global minimum point is hopefully included. If only one reduced problem was solved given a $\upsilon$, it would not be possible to infer about the quality of the solution, and the calculated gap for that particular instance could not represent the global domain of the full problem, potentially leading to an overestimation.

For the global optimization task, Algorithm 1 requires as parameters $\upsilon_{o_{min}}$, $\upsilon_{o_{\delta}}$, and $\upsilon_{o_{max}}$ for the global optimization task. Naturally, $\upsilon_{o_{min}}\geq\upsilon_{f_{max}}$, and it is reasonable to consider $\upsilon_{o_{\delta}}>\upsilon_{f_{\delta}}$. By proper adjustment of the previous parameter, in best case scenario, the full Algorithm is concluded for $k_{f}=1$ and $k_{o}=2$.

Although for every iteration the maximum required gap $\epsilon$ is equally fixed, the equivalent calculated gap, having as reference the full-size domain, varies. Larger values of $\upsilon_{k_{o}}$ lead to equal or lower values of $\kappa_{k_{o}}$. This causes that in general, $\tau_{k_{o}}$ is also lower, until the ideal reduced search space is found, when equal values of $\tau_{k_{o}}$ should be obtained.

Therefore, after the termination of the algorithm, a gap updating procedure is performed based on the last calculated value of $\kappa_{k_{o}}$, to recalculate the relative difference for all previous iterations respect to this value (line 36). Let the recalculated gap in the global iteration $k_{g}$, including the feasibility and global optimization problems, be $\epsilon_{k_{g}}$. In this sense, an evolution of the gap in function of the iterations is available, providing further insights and the sense of convergence, as the objective value decreases monotonically.

V Computational Experiments

The following experiments have been carried out on an Intel Core i7-6600U CPU running at 2.50 GHz and with 16 GB of RAM. The chosen MILP solver is the branch-and-cut solver implemented in IBM ILOG CPLEX Optimization Studio V12.7.1 [29].

In the Section V-A a sensitivity analysis of the parameters $\upsilon_{o_{min}}$, and $\upsilon_{o_{\delta}}$ is presented. This is achieved by applying the proposed methodology of Section IV, on two OWFs using several sets of parameters.

In Section V-B, the proposed method is benchmarked against the results obtained through a different approach published in the scientific literature [9]. With this aim, the same testbed is employed, while assuming the same considerations, such as, objective (8)-(11), and constraints embodied by the equations (14)-(23).

V-A Algorithm’s Parameters Sensitivity Analysis

The two real OWFs West of Duddon Sands (WDS) and Thanet (TH) are used for the parameter sensitivity analysis. The information regarding these OWFs is provided in [30].

Table I displays the main parameters for the sensitivity analysis. In the case of WDS, the objective function is defined by a combined total economic cost including the initial investment, and the total electrical power losses of the cable layout (IP). Whilst for TH, the target is to minimize the initial investment of the cable layout, defined as the cables capital expenditures, and installation costs.

The difference between objective functions is not relevant for the purposes of the sensitivity analysis, but rather, to show the capability of the method to support both cases. Economic parameters for the discounted cash flows calculations are also presented.

Both instances are challenging to solve given their large size (with more than $100$ WTs each), the large number of cables sizes considered (up to three cables types), and maximum capacity $U=\max{\bm{U}}$, spanning from $10$ to $13$ (see Table I). High-level information such as WT power and number ($P_{n}$ and $n_{w}$), OSSs number ($n_{o}$), and the maximum number of allowed feeders connected to the OSSs ($\phi$) is available. The limit $\phi$ provides a hard binding constraint. It is fixed to a practical value usually considered by OWFs developers.

Finally, the rest of parameters for the optimization are displayed. The collection systems voltage level ($V_{n}$) is the traditionally used, and two different set of cables are considered. The electrical and thermal information of the cables is available in [31].

The capacity set $\bm{U}$ is determined for each instance, similarly, for the total cost per kilometer, including capital expenditures ($\bm{C_{c}}$), and installation costs $\bm{C_{p}}$. The cables $\left\{1,2,3\right\}$ have cross-sections of $240\text{\,}\mathrm{m}\mathrm{m}^{2}$, $500\text{\,}\mathrm{m}\mathrm{m}^{2}$, and $1000\text{\,}\mathrm{m}\mathrm{m}^{2}$, respectively.

The set represented by $\left\{4,5\right\}$ is used for TH, and it has associated electrical and economic parameters matching those considered in the benchmark work [9], while neglecting the capacitive currents (see Section V-B for benchmark analysis).

The whole framework compacted in Algorithm 1 intends to find an approximation of the minimum search space of a given instance. It is presumed the finding the model best feasible point (or near to it) for a given required maximum gap, having as reference the full (global) problem size.

Each problem instance is composed by $k_{g}=k_{f}+k_{o}$ sub-instances (or general iterations), which in turn include those iterations required for solving the feasibility ($k_{f}$), and global optimization serial problems ($k_{o}$); this implies that, in the best-case scenario, a maximum number of iterations equal to $k_{f}=1$ and $k_{o}=2$ should be enough to solve any problem instance.

However, factors such as the adequate setting of the algorithm’s parameters and the specific spatial characteristics influence the number of global iterations $k_{g}$, and the circumvention of local minima.

The avoidance of local minima is crucial, as it has a lifetime impact over a project. Successful implementations depend on the chosen values of $\upsilon_{o_{min}}$, and $\upsilon_{o_{\delta}}$.

The sensitivity analysis elaborates on how the iterative reduction of the search space through Algorithm 1 can lead to sub-optimal solutions.

In pursuance of this, two values associated to the initial value parameter, $\upsilon_{o_{min}}$, are examined: $\{6,15\}$; while the increase steps parameter, $\upsilon_{o_{\delta}}$, is set between $\{2,5,10\}$. All possible combinations of both parameters are applied to WDS and TH OWFs through the Algorithm 1.

The results are graphically shown in Figure 1 presenting quality of the solution (ordinate) versus computing time (abscissa). Both are normalized to the best corresponding case. Each case is depicted in the legend with [$\upsilon_{o_{min}}$,$\upsilon_{o_{\delta}}$], i.e. [6,2] and so on.

The results show that for the WDS OWF, the solution for $\upsilon_{o_{min}}=6$, and $\upsilon_{o_{\delta}}=2$ is nearly $2\%$ more expensive than the other cases (see Figure 1a, case [6,2]). For the rest of parameters sets, the obtained solutions are, in practical terms, the same. The costs differences among them are due to the inclusion of total electrical power losses, which are sensitive to small changes of cables lengths. The computing times varies almost in a positively correlated fashion with the values of these parameters, with the fastest case being also the most expensive.

This could be explained by the particular spatial characteristics of WDS, which impact the number of minimum candidate arcs to cover the global minimum. The optimized layout shown in Figure 2a ($\upsilon_{o_{min}}=15$ and $\upsilon_{o_{\delta}}=5$) evidentiate the non-uniform distribution of WTs in the plane, with empty areas around WTs number 43 and number 44, for instance. This can be interpreted as a higher degree of freedom in the design, having the possibility to interconnect WTs located further away from each other.

One should note that the optimized layout has a connection between WTs number 72 and 43, the latter being out of the top-15 closest WTs for the former (in fact, is the 17th closest), hence, making it necessary to increase $\upsilon_{o}$ to $20$, to include this connection.

Returning to the case of $\upsilon_{o_{min}}=6$, and $\upsilon_{o_{\delta}}=2$, the step increase value is not enough to provide a significantly broader search space to improve the solution quality.

The results suggest that the values of $\upsilon_{o_{\delta}}=5$, and $\upsilon_{o_{\delta}}=10$, result in an improvement of the solution in each iteration, and confirm in the last iteration the covering of minimum candidate arcs set. A value of $\upsilon_{o_{\delta}}=2$ in this case provides the same output, but evidently with potentially higher chances of falling into suboptimal points for other OWFs.

For the TH OWF, all the cases generate the same result, as shown in Figure 1b. In contrast to WDS OWF, Thanet presents a regular (grid-based) micrositing layout (see Figure 2b), hence, less degree of flexibility, which translates in an optimum point represented by WTs each connected to maximum the 6th closest WT.

In order to increase the likelihood for getting a balance between solution quality, computing time, and memory requirements, a heuristic rule of considering $\upsilon_{o_{min}}=15$, and $\upsilon_{o_{\delta}}=5$ is implemented in this manuscript.

Certainly, every single OWF should be individually analysed to come up with tailored parameters, but results point out that these settings may lead to adequate terminations in most of the problems.

Table II summarizes in-detail the results of the experiments using the proposed framework with $\upsilon_{f_{min}}=5$, $\upsilon_{f_{\delta}}=1$, $\upsilon_{f_{max}}=15$, $\upsilon_{o_{min}}=15$, $\upsilon_{o_{\delta}}=5$, and $\upsilon_{o_{max}}=50$.

For each of the sub-instances, the output CT1 includes the total processing time in order to generate the Constraints (14) to (23), and additionally, the solution of the independent sub-problems as defined in the model from (8) to (11). Note that in the case of the feasibility problems the latest procedure is circumvented by fixing the cost coefficients equal to zero.

Likewise, CT2 is the computing time to solve the main mathematical model (from (13) to (23)) for a given maximum gap, in the case of WDS $\epsilon=0.5\%$, and for TH $\epsilon=0.3\%$. An additional experiment with $\epsilon=0.5\%$ for TH has been run, with a duration of around an hour, finding a solution only $0.07\%$ more expensive than the one presented in Table II, with a calculation time approx. $50\%$ smaller. In principle, any imposed maximum gap is supported, at the potential expense of a very steep increase on computing time.

The column Obj presents the objective value of the best feasible point obtained under the explained conditions, after the termination of each sub-instance calculation. Finally, the columns $\epsilon_{k_{f}}$/$\epsilon_{k_{o}}$, and $\epsilon_{k_{g}}$ gather the calculated gap of the best feasible point, and the recalculated gap for each sub-instances, respectively, after the Algorithm 1 terminates as explained in Section IV-B.

By means of this strategy, is surmised the delimitation of a reduced search space representative of the global problem, including the global minimum point.

In the particular case of WDS, the results say that an unnecessary continuation to a fifth cycle, i.e. $k_{g}=5$, was avoided due to the warm start strategy. The objective value does not change between $k_{o}=2$ and $k_{o}=3$, this being a clear (but not mathematically conclusive) signal of the successful convergence of the method.

Furthermore, an overestimation of $\epsilon_{k_{o}}$ is seen for $v_{o}=15$ when referenced to the approximated minimum search space, as $\epsilon_{k_{g}}=2.17\%>0.5\%$. It should be noted that the first feasible point is found in almost $45$ seconds, and from there, the best feasible point is found in very reasonable computing time. A reduction of the gap from $38.62\%$ to $0.50\%$ is obtained in only $1$ h:$53$ min (one hour and $53$ minutes).

In contrast to WDS, only three global iterations are required for TH, mainly due to the more uniform distribution of WTs, as shown in the Figure 2b. For TH OWF, a feasible point is obtained in $17$ s, and in $1$ h:$47$ min the gap is improved from $33.60\%$ to $0.30\%$.

The proposed procedure can be impacted by multiple valid solutions around the required maximum gap. Nevertheless, this would cause only false continuations of the Algorithm 1, rather than affect the solution quality. The previous point is partially palliated by the likely decrease of computing time in subsequent iterations, given the warm start point provided from the previous step. If the objective value was used as criterion to stop the algorithm, a maximum ratio of objective change per subsequent iterations should be set, which would be open to different assessment criteria, and potentially could lead to false terminations.

For both WDS and TH, the values $\epsilon_{k_{g}}$ reported in Table II for the last iteration, are representative of the full problem; fact validated when solving them without reduced search space (full arcs set), provided the best feasible point available.

V-B Benchmarking

A testbed of real-world cases, presented in Table III, is employed. The OWFs names are given in the Acronyms definition. Results comparison are available in Table IV.

This testbed has been mostly extracted from [9], by selecting the most challenging instances. In [9] a significantly different approach (different mathematical formulation, for instance) from the one proposed in this article, consisting in a matheuristic model, has been designed, implemented, and tested.

The results in [9] where obtained with computational resources similar to those used in this manuscript (IX CPU X5550 running at 2.67GHz, CPLEX 12.6).

In order to provide a fair comparison, all practical and technical constraints are conceptually equivalent, while the objective function is the initial investment, as losses are computed differently. Likewise, capacitive currents have been neglected as they are not included in [9].

The comparison results are presented in Table IV (instances 1-9). A small (O), two large (HR1 and DT), and a very large OWF (TH) are studied. Each instance is defined by a OWF, and a set of cables available $\bm{T}$.

Three aspects are compared: (i) solution quality, (ii) calculated gap, and (iii) computing time; each of them are directly compared by inspecting the columns Diff.[M€ ], Diff.[%], and Diff.[min], respectively. In all metrics, positive value means better performance for the method proposed in this article.

Regarding solution quality, it can be seen that, for all instances, the obtained solutions are equal or lower than in the benchmark work. For instance nr.7, around € 550,000, representing around $1\%$ of the total cost, are saved. Particularly for this instance, the gap is improved from $7.36\%$ to $0.01\%$, while simultaneously reducing the computing time.

When the calculated gap in both methods is lower than $0.01\%$, the objective values are essentially the same. This validates that the dual values between the two models are also equivalent.

In any instance the proposed method provides equally tight or even tighter solutions. The gap values reported in [9] have been recalculated in this manuscript using the best feasible point as reference (instead of the best dual value obtained in $24$ h), to make them comparable with the proposed approach.

Finally, in almost all the instances, the computing time is shorter, with exception of instance nr.8, where a considerable difference of $622$ min is observed.

It is important to clarify that the reported computing times of the proposed method are for the whole running of Algorithm 1, including in all instances the final iteration $k_{o}=2$, necessary to confirm finding the global point; in contrast, in the benchmark work, the reported time is when the best feasible point has been found. Similarly, as mentioned before, $24$ h is the time limit to obtain the best dual value.

For all the instances, the proposed method calculates feasible points in less than $40$ s, with recalculated gap of maximum $41\%$.

Besides the benchmark cases, an extra instance - not implemented in [9] - is included to show the method’s applicability to OWFs with multiple OSSs. LA OWF is the second largest project (measured by installed power) under operation, surpassed only by Walney Extension OWF (although with less number of WTs therefore potentially easier to solve). For this instance, an initial feasible point is obtained in only $3$ min, and the best feasible point is calculated in $12$ h:$28$ min, with the gap being improved by $46.87\%$.

VI Conclusion

The proposed method provides a global optimization approach to solve the cable layout of OWFs collection systems.

The main novelties of this manuscript are: (i) proposition of a model, able to provide very good solutions in very reasonable computing times, (ii) possibility to provide very tight solution quality certificates, (iii) integration of realistic and high-fidelity models to calculate total electrical power losses, and capacitive currents in the collection systems.

An algorithmic framework for reducing the search space iteratively is the main technique used. The objective function supports the total economic costs, including initial investment, and lost revenues due to total electrical power losses in the project lifetime.

The proposed methodology has been benchmarked against a state-of-the art method with significantly different approach (different MILP model, and application of up to four distinctive heuristics); with all practical and technical constraints conceptually equivalent.

Ten real-world problem instances have been considered in the benchmark. The numerical results indicate that (i) the proposed algorithm provides, in general, at least equally good solutions, and in some cases, sizeable cheaper ones than the benchmark work, (ii) tighter gaps are calculated, in shorter computing times.

The proposed algorithm also performs satisfactorily for large OWFs with multiple OSSs, where the clustering is intrinsically defined in the mathematical formulation. It does not require predecessor algorithms to group WTs into OSSs, avoiding in this way artificially biased solutions.

Acknowledgment

This research has received funding from the Baltic InteGrid Project (www.baltic-integrid.eu/).

The authors thank Daniel Hermosilla Minguijón for his practical inputs throughout the development of the computational experiments, and Martina Fischetti, Lead Engineer - Operational Research in Vattenfall, for providing the data to perform the benchmark analysis.