Number of power lines in operation
Number of buses in the system
Number of loading scenarios representing given time frame
Number of scenarios representing planning horizon
Number of time frames representing horizon
Index of a decision point
Index of a scenario at time frame t
Occurrence probability of a scenarioat time frame
Vector of initial line inductances
Vector of maximum (minimum) active power generator outputs
Vector of maximum (minimum) reactive power generator outputs
Vector of active (reactive) power demands
Vector of line apparent power limits
Vector of maximum (minimum) allowed voltages
Cost per Ohm of a series FACTS device
Cost per MVAr of a shunt FACTS device
- Optimization variables (operational, scenario dependent):
Vector of bus voltage magnitudes
Vector of bus voltage angles
Vector of generator active power injections
Vector of generator reactive power injections
Vector of line inductances modified by SC devices
Vector of series FACTS settings
Vector of shunt FACTS settings
- Optimization variables (investment, scenario independent):
Vector of series FACTS capacities built at decision point t
Vector of shunt FACTS capacities built at decision point t
Energy market deregulation and massive installation of renewables are among significant stress factors forcing transmission systems across the world to operate closer to their limits. When a transmission system is constrained building new lines seems a natural remedy [1, 2, 3, 4]. However, this option is costly and severely limited in many countries due to social and environmental concerns, hence rising the question: if the system can be upgraded creatively and gracefully by installing Flexible Alternating Current Transmission System (FACTS) devices which are both less expensive and whose spatial footprint is much smaller?
This question was addressed in the literature [5, 6, 7, 8, 9] suggesting that FACTS devices may indeed be effective in increasing transmission capacity and keeping the system safe and operational as the demand growth. However, it also became clear from the studies that this sustainable upgrade needs to be creative – types, locations and capacities of the newly introduced FACTS devices must be chosen carefully in order to exploit their benefits. Several objectives, including decreasing operational [7, 10, 11] and investment costs [5, 9], reducing transmission losses [6, 10], increasing power system loadability and managing system congestion [12, 10, 13, 14], reducing load curtailment  and improving voltage profile  and voltage stability index [8, 14], have been considered. Formally, these problems were stated as mixed Integer Nonlinear Programming (MINLP) optimizations which exact solutions are limited to only very small (typically not exceeding tens of nodes) power system models. Sensitivity analysis [15, 16], Mixed Integer Linear Programming (MILP) 18, 19, 20, 21] were suggested to resolve the MINLP formulations approximately. The sensitivity based methods consist in computing a few indicators to identify the most critical lines which are thus suggested as good option for FACTS placement. However, this ad-hoc methodology is obviously not optimal and it also does not offer any suggestion on how to size the devices. Genetic algorithms aim at finding optimal solutions but come with unacceptably high computational cost. Relaxation and approximation techniques were suggested to convert MINLP to MILP . Approximation techniques are computationally less expensive as they use the simplified line flow based model  or DC power flow model . However, such approximations have serious limitations especially for planning installation of shunt FACTS devices.
One significant complication in resolving MINLP is related to combinatorial explosion in the number of choices for possible FACTS placements. It is thus desirable to enforce sparcity of the installation. A principal way to achieve sparse placement in a computationally feasible way was suggested in [23, 24]. An alternative way to achieve sparse placement was also discussed in .
Uncertainty of the economic growth is another notable obstacle for efficient and practical implementation of this and other related planning problems. A solution for resolving this problem was suggested  where the uncertainty was modeled via exogenously described multiple operational samples of loads and related nonlinear AC power flow solutions was embedded in the large-scale optimization explicitly.
In this manuscript we extend the approach of  to resolve the last remaining issue in the optimal multiple-scenario aware, AC-based and sparse placement and sizing of the FACTS devices in a large transmission system. Here, we choose to represent future not in one time step, as was done in , but in multiple time steps. In other words here we complete general description of the framework started in [23, 24, 25] and propose a comprehensive resolution for finding optimal locations of FACTS devices in a large transmission system by preparing the system for future loading gradually through multi-stage, properly paced investments. Main highlights of our comprehensive approach are as follows:
Planning horizon is represented by multiple decision points (multiple time frames). At each new time frame a set of new FACTS devices can be installed, and they are assumed available for operations immediately such that respective operational values do not exceed the installed capacities. Therefore, installation of FACTS devices can be paced. In this manuscript we work with a finite number, , of the time-horizon sub-intervals. Notice, however, that extension of the approach to the case of a receding horizon, thus accounting at each step for updated forecast, is straightforward. (We plan to conduct extensive discussion of the receding horizon experiments in a future journal version of the manuscript.)
Future operational conditions are represented through multiple loading scenarios and associated probabilities broken into time frames. Our framework is set in the way that the scenarios are stated as an exogeneous input, which allows us to separate the problem of scenario generation from the intrinsic optimization details. Given the exogeneously prescribed scenarios, optimal installation of FACTS is resolved within the optimization framework by accounting for both investment variables and operational variables, characterizing installation decisions and operational implementations (per scenario) respectively. It is important to stress that an optimal constraints (i.e. feasible) all the scenarios. (This is in a contrast with the worst case planning approach.)
Both capital and operational expenditures are optimized simultaneously. To the best of our knowledge no prior works have considered optimizing them at the same time. But for practical planning horizons operational cost is much bigger than the cost of FACTS installation. Thus relatively small additional investment allows to save a significant amount of money by reducing congestion additionally to resolving infeasibility of particular load scenarios.
A novel optimization iterative heuristics is developed, which is a combination of analytic linearization of non-linear constraints, a solution of Quadratic Programming (QP) or Linear Programming (LP) (depending on generation cost) problem for finding of investment variables and operational settings for all scenarios and Alternating Current (AC) Power Flow (PF) resolution (for each scenario) to update previously found states.
The developed heuristics for finding optimal locations of FACTS devices considering multiple loading scenarios and multiple time frames can be applied to large power systems. In other words, the developed approach is scalable. To the best of our knowledge the system considered in the literature addressing optimal locations of FACTS placement consists of a maximum of 1228 buses . The results of the proposed methodology is demonstrated on 2736-bus Polish system, thereby proving its scalability. Moreover, algorithm provides upper bound solutions of the objective function with the gap less than 0.1
The material in the manuscript is organized as follows. Our basic optimization model is introduced in Section II. The solution algorithm is described in Section III. Section IV is at the core of the manuscript - it describes in details our experiments. We conclude and discuss path forward in Section V. Appendixes describe constructions and steps needed to run our algorithm. (This material is largely a repetition from  reproduced here for completeness.)
Ii Optimization model
This section describes our optimization framework for operations-aware installation of FACTS devices taking into account multiple future decision points (or multiple time intervals).
Assume that the planning time horizon is , is the number of time intervals, is number of given loading configurations (scenarios) per each time frame (the number of scenarios per time frame may also vary with the time frame). In this setting we aim to place and size the Series Compensation (SC) and Static Var Compensation (SVC) devices, where an SC device, installed at a line, modifies inductance of the line (thus allowing to reroute apparent power), while an SVC device, installed at a node, injects or consumers reactive power at the node thus helping to balance the voltage locally.
Fig. 1 illustrates the setting. Since scenarios are generated within each time interval independently, the total number of paths accounted for within our optimization formulation is , where a path is a sequence of scenarios (each per time interval). Notice that even though the number of paths is exponential in , the total number of the operational constraints in the optimization formulation, , scales linearly in .
The overall problem is to minimize a combination of the sum (over the time intervals) of the investment cost and the sum of operational costs over all the scenarios taking into account (a) operational constraints for every scenarios (per time interval) and (b) investment constraints requiring that the operational variables (for every scenario per tiem interval) do not exceed the respective installed capacities. Operational settings can be different for different scenarios but installed capacities of the devices are the same for all the scenarios representing given time interval.
Mathematically the optimization problem is stated as follows:
where labels the scenarios; upper index labels the time intervals, ; and denotes the set of nodes and the set of (undirected) edges, of the grid-graph, where a node can be of the load type, , or of the generator type, .
The objective function in (II) consists of three terms. The first two, sparsity promoting terms [23, 24], express the capital investment costs of the installation of the two types of FACTS devices (investment can be performed at each decision point ). The third term stands for the operational cost in which the summation is over all the scenarios for each time frame and over the time frames ( is a shortcut for, ) accounting for respective occurrence probability multiplied by the number of years (service period). Therefore, the optimization (II) is nothing but an operational aware planning.
) represent total available capacity at the decision moment. (5) and (6) ensures that the already installed capacities are inherited in the future time frames. (7) bounds actual line inductances, which are adjusted according to the operational value of the installed series compensation for each scenario, within their respective installed capacities (represented by (14)). (8) and (9) represent active and reactive power balances at each bus of the network. Components of the vectors () and () are assumed equal to zero at the buses containing, respectively, no generators or loads. (10) and (11) represent the net active () and reactive () power injections at the system buses. The term expresses shunt compensation by SVC adjusted to a scenario bounded by the respective installed capacities (represented by (15)). Active and reactive power generation limits are set by (12) and (13). Voltage and thermal line flow constraints are represented by (16) and (17). and stand, respectively, for the apparent power flows from to and to from along the line . One also accounts for the budget constraint per time step, (18), and/or for the maximum built capacity constraint per time step, (19).
The main challenge in resolving the optimization is related to nonlinearity of the Power Flow relations (10), (11) and also to nonlinearity of the line thermal limits (17). Available non-linear solvers, such as IPOPT, are not effective in resolving the nonlinearities for large systems efficiently. This has motivated us to develop an heuristic algorithm consisting in sequential linearization of the nonlinear constrains discussed in the following Section (and, specifically, in Subsection III-A.
Iii The Algorithm
This Section describes the algorithm which allows us to resolve efficiently (and in spite of its complexity) the optimization problem just stated. Our algorithm consists of the following steps:
Scenarios are generated for each time frame according to the methodology suggested and described in details in 
. Briefly, one picks the base case, re-scale it for different time frame (taking into account the economic growth), and then introduce fluctuations around the re-scaled solutions to represent the forecasted load uncertainty. The fluctuations are chosen to be Gaussian with the standard deviation proportional to the mean.
Generation is initialized (for each load scenario) according to scheme explained in Appendix C.
If some of the constraints (7)-(17) are violated the initial state of the system is outside of the feasible domain defined by them. The non-linear constraints (10), (11) and (17) are linearized around the current state. This allows to construct current linearized version of the non-linear optimization problem (II)-(17).
The resulting linearized problem is solved by QP (or LP, depends on generation cost functions) using one of the available algorithms of the CPLEX solver .
AC power flow (AC-PF) is solved to update the state obtained at the previous step. This step is needed to prepare a feasible solution for the next iteration.
Steps 2-5 are repeated till either no constraints remain violated or the target precision is reached or the maximum allowed number of iterations is reached.
It is important to emphasize that, by construction, the algorithm maintains a feasible physical states at each iteration of its main loop including linearization, solution of the current QP optimization and back projection to the non-linear PF equations (achieved through the AC PF step).
Below we present details of the main steps of the algorithm.
The Power Flow (PF) equality constraints (10), (11) and the thermal limit constraints (17) are nonlinear. We choose to add to the optimization formulation auxiliary variables – active and reactive power flows expressed via voltages, phases and system parameters (see Appendix A for modeling details). This allows to localize non-linearities in local relations between the line power flows and respective voltages and phases at the two ends of the lines. Details of this technical trick, aimed at improving performance of the CPLEX solver, are as follows.
Introduce the following operational state.
and substitute the left hand side of Eq. (17) by its Taylor expansion around the current state, ,
thus arriving at the two equations representing a line
where the newly introduced auxiliary variables should also be substituted by the respective linearized expressions
Three comments/clarifications are in order. First, notice that the operational variables are adjusted independently for each scenario, thus enabling devices’ efficient utilization. Second, to manage possible degeneracy of the resulting system of linear constraints and limit the change of reactive flows at every we add to Eqs. (II)-(19) the following soft constraints for reactive power dispatch at each QP/LP step of the procedure
(In the case of LP, when the degeneracy is stronger, we also add similar constraints imposed on the active power.) Finally, third, to speed up the QP/LP computations one uses a cutting plane (constraint management) procedure. We split the whole set of constraints (17) into “active” and “inactive” sets including the constraints which were overloaded and, respectively, not overloaded, at the current state (of the previous iteration) or at any of the preceding steps. Only active constraints are explicitly accounted for in the optimization, while the validity of the inactive constraints is verified post-factum and the active/passive split is updated at every LP/QP step.
Iii-B QP/LP implementation
Standard CPLEX solver is called at each QP/LP step which outputs operational variables for each scenario along with investment variables for each time frame, and .
Iii-C AC-PF feasibility
The QP/LP step is followed by the AC-PF step, which is needed to maintain the AC PF feasibility destroyed by the linearization. Overall, combination of the QP/LP and AC-PF steps allow to maintain solution and resolve contingencies of the system simultaneously and gracefully.
Iv Case studies
The AC PF and optimization algorithms are implemented in Julia/JuMP. (See  and references therein.) QP/LP optimizations (called at internal steps of our algorithm) are resolved by CPLEX . When possible we utilize IPOPT , called from JuMP, to solve the optimization problem (in its original, nonlinear formulation). The (brute-force) IPOPT solution is computationally expansive and it is used as a ground truth (to validate our heuristic algorithm). Computational performance of the algorithms is analyzed on a Macbook Pro laptop (Core i7 3.3 GHz (2 Cores), 16 Gb of RAM).
Iv-a Algorithm validation
The algorithm is validated by comparison with the IPOPT solution. Since IPOPT is not able to resolve the 2376 bus-large Polish model (even with a single scenario) we perform the initial validation study on the 30 bus IEEE model. (Both the Polish model and the IEEE 30 bus model are available within the MathPower package .) Table I presents results of the (IPOPT vs our heuristics) comparison, where the planning horizon is taken to be year, , and for each time interval. Budget constraint per time interval is set to .
IPOPT objective is . Main algorithm objective is . Our main algorithm gives upper bounds solutions of the objective function with the gap less than 0.1.
Iv-B Scalability analysis
Our next step (after completion of the aforementioned validation study) was to perform a comparative analysis of the algorithms’ (computational) scalability.
First, we fix the number of scenarios ( per time interval) and study dependence on the number of time intervals. The results are shown in Fig. 2. Then, we consider one time interval and study dependence on the number of scenarios. These results are shown in Fig. 3. Both tests are done still on the 30 bus model with the quadratic generation cost (de-fault in the Mathpower package). In both cases we compare performance of the brute-force IPOPT solver with performance of our main algorithm and of the main algorithm reinforced by the cutting plane.
Comparing performance of the IPOPT in the two setting one observes a drastic difference. In the case of Fig. 2 IPOPT shows a surprisingly good performance, outperforming in speed both of our algorithms. (We relate this good performance of the IPOPT to using the primal simplex option within the IPOPT.) However the situation is reversed in the case of Fig. 3 where, moreover, performance of the IPOPT degrades exponentially with the number of scenarios. Our main algorithm and the reinforced (by cutting plane) algorithm show similar scaling performance in both cases (with the reinforced algorithm performing slightly better). Notice that juxtaposing here our algorithms against each other has a sense because of the comparable number of constraints contributing the two optimization settings.
Moving to the scaling analysis of the Polish model, one first of all note that in this case the IPOPT fails to converge. To illustrate performance of our reinforced algorithm we focus on analyzing dependence on the number of time intervals. The results of our scaling experiments with the Polish model are shown in Fig. (4), where the dashed line show a (rather satisfactory) quadratic match (for dependence of the overall computational time on the number of intervals.) In this case we use the Primal Simplex CPLEX solver at each LP step of our algorithm. (This is LP and not QP, as in the 30 bus model, because the de-fault generation cost is linear in the Polish model of Mathpower.)
We conclude this Section with a number of preliminary, and not yet fully conclusive but calling for further investigation, remarks. First of all, we have observed that developing efficient computational strategy for our linearization algorithm/heuristics became the task which is rather sensitive to the functional form of the generation cost and the choice of variables. If the cost is linear (in generated power) introducing auxiliary (line flow) variables and using the Primal Simplex solver is the winning strategy. However, the same approach in the case of quadratic cost (QP step replacing LP step) leads to slower convergence for the Primal Simplex algorithm while the barrier algorithm fails to converge at all.
Iv-C Gradual investments to resolve congestion
We apply our newly developed algorithm to study effect of the gradual investment, available only within the multi-time period framework, on the overall cost. We study the Polish model in the case of a single scenario with the optimization horizon of one year broken in periods. The (single) loading scenario, chosen to be stressed but still feasible (it is only away from the boundary of the AC OPF infeasibility - see Appendix B for details), stays the same over time. The congestion cost of the initial loading scenario (yet no investments in FACTS) is 17000 . The investment budget is limited to per (one month) time interval. The results of optimal investment generated by our reinforced algorithm are illustrated in Figs. 5,6,7. We observe that only SC devices were installed at lines, of which only two would be overloaded (if the line limits are, first, ignored while solving AC OPF and then checked for the overload). (Lines which are both overloaded, and thus contained in the active set of our cutting plance algorithm, and which are also selected for optimal SC installation are shown blue in Fig. 5. Lines which are shown green were not overloaded but chosen for SC installation. Lines which are shown red were overloaded but were not chosen for SC installation.)
We observe that the optimal installation is gradual. Moreover, all (constrained) available money are spent at each (time interval) decision. Fig. 6 shows how the congestion reduces with time, thus leading to reduction of the operational cost (blue bars) as time progresses. Orange line marks result of the AC OPF before investments start. Red line marks result for the (initial, i.e. before investments) AC OPF with the thermal limits ignored.
Fig. 7 shows that distribution of investments over lines and time period is nontrivial, therefore utilizing the newly available SC-capacities with other operational degrees of freedom.
We also show in Fig. 6 and Fig. 7 the result of optimal investment when the entire year (time horizon) is considered as one time interval. We observe that in the latter case the entire installation budget is used immediately such that the final solutions in the case one and intervals are the same and equal to the available budget. Notice, however, that making one investment upfront for the entire year, as opposed to breaking it into , periods is preferable because the overall (integrated over the year) cost of generation is reduced.
In this manuscript, we propose a new optimization framework for algorithmic resolution of placement and sizing of FACTS devices in a large transmission grids. Our algorithmic solution of the problem takes into account non-linear power flow equations. The most important features of the newly developed algorithms include, scalability, allowing to resolve congestion over practical (thousands of buses) size transmission systems, and also the ability to resolve multiple scenarios and over multiple time intervals simultaneously. The optimization can be considered as generalizing standard AC-OPF over multiple scenarios and multiple time intervals with an added cost of installation. This optimization accounts for both installation and operations, thus allowing the installed devices to adjust operationally to a particular scenario within the bounds set by the installed capacity.
Many interesting cases were analyzed experimentally. In all the cases considered the output (optimal placement) is spatially sparse also showing strong non-locality (in the sense that placement of a device may resolve congestion in a distant region of the grid). It is also observed that under highly loaded conditions FACTS devices are beneficial in reducing the total cost of generation. Optimal installation of the devices helps to resolve infeasibilities that are projected to become even more severe in the future.
Main technical achievement reported in this paper is the development of an efficient heuristics for solving the non-linear, non-convex and multi-time-interval optimization. The developed algorithm builds a convergent sequence of convex optimizations with linear constraints. Each constraint is represented explicitly through exact analytical linearization of the original nonlinear constraints (e.g. representing power flows and apparent power line limits) over all the degrees of freedom (including FACTS corrections) around the current operational point for particular loading scenario over particular time interval. In order to represent uncertainty in the projected growth of the system (loads) a custom scenario sampling split over multiple time intervals is introduced. Practicality of our approach for resolving the problem of investment (new installation) planning is illustrated on the IEEE 30-bus model and 2736-bus Polish model. It is evident from the experimental results that the approach is capable of both improving the system’s economy (reduce congestion price and generation cost) and also of resolving feasibility issues by introducing additional degrees of freedom (associated with the newly installed FACTS devices).
The work at LANL was supported by funding from the U.S. Department of Energy’s Office of Electricity as part of the DOE Grid Modernization Initiative.
Appendix A Transmission system modeling
Figure 8 shows a pictorial representation of model of a transmission line.
The model is parametrized by series impedance , total charging susceptance , transformation ratio and shift angle . A transmission line has two ends: sending (often called “from” end) and receiving (often called “to” end).
Explicit expressions for the apparent powers injected at “from” and “to” ends of a line in terms of voltages and phases are:
where and . Note that not only , but also expressions for and are non-symmetric with respect to the change of and indices. The symmetry is broken as the transformer is positioned at “from” side of the line.
Generators and loads:
Generators and loads are modeled as apparent power injections or consumptions:
Generation cost is a quadratic function of active power generation.
FACTS devices are described by two variables - first is capacity, second is setting.
Series Compensation (SC) device capacity and setting for the line :
Static Var Compensation (SVC) device capacity and setting for the load node :
Appendix B Generation of scenarios
The method of scenario generation/sampling is used to include the uncertainty related to system load for the planning period. Standard power system load growth over the time horizon is modeled via the LD curve. The base LD curve is illustrated in Figure 9. In our simulations we usually have just a single base case loading profile for a power case. In order to define where it is situated on the corresponding LD curve we define loading level to represent condition which is from ACOPF infeasibility by uniform load rescaling (we call it top point conditions).
We use the base LD curve, first, to generate LD curvess for consecutive years, re-scaling the base LD curve by the load growth factor of a year. Second, each early LD curve is split into piece-wise-constant parts. (We choose in our experiments.) Finally, each piece of an LD curve is used to generate a scenarios according to a random (thus called sampling) procedure described below. This scheme of scenario generation/sampling models variations in the distribution of loads thus simulating power system behavior during an extended period of time in the future.
We assume (and this assumption is maintained in all of our experimental tests) that each of the generated (sampled) load scenario is ACOPF feasible when the line constraints are ignored. (In other words, we consider the setting when there is enough of generation capacity even for the stressed cases.) Depending on the sampled scenario, 3 situations may arise.
ACOPF is feasible and congestion price is zero (low loading level).
ACOPF is feasible and congestion price is positive (higher loading level representing peak conditions).
ACOPF is infeasible due to either congestion of lines and/or voltage constraints but the system has enough generation capacity. ACOPF without apparent power limits on lines (and without voltage constraints if infeasible) is feasible (overloaded conditions which are possible in the future).
The aim of planning installation of FACTS devices at the right locations with their corresponding capacities is to reduce generation cost for point 2, and to improve or extend feasibility domain of the system for the point 3. Extra years of service can hence be added to the existing grid by making it more flexible, thereby allowing to delay investments into new lines and generators.
B-a Scenarios sampling for each segment
The loading level for a segment is represented by:
Future loading configurations are obtained from the base case by re-scaling all active and reactive loads by uniformly. The resulting vector of loads for a segment is thus given by:
Loading configurations are generated, for each segment and each , through modification of initial . It is done by adding Gaussian correction to each load with zero expected value and a respective standard deviation:
where, is given by:
The choice of parameters used in our experimental test to sample the scenarios is described in Table II.
Congestion Analysis Correction
If we study a case where for given load configuration standard AC-OPF outputs solution which is not congested, i.e. solution for which each constraint (on line flows or voltages) is satisfied with a margin, then this scenario does not require any FACTS installations. If the whole segment (from the procedure described in the preceding Subsection) is of this “zero-congestion” type, then obviously do not need to generate many samples representing the segment. Instead, we pick one re-scaled base scenario to represent the whole segment.
Appendix C Initializing the generation profile
Generation capacity is assumed to be large enough, i.e. not limiting, for the loading level considered. The initial profile of the generation for each load scenario have to be determined to run the algorithm. The initial generation profile is derived following the following two steps. 1) Re-solving ACOPF with the thermal limits ignored. 2) Setting up proportional response of the generators. The second step includes a) search for the smallest load re-scaling factor lowering the load and thus making the resulting case feasible; b) resolving ACOPF with this new re-scaled loading; c) increasing the generation and load proportionally to match the initial system loading (voltages are kept equal to these provided by the ACOPF solution); finally, d) solving ACPF to obtain generation matching the initial load.
-  N. Alguacil, A. L. Motto, and A. J. Conejo, “Transmission expansion planning: a mixed-integer lp approach,” IEEE Transactions on Power Systems, vol. 18, no. 3, pp. 1070–1077, Aug 2003.
-  R. Bent, G. Toole, and A. Berscheid, “Transmission network expansion planning with complex power flow models,” IEEE Transactions on Power Systems, vol. 27, no. 2, pp. 904–912, 2012.
-  A. H. Fuldner, “Upgrading transmission capacity for wholesale electric power trade,” :http://www.eia.doe.gov/cneaf/pubs_html/feat_trans_capacity/w_sale.html, accessed: 2016-11-16.
-  H. Nagarajan, R. Bent, P. Van Hentenryck, S. Backhaus, and E. Yamangil, “Resilient Transmission Grid Design: AC Relaxation vs. DC approximation,” arxiv:1703.05893, 2017.
-  K. Aoki, M. Fan, and A. Nishikori, “Optimal VAR planning by approximation method for recursive mixed-integer linear programming,” IEEE Trans. Power Syst., vol. 3, no. 4, pp. 1741–1747, Nov. 1988.
-  Y. T. Hsiao, C. C. Liu, H. D. Chiang, and Y. L. Chen, “A new approach for optimal VAR sources planning in large scale electric power systems,” IEEE Trans. Power Syst., vol. 8, no. 3, pp. 988–996, Aug. 1993.
-  D. Chattopadhyay, K. Bhattacharya, and J. Parikh, “Optimal reactive power planning and its spot-pricing: an integrated approach,” IEEE Trans. Power Syst., vol. 10, no. 4, pp. 2014–2020, Nov. 1995.
-  F. Dong, B. H. Chowdhury, M. L. Crow, and L. Acar, “Improving voltage stability by reactive power reserve management,” IEEE Trans. Power Syst., vol. 20, no. 1, pp. 338–345, Feb. 2005.
-  G. Y. Yang, G. Hovland, R. Majumder, and Z. Y. Dong, “TCSC allocation based on line flow based equations via mixed-integer programming,” IEEE Trans. Power Syst., vol. 22, no. 4, pp. 2262–2269, Nov. 2007.
-  A. Lashkar Ara, A. Kazemi, and S. A. Nabavi Niaki, “Multiobjective optimal location of FACTS shunt-series controllers for power system operation planning,” IEEE Trans. Power Del., vol. 27, no. 2, pp. 481–490, Dec. 2011.
-  S. J. Galloway, I. M. Elders, G. M. Burt, and B. Sookananta, “Optimal flexible alternative current transmission system device allocation under system fluctuations due to demand and renewable generation,” IET Gen. Trans. Distrib., vol. 4, no. 6, pp. 725–735, June 2010.
-  R. Minguez, F. Milano, R. Zarate-Miano, and A. J. Conejo, “Optimal network placement of SVC devices,” IEEE Trans. Power Syst., vol. 22, no. 4, pp. 1851–1860, Nov. 2007.
-  E. Ghahremani and I. Kamwa, “Optimal placement of multiple-type FACTS devices to maximize power system loadability using a generic graphical user interface,” IEEE Trans. Power Syst., vol. 28, no. 2, pp. 764–778, Aug. 2013.
-  R. S. Wibowo, N. Yorino, M. Eghbal, Y. Zoka, and Y. Sasaki, “FACTS devices allocation with control coordination considering congestion relief and voltage stability,” IEEE Trans. Power Syst., vol. 26, no. 4, pp. 2302–2310, Nov 2011.
-  J. G. Singh, P. Guha Thakurta, and L. Soder, “Load curtailment minimization by optimal placements of SVC/STATCOM,” Int. Trans. Elect. Energy Syst., vol. 25, no. 11, pp. 2769–2780, Sept. 2014.
E. A. Leonidaki, D. P. Georgiadis, and N. D. Hatziargyriou, “Decision trees for determination of optimal location and rate of series compensation to increase power system loading margin,”IEEE Trans. Power Syst., vol. 21, no. 3, pp. 1303–1310, Jul. 2006.
-  F. G. M. Lima, F. D. Galiana, I. Kockar, and J. Munoz, “Phase shifter placement in large-scale systems via mixed integer linear programming,” IEEE Trans. Power Syst., vol. 18, no. 3, pp. 1029–1034, Aug. 2003.
-  S. Rahimzadeh, M. Tavakoli, and B. A. H. Viki, “Simultaneous application of multi-type FACTS devices to the restructured environment: achieving both optimal number and location,” IET Gen. Trans. Distrib., vol. 4, no. 3, pp. 349–362, Mar. 2010.
-  S. Gerbex, R. Cherkaoui, and A. J. Germond, “Optimal location of multi-type FACTS devices in a power system by means of genetic algorithms,” IEEE Trans. Power Syst., vol. 16, no. 3, pp. 537–544, Aug. 2001.
-  P. Paterni, S. Vitet, M. Bena, and A. Yokoyama, “Optimal location of phase shifters in the French network by genetic algorithm,” IEEE Trans. Power Syst., vol. 14, no. 1, pp. 37–42, Feb. 1999.
-  L. Ippolito and P. Siano, “Selection of optimal number and location of thyristor-controlled phase shifters using genetic based algorithms,” IEE Proc. Gen. Trans. Distrib., vol. 151, no. 5, pp. 630–637, Sept. 2004.
-  C. Duan, W. Fang, L. Jiang, and S. Niu, “FACTS devices allocation via sparse optimization,” IEEE Trans. Power Syst., vol. 31, no. 2, pp. 1308–1319, Mar. 2016.
-  V. Frolov, S. Backhaus, and M. Chertkov, “Efficient algorithm for locating and sizing series compensation devices in large power transmission grids: I. model implementation,” New Journal of Physics, vol. 16, no. 10, 2014.
-  ——, “Efficient algorithm for locating and sizing series compensation devices in large power transmission grids: II. solutions and applications,” New Journal of Physics, vol. 16, no. 10, 2014.
-  V. Frolov, P. Guha Thakurta, S. Backhaus, J. Bialek, and M. Chertkov, “Optimal Placement and Sizing of FACTS Devices to Delay Transmission Expansion,” arxiv:, 2016.
-  “ILOG CPLEX optimization studio, http://www-03.ibm.com/software/products/ru/ibmilogcpleoptistud.”
-  “[online] julia/jump package for steady-state power network optimization,” https://github.com/lanl-ansi/PowerModels.jl, accessed: 2017-06-17.
-  “[online] interior point optimizer,” https://projects.coin-or.org/Ipopt, accessed: 2017-06-18.
-  R. D. Zimmerman, C. E. Murillo-Sanchez, and R. J. Thomas, “MATPOWER: Steady-State Operations, Planning and Analysis Tools for Power Systems Research and Education,” IEEE Trans. Power Syst., vol. 26, no. 1, pp. 12–19, Feb 2011.