A hybrid Benders decomposition and bees algorithm matheuristic approach to transmission expansion planning considering energy storage

TL;DR

This paper presents a new hybrid algorithm combining Benders decomposition and the Bees Algorithm, designed for large-scale transmission expansion and energy storage planning, offering fast, high-quality solutions with potential for optimality guarantees.

Contribution

The paper introduces the Bee-Benders Hybrid Algorithm (BBHA), a novel hybrid matheuristic that enhances solution speed and quality for complex MILP problems, especially in energy network planning.

Findings

BBHA performs at least as well as individual methods.

BBHA significantly outperforms standalone approaches.

Algorithm is highly effective for large-scale energy planning problems.

Abstract

This paper introduces a novel hybrid optimisation algorithm that combines elements of both metaheuristic search and integer programming. This new matheuristic combines elements of Benders decomposition and the Bees Algorithm, to create the Bee-Benders Hybrid Algorithm (BBHA) which retains many of the advantages both of the methods. Specifically it is designed to be easily parallelizable, to produce good solutions quickly while still retaining a guarantee of optimality when run for a sufficiently long time. The algorithm is tested using a transmission network expansion and energy storage planning model, a challenging and very large scale mixed integer linear programming problem. Transmission network planning problems are already difficult on their own. When including the planning for storage systems in the network, the variation of demand over time has to be taken into account significantly increasing the size and difficulty of the optimization problem. The BBHA is shown to be highly effective hybrid matheuristic algorithm that performs at least as well as either Benders decomposition or the Bees Algorithm where these are effective on their own, and significantly improves upon the individual approaches where neither component part has a pronounced advantage. While the paper demonstrates the effectiveness in terms of the concrete electricity network planning problem, the algorithm could be readily applied to any mixed integer linear program, and is expected to work particularly well whenever this has a structure that is amenable to Benders decomposition.

Tables

Table 1. Table 1: Parameters for the BBHA.

Name	Description	Default value
ne	number of elite sites	1
nb	number of best sites	2
nre	recruited bees for elite sites	10
nrb	recruited bees for remaining best sites	5
ngh	maximum size of neighbourhood for local search	8

Table 2. Table 2: Tuning results for IEEE 25-bus test system. Note: We use a composite trapezoidal rule to integrate along the time axis and then rescale (Scaled Trapz) which gives a measure of average solution quality over the half hour duration of the run. Lower is better.

Scenario	ne	nb	nre	nrb	Objective	iterations	Scaled Trapz
					(US$103)
Long peak	$2$	$3$	$10$	$5$	$44,703$	$93$	$0.26$
Long peak	$1$	$4$	$10$	$5$	$44,703$	$92$	$0.30$
Long peak	$1$	$2$	$30$	$10$	$43,812$	$56$	$0.32$
Long peak	$1$	$4$	$20$	$15$	$44,703$	$39$	$0.37$
Long peak	$1$	$2$	$10$	$5$	$44,703$	$133$	$0.40$
Long peak	$1$	$3$	$30$	$25$	$44,703$	$32$	$0.41$
Long peak	$1$	$2$	$20$	$15$	$55,282$	$67$	$0.41$
Long peak	$1$	$2$	$30$	$15$	$43,812$	$50$	$0.44$
Long peak	$1$	$3$	$20$	$10$	$44,703$	$57$	$0.45$
Long peak	$3$	$4$	$20$	$10$	$45,378$	$36$	$0.46$
Long peak	$1$	$3$	$30$	$20$	$44,703$	$36$	$0.47$
Long peak	$1$	$3$	$10$	$5$	$45,378$	$88$	$0.50$
Long peak	$2$	$3$	$20$	$5$	$44,479$	$48$	$0.50$
Long peak	$2$	$3$	$20$	$10$	$43,812$	$52$	$0.50$
Long peak	$3$	$4$	$20$	$5$	$46,364$	$30$	$0.54$
Long peak	$3$	$4$	$30$	$10$	$56,454$	$22$	$0.56$
Long peak	$1$	$3$	$20$	$15$	$55,088$	$53$	$0.58$
Long peak	$1$	$2$	$20$	$10$	$55,088$	$63$	$0.58$
Long peak	$2$	$3$	$30$	$20$	$44,709$	$30$	$0.60$
Long peak	$1$	$3$	$30$	$15$	$44,479$	$44$	$0.60$
Long peak	$3$	$4$	$30$	$20$	$51,489$	$24$	$0.62$
Long peak	$3$	$4$	$10$	$5$	$43,812$	$69$	$0.62$
Long peak	$1$	$2$	$30$	$20$	$44,703$	$50$	$0.62$
Long peak	$2$	$3$	$30$	$15$	$44,703$	$32$	$0.63$
Long peak	$3$	$4$	$30$	$15$	$55,120$	$24$	$0.64$
Long peak	$1$	$4$	$30$	$10$	$48,880$	$40$	$0.64$
Long peak	$1$	$4$	$30$	$15$	$51,785$	$35$	$0.65$
Long peak	$3$	$4$	$20$	$15$	$43,812$	$30$	$0.66$
Long peak	$1$	$3$	$30$	$10$	$46,045$	$47$	$0.66$
Long peak	$2$	$3$	$30$	$10$	$55,282$	$35$	$0.68$
Long peak	$2$	$3$	$20$	$15$	$60,263$	$34$	$0.69$
Long peak	$1$	$4$	$30$	$25$	$50,092$	$22$	$0.76$
Long peak	$1$	$4$	$20$	$10$	$45,414$	$38$	$0.77$
Long peak	$1$	$4$	$30$	$20$	$52,760$	$20$	$1.00$

Table 3. Table 3: Results for the BBHA compared with the basic Bee Algorithm and Benders Decompostion on their own. The BBHA & Bee algorithm were run three times with different parameter settings for each instance.

Network	Scenario	Params	BBHA worst	BBHA mean	BBHA best	Bee worst	Bee mean	Bee best	Benders
			(US$103)	(US$103)	(US$103)	(US$103)	(US$103)	(US$103)	(US$103)
46-bus	Long peak	1 2 10 5	$121,394.82$	$107,205.36$	$100,110.63$	$216,458.67$	$170,644.48$	$144,749.27$	$111,840.23$
46-bus	Long peak	1 2 30 10	$127,617.19$	$111,251.96$	$100,110.63$	$142,615.21$	$128,806.83$	$115,754.15$	$111,840.23$
46-bus	Long peak	2 3 10 5	$119,453.01$	$113,597.68$	$110,321.02$	$241,328.47$	$201,617.23$	$180,162.04$	$111,840.23$
46-bus	Short peak	1 2 10 5	$72,355.41$	$72,355.41$	$72,355.41$	$138,245.85$	$113,529.47$	$98,384.70$	$72,355.41$
46-bus	Short peak	1 2 30 10	$72,355.41$	$72,355.41$	$72,355.41$	$116,960.58$	$108,905.22$	$102,433.18$	$72,355.41$
46-bus	Short peak	2 3 10 5	$72,355.41$	$72,355.41$	$72,355.41$	$165,592.33$	$145,689.30$	$125,070.37$	$72,355.41$
46-bus	SGSC summer	1 2 10 5	$46,434.71$	$46,434.71$	$46,434.71$	$72,356.25$	$65,746.02$	$52,702.91$	$46,434.71$
46-bus	SGSC summer	1 2 30 10	$46,434.71$	$46,434.71$	$46,434.71$	$59,616.67$	$55,531.41$	$48,323.52$	$46,434.71$
46-bus	SGSC summer	2 3 10 5	$46,434.71$	$46,434.71$	$46,434.71$	$108,743.60$	$93,496.44$	$78,362.36$	$46,434.71$
46-bus	SGSC winter	1 2 10 5	$59,952.72$	$59,952.72$	$59,952.72$	$94,794.21$	$94,447.67$	$93,841.26$	$59,952.72$
46-bus	SGSC winter	1 2 30 10	$59,952.72$	$59,952.72$	$59,952.72$	$85,922.72$	$77,248.63$	$68,314.38$	$59,952.72$
46-bus	SGSC winter	2 3 10 5	$59,952.72$	$59,952.72$	$59,952.72$	$118,503.15$	$101,681.93$	$83,985.08$	$59,952.72$
93-bus	Long peak	1 2 10 5	$1,834.13$	$1,743.23$	$1,581.22$	$1,977.77$	$1,705.92$	$1,429.31$	$9,537.89$
93-bus	Long peak	1 2 30 10	$2,434.92$	$2,120.86$	$1,891.62$	$1,927.81$	$1,723.29$	$1,550.24$	$9,537.89$
93-bus	Long peak	2 3 10 5	$2,391.34$	$2,258.85$	$2,137.23$	$1,882.27$	$1,832.58$	$1,740.81$	$9,537.89$
93-bus	Short peak	1 2 10 5	$1,110.53$	$758.89$	$579.23$	$1,561.26$	$1,526.77$	$1,490.24$	$2,181.05$
93-bus	Short peak	1 2 30 10	$960.64$	$835.98$	$704.98$	$1,707.68$	$1,539.22$	$1,286.71$	$2,181.05$
93-bus	Short peak	2 3 10 5	$1,425.75$	$1,066.31$	$829.03$	$1,792.89$	$1,674.09$	$1,582.94$	$2,181.05$
93-bus	SGSC summer	1 2 10 5	$1,076.99$	$928.03$	$848.93$	$1,519.32$	$1,403.26$	$1,266.19$	$2,592.20$
93-bus	SGSC summer	1 2 30 10	$1,189.71$	$1,039.27$	$828.64$	$1,524.21$	$1,421.82$	$1,305.48$	$2,592.20$
93-bus	SGSC summer	2 3 10 5	$1,093.63$	$1,078.31$	$1,059.19$	$1,733.61$	$1,639.30$	$1,530.99$	$2,592.20$
93-bus	SGSC winter	1 2 10 5	$1,444.57$	$1,097.66$	$897.01$	$1,507.05$	$1,456.26$	$1,360.55$	$1,077.04$
93-bus	SGSC winter	1 2 30 10	$904.47$	$812.52$	$738.34$	$1,636.75$	$1,558.55$	$1,474.41$	$1,077.04$
93-bus	SGSC winter	2 3 10 5	$1,501.19$	$1,302.51$	$1,143.90$	$1,853.96$	$1,764.57$	$1,704.73$	$1,077.04$

Equations

$$\sum _{i\in Z}{x}_i+\sum _{i\in N\setminus Z}(1-x_i)\le ngh.$$

$$z=\sum _{(i,j)}{c}_{ij}{y}_{ij}^p+v$$

$$y_{ij}^p\ge y_{ij}^{p+1}\forall \left(i,j\right)\in {\Omega }_c,\forall p\in \left\{1\dots {\bar{n}}_{ij}-1\right\}$$

$$v$$

$$y_{ij}^p$$

$$v=\sum _{k\in \Gamma }{b}_k{x}_k+\sum _{t\in \Psi }\sum _{k\in \Gamma }{\alpha }_{tk}{r}_{tk}$$

$$\zeta +g_{tk}+r_{tk}-{\beta }_{tk}=d_{tk}\forall t\in \Psi ,\forall k\in \Gamma$$

$$\zeta =\sum _{\left(i,k\right)\in {\Omega }_0}{f}_{tik}^0-\sum _{\left(k,j\right)\in {\Omega }_0}{f}_{tkj}^0+\stackrel{{\bar{n}}_{ij}}{\sum _{p=1}}\sum _{\left(i,k\right)\in {\Omega }_c}{f}_{tik}^p-\stackrel{{\bar{n}}_{ij}}{\sum _{p=1}}\sum _{\left(k,j\right)\in {\Omega }_c}{f}_{tkj}^p$$

$$\begin{array}{r}{f}_{tij}^0-{\gamma }_{ij}{n}_{ij}^0\left({\theta }_{ti}-{\theta }_{tj}\right)=0\forall t\in \Psi ,\forall \left(i,j\right)\in {\Omega }_0\end{array}$$

$$\begin{array}{r}|f_{tij}^p-{\gamma }_{ij}\left({\theta }_{ti}-{\theta }_{tj}\right)|\le M_{ij}(1-{\hat{y}}_{ij}^p)\forall t\in \Psi ,\forall \left(i,j\right)\in {\Omega }_c,\forall p\in \left\{1\dots {\bar{n}}_{ij}\right\}\end{array}$$

$$\left|f_{tij}^0\right|$$

$$\left|f_{tij}^p\right|$$

$$l_{1k}=l_{Tk}+{\beta }_{1k}\forall k\in \Gamma$$

$$l_{tk}=l_{t-1,k}+{\beta }_{tk}\forall t\in \Psi ,\forall k\in \Gamma$$

$$0\le l_{tk}\le x_k\forall t\in \Psi ,\forall k\in \Gamma$$

$$0\le x_k\le {\bar{x}}_k\forall k\in \Gamma$$

$$0\le g_{tk}\le {\bar{g}}_k\forall t\in \Psi ,\forall k\in \Gamma$$

$$0\le r_{tk}\le d_{tk}\forall t\in \Psi ,\forall k\in \Gamma$$

$$f_{tij}^0,f_{tij}^p,{\beta }_{tk},{\theta }_{tk}\mathrm{unbounded}$$

$$\begin{array}{r}v-\sum _{t\in \Psi }\sum _{\left(i,j\right)\in {\Omega }_c}\left[{\pi }_{f_{tij}^{+p}}{y}_{ij}^p\overline{f_{ij}}+{\pi }_{f_{tij}^{-p}}{y}_{ij}^p\overline{f_{ij}}\right]-\sum _{t\in \Psi }\sum _{\left(i,j\right)\in {\Omega }_c}\left({\pi }_{{\gamma }_{tij}^{+p}}+{\pi }_{{\gamma }_{tij}^{-p}}\right)\left(M_{ij}(1-y_{ij}^p)\right)\ge \\ \sum _{t\in \Psi }\sum _{k\in \Gamma }{d}_{tk}{\pi }_{d_{tk}}+\sum _{t\in \Psi }\sum _{\left(i,j\right)\in {\Omega }_0}\left[{\pi }_{f_{tij}^{+0}}{n}_{ij}^0\overline{f_{ij}}+{\pi }_{f_{tij}^{-0}}{n}_{ij}^0\overline{f_{ij}}\right]+\sum _{t\in \Psi }\sum _{k\in \Gamma }\overline{g_k}{\pi }_{g_{tk}}+\sum _{t\in \Psi }\sum _{k\in \Gamma }{d}_{tk}{\pi }_{r_{tk}}+\sum _{k\in \Gamma }\overline{x_k}{\pi }_{x_k}\end{array}$$

Full text

Index terms— hybrid heuristic, optimization, power transmission, energy storage

1 Introduction

The need to solve large scale mixed integer programming problems arises in many applications. In this paper the planning of electricity networks to cope with renewable generation has been used to both motivate the need for the proposed new algorithm and to evaluate its effectiveness.

Integrating renewable energy generation, especially variable generators such as wind and solar, into the electrical transmission network is a considerable design challenge currently facing network planners. For example, a recent blackout in South Australia saw 315MW of wind generation disconnect from the grid amid voltage dips and loss of load [7], recently a 100 mega-watt battery was installed by Tesla to prevent this from happening again [24]. Correspondingly, there has been a renewed interest in electricity network planning problems [20]. One such problem is the transmission expansion planning problem (TEP). In this problem the objective of the planning is minimize the investment and operational costs of the network while meeting a set of operational constraints, for example, generation, demand, geographical, and environmental constraints [27]. Transmission lines require a huge initial investment but have long life whereas storage systems have a short life but can be installed very quickly and be expanded gradually, [39].

Hydro is the most common type of transmission scale Energy Storage System (ESS), but its feasibility is determined by climate, geography, and environmental constraints. Batteries have also been successfully deployed to smooth the 5 minute ramp rate of a wind farm [50] and given their high power and energy capacities compressed air technologies remain viable but expensive [52]. The transmission network expansion and energy storage planning (TESP) model considers a generic ESS, the primary purposes of which are transmission upgrade deferral and demand shifting. Transmission upgrade deferral occurs when need for additional or larger capacity transmission lines is avoided, in this case by using ESS located near sites of generation or demand, to store energy and release it at a steady rate over time. Subtly different, demand shifting entails using stored energy generated in a prior time interval to meet demand in the current time interval. Storage facilities are also considered as decision makers dealing with risks arising from the costs, environmental impact and supply issues [44].

The TEP and related problems are often modeled as a mixed integer nonlinear program (MINLP), or in a correspondent disjunctive mixed integer linear programming (MIP) form. An overview of the standard models and test systems is given in [41].

Advances in commercial solver technology mean that simpler linear models of small networks can generally be solved to optimality within a few minutes. However, a considerable body of research is dedicated to solving larger or more complicated instances. Novel approaches to these problems include branch and bound with a GRASP meta-heuristic [2], Projection-Adapted Cross Entropy [17], and particle swarm optimization [1]. Often TEP problems can be decomposed into investment and operational subproblems. Benders decomposition with alternately continuous or discrete decision variables in the master (investment) problem, and DC approximation or transportation operational subproblems is investigated in [34]. Fencing constraints and additional constraints on new paths are shown to substantially reduce the number of iterations when added to the master problem [21], and adding Gomory cuts evaluated from the master problem to a Benders decomposition of a linear disjunctive MIP model is shown to result in significant CPU time savings [5]. More recently, local branching is used to accelerate the Benders decomposition of a TEP problem in [13]. A useful survey of the literature is given in [46].

Where the modeling becomes more complex and computationally demanding, meta-heuristic approaches have been shown to produce good results. If the transmission expansion planner is concerned only with determining a final network plan, the planning is considered static, whereas dynamic planning involves the determination of one or more intermediate plans over multiple periods. A specialized genetic algorithm is shown to produce good solutions for coordinated, multistage planning problems[16]. A Differential evolution algorithm is used to solve a similar problem in [47].

Another common complication to the planning is (n-1) redundancy. Systems operating under this scheme must not shed load if a single component, in this context a circuit, fails. An adaption of the Chu-Beasley[8] genetic algorithm was used to solve a TEP model with (n-1) security constraints in [12].

Some recent works have considered incorporating energy storage in transmission networks, however the time dimension is largely ignored allowing the ESS to behave as an alternative type of generation [22]. A pair of linear programming models that take into account both variable and dispatchable generation, as well as energy storages are compared in [9].

Locating and sizing small scale ESS in distribution networks has been approached using a modified particle swarm optimization (PSO) to optimize a multi-period design problem[43], and using a genetic algorithm combined with simulated annealing to plan a low voltage network with high solar photovoltaic generation [11].

For solving large scale optimisation problems such as encountered in network planning, there are two categories of methods that are commonly used: (1) meta-heuristics that aim to get good solutions quickly but have no guarantee of optimality, and (2) exact methods such as integer programming and constraint programming that at least in principle produce optimal solution but may take too long to run in practice. Over the last decade there has been an increasing interest in hybrid methods that combined elements taken from both worlds. While a complete review of such methods is beyond the scope of this paper, a few examples show the variety of approaches tried in this category. [38] developed a method they called MetaBoosting to use metaheuristics to improve column generation or cut separation in an integer programming framework. Using Ant Colony Optimisation with Constraint Propagation has been shown to be effective in [48]. Lagrangian relaxation has been combined with Particle Swarm Optimisation [15, 19] or Ant Colony Optimisation [49]. Genetic programming has been combined with MIP though largely as a high-level tuning mechanism for the large number of parameters available in a MIP solver [26]. The rise of these types of methods that, in principle, can operate on any integer linear programming problem, has led to the term matheuristics being used to describe such methods that combine mathematical programming with heuristics [6].

Of particular interest from the point of view of this paper are methods based on combining Benders decomposition with metaheuristic algorithms. Benders decomposition breaks large MIPs into a master problem and one or more subproblems. The subproblems are used to evaluate and test the feasibility of solutions proposed by an optimisation process solving the master problem. Information from the subproblems is also fed back to the master in terms of additional constraints (cuts) generated from the subproblem solutions. [37] created a method based on Benders decomposition where the master problem is always solved with a Genetic Algorithm (GA) while the subproblem is solved using Linear Programming (LP) as in the standard Benders approach. When tested on general MIPs this proved to be more effective than using the Benders method on its own but often less effective than simply using the CPLEX solver without any decomposition. A little earlier [45] independently developed a nearly identical GA-Benders hybrid and applied it to a problem in power generation expansion planning, a similar problem to the one considered here but only looking at power generation investment without any consideration of network transmission or energy storage.

In this paper, we present a hybrid exact/meta-heuristic algorithm that melds Benders decomposition and a Bees Algorithm (BA)[36] inspired approach. Unlike previous proposed matheuristics, this method retains the ability to generate provably optimal solutions. The ideas of multiple bees that have different functions, scouts and workers, is used to balance diversification and intensification in the solution of the master problem. Multiple parallel optimisation processes are used to speed up the search, that learn from each other not just through the exchange of good feasible solutions but also by sharing dual (cut) information obtained when solving Benders subproblems. Using the transmission network expansion and energy storage planning model (TESP) to test the model, we show the Bee-Benders hybrid algorithm (BBHA) to be an effective hybrid algorithm that exhibits equivalent performance to its component parts in the segments of the problem domain where those parts are strongest, and significantly improves upon the individual approaches where neither component part has a pronounced advantage.

The rest of this paper is organized as follows. The Bee-Benders algorithm is introduced in Section 2. A MIP formulation of the TEP with storage model is given in Section 3. Numerical results, in which the algorithm is evaluated using the Brazilian 46-bus and Colombian 93-bus test systems are discussed in Section 4. We conclude in Section 5.

2 The Bee-Benders Hybrid Algorithm

2.1 The Bees Algorithm

Here we present a hybrid exact/meta-heuristic algorithm that combines Benders decomposition with an approach inspired by the Bees Algorithm. There are many variants of optimisation metaheuristics inspired by the behaviour of bees (see for example [23] for a review of one of the alternatives, the Artificial Bee Colony optimisation). Here we will follow the Bees Algorithm as proposed by [36] and [35].

In the most basic form, the algorithm comprises two phases: global search, and local search. A pseudo-code description of this metaheuristic has been provided in Algorithm 1. Each solution in the solution space is referred to as a flower in the terminology of the Bees Algorithm, and the local neighbourhood of a solution is called a flower patch. In the initialisation phase, “scout” bees leave the hive and fly to a random flower. The fitness of the flower is evaluated and the scout bees return to the hive. During the local search phase, the scouts who discover the $ne$ elite and the $nb$ best flowers (solutions) recruit “worker” bees to explore their respective flower patches, that is, the flowers in the neighbourhood of those the scouts discovered. Recruited worker bees fly to a random flower within the flower patch and evaluate its fitness. The fittest flower from the elite and best flower patches are combined with the fittest new flowers discovered scouts during the subsequent global search phase to produce a new pool of elite and best solutions for further local exploration. Stopping conditions may include time, the number of iteration, or a test for convergence.

This can be thought as a multi-start local search algorithm which always works on a subset of best known solutions, with more effort expended on the local search in the neighbourhood of the elite solutions than the remaining solutions. It should be noted that this algorithm can be parallelized in a fairly straight forward manner by carrying out the search in each flower patch (neighbourhood of an elite or best solution) in parallel. Also the “scouting”, that is, generation of new random solutions, can be carried out independently.

The BA has been applied to numerous combinatorial optimization problems such as the generalized assignment problem[32] and machine scheduling [35], and as also shown value for applied industrial applications such as crack detection of beam-type structures [31].

We have chosen to combine the BA with Benders decomposition because it has been shown to perform at least as well as standard evolutionary approaches, to be less sensitive to tuning parameters than other swarm approaches such as PSO, and yet retains an extraordinary simplicity of implementation [36].

2.2 Benders Decomposition

Benders decomposition is a technique that allows a large, intractable problem, such as the TESP model described in Section 3, to be divided into more tractable component parts [3]. The first part is called the master problem and consists of a MIP that includes all of the integer variables and any applicable continuous variables. The second part consists of one or more subproblems, that collectively contain the remaining continuous variables [18]. Typically both master and subproblem(s) have not only significantly fewer variables than the original problem but also far fewer constraints, making these much easier to solve using a MIP solver. The master problem is solved to yield a candidate solution which is used to fix the complicating variables that would otherwise be present in the subproblem. The dual solution of the subproblem is used to produce a feasibility or optimality cut to be added to the master, and the master problem is solved again. This iterative procedure continues until it no further cuts are necessary. Note that in each iteration a candidate solution is evaluated to determine if it is both feasible and better than any seen previously. The addition of the cuts ensures that the master problem does not repeatedly generate the same candidate solution. For problems where only optimality cuts are possible, as is the case in this paper, the collection of cuts represent a piecewise linear approximation to the objective function. This approximation is iteratively refined and improved as more cuts are added.

Benders decomposition has been applied to numerous optimization problems such as the fixed charge network design problem[10], the unit commitment problem[28], the network-constrained unit commitment problem [51], and the scheduling of crude oil in an oil refinery [42]. It is also proven quite effective on multi-stage stochastic energy planning problems [33, 40] and CCHP-microgrid operation involving battery storage [30].

2.3 The hybrid method

The Bee-Benders hybrid algorithm (BBHA) is a hybrid of Benders decomposition and a local search phase that is largely based on the Bees Algorithm. The algorithm operates on large MIP that has been decomposed in a manner suitable for Benders decomposition. The master problem contains binary variables representing certain investment decisions, and an LP subproblem containing continuous variables and largely operational constraints. The particulars of the mathematical model detailed in Section 3.

2.4 The BBHA in detail

As with the BA, the algorithm comprises global search and local search phases. The global search phase commences as a conventional Benders decomposition using a “single tree” master approach. Lazy constraint callbacks are used to separate Benders cuts, as opposed to solving the master problem to optimality at each iteration. This single tree branch and bound search fulfills the role of the “scout” bee in the BA algorithm, ensuring that eventually the whole solution space is searched and the method can never remain stuck in a local optimum. Meanwhile, an initial set of random solutions is generated for exploration during the local search phase.

During the local search phase, “worker” bees (henceforth known as “workers”) explore the local neighbourhood (subsequently referred to as a “site”) of each solution, by estimating the fitness a subset of solutions using a matheuristic based on the set of known Benders cuts. The most heuristically promising solution discovered at the site is selected for full evaluation of the LP.

As with the BA, the fittest solution from both the elite and best sites are combined with the incumbent solution of the Benders decomposition to produce a new pool of elite and best solutions for further local search. The algorithm iterates in this way until stopping condition is met, or the Benders decomposition finds and proves the optimal solution. A pseudo-code description of the BBHA is given in Algorithm 2. It should be noted that while the algorithm description is essentially serial, the intention is for each “bee” to be executing as a parallel process that is carrying out its search independent of the other processes. The BBHA is discussed in greater detail below.

2.4.1 Initialization

The algorithm is initialized with a population of $nre+nrb$ workers, which are uniformly randomly distributed over the solution space. The fitness of each solution is evaluated by solving the LP subproblem. Each LP subproblem produces a Benders cut which is stored in a pool of cuts shared by the Benders decomposition. The fitness scores are ranked and the $nb$ best “flower patches” are selected for neighbourhood search. The algorithm enters the main loop.

Simultaneously, the algorithm commences solving the Benders decomposition using the “single tree” master problem approach: Lazy constraint callbacks are used to solve the LP subproblem and separate the cuts. This means that the master problem need only be solved to optimality once as opposed to once per iteration. Any generated cuts are add to the shared pool of cuts.

2.4.2 The main loop

The main loop consists of two main phases: neighbourhood search and cut sharing. The neighbourhood search is carried out by each process independently, while the cut sharing represents a communication or synchronisation step between the processes.

2.4.3 Neighbourhood search

Each iteration, the workers that discovered the $ne$ elite solutions each recruit $nre$ workers for neighbourhood search. Likewise, the workers who discovered the remaining $nb-ne$ best solutions each recruit $nrb$ workers for neighbourhood search. The $ns-nb$ workers who failed to find a best solution rejoin the pool of workers.

Neighbourhood search at a given site is performed by each worker producing a pool of candidate solutions using a Hamming distance function which randomly selects at most $ngh$ binary variables to alter. When solving arbitrary MIPs with binary variables $x_{i}$ for $i\in N$ this is equivalent to imposing the constraint

[TABLE]

That is, of all the variables $x_{i},\ i\in Z$ which have $x_{i}=0$ in the current solution and all the remaining variables that have $x_{i}=1$, only $ngh$ may change their value. In our application a more specialised neighbourhood move can be defined based on the structure of the problem. We have binary variables that represent a unit increment in transmission capacity between two locations (a right of way). Each right of way (analogous to a set of edges $ij$ on a multigraph with nodes $i$ and $j$) has $p$ binary variables denoting the installation of a equivalent line. This means that individually installing the $1^{st}$ line is equivalent to installing the $2^{nd}\dots p^{th}$ line. Clearly it is undesirable for the Hamming distance function to randomly replace the installation of one line on a right of way with another. For this reason the function operates on groups of binary variables representing a single right of way. If two or more changes are made to a right of way they are directionally consistent i.e. if the first change added a line, subsequent changes will also add a line until the maximum of $p$ lines are installed.

The fitness of each candidate solution in the pool is estimated using the matheuristic given in Algorithm 3. Here the cost of the master problem is calculated from the candidate solution. If the cost exceeds the incumbent fitness value the evaluation stops. Otherwise, the shared set of Benders cuts are used to estimate the cost of the subproblem. The piecewise linear approximation of the original problem objective is of the form $c^{T}y+\max_{i\in Cuts}\{r_{i}-B_{i}y\}$. Here $y$ are the master problem (binary) decision variables with cost vector $c$. Each Benders cut takes the form $\hat{v}\geq r_{i}-B_{i}y$ where $B_{i}$ is a row vector and $r_{i}$ a constant. Collectively these cuts enforce that $\hat{v}\geq\max_{i\in Cuts}\{r_{i}-B_{i}y\}$ with minimisation of $\hat{v}$ ensuring that at least one of the inequalities holds with equality. Hence if $v=r-By$, the estimated fitness of the candidate solution comprises the cost of the master problem plus the maximum value of the vector $v$.

Each worker then solves the LP subproblem for the most promising heuristically determined solution in their solution pool, and the generated Benders cuts are added to the shared pool of known cuts.

The fitness scores of the solutions found by the workers are combined with the incumbent solution of the Benders decomposition and are ranked from best to worst. The $nb$ best solutions are selected for neighbourhood search during the next iteration.

2.4.4 Cut sharing

At the conclusion of the neighbourhood search phase any Benders cuts produced by the Benders decomposition are added to the shared pool of cuts. Any cuts produced by the workers are likewise made available to the Benders decomposition, and may be added to the pool of cuts managed by CPLEX during a subsequent execution of the lazy constraint callback. The effect of this cut sharing is that each of the bees has a more accurate approximation of the objective function. With this approximation the scout bee (branch and bound tree) avoids searching any solution that is not at least potentially better than the best found so far. The worker bees use the approximation to quickly evaluate solutions in the site to identify the most promising solution for which a full LP is solved.

2.4.5 Termination

The algorithm may terminate in several ways: After $n^{max}$ iterations, $t^{max}$ seconds, or if the Benders decomposition identifies and proves the optimal solution. Note that even if the time or iteration limit does not allow a provably optimal solution to be found, we are still able to extract a lower bound from the branch and bound tree that the scout bee was searching. Thus even in the case where the method is only a heuristic, we have an estimate of the maximum gap to the globally optimal solution.

3 Mathematical model of TESP

We consider a electricity network consisting of nodes and arcs (referred to as “rights of way”). The ability to generate power and demand for power occurs at the nodes over a number of time periods. The objective of the complete TESP model is the is to minimize the investment cost of expanding the transmission network while simultaneously minimizing a penalty for load curtailment at nodes with net demand. A discrete number of new or reinforcing circuits may be installed on each right of way, and the location and size of any ESS at the nodes are determined.

Cyclic discrete time is used to model the period of operation, and therefore the state of any installed ESS in the last time interval must be identical to the state in the initial time interval. This might model the typical power use over a day with the end of one day matching the start of the next. Generation is re-dispatchable and demand may vary between time intervals. Despite the introduction of time to the model, the planning is static, and only a single final expansion plan is produced. More complex models, for example with multiple scenarios for demand and renewable power generation capacity, are possible but not considered in this paper.

The model determines the network expansion plan, and operational characteristics such as the amount of energy stored in the ESS, the network flows, and the phase angles at each bus for each time interval. As with other variants of the disjunctive TEP, power flows are modeled using a DC approximation [25, p.36].

The mathematical model presented here, as well as alternative modelling approaches in the literature, is discussed in detail in [29]. As such, only an abridged discussion of the decomposed model follows. The key point to note here is that the main integer (binary) variables to be determined relate to the rights of way to be installed, while a very large number of continuous variables and associated constraints have to be considered to determine the optimal generation, storage and power flows for any choice of network expansion.

The following notation will be used throughout this paper to define the TESP:

*Sets

$\Gamma$ the set of indices for buses;

$\Omega_{0}$ the set of rights of way for existing circuits;

$\Omega_{c}$ the set of rights of way for candidate circuits;

$\Psi$ the set of uniform time intervals $\left\{1,2,\ldots,T\right\}$;

*Parameters

$\alpha_{tk}$ cost of curtailment at time $t$ at bus $k$;

$b_{k}$ cost of installing storage at bus $k$;

$c_{ij}$ cost of installing a circuit on right of way $ij$;

$d_{tk}$ demand at time $t$ at bus $k$;

$\bar{f}_{ij}$ maximum possible power flow on right of way $ij$;

$\bar{g}_{k}$ maximum possible generation at bus $k$;

$\gamma_{ij}$ susceptance of circuits installed on right of way $ij$;

$M_{ij}$ the disjunctive parameter for right of way $ij$

$n^{0}_{ij}$ number of existing circuits on right of way $ij$;

$\bar{n}_{ij}$ maximum number of installable circuits on right of way $ij$;

$\bar{x}_{k}$ maximum installable storage capacity at bus $k$;

$\hat{y}^{p}_{ij}$ binary parameter denoting installation of the $p$th candidate circuit on right of way $ij$;

*Decision variables

$\beta_{tk}$ power flow to storage at bus $k$ at time $t$;

$g_{tk}$ generation at time $t$ at bus $k$;

$f^{0}_{tij}$ power flow for existing circuits at time $t$ on right of way $ij$;

$f^{p}_{tij}$ power flow for the $p$th candidate circuit at time $t$ on right of way $ij$;

$l_{tk}$ level of storage at bus $k$ at time $t$;

$r_{tk}$ demand curtailment at time $t$ at bus $k$ ;

$\theta_{tk}$ phase angle at time $t$ at bus $k$;

$x_{k}$ storage capacity installed at bus $k$;

$y^{p}_{ij}$ binary variable denoting installation of the $p$th candidate circuit on right of way $ij$;

$v$ estimate of the contribution of the subproblem to the objective function of the master problem

3.1 The master problem

The objective of the master problem is to minimize the function

[TABLE]

where $c_{ij}$ is cost of installing a line on right of way $ij$ and $y_{ij}^{p}$ is a binary variable denoting the installation of the $p$th candidate line on $ij$. The estimated contribution of the subproblem to the objective function is given by $v$.

The following constraints are necessary to the master problem:

Symmetry breaking constraints

[TABLE]

The lexicographical constraint (3.2) eliminates the symmetry of the binary decision variables by mandating the order of installation of parallel circuits be arbitrary.

Other constraints

[TABLE]

3.2 The subproblem

Given a set of new circuit installations determined by the master problem, the subproblem determines the cost of any installed ESS, and a penalty for load curtailment.

The objective of the subproblem is to minimize the function

[TABLE]

where $b_{k}$ is the fixed cost of installing $x_{k}$ MW of storage at bus $k$, and $\alpha{tk}$ the cost of curtailing $r_{tk}$ in each time interval $t$. It is assumed that the variable operating cost of ESS is low relative to fixed costs, and these are therefore omitted from the objective function.

The following technical constraints govern the operation of the network:

Nodal balance and power flow

[TABLE]

where

[TABLE]

Nodal balance i.e. Kirchhoff’s current law is ensured for each time interval by constraint (3.6).

Power flows are modeled using a DC approximation requiring that the phase angle at each bus be determined for each time interval:

[TABLE]

[TABLE]

Kirchhoff’s voltage law is implemented for existing circuits by (3.8), and for candidate circuits by (3.9). Note that absolute values are given to simplify the notation, however in practice these are readily expanded into pairs of ranged linear constraints. The disjunctive parameter $M_{ij}$ must be large enough that it does not limit the difference in phase angles of buses $i$ and $j$. Minimal values of $M_{ij}$ may be calculated by following the procedure given in [5].

[TABLE]

Constraint (3.10) and constraint (3.11) enforce nominal thermal limits on existing and candidate circuits respectively.

Storage level and charge/discharge limits

[TABLE]

[TABLE]

The set of time intervals $\Psi$ is assumed to be cyclic to allow the operation of the storage throughout the desired time period, for example, a typical day. As such, the storage level at the end of the day is required to match the initial storage state. The “wrap around” constraint (3.12) implements this requirement. For all other time intervals the storage level is determined by (3.13).

[TABLE]

[TABLE]

Constraint (3.15) establishes bounds on the installable storage capacity at bus $k$, while constraint (3.14) ensures the stored energy does not exceed the installed capacity.

Generation bounds

[TABLE]

Constraint (3.16) imposes bounds on generator re-dispatch.

Curtailment bounds

[TABLE]

Load curtailment at any bus $k$ cannot exceed demand during the same time interval $t$..

Other constraints

[TABLE]

3.3 Optimality cut

As noted above, load curtailment is permitted at any bus during any time interval so long as it does not exceed demand at that bus during the same time period. Therefore, the dual of the subproblem remains bounded for any feasible solution to the master problem. Accordingly, we need only consider the following optimality cut:

Let the dual variables $\pi_{d_{tk}}$ correspond to constraint (3.6), $\pi_{\gamma_{tij}}$ to constraint (3.8), $\pi_{\gamma^{+p}_{tij}}$ and $\pi_{\gamma^{-p}_{tij}}$ to constraint (3.9), $\pi_{f^{+0}_{tij}}$ and $\pi_{f^{-0}_{tij}}$ to (3.10), and $\pi_{f^{+p}_{tij}}$ and $\pi_{f^{-p}_{tij}}$ to (3.11). The dual variables $\pi_{s_{tk}}$ correspond to constraints(3.12) and (3.13), and $\pi_{\bar{l}_{k}}$ to (3.14). Lastly, let the dual variables $\pi_{g_{tk}}$, $\pi_{r_{tk}}$, and $\pi_{x{k}}$ correspond to the bounds (3.15 - 3.17) respectively.

The optimality cut is therefore

[TABLE]

3.4 Limitations

The relative simplicity of this TESP formulation comes at the cost of addressing certain features of a real world electrical transmission network. The most obvious limitation is that power flows are modelled using a DC approximation to the AC power flow of most transmission networks.

It is also assumed that the variable operating cost of ESS is negligible, at least relative to the fixed cost of installing and operating the ESS over its lifetime, and that fixed costs increase linearly with capacity. Furthermore, power flow to and from ESS is limited only be the total capacity and current level of the storage. The model allows that the storage completely charge or discharge within a single time interval. Furthermore, the model assumes 100% efficiency for storage and losses are not considered.

The model also assumes generator re-dispatch does not incur cost, and that generators are not subject to technical constraints such as generation ramp rates.

While it is possible to address these and other limitations of the model with additional variables and constraints, these come at the cost of significant complexity in both notation and implementation. Here we have sought to balance to the realism of the modelling with the intent to use the model simply to demonstrate the use of the algorithmic approach.

4 Numerical results

In each of the numerical experiments described in this section the model is implemented in Python 3.4.3 and, where appropriate, makes use of the Python library for IBM ILOG CPLEX 12.6.3. Parallelization is achieved using multiple processes, not threading. The Benders decomposition is implemented with a “single tree” master using lazy constraint callbacks. Preprocessing is disabled by default, and while the LP solver may take advantage of multi-threading, the branch and cut is single threaded.

4.1 Parameter tuning

There are a number of parameters to the BBHA algorithm which may be tuned to find a set of default values that empirically demonstrate good performance. The tuneable parameters are given in Table 1.

The IEEE-25 bus test system is used to benchmark combinations of parameters presumed likely to perform well. A schematic and tabulated data are available in [14]. The system has 25 buses and 36 rights of way with a total demand of 2750 MW. Without storage, and permitting a maximum of 4 new or reinforcing circuits on each right of way, the cost of the optimal expansion plan is US$107.7 million.

While it would be preferable to incorporate real world storage costs into the model, the cost per MW of long term energy storage technology is currently high enough to prevent the installation of any storage in the test systems discussed in this paper. Therefore, an arbitrary cost coefficient of US$2000MW/h is used for each network to ensure storage is installed.

Under the long peak scenario shown in Figure 3 the cost of the optimal expansion plan is US$43.8 million. This result is the benchmark objective for the parameter tuning.

In this tuning exercise, 34 sets of parameters are compared over the first 1800 seconds (30 minutes) of the optimization. The results are given in Table 2. We use a composite trapezoidal rule to integrate along the time axis and then rescale against the worst (largest) integral (scaled trapz). The parameter sets are ranked then ranked. The best result is that with lowest value. Only 5 sets of parameters ([ne: 1, nb: 2, nre: 30, nb: 10], [ne: 1, nb: 2, nre: 30, nb: 15], [ne: 2, nb: 3, nre: 20, nb: 10], [ne: 3, nb: 4, nre: 10, nb: 5], and [ne: 3, nb: 4, nre: 20, nb: 15]) find the optimal solution within the 30 minute window. The timeseries of the incumbent value of the best of these parameter sets is plotted with the best and worst parameter sets in Figure 1.

In general, parameters sets with a relatively modest number of workers and associated high number of iterations appear to do well. An exception is the parameter set (ne: 1, nb: 4, nre: 20, nb: 15) which requires only 39 iterations to match the best sub-optimal objective function value. This is explained by the proximity of the 3 elite search neighbourhoods, and subsequent thorough exploration of the combined neighbourhood. Other similar parameter sets that match this objective function value by the end of the optimization do not converge as quickly, as evidenced by their larger scaled trapz scores.

The Benders scout ensures that the BBHA is guaranteed to find the exact optimal solution to the problem given sufficient time to run to completion. Of course this may take a significant amount of time. The objective of the BBHA is to discover high quality solutions quickly, and as such we favour parameter sets which rapidly converge to such solutions in the case studies that follow. We explore the performance of three sets of parameters:

• [ne: 2, nb: 3, nre: 10, nrb: 5]: the parameter set with the smallest scaled trapz measure.
• [ne: 1, nb: 2, nre: 30, nrb: 10]: the parameter set that converges to the optimal solution the fastest using the scaled trapz measure.
• [ne: 1, nb: 2, nre: 10, nrb: 5]: the parameter set with the largest number of iterations.

The final parameter to consider is the size of the neighbourhood for local search $ngh$, described in Subsection 2.4.3. A histogram showing the distribution of the hamming distance over the range 1-10 required to produce the best improved solution of nearly 14000 workers is shown in Figure 2. A value of 2 accounts for the largest number of improved solutions. This is perhaps unsurprising as it reflects the somewhat routine circuit swap in which one circuit is deselected and another selected. Given that the long tail of larger hamming distances typically resulted in improved solutions only at the beginning of the optimization the value of $ngh$ was reduced to 8 for the case studies.

4.2 Case study: 46-bus network

Representing the southern part of the Brazilian transmission network, the 46-bus test system comprises 46 buses and 79 rights of way. Total demand in the network is 6880MW. Tabulated data are provided in [21]. The investment cost of the optimal expansion plan without ESS is US$154.42 million.

In this case study we allow the installation of a maximum of 5 new or reinforcing circuits on each right of way. Storage may be installed at any bus at an arbitrary cost of US$2000MW/h.

As the amount of storage installed depends upon the demand scenario under which it is operated, four demand scenarios are considered. The short peak and long peak scenarios are described in [29], and the Smart Grid, Smart City (SGSC) residential winter and summer scenarios are generic load profiles taken from [4]. Each scenario describes a 24 hour period with a 30 minute time step, and is shown in Figure 3.

Each scenario is optimized $N=5$ times for both the BBHA and Bees algorithm, and once using Benders decomposition which as a deterministic method exhibits little variance. Each optimization is limited to 4 hours. Tabulated results are given in Table 3.

The parameter set [ne: 1, nb: 2, nre: 10, nrb: 5] typically matches or exceeds the mean performance of the other parameter sets under investigation for the BBHA, whereas the parameter set [ne: 1, nb: 2, nre: 30, nrb: 10] exhibits better performance for the BA. The BBHA finds the optimal solution for the short peak and SGSC summer and SGSC winter demand scenarios, and the Benders scout is able to prove optimality. This is also true of the Benders decomposition run. The range of incumbent solution values over time are shown in Figures 5, 6 & 7 for the short peak, SGSC summer and SGSC winter scenarios respectively.

For the long peak scenario the best BBHA runs find the optimal solution, but optimality is not proven. However, it is possible to use the pool of generated cuts to prove a lower bound if necessary. A plot of the range of incumbent solution values over the duration of the optimization is given in Figure 4.

4.3 Case study: 93-bus network

The Colombian 93-bus network is a medium complexity transmission network with 93 buses and 155 possible rights-of-way. The planning horizon includes 3 discrete stages making this test system useful for testing multi-stage optimization techniques [16]. In this case study will consider only the total demand of 14559 MW in the final stage of the planning horizon.

A maximum of 4 new or reinforcing circuits is permitted to be installed on each right of way. As with the previous case storage may be install at any bus at an arbitrary cost of US$2000MW/h. Network expansion plans are optimized for the long peak, short peak, SGSC summer, and SGSC winter scenarios over a 4 hour period. Tabulated results are included Table 3.

For this test system the parameter set [ne: 1, nb: 2, nre: 10, nrb: 5] exhibits consistently good performance for both the BBHA and BA. The BA achieves a lower mean for the long peak scenario. The Benders decomposition tends to lag behind both approaches for all scenarios except the SGSC winter demand profile, as shown in Figures 12 & 13.

4.4 Discussion

The BBHA exhibits the essential characteristics of a hybrid optimization method. Where the problem is readily solved by one of the component optimization methods the BBHA performs comparably at minimum. Where each component optimization method performs similarly on a given problem, the hybrid approach exceeds this individual performance. In short, the whole is greater than the sum of its parts.

Figure 5 shows the 46-bus test system with the short peak demand scenario, a problem known to amenable to Benders decomposition. The BBHA performs comparably to the Benders decomposition, and in most runs discovers the optimal solution earlier. This can be observed by the incremental improvements to the incumbent value over the first 1000 seconds. Both methods are able to prove optimality within the time limit, however in this case the BBHA takes longer. For the SGSC summer and SGSC winter scenarios shown in Figures 6 and 7 respectively, the BBHA not only discovers the optimal solution heuristically well in advance of the Benders decomposition, but is also able to prove optimality prior.

In the case of the 46-bus test system and long peak scenario shown in Figure 4, the best BBHA run discovers the optimal solution but is not able to prove it optimal within the time limit. The mean solution is slightly better than the incumbent of the Benders decomposition, however the worst BBHA solution is 109% of the Benders decomposition incumbent.

Where the problem favours the new approach, such as for the 93-bus problem under the long peak scenario shown in Figure 8, there is little discrepancy between the ranges of the BBHA and the BA although the BA has better mean performance. Note: The Benders decomposition incumbent value does not fall within the plotted range.

Like any other hybrid approach the BBHA is a compromise. A straight Benders decomposition implementation running on the same computing infrastructure will evaluate more of the search tree than the BBHA scout. Likewise, without the continuously running scout or the trade off between producing Benders cuts and local search the BA approach can dedicate more cores to evaluating candidate solutions. As a result the straight Benders approach tends to outperform the BBHA in terms of the lower bounds it produces. However since these are generally of much less interest than finding good feasible solutions, and because the Bees algorithm is unable to produce any lower bounds, we have restricted ourselves to only compare upper bound solution quality in this analysis.

However, empirically we have found the benefits of cut sharing largely negate any compromise. In the first instance, the cuts generated by the Benders scout improve the heuristic estimate of the fitness of the candidate solutions in the worker solution pool. Likewise, the cuts generated in parallel by the elite workers are typically in the neighbourhood of the incumbent solution they prove useful to the Benders decomposition. Perhaps the clearest example of this is shown in Figure 7. Here, by sharing cuts between workers and the simultaneous Benders decomposition, the BBHA is able to prove the optimal solution faster than Benders decomposition alone, even though the Benders decomposition has a resource advantage on the compute infrastructure. The effect is also evident to a lesser extent in Figure 6.

Figure 9 shows the random effects of local search with few workers, and the responsiveness of the BBHA to “good” Benders cuts. These runs display a very large range of incumbent values because the worst of the runs was unable to fully exploit the cut sharing. We observe similar random effect in Figure 11. Here the BBHA shows sensitivity to the parameters, and by increasing the workers available for local search (parameter set: [ne: 1, nb: 2, nre: 30, nre: 10]) more of the search space is evaluated each iteration and we observe less variance and broadly better. Results across the 5 sample runs shown in Figure 10.

Increasing the number of workers also significantly improves the optimization of the 93-bus test system under the SGSC winter demand scenario shown in Figures 12 (fewer workers) and 13 (more workers). In this case the entire range of BBHA objective values improve upon that of the Benders decomposition by the end of the optimization.

As noted in Section 4.3, although the BBHA is a hybrid matheuristic optimization technique, the use of Benders decomposition ensures that the solution can be proven optimal if the algorithm is allowed to run for sufficient time. If not allowed to run until completion, the library of cuts may be used to produce a valid lower bound.

5 Conclusion

In this paper, we introduced a hybrid exact/meta-heuristic algorithm that combined Benders decomposition and a Bees algorithm inspired approach. To the best of our knowledge this is the first such matheuristic based on Benders decomposition and the Bees algorithm.

The BBHA approach was demonstrated using a transmission network expansion and energy storage planning model that is known to become more tractable when decomposed into investment and operational subproblems.

The approach as been shown to combine the best performance of its component parts in the segments of the problem domain where those parts excel, and to improve upon the individual approaches where neither shows a substantial advantage.

As the BBHA is general in nature and does not require any special problem structure beyond the decomposition required by the Benders decomposition method. The approach may be applied to any general decomposable mixed integer programming problem.

Future work will include a stochastic BBHA variant with multiple probabilistically weighted subproblems. For problems such as the TESP which rapidly become intractable as complexity and realism increases, a highly parallelized stochastic algorithm is expected to significantly advance solution quality.

Acknowledgment

Melih Ozlen is supported by the Australian Research Council under the Discovery Projects funding scheme (project DP140104246).