*3.2. Constraints*

The investments in new transformers (Equation (12)), wind units (Equation (13)), and PV modules (Equation (14)) are limited per node n.

$$0 \le \sum\_{t \in \Omega^{\Gamma}} \mathcal{Y}\_t^{\text{SS}; \mathfrak{n}} \le \overline{\mathcal{S}^{\text{NEW}; \mathfrak{n}}} / \overline{\mathcal{S}^{\text{SS}}}; \ \forall \mathfrak{n} \in \Omega^{\text{SS}} \tag{12}$$

$$0 \le \sum\_{t \in \Omega^{\Gamma}} Y\_t^{\text{ord}, n} \le \overline{Y\_n^{\text{ord}}}; \forall n \in \Omega^L \tag{13}$$

$$0 \le \sum\_{t \in \Omega^{\Gamma}} \mathsf{Y}\_t^{\mathrm{pv}, n} \le \overline{\mathsf{Y}\_n^{\mathrm{pv}}}; \ \forall n \in \Omega^L. \tag{14}$$

The maximum value of the power allowed to install at each bus n is limited by Equation (15).

$$\overline{P^{m\text{od}}} \ge \sum\_{t \in \Omega^T} \left( \overline{P^{m\text{od}}} \, \, \, \mathcal{Y}\_t^{\text{od},\mu} + \, \, \overline{P^{p\text{v}}} \, \, \mathcal{Y}\_t^{\text{ev},\mu} \right); \, \forall n \in \Omega^L. \tag{15}$$

The investments are limited by Equations (16) and (17). Equation (16) refers to the annual investment payment limit that describes the budget available for each investment period, while the portfolio investment (Equation (17)) is related to the amount of money that is available for investment in the long term.

$$c\dot{c}\_{l} \le c\dot{l}^{\text{bgt}};\ \forall t.\tag{16}$$

$$\sum\_{t \in \Omega^{\rm T}} \beta^t \left[ \sum\_{n \in \Omega^{\rm SS}} c i^{\rm SS} Y\_t^{\rm SS, n} + \sum\_{n \in \Omega^{\rm L}} \left( c i^{\rm pv} Y\_t^{\rm pv, n} + c i^{\rm rad} Y\_t^{\rm wd, n} \right) \right] \le c I\_{\rm LC}^{\rm ggt}. \tag{17}$$

The power flow constraints (Equations (18)–(41)) refer to both active and reactive power equations. The distribution system is composed of two types of buses (load and substation buses) where each kind of bus has a different load flow expression.

$$\begin{split} \sum\_{n \in \Omega^{\rm N}} \left( P\_{t,k,\omega}^{+,\rm n,m} - P\_{t,k,\omega}^{-,\rm n,m} \right) &+ P\_{t,k,\omega}^{\rm n\varepsilon,m} + P\_{t,k,\omega}^{\rm nd} + P\_{t,k,\omega}^{\rm pw,m} - \sum\_{n \in \Omega^{\rm N}} \left( P\_{t,k,\omega}^{+,\rm n,m} - P\_{t,k,\omega}^{-,\rm m,n} + R^{\rm m,n} I\_{t,k,\omega}^{\rm aqr,m} \right) \\ &= f\_{t} f\_{k,\omega}^{\rm fd} P\_{t,\omega}^{\rm kd,m}; \forall \{m \in \Omega^{\rm L}\_{-}, t, k, \omega\} \end{split} \tag{18}$$

$$\sum\_{n \in \Omega^{N}} \left( P^{+,n,m}\_{t,k,\omega} - P^{-,n,m}\_{t,k,\omega} \right) - \left( P^{+,m,n}\_{t,k,\omega} - P^{-,m,n}\_{t,k,\omega} + R^{m,n} I^{\text{op},m,n}\_{t,k,\omega} \right) = f\_{t} f^{kl}\_{k,\omega} P^{l,m}, \forall \left( m \in \Omega^{\text{SS}}, t, k, \omega \right) \tag{19}$$
 
$$\sum\_{n \in \Omega^{N}} \left( Q^{+,n,m}\_{t,k,\omega} - Q^{-,n,m}\_{t,k,\omega} \right) + Q^{\text{os},m}\_{t,k,\omega} + Q^{\text{un},m}\_{t,k,\omega} + Q^{\text{up},m}\_{t,k,\omega} - \left( Q^{+,m,n}\_{t,k,\omega} - Q^{-,m,n}\_{t,k,\omega} + X^{\text{un}}\_{t} \right)\_{t,k,\omega}^{\text{op},m,n} \tag{20}$$
 
$$= f\_{t} f^{kl}\_{k,\omega} Q^{l,m} 0; \forall \left( m \in \Omega^{\text{}}, t, k, \omega \right) \tag{21}$$
 
$$\sum\_{n \in \Omega^{\text{N}}} \left( Q^{+,n,m}\_{t,k,\omega} - Q^{-,n,m}\_{t,k,\omega} \right) + Q^{\text{os},m}\_{t,k,\omega} - \left( Q^{+,m,n}\_{t,k,\omega} - Q^{-,m,n}\_{t,k,\omega} + X^{\text{un}}\_{t} \right)\_{t,k,\omega}^{\text{op},m,n} \tag{21}$$
 
$$= f\_{t} f^{kl}\_{k,\omega} Q^{l,m}; \forall \left( m \in \Omega^{\text{SS}}, t, k, \omega \right). \tag{22}$$

The voltage is related to the different electrical values (Equation (22)) and is limited by the maximum and minimum values (Equation (23)).

$$\begin{aligned} \left(V\_{t,k,\nu}^{\mathrm{sqr},m} - \mathbb{Z}^{m,n}\mathbf{1}\_{t,k,\nu}^{\mathrm{sqr},m,n} - V\_{t,k,\nu}^{\mathrm{sqr},n} - \mathbb{2}\left(\mathbb{R}^{m,n}\begin{pmatrix} P\_{t,k,\nu}^{+,m,n} - P\_{t,k,\nu}^{-,m,n} \\ t,k,\nu \end{pmatrix} + X^{m,\nu} \begin{pmatrix} Q\_{t,k,\nu}^{+,m,n} - Q\_{t,k,\nu}^{-,m,n} \\ t,k,\nu \end{pmatrix}\right) \\ = 0; \;\forall \left(m,n \in \mathbb{Z}^{N}, t, k, \nu\right) \end{aligned} \tag{22}$$

$$
\underline{V}^2 \le V\_{t,k;w}^{sqr,m} \le \overline{V}^2; \ \forall \left(m \in \Omega^N, t, k, w\right). \tag{23}
$$

Equation (24) represents the maximum current permitted through the distribution feeders. These limits are applied to the active (25) and (26) and reactive power Equations (27) and (28). Equations (29) and (30) are necessary to avoid inverse flows.

$$0 \le I\_{t,k,\omega}^{\kappa qr,m,n} \le \overline{I^{m,n}}^2; \ \forall \left(m, n \in \Omega^N, t, k\_\prime\prime\,\omega\right) \tag{24}$$

$$P\_{t,k,\omega}^{+,m,n} \le V^{mnm} \overline{I^{m,n}} \, \mathcal{Y}\_{t,k,\omega}^{P+,m,n}; \; \forall \left(m, n \in \Omega^N, t, k, \omega\right) \tag{25}$$

$$\mathbb{P}\,P\_{t,k,\omega}^{-,m,n} \le V^{m\,\mathrm{nm}} \,\overline{I^{m,n}} \,\mathrm{Y}\_{t,k,\omega}^{P-,m,n}; \;\forall \left(m,n \in \Omega^N, t, k, \omega\right) \tag{26}$$

$$\mathbb{Q}\begin{array}{c}Q^{+,m,n}\_{t,k,\omega}\leq V^{\text{nom}}\ \overline{I^{m,n}}\ \mathcal{Y}^{Q+,m,n}\_{t,k,\omega};\ \forall \left(m,n\in\Omega^{N},t,k,\omega\right)\end{array}\tag{27}$$

$$\mathbb{Q}\,\mathrm{Q}^{-,\mathrm{m},n}\_{t,k,\omega} \le V^{\mathrm{nom}}\,\overline{I^{\mathrm{m},n}}\,\mathrm{Y}^{Q-,\mathrm{m},n}\_{t,k,\omega};\,\forall \Big(m,n \in \Omega^N, \mathfrak{t}, k, \omega\Big)\tag{28}$$

$$\mathcal{Y}\_{t,k,\omega}^{P+,\mathfrak{m},\mathfrak{n}} + \mathcal{Y}\_{t,k,\omega}^{P-,\mathfrak{m},\mathfrak{n}} \le 1; \ \forall \left(m, n \in \mathcal{Q}^N, t, k, \omega\right) \tag{29}$$

$$\mathcal{Y}\_{\cdot}^{Q+,\mathfrak{m},\mathfrak{n}} \pm \mathcal{Y}\_{\cdot}^{Q-,\mathfrak{m},\mathfrak{n}} < 1 \cdot \mathcal{Y}(m, n \in \mathcal{Q}^N, t \mid k \mid\_{\mathcal{Q}}) \tag{30}$$

$$Y\_{t,k,\omega}^{Q+,m,n} + Y\_{t,k,\omega}^{Q-,m,n} \le 1; \ \forall \left(m, n \in \Omega^N, t, k, \omega\right) \tag{30}$$

The power flow has been linearized using Equations (31) to (41) as explained in [18] and [19].

$$\Psi\_{t,k,\omega}^{P+,n,\mathfrak{m}} \in \{0,1\}; \ \mathfrak{d}\left(m, n \in \Omega^N, t, k, \omega\right) \tag{31}$$

$$\mathcal{Y}\_{t,k,\omega}^{P-\text{;n,m}} \in \{0,1\}; \ \forall \left(m, n \in \mathcal{Q}^N, t, k, \omega\right) \tag{32}$$

$$\{Y\_{t,k,\omega}^{Q+;n,m} \in \{0,1\}; \ \forall \{m,n \in \mathcal{Q}^N, t, k, \omega\} \tag{33}$$

$$\{Y\_{t,k,\omega}^{Q-\mu,m} \in \{0,1\}; \ \forall \{m,n \in \mathcal{Q}^N, t, k, \omega\} \tag{34}$$

$$0 \le \Delta P\_{t,k,\omega}^{n,m,r} \le \Delta S\_{t,k,\omega}^{m,n,r}; \ \forall (m, n \in \Omega^N, r \in \Omega^R, t, k, \omega) \tag{35}$$

$$0 \le \Delta Q\_{t,k,\omega}^{n,m,r} \le \Delta S\_{t,k,\omega}^{m,n,r}; \ \forall \left(m, n \in \Omega^N, r \in \Omega^R, t, k, \omega\right) \tag{36}$$

$$\sum\_{r \in \Omega^{\mathbb{R}}} \left( m\_{t,k,\omega}^{m,n,r} \, \Delta P\_{t,k,\omega}^{m,n,r} \right) + \sum\_{r \in \Omega^{\mathbb{R}}} \left( m\_{t,k,\omega}^{m,n,r} \, \Delta Q\_{t,k,\omega}^{m,n,r} \right) = V^{mom2} I\_{t,k,\omega}^{spr,n,r}; \; \forall \left( m,n \in \Omega^{N}, t, k, \omega \right) \tag{37}$$

$$P\_{t,k,\mu}^{+,\text{;m,n}} + P\_{t,k,\mu}^{-,\text{;m,n}} = \sum\_{r \in \Omega^R} \Delta P\_{t,k,\mu}^{\text{;n,n},r}; \forall \left(m, n \in \Omega^N, t, k, \omega\right) \tag{38}$$

$$\Delta Q\_{t,k,\omega}^{+,m,n} + Q\_{t,k,\omega}^{-,m,n} = \sum\_{r \in \Omega^R} \Delta Q\_{t,k,\omega}^{m,n,r}; \forall \left(m, n \in \Omega^N, t, k, \omega\right) \tag{39}$$

$$m\_{m,n,r,t,k,\omega} = (2r-1)\Delta S\_{m,n,r,t,k,\omega}; \; \forall \left(m, n \in \Omega^N, r \in \Omega^R, t, k, \omega\right) \tag{40}$$

$$
\Delta S\_{m,n,r,t,k,\omega} = \frac{\left(V\_{nom}\,\,\overline{l\_{m,n}}\right)}{\overline{R}};\ \forall \left(m,n\in\Omega^N, r\in\Omega^R, t, k, \omega\right). \tag{41}
$$

The not supplied active and reactive power must be less than the active and reactive power demand, respectively. Note that the demand is updated each period.

0 ≤ *Pns*,*<sup>n</sup> <sup>t</sup>*,*k*,<sup>ω</sup> <sup>≤</sup> *ft <sup>f</sup> ld <sup>k</sup>*,<sup>ω</sup>*Pld*,*n*; <sup>∀</sup> - *n* ∈ Ω*L*, *t*, *k*, ω (42)

$$0 \le \mathcal{Q}\_{t,k,\omega}^{\rm us,n} \le f\_t \mathcal{f}\_{k,\omega}^{\rm ld} \mathcal{Q}^{\rm ld,n};\ \forall \left(n \in \Omega^{\rm l}, t, k, \omega\right). \tag{43}$$

The substation available power is subject to yearly updates (44). New investments are made considering the already existing power and the new possible power of the candidate transformers (45), (46).

$$ST\_t^{SS; \mathfrak{n}} = S^{SS; \mathfrak{n}} + S\_t^{NE; \mathfrak{n}}; \; \forall \left(n \in \Omega^{SS}, t\right) \tag{44}$$

$$S\_t^{NEW;n} = Y\_t^{SS\mu} \,\overline{S^{SS}} / S^{base}; \; \forall \left(n \in \Omega^{SS}, t = 1\right) \tag{45}$$

$$S\_t^{NEN\mu} = S\_{t-1}^{NEN\mu} + Y\_t^{SS\mu} \overline{S^{SS}} / S^{\text{base}}; \ \forall \left(n \in \Omega^{SS}, t > 1\right). \tag{46}$$

The renewable power (wind turbines and PV modules) available is updated every year. Note that the vector of the candidate nodes appears in Equations (47) to (50).

$$P T\_t^{\text{ord},n} = \overline{P^{\text{ord}}} \ Y\_t^{\text{ord},n} \ C^{\text{ord};n} ; \ \mathsf{V} \left( n \in \Omega^L, t = 1 \right) \tag{47}$$

$$PT\_t^{\text{wd,n}} = \overline{P^{\text{wd}}} \ Y\_t^{\text{wd,n}} \text{ C}^{\text{wd,n}} + PT\_{t-1}^{\text{wd,n}}; \ \forall \left(n \in \Omega^L, t > 1\right) \tag{48}$$

$$PT\_t^{\text{pv},n} = \overline{P^{\text{pv}}} \ Y\_t^{\text{pv},n} \ C^{\text{pv},n}; \ \forall \left(n \in \Omega^L, t = 1\right) \tag{49}$$

$$PT\_t^{\text{pv},n} = \overline{P^{\text{pv}}} \ Y\_t^{\text{pv},n} \ C^{\text{pv},n} + PT\_{n,t-1}^{\text{pv},n}; \ \forall \left(n \in \Omega^L, t > 1\right). \tag{50}$$

The maximum amount of generation that is available is limited by the generation levels of each of the technologies available (Equations (51) and (52)).

$$0 \le P\_{t,k,\omega}^{\text{und},\text{n}} \le f\_{k,\omega}^{\text{nd}} P T\_t^{\text{ud},\text{n}}; \ \forall \left(n \in \Omega^{L}, t, k, \omega\right) \tag{51}$$

$$0 \le \mathcal{P}\_{t,k,\mu}^{\mathrm{pv},\mathrm{u}} \le f\_{k,\mu}^{\mathrm{pv}} P T\_t^{\mathrm{pv},\mathrm{u}};\ \forall \left(n \in \mathcal{Q}^L, t, k, \omega\right). \tag{52}$$

The active power injected into the distribution system must keep a constant power factor (Equation (53)). In this way the reactive power is also limited (Equation (54)).

$$P\_{t,k,\omega}^{\rm SS,n} \le S T\_t^{\rm SS,n} / \sqrt{\left(1 + \tan(\eta^{\rm SS})^2\right)}; \ \forall \left(n \in \mathcal{Q}^{\rm SS}, t, k, \omega\right) \tag{53}$$

$$Q\_{t,k,\omega}^{\text{SS},n} \le \tan(qr^{\text{SS}}) P\_{t,k,\omega'}^{\text{SS},n}; \ \forall \left(n \in \Omega^{\text{SS}}, t, k, \omega\right). \tag{54}$$

The reactive renewable power that is injected in the distribution system must have a constant power factor (Equations (55) and (56)).

$$0 \le \mathcal{Q}\_{t,k,\omega}^{\text{ind},n} \le P\_{t,k,\omega}^{\text{snd},n} \tan(\phi^{\text{rad}}); \ \forall \left(n \in \Omega^L, t, k\right) \tag{55}$$

$$0 \le \mathcal{Q}\_{t,k,\omega}^{pv,n} \le \mathcal{P}\_{t,k,\omega}^{pv,n} \tan(\varphi^{pv});\ \forall \left(n \in \mathcal{Q}^L, t, k\right). \tag{56}$$

#### **4. Solution Procedure: Benders' Decomposition**

Model (1)–(56) is set as a MILP problem whose size depends on the number of time blocks, scenarios, and years considered in the study. In order to obtain informed expansion decisions, a large number of time blocks and scenarios are needed to obtain a good representation of the loads and renewable power outputs, as well as a large number of years. However, in such a case the MILP problem (1)–(56) may become computationally intractable. Nevertheless, note that if variables *YSS*;*<sup>n</sup> <sup>t</sup>* , *<sup>Y</sup>pv*;*<sup>n</sup> <sup>t</sup>* , and *Ywd*;*<sup>n</sup> <sup>t</sup>* , i.e., the expansion decision variables, are fixed to given values, then, problem (1)–(56) can be decomposed into a set of smaller problems, one for each time block, scenario, and year. This allows us to use Benders' decomposition to solve the problem.

Model (1)–(56) uses binary variables for the linearization of power flow equations. Since Benders' decomposition cannot have binary variables in the sub-problem, a procedure to fix the values of these variables must be used. Further details are shown in step 1.

In the following sections superindex (υ) indicates the values in the υ-th Benders' iteration. The proposed Benders' decomposition algorithm comprises the steps below:

1. Step 0. Initialization: Initialize the iteration counter, υ = 1, the upper bound of the objective function, z(1) up = ∞ and its lower bound, *z* (1) *down* <sup>=</sup> −∞ and solve the following problem:

$$\begin{array}{c} \underset{\begin{subarray}{c} \mathbf{Y}^{SS;\mu}, \mathbf{Y}^{\text{pv:pu}}, \mathbf{Y}^{\text{und};\mu}, \alpha \end{array}}{\text{minimize}} \quad \underset{\begin{subarray}{c} t \in \Omega^{\top} \end{subarray}}{\text{minimize}} \beta\_{t} \\ \text{c.t.} \tag{57}$$

subject to constraints (12)–(17) and

$$
\alpha \ge \alpha^{down}.\tag{58}
$$

This problem has the trivial solution: α(1) = α*down*., *YSS*;*<sup>n</sup> <sup>t</sup>* (1), *<sup>Y</sup>pv*;*<sup>n</sup> <sup>t</sup>* (1), *<sup>Y</sup>wd*;*<sup>n</sup> <sup>t</sup>* (1) <sup>=</sup> 0, <sup>∀</sup>(*t*, *<sup>n</sup>*).

After problem (57), (58), (12)–(17)above is solved, the lower bound of the optimal value of the objective function (1) is updated:

$$z\_{\text{down}}^{(1)} = \sum\_{t \in \mathcal{O}^T} \beta\_t \text{ci}\_t^{(1)} + a^{(1)}.\tag{59}$$

2. Step 1. Sub-problem solution: The variables *YSS*;*<sup>n</sup> <sup>t</sup>* , *<sup>Y</sup>pv*;*<sup>n</sup> <sup>t</sup>* , and *<sup>Y</sup>wd*;*<sup>n</sup> <sup>t</sup>* , <sup>∀</sup>(*t*, *<sup>n</sup>*) are set, to their optimal values from the previous step and solve the following sub-problems, one for each year *t*, time block *<sup>k</sup>*, and scenario <sup>ω</sup>. <sup>⎧</sup>

$$\begin{cases} \underset{\Xi^1\_{t,k,\omega}}{\min} \quad \beta\_t \mathsf{N}^{\mathbb{I}}\_k \; \gamma\_{k,\omega} \mathsf{com}\_{t,k,\omega} \end{cases} \tag{60}$$

subject to constraints (7)–(11), (18)–(56) and

$$\mathcal{Y}\_t^{\text{SS;n}} = \mathcal{Y}\_t^{\text{SS;n}(\nu)} ; \forall n \in \mathcal{Q}^{\text{SS}} \tag{61}$$

$$\mathcal{Y}\_t^{\mathbb{M};\mathbb{M}} = \mathcal{Y}\_t^{\mathbb{M};\mathbb{M}(\mathbb{U})} ; \mathsf{\forall}{n} \in \Omega^L \tag{62}$$

$$\left| Y\_t^{\text{ord};n} = Y\_t^{\text{ord};n} ; \forall n \in \Omega^L \right| ; \; \forall (t, k, \omega) \tag{63}$$

where:

$$
\begin{bmatrix}
\boldsymbol{\Sigma}\_{\boldsymbol{\upbeta},\boldsymbol{\omega},\boldsymbol{\upmu}}^{1} \\
\Gamma\_{\boldsymbol{\upbeta},\boldsymbol{\upmu},\boldsymbol{\upmu}}^{\text{spr},\boldsymbol{\upmu}},\boldsymbol{P}\_{\boldsymbol{\upbeta},\boldsymbol{\upmu},\boldsymbol{\upmu}}^{\text{spr},\boldsymbol{\upmu}},\boldsymbol{P}\_{\boldsymbol{\upbeta},\boldsymbol{\upmu},\boldsymbol{\upmu}}^{\text{spr},\boldsymbol{\upmu}},\boldsymbol{P}\_{\boldsymbol{\upbeta},\boldsymbol{\upmu},\boldsymbol{\upmu}}^{\text{sp},\boldsymbol{\upmu},\boldsymbol{\upmu}},\boldsymbol{P}\_{\boldsymbol{\upbeta},\boldsymbol{\upmu},\boldsymbol{\upmu}}^{\text{sp},\boldsymbol{\upmu},\boldsymbol{\upmu}},\boldsymbol{Q}\_{\boldsymbol{\upbeta},\boldsymbol{\upmu}}^{\text{sp},\boldsymbol{\upmu}},\boldsymbol{Q}\_{\boldsymbol{\upbeta},\boldsymbol{\upmu}}^{\text{sp},\boldsymbol{\upmu}},\boldsymbol{Q}\_{\boldsymbol{\upbeta},\boldsymbol{\upmu},\boldsymbol{\upmu}}^{\text{sp},\boldsymbol{\upmu}},\boldsymbol{Q}\_{\boldsymbol{\upbeta},\boldsymbol{\upmu},\boldsymbol{\upmu}}^{\text{sp},\boldsymbol{\upmu}},\boldsymbol{Q}\_{\boldsymbol{\upbeta},\boldsymbol{\upmu},\boldsymbol{\upmu}}^{\text{sp},\boldsymbol{\upmu}},\boldsymbol{Q}\_{\boldsymbol{\upbeta},\boldsymbol{\upmu}}^{\text{sp},\boldsymbol{\upmu}},\boldsymbol{Q}\_{\boldsymbol{\upbeta},\boldsymbol{\upmu}}^{\text{sp},\boldsymbol{\upmu}},\boldsymbol{Q}\_{\boldsymbol{\upbeta},\boldsymbol{\upmu}}^{\text{sp},\boldsymbol{\upmu}},\boldsymbol{V}\_{\boldsymbol{\upbeta},\boldsymbol{\upmu}}^{\text{sp},\boldsymbol{\upmu}},\boldsymbol{V}\_{\boldsymbol{\upbeta},\boldsymbol{\upmu}}^{\text{sp},\boldsymbol{\upmu}},\boldsymbol{V}\_{\boldsymbol{\upbeta},\boldsymbol{\upmu}}^{\text{sp},\boldsymbol{\upphi}},\$$

Each instance of the sub-problem, one for each year *t*, time block *k*, and scenario ω, is solved twice. The first simulation solves problem (60)–(63), (7)–(11), (18)–(56) as a MILP problem, to obtain the values of *YP*+,*m*,*<sup>n</sup> <sup>t</sup>*,*k*,<sup>ω</sup> , *<sup>Y</sup>P*−,*m*,*<sup>n</sup> <sup>t</sup>*,*k*,<sup>ω</sup> , *<sup>Y</sup>Q*+,*m*,*<sup>n</sup> <sup>t</sup>*,*k*,<sup>ω</sup> , and *<sup>Y</sup>Q*−;*n*,*<sup>m</sup> <sup>t</sup>*,*k*,<sup>ω</sup> , which are the sub-problem's binary variables. Then, the second run uses the optimal values of variables *YP*+,*m*,*<sup>n</sup> <sup>t</sup>*,*k*,<sup>ω</sup> , *<sup>Y</sup>P*−,*m*,*<sup>n</sup> <sup>t</sup>*,*k*,<sup>ω</sup> , *<sup>Y</sup>Q*+,*m*,*<sup>n</sup> <sup>t</sup>*,*k*,<sup>ω</sup> , and *<sup>Y</sup>Q*−;*n*,*<sup>m</sup> <sup>t</sup>*,*k*,<sup>ω</sup> , obtained in the first simulation, as parameters for solving sub-problem (64)–(67), (7)–(11), (18)–(56) i.e., in this case we solve sub-problem (64)–(67), (7)–(11), (18)–(56) as the following LP problem.

$$\begin{cases} \underset{\Xi^2\_{t,k,\omega}}{\text{minimize}} & \beta\_l N\_k^{\text{lt}} \; \gamma\_{k,\omega} \text{com}\_{t,k,\omega} \\ & \sum\_{k=1}^{2} \beta\_l \; \text{as}\_{t,k,\omega} \end{cases} \tag{64}$$

subject to constraints (7)–(11), (18)–(56) and

$$\mathcal{Y}\_t^{SS; \mu} = \mathcal{Y}\_t^{SS; \mathfrak{n}(v)} \; : \; \lambda\_{1, \mathfrak{n}, t, k, \omega} ; \; \forall n \in \mathcal{Q}^{SS} \tag{65}$$

$$\mathcal{Y}\_t^{\text{pv};\mu} = \mathcal{Y}\_t^{\text{pv};n(\nu)} : \ \lambda\_{2,n,t,k,\mu} ; \ \forall n \in \mathcal{Q}^L \tag{66}$$

$$Y\_t^{\text{adv};n} = Y\_t^{\text{adv};n(\nu)} \;:\; \lambda\_{3,n,t,k,\omega}; \; \forall n \in \Omega^L \big| \; \forall \; (t,k,\omega). \tag{67}$$

where:

$$
\Sigma\_{\mathbf{t},\mathbf{k},\omega}^{2} = \left\{ \Gamma\_{\mathbf{t},\mathbf{k},\omega}^{\mathrm{spr},\mathbf{n}}, \Gamma\_{\mathbf{t},\mathbf{k},\omega'}^{\mathrm{usp},\mathbf{n}}, \Gamma\_{\mathbf{t},\mathbf{k},\omega'}^{\mathrm{pv},\mathbf{n}}, \Gamma\_{\mathbf{t},\mathbf{k},\omega'}^{\mathrm{S}\mathbf{S}\mathbf{n}}, \Gamma\_{\mathbf{t},\mathbf{k},\omega'}^{\mathrm{S}\mathbf{S}\mathbf{n}}, \zeta\_{\mathbf{t},\mathbf{k},\omega'}^{\mathrm{S}\mathbf{S}\mathbf{n}}, \zeta\_{\mathbf{t},\mathbf{k},\omega'}^{\mathrm{pv},\mathbf{n}}, \zeta\_{\mathbf{t},\mathbf{k},\omega'}^{\mathrm{pv},\mathbf{n}}, \Gamma\_{\mathbf{t},\mathbf{k},\mathbf{n}'}^{\mathrm{spr},\mathbf{n}}, \Lambda\_{\mathbf{t},\mathbf{k},\omega}^{\mathrm{p},\mathbf{n},\mathbf{r}} \right\}.
$$

The outputs of sub-problem (60)–(63), (7)–(11), (18)–(56) are variables of set Ξ<sup>1</sup> *<sup>t</sup>*,*k*,ω, and the outputs of sub-problem (64)–(67), (7)–(11), (18)–(56) are variables of set Ξ<sup>2</sup> *<sup>t</sup>*,*k*,ω, as well as dual variables λ1,*n*,*t*,*k*,ω, λ2,*n*,*t*,*k*,ω, and λ3,*n*,*t*,*k*,ω. Note that the sensitivities used for Benders' cuts in (65), (66) and (67) are only obtained after the second run.

A sub-problem (64)–(67), (7)–(11), (18)–(56) is solved for each year *t*, time block *k*, and scenario ω. In order to formulate the Benders' cuts in the Master problem (see Step 3), it is necessary to define and compute parameters λ˜ 1,*n*,*t*, λ˜ 2,*n*,*t*, and λ˜ 3,*n*,*<sup>t</sup>* as follows:

$$\bar{\lambda}\_{1,n,t} = \sum\_{k \in \mathcal{Q}^{\mathcal{K}}} \sum\_{\alpha \in \Psi\_k^{\omega}} \lambda\_{1,n,t,k,\omega} \cdot \forall \left(n \in \mathcal{Q}^{\mathcal{S}S}, t\right) \tag{68}$$

$$\tilde{\lambda}\_{2,n,t} = \sum\_{k \in \Omega^K} \sum\_{\alpha \in \Psi\_k^{\omega}} \lambda\_{2,n,t,k,\omega} \cdot \mathbb{V}\left(n \in \Omega^L, t\right) \tag{69}$$

$$\lambda\_{3,n,t} = \sum\_{k \in \Omega^K} \sum\_{\omega \circ \Psi\_k^{\omega}} \lambda\_{3,n,t,k,\omega}; \forall \{n \in \Omega^L, t\}. \tag{70}$$

Finally, the upper bound of the optimal value of the objective function (1) is updated as:

$$z\_{up}^{(\nu)} = \sum\_{t \in \Omega^{\Gamma}} \beta\_t c i\_t^{(\nu)} + \sum\_{t \in \Omega^{\Gamma}} \sum\_{k \in \Omega^K} \sum\_{\alpha \in \Psi\_k^{\omega}} \beta\_t \mathbf{N}\_k^{\text{th}} \boldsymbol{\gamma}\_{k,\omega} c o m\_{t,k,\omega}^{\text{(\nu)}} \,. \tag{71}$$


$$\chi\_{\mathbf{Y}\_t^{SS;\mathbf{u}},\mathbf{Y}\_t^{\text{pv};\mathbf{u}},\mathbf{Y}\_t^{\text{adv};\mathbf{u}},\alpha} \quad \sum\_{t \in \Omega^T} \beta\_t c i\_t + \alpha \tag{72}$$

subject to constraints (12)–(17) and

$$a \ge a^{down} \tag{73}$$

$$\begin{split} \sum\_{\boldsymbol{\theta},\boldsymbol{\lambda},\boldsymbol{\omega}} & \sum\_{\boldsymbol{\lambda},\boldsymbol{\omega}} \beta\_{\boldsymbol{\lambda}} \mathbf{N}\_{\boldsymbol{\lambda}}^{\text{H}} \boldsymbol{\gamma}\_{\boldsymbol{\lambda},\boldsymbol{\omega}} \mathbf{com}\_{\boldsymbol{\lambda},\boldsymbol{\lambda},\boldsymbol{\omega}} (\boldsymbol{v}) \\ & & + \sum\_{\boldsymbol{\ell} \in \mathcal{LI}^{\text{T}}} \left( \sum\_{\boldsymbol{\sigma} \in \mathbb{Z}^{\text{S} \times \mathbb{Z}}} \widetilde{\lambda}\_{1,\boldsymbol{\mu},\boldsymbol{\ell}}^{(\text{h})} \left( \mathbf{Y}\_{\boldsymbol{\ell}}^{\text{SS;\mu}} - \mathbf{Y}\_{\boldsymbol{\ell}}^{\text{SS;\mu} \times \boldsymbol{\uprho}} \right) \right) \\ & & + \sum\_{\boldsymbol{\kappa} \in \mathbb{Z}^{\text{L}}} \left[ \widetilde{\lambda}\_{2,\boldsymbol{\mu},\boldsymbol{\ell}}^{(\text{h})} \left( \mathbf{Y}\_{\boldsymbol{\ell}}^{\text{R;\mu}\text{H}} - \mathbf{Y}\_{\boldsymbol{\ell}}^{\text{R;\mu} \times \boldsymbol{\uprho}} \right) + \widetilde{\lambda}\_{3,\boldsymbol{\mu},\boldsymbol{\ell}}^{(\text{h})} \left( \mathbf{Y}\_{\boldsymbol{\ell}}^{\text{R;\mu}\text{H}} - \mathbf{Y}\_{\boldsymbol{\ell}}^{\text{R;\mu} \times \boldsymbol{\uprho}} \right) \right] \right) \leq \boldsymbol{\alpha}; \ \forall \boldsymbol{h} \\ & & = 1,\ldots,\boldsymbol{\nu} - 1 \end{split} \tag{74}$$

where *cit* is defined in Equations (5) and (6) and their parameters in Equations (2)–(4).

Each constraint in (74) is known as a Benders' cut [20]. Benders' cuts approximate objective function (1) from below as a function of the investment variables. Notice that at every iteration the size of master problem (72)–(74), (12)–(17) increases since a new constraint is added.

Next, the lower bound of the objective function is updated (1):

$$z\_{\text{down}}^{(\upsilon)} = \sum\_{t \in \Omega^{\mathbb{T}}} \beta\_t \text{c}i\_t^{(\upsilon)} + \alpha^{(\upsilon)}.\tag{75}$$

Then, the algorithm continues at Step 1.

For clarity purposes, the flowchart of the algorithm is shown in Figure 1.

**Figure 1.** Benders' decomposition flowchart.

Note that the Benders' decomposition algorithm described in this section can solve each sub-problem separately, and hence, the computational time of the proposed algorithm grows linearly with the number of scenarios, time blocks, and years, being its scalability excellent.
