Updated results for ORGDT
%\footnotesize %\small % ``The gang" Harsha Nagarajan$^{\dag
As illustrated in recent years (Superstorm Sandy, the Northeast Ice Storm of 1998, etc.), extreme weather events pose an enormous threat to the electric power transmission systems and the associated socio-economic systems that depend on reliable delivery of electric power. Besides inevitable malfunction of power grid components, deliberate malicious attacks can cause high risks to the service. These threats motivate the need for approaches and methods that improve the resilience of power systems. In this paper, we develop a model and tractable methods for optimizing the upgrade of transmission systems through a combination of hardening existing components, adding redundant lines, switches, generators, and FACTS and phase-shifting devices. While many of these controllable components are included in traditional design (expansion planning) problems, we uniquely assess their benefits from a resiliency point of view. More importantly, perhaps, we evaluate the suitability of using state-of-the-art AC power flow relaxations versus the common DC approximation in resilience improvement studies. The resiliency model and algorithms are tested on a modified version of the RTS-96 (single area) system.
Algo.tex
%!TEX root = ../pscc2016_orgdt.tex

%+++++++++++++++++++++++++++;
%  Section: Algorithms      ;
%+++++++++++++++++++++++++++;
\section{Algorithms}
\label{Sec:algo}
In this section we discuss algorithms we use to solve the ORGDT. The ORGDT is a two-stage Mixed-Integer Quadratically Constrained Program (MIQCP) with a block diagonal structure. 
%that includes coupling variables between the blocks.
In order to exploit this structure, we generalize the  scenario-based decomposition (SBD) techniques of \cite{yamangil2015resilient} to solve the ORGDT. In the remainder of this paper, let $\ORGDT(\Disasters^\prime)$ denote the ORGDT with the scenario set $\Disasters^\prime \subseteq \Disasters$ and $\sigma$ denote the vector of construction variables $\lineVariable_{ij}$, $\switchVariable_{ij}$, $\hardenVariable_{ij}$,$\gamma_{ij}$, $\delta_{ij}$ for all $ij \in \Edges$ and $\facilityVariable_{i}, zp_i$ for all $i \in \Nodes$. 
%\vspace{-.5cm}
%\subsection{Exact Algorithm (Scenario-Based Decomposition (SBD))}
%\label{subsec:SBD}
%Decomposition is often used for solving two-stage stochastic MIPs \cite{VW:10}, and it 
SBD can be applied to ORGDT after the following key observation:

%\vspace{-.2cm}
\begin{obs}
The second stage variables do not appear in the objective function. Therefore any optimal first stage solution based on a subset of the second stage subproblems that is feasible for the remaining scenarios, is an optimal solution for the original problem.
\end{obs}


%Based on this observation, we can apply SBD to solve ORGDT as outlined in 

The SBD approach is outlined in Algorithm \ref{algo:sbd}. At a high level, Algorithm \ref{algo:sbd} solves problems with iteratively larger sets of scenarios until a solution is obtained that is feasible for all scenarios. The algorithm takes as input the set of scenarios and an initial scenario to consider, $\Disasters^\prime$.  Line 2 solves the ORGDT on $\Disasters^\prime$,
where $\ORGDT(\Disasters^\prime)$ and $\sigma^*$ are used to denote the problem and solution respectively.  Line  3 then evaluates   $\sigma^*$ on the remaining scenarios in $\Disasters \setminus \Disasters^\prime$.
The function $l:P^\prime(s,\sigma^*) \rightarrow \mathbb{R}_+$, is an infeasibility measure that is 0 if the problem is feasible for scenario $s$, positive otherwise. This function is implemented by solving the feasibility subproblem reliability constraints, i.e. total and critical demands are satisfied.
We use a version of (\ref{eqn:primal_r}) where $M=1$ and the non convex constraints are replaced with their convex counterparts.
%It measures the gap between the delivered and the required demand (right hand side of constraints \ref{eqn:p_cr}-\ref{eqn:q_ncr}). 
This function prices the current solution over $s\in \Disasters\setminus \Disasters^\prime$. If all prices are 0, then the algorithm terminates with solution $\sigma^*$ (lines 4-5).  Otherwise, the algorithm adds the scenario with the worst infeasibility measure to $\Disasters^\prime$ (line 7).

\begin{algorithm}[!ht]
%\SetLine
\small
\SetKwInOut{Input}{input}
\Input{A set of disasters $\Disasters$ and let $\Disasters^\prime = \Disasters_0$\;}
\While{$\Disasters\setminus \Disasters^\prime \neq \emptyset$}{
	$\sigma^* \leftarrow$ Solve $\ORGDT(\Disasters^\prime)$\;
%	Let $\bar{x}$ denote the optimum first stage solution\;
%	Let $I\leftarrow \emptyset$\;
	$I \leftarrow \left<s_1, s_2 \ldots s_{|\Disasters\setminus \Disasters^\prime|} \right> s\in \Disasters\setminus \Disasters^\prime :$ $l(\ORGDT^\prime(s_i,\sigma^*)) \ge l(\ORGDT^\prime(s_{i+1},\sigma^*))$\;
	%\For{$s\in S\setminus S^\prime$}{
	%	$I\rightarrow push\_back (l(P^\prime(s,\bar{x})))$\;
%	}
	\eIf{$l(\ORGDT^\prime(I(0),\sigma^*)) \leq 0$}{
		\Return{$\sigma^*$}\;
	}{
		$\Disasters^\prime \leftarrow \Disasters^\prime \cup  I(0)$\;
	}
\Return{$\sigma^*$}
}
\caption{Scenario-Based Decomposition}
\label{algo:sbd}
\end{algorithm}

\begin{rmk}
We observed that the LP-relaxation for the ORGDT is very loose. To overcome this issue, we augmented every iteration of Algorithm \ref{algo:sbd} with the previous optimum objective value as a lower bound for the current iteration.
\end{rmk}

\noindent
{\bf Cutting-plane algorithm for quadratic constraints:}
%For the sake of completeness, we invoke the following quadratic and SOC constraints from Sec. \ref{Sec:formulation}. 
%{\fontsize{9}{8}\selectfont
%\begin{subequations}
%\label{eqn:SOCs}
%\begin{align}
%& {p^s}^{2}_{ij} + {q^s}^{2}_{ij} \leq \loss_{ij} {\hat{v}^s}_{i},  \label{eqn:soc1} \\
%&{p^s_{ij}}^2+{q^s_{ij}}^2 \leq \tilde{x}^s_{ij} T^2_{ij},  \ {p^s_{ji}}^2+{q^s_{ji}}^2 \leq \tilde{x}^s_{ij} T^2_{ij},  \label{eqn:soc3} \\
%& {\voltageVariable}_i^2 \leq \widehat{\voltageVariable}_{i}, \ \mathrm{Constraint \ \eqref{eqn:cos_qc}} \label{eqn:soc5}
%\end{align}
%\end{subequations}}
%\vspace{-0.3cm}
Even though optimization theory guarantees that the set of convex inequalities in ORGDT
%\eqref{eqn:SOCs} 
can be solved efficiently, several numerical experiments demonstrated that it was challenging to solve even moderately sized problems using
state-of-the-art MIQCP solvers (CPLEX/Gurobi). Either the  solver convergence was very slow or it terminated with ``numerical trouble". In order to circumvent this issue, we adopted a cutting-plane approach.
Let $\mathcal{P}^{miq}$ be a general MIQCP which contains the following quadratic constraint (rotated second-order cone):
%, representing the generalization of constraints \eqref{eqn:SOCs}: 
\begin{equation}
\label{eqn:SOC_gen}
\mathbf{x}^T\mathbf{x} \leq yz, \ y\geq 0,z\geq 0, \ (\mathbf{x},y,z) \in \mathbb{R}^{|\mathbf{x}|+2}
\end{equation}

Algorithm \ref{algo:cp} outlines the cutting plane procedure to solve $\mathcal{P}^{miq}$. The key idea behind this algorithm is to solve a finite sequence of MILPs to obtain an optimal solution for the original MIQCP. Let $\mathcal{P}^{mil}$, a MILP, represent the original problem without the quadratic constraints. The solution in step 3 of Algorithm \ref{algo:cp} is then guaranteed to be a lower bound to $\mathcal{P}^{miq}$. This lower bound is tightened further for every infeasible solution (step 4) by adding a valid cutting plane (step 5). This cutting plane is valid as it is an outer-approximation of the original feasible convex set. This procedure is repeated until the solution obtained is feasible (satisfies step 4), and hence optimal, for $\mathcal{P}^{miq}$.  

\begin{algorithm}[ht]
\small
\caption{Cutting plane algorithm for MIQCPs}
\label{algo:cp}
\textbf{Notation}: $f(\mathbf{x},y) = \frac{\mathbf{x}^T\mathbf{x}}{y}$ \\
\textbf{Input}: {$\mathcal{P}^{mil}$, $tol>0$, $\mathcal{L}=\varnothing$} \\
$\hat{\mathbf{x}},\hat{y},\hat{z}$ $\leftarrow$ Solve $\mathcal{P}^{mil}(\mathbf{x},y,z)$ \\
\While{$\hat{\mathbf{x}}^T\hat{\mathbf{x}} > \hat{y}\hat{z} + tol$}{
	Augment $\mathcal{L}$ with the following cut: 
	$f(\hat{\mathbf{x}},\hat{y}) +\sum_{i=1}^{|\mathbf{x}|} \frac{\partial f(\mathbf{x},y)}{\partial x_i}(x_i-\hat{x}_i) + \frac{\partial f(\mathbf{x},y)}{\partial y}(y-\hat{y}) \leq z$ \\
% \Return{$\sigma^*$}
$\hat{\mathbf{x}},\hat{y},\hat{z}$ $\leftarrow$ Solve $\left(\mathcal{P}^{mil}(\mathbf{x},y,z) \cup \mathcal{L}\right)$
}
\end{algorithm}

%Numerical experiments suggest that this approach is more stable and faster (mostly) than the conventional SOC solvers. 

% The basic idea of the cutting-plane approach is to successively solve tighter polyhedral outer approximations of the original set until a solution is found that is feasible to the original set. 


%\noindent
%{\bf Improving SBD with a heuristic solution} SBD has one weakness, that is for each iteration, the optimum solution found in the previous iteration is infeasible. Therefore, we do not have an initial feasible solution to the master problem at each iteration. To  this issue, we can heuristically solve the master problem to get an initial feasible solution and warm start the master problem from this solution. The algorithm proceeds exactly the same as Algorithm 1, except that $Solve(\ORGDT(S^\prime))$ is replaced by $Solve(\ORGDT(\Disasters^\prime)) \textit{ starting from } Heuristic(\ORGDT(\Disasters^\prime))$ in line 2, where $Heuristic$ defines an arbitrary heuristic. 

%\subsection{Heuristic Algorithms}


%\vspace{-.9cm}

%%% RBENT, since we are not currently using this anywhere, I commented it out... we'll want it in later papers

%\noindent
%{\bf Construction Heuristic (Greedy)} A computationally efficient way of generating feasible solutions to the ORGDT relaxes the coupling first %stage variables and solves each scenario $s \in \Disasters$ individually.
%The solutions are combined by taking the maximum of each construction variable 
% (${\cal{X}} = \lineVariable \cup \switchVariable \cup \hardenVariable \cup \microgridVariable \cup \facilityVariable$)
%over all scenarios (Algorithm 2).
%The switch construction cost is determined by switches that are needed to reduce the network into a tree for every scenario (line 4). Although %the Greedy Algorithm is simple and fast, it rarely results in an optimal investment decision. However, it is representative of the types of %heuristics used by the industry: see Reference \cite{Munoz2014} for a survey.
%\vspace{-.4cm}

%\begin{algorithm}[!ht]
%\SetLine
%\small
%\SetKwInOut{Input}{input}
%\Input{A set of disasters $\Disasters$\;}
%\For{$s\in \Disasters$}{
%	$\sigma^s \leftarrow Solve(P^\prime(s))$\;
%}
%$\sigma^*(x) = \max\{\sigma^s(x) | \forall s\in \Disasters\},\; \forall x \in {\cal{X}}$\;
%Update $\sigma^*(x_i)$ with switches to preserve feasibility\;
%\Return{$\sigma^*$}
%\caption{Greedy}
%\label{greedy}
%\end{algorithm}

%\vspace{-.9cm}
% \noindent
%{\bf Improvement Heuristic (Variable Neighborhood Search (VNS))} To overcome the computational limitations of exact methods,
%we use a VNS heuristic.
%We can describe the VNS in the following high level form: Given an incumbent to the original problem, the algorithm fixes a subset of first stage variables to their incumbent value and searches the remaining variables for a better solution. If we can find a better solution, the incumbent is reset by the better solution and we restart the search from the new incumbent, else, we increase the neighborhood size (number of unfixed variables). If we reach a maximum number of iterations with no improvement, the algorithm restarts the search with a random subset of variables fixed to their incumbent values. The details of the algorithm can be found at \cite{yamangil2015resilient} (this is a modification of the general VNS for an arbitrary MIP proposed by \cite{LHMU:10}).

%It is important to note that, if all the first stage variables are fixed, the problem decomposes into $|S|$ separate problems that are easily solved and provide heuristic justification for searching the neighborhood of only the first stage variables.

%Algorithm 3 describes the VNS procedure. Line  1 computes the solution to the LP relaxation of the ORGDT, $(\sigma^{LP})$. Line 4 counts the number of variable assignments that are different between the solution to LP relaxation $(\sigma^{LP})$ and the best known solution $\sigma^*$ ($\sigma(x)$ denotes the variable assignment of $x$ in solution $\sigma$). Line 5 orders the variables of $\cal{X}$ by the difference between their assignments in $\sigma^*$ and $\sigma^{LP}$.  Heuristically, those variables whose assignments are furthest from their LP assignment represent good opportunities to improve $\sigma^*$. The algorithm updates the rate at which the neighborhood size is increased ($step$) based on whether or not the algorithm is in a restart situation (lines 8 and 11).  If the algorithm is in a restart, the ordering of the variables is also randomized (line 9).  Line 13 computes the best solution in the neighborhood of $\sigma$ where the first $k$ elements of $J$ are fixed.  If the resulting solution is better, then the algorithm proceeds with a new $\sigma^*$ (lines 15-18)--$f$ is used as shorthand for Eq.~\ref{eqn:obj}. Otherwise, the size of the neighborhood is increased (lines 20-23). The iterations terminate when the maximum number of restarts is reached (line 2), the maximum number of neighborhood re-sizings is reached (line 12), or a time limit is reached. In this paper, $\textsc{maxRestarts} = 10$, $\textsc{maxIterations} = 4$, $\textsc{maxTime} = 48$ CPU hours, and $d = 2$.
%\vspace{-.4cm}

%\vspace{-.9cm}
%\noindent
%{\bf SBD-VNS Heuristic for ORGDT} In this algorithm we combine the strengths of VNS with SBD. SBD provides an exact solution to ORGDT while VNS generates high quality solutions with limited CPU time. Understanding that the master problem takes a substantial amount of the computation time, it is natural to create a hybrid algorithm. The hybrid algorithm proceeds exactly the same as Algorithm 1, except that the exact solver for $Solve(\ORGDT(\Disasters^\prime))$ is replaced by $VNS(\ORGDT(\Disasters^\prime))$ in line 2.

%\begin{algorithm}[!ht]
%\SetLine
%\small
%\SetKwInOut{Input}{input}
%\While{$S\setminus S^\prime \neq \emptyset$}{
%	$\sigma^* \leftarrow$ Solve $P(S^\prime)$\ {\bf using VNS}\;
%	$I \leftarrow \left<s_1, s_2 \ldots s_{|S\setminus S^\prime|} \right> s\in S\setminus S^\prime :$ $l(P^\prime(s_i,\sigma^*)) \ge l(P^\prime(s_{i+1},\sigma^*))$\;
	%\For{$s\in S\setminus S^\prime$}{
	%	$I\rightarrow push\_back (l(P^\prime(s,\bar{x})))$\;
%	}
%	\eIf{$l(P^\prime(I(0),\sigma^*)) \leq 0$}{
%		\Return{$\sigma^*$}\;
%	}{
%		$S^\prime \leftarrow S^\prime \cup  I(0)$\;
%	}
%\Return{$\sigma^*$}
%}
%\caption{SBVNDS}
%\label{benders}
%\end{algorithm}

% TODO: DID WE TEST BENDERS?
%We also tested other decomposition strategies such as Benders and  Dantzig-Wolfe, however, their performance was tempered by the ORDGDP  structure. The ORDGDP has MIP formulations at both stages of the problem and does not contain optimality conditions in the second stage (only feasability conditions).  These approaches rarely out performed the commerical MIP solver.

%{\color{red}....I can't follow this paragraph....}  Therefore we do can apply the following simple algorithm to solve the Problem \ref{fig:optimization_model}, where $P(S^\prime)$ denote the problem induced by set $S^\prime \subseteq S$ and $P^\prime(s,\bar{x})$ denote the subproblem corresponding to the scenario $s\in S$ induced by a first stage solution $\bar{x}$. We also assume $l:P^\prime(s,\bar{x}) \rightarrow \mathbb{R}_+$, an infeasibility measure is available. The infeasibility measure is equivalent to 0 if the problem is feasible, positive otherwise. This can be easily implemented as maximizing the total or critical demand satisfied, or the gap between the delivered and the required demand, the right hand side of the constraints \ref{eqn:critical_constraint} and \ref{eqn:loadserve_constraint}, respectively. This is essentially pricing the current solution over the subproblems corresponding to every scenario $s\in S\setminus S^\prime$ that is currently out of the master problem.

%This scenario-based decomposition shares a number of key features with Benders decomposition \cite{VW:10}.  Both Benders and SBD solve the subproblems based on iterative solves to the first-stage problem.  Benders will typically add a single constraint in the form of cut that represents a facet of the subproblem, whereas in SBD we add the entire polyhedron associated with the sub-problem.
%Benders is a successful approach on similar expansion, commodity flow, and interdiction problems, but is unsuccessful here due to the non-convex (discrete) nature of the second stage. Hence our generalization to SBD.
%\vspace{-.4cm}

Conclusions.tex
%!TEX root = ../pscc2016_orgdt.tex

%+++++++++++++++++++++++++++++++;
%  Section: Conclusions         ;
%+++++++++++++++++++++++++++++++;

\vspace{-.2cm}

\section{Conclusions}
\label{Sec:conc}

In this paper, we formulated, modeled, and developed a model of ORGDT. 
Our contributions include a comparison of solutions to the ORGDT with the DC approximation and the QC relaxation that demonstrates the necessity of AC power flow modeling in resilience studies. We also developed an approach for recovering fully AC feasible solutions that indicate that solutions based on the QC relaxation tend to be tight in practice. Finally, we have also shown that a small set of strategically located PST and FACTS devices as design options can reduce the need for other design considerations like line hardening or distributed generation.

%exact and heuristic algorithms for solving the ORGDT that exploit the decomposable structure of the problem. We have also provided an analysis that demonstrates one of the key challenges in resilient design is handling constraints associated with enforcing stability-related phase angle differences on power lines. This constraint drives costs orders of magnitude higher, however, we have shown that these costs can be managed by introducing technologies, like FACTS, for resilience purposes.   

There remain a number of interesting future directions for work in this area. First, it will be important to scale the approaches of the paper to larger, more realistic transmission grids. An important idea in this area is to limit the upgrades to those that are deemed practical by subject matter experts. Second, it will be important to introduce tighter and improved models of FACTS and PST devices (such as continuous set points). Finally, it will be interesting to consider= approaches for using FACTS and PST devices to improve the restoration process of transmission grids, as seen in \cite{Mak2014}.

%Third, we are also working on improving our approach for generating primal feasible solutions from the relaxed solution. Note: I think being able to solve p_0 on the majority of the instances so far demonstrates this 

% \section*{Acknowledgements} This work was supported by the Defense Threat Reduction Agency  project \textit{Robust Network Interdiction Under Uncertainty Advancing Knowledge of Network Theory for Understanding Robustness} and the Center for Nonlinear Studies at LANL.



    % It can be observed that the savings range from 32.5 to 39.0\%, which is a substantial decrease in the upgrade cost (roughly \$60M). At lower values of $\theta^u$, since the networks are more congested, one must consider resiliency options like building new lines or hardening the congested portions of the network to increase throughput under extreme events. However, here we demonstrate that a small set of strategically placed, not-so-expensive FACTS devices can also accomplish the same objective (Fig. \ref{Fig:15deg}), which is very encouraging from an engineering perspective. 


Formulation.tex
%!TEX root = ../pscc2016_orgdt.tex
%+++++++++++++++++++++++++++++++++++++++++;
%  Section: Mathematical formulation      ;
%+++++++++++++++++++++++++++++++++++++++++;
\section{ORGDT Optimization Model}
\label{Sec:formulation}

% Insert Nomenclature
{\fontsize{8.3}{6}\selectfont
\input{Nomenclature}
\printnomenclature
}

We formulate the ORGDT as a two-stage mixed-integer nonlinear program, $\ORGDT$. The first-stage variables (without superscript $s$) specify the new infrastructure enhancements and the second-stage variables (with superscript $s$) describe the operation of invested infrastructure to serve loads for each damage scenario $s \in \Disasters$. 

{\fontsize{8.7}{8}\selectfont
\begin{subequations}
\begin{flalign}
\ORGDT := \min & \sum_{ij\in \Edges} \left( c^{x}_{ij} x_{ij} + c^{\tau}_{ij} \tau_{ij} + c^{t}_{ij} t_{ij} + c^{\delta}_{ij}\delta_{ij} + c^{\gamma}_{ij}\gamma_{ij}\right)  \nonumber\\
& + \sum_{i \in \Nodes}\left(c^{u}_{i} u_{i} + c^{zp}_{i} zp_{i} \right)  \label{eqn:obj}  \\
\text{s.t.}\quad & x^s_{ij} \leq x_{ij} \quad \quad \forall ij\in \Damages_s, s\in \Disasters  \label{eqn:assign:1a} \\
& x^s_{ij} = x_{ij} \quad \quad \forall ij\notin \Damages_s, s\in \Disasters  \label{eqn:assign:1b} \\
& \tau^s_{ij} \leq \tau_{ij}, \  t^s_{ij} \leq t_{ij}   \quad \quad \forall ij\in \Edges, s\in \Disasters  \label{eqn:assign:2} \\
& \delta^s_{ij} \leq \delta_{ij}, \  \gamma^s_{ij} \leq \gamma_{ij}   \quad \quad \forall ij\in \Edges, s\in \Disasters  \label{eqn:assign:3} \\
% & t^s_{ij} \leq t_{ij} \quad \quad \forall ij\in \Edges, s\in \Disasters  \label{eqn:assign:3} \\
& zp^s_{i} \leq zp_{i}, \ zq^s_{i} \leq zq_{i} \quad \forall i\in \Nodes, s\in \Disasters  \label{eqn:assign:4a} \\
& 0 \leq zp_{i} \leq zp^u_{i} u_{i}, \ |zq_{i}| \leq zq^u_{i} u_{i} \quad \forall i\in \Nodes   \label{eqn:assign:6a} \\
& zp_{i} \geq 2|zq_{i}|\quad \forall i\in \Nodes   \label{eqn:assign:6b} \\
& (\mathbf{x}^s,\boldsymbol{\tau}^s,\mathbf{t}^s,\boldsymbol{\delta}^s, \boldsymbol{\gamma}^s, \mathbf{zp}^s, \mathbf{zq}^s,\mathbf{u}) \in \SubProblem \quad \forall s\in \Disasters \label{eqn:feasible_network} \\
& x_{ij},\tau_{ij}, t_{ij},\delta_{ij},\gamma_{ij},u_{i} \in \{0,1\} \quad   \forall ij\in \Edges, i\in \Nodes \label{eqn:discrete_2}
% & u_{i} \in \{0,1\} \quad \quad  \forall i\in \Nodes  \label{eqn:discrete_3}
\end{flalign}
\label{eq:master_AC}
\end{subequations}}
\noindent
In $\ORGDT$, Eq. \eqref{eqn:obj} minimizes the total upgrade cost, which includes new lines, switches, hardening, FACTS, phase shifters, generators and real power capacity, respectively. Constraints \eqref{eqn:assign:1a}-\eqref{eqn:assign:4a} link the first-stage (construction) decisions with second-stage variables in $\SubProblem$.
Constraints \eqref{eqn:assign:6a} and \eqref{eqn:assign:6b} model generation 
capacity upgrades. Without loss of generality, reactive power capacity is limited to half of the real power capacity.
%facility location options such that the purchased real capacity is always greater than twice the reactive capacity. 
Constraint \eqref{eqn:assign:6b} ensures that the generators are not built purely for reactive power support. Constraint \eqref{eqn:feasible_network} states that the mixed-integer vector $(\mathbf{x}^s,\boldsymbol{\tau}^s,\mathbf{t}^s,\boldsymbol{\delta}^s, \boldsymbol{\gamma}^s, \mathbf{zp}^s, \mathbf{zq}^s,\mathbf{u}) \in \SubProblem$ is an AC feasible transmission network for scenario $s$ subject to the constraints of $\SubProblem$. 

The constraints of $\SubProblem$ describe the AC-power flow equations and budget constraints on resiliency options. 
%The two fundamental physics laws that govern the network flows are a) 
Ohm's law is modeled in constraints \eqref{eqn:acopf_real_ij}-\eqref{eqn:acopf_reac_ji}. For a given topology of the network. $\tilde{x}^s$, flow on line $ij$ is forced to zero when $\tilde{x}^s_{ij}=0$. %and b) 
Kirchhoff's current law (flow conservation) is given in constraints \eqref{eqn:balance1} and \eqref{eqn:balance2}.
%, capturing the asymmetry of AC power flow. 
Constraints \eqref{eqn:acopf_loss1}-\eqref{eqn:acopf_loss3} model the connection between current magnitude ($l^s_{ij}$) and the power loss equations as discussed in \cite{baran1989optimal}.
Operational constraints \eqref{eqn:thermalij}-\eqref{eqn:v_bnd} represent the thermal limits, phase angle difference limits %(worst-case bound is captured if $\tilde{x}^s_{ij}=0$) 
and voltage bounds, respectively. Constraint \eqref{eqn:damage} models the damaged lines of the scenario $s\in S$, i.e., a line is inoperable when damaged and unhardened. Constraint \eqref{eqn:lineminusswitch} defines the topology for scenario $s$. A switch is available for operation only if the line is active. Constraints \eqref{eqn:mg_a} bound the real and reactive power capacities that can be built (if $u_i=1$) at bus $i\in \Nodes$. Constraints in \eqref{eqn:gp_a} and \eqref{eqn:gq_a} represent total real and reactive power capacities including the existing generators of the network. Finally, constraint \eqref{eqn:shed} captures the continuous power shedding subject to the resiliency constraints in \eqref{eqn:p_cr}-\eqref{eqn:q_ncr}. 
User-defined resiliency parameters ($lp_{cr},lq_{cr},lp_{ncr},lq_{ncr}$) ensure that a minimum fraction of critical and non-critical loads is served during every damage scenario $s \in \Disasters$. 
\vspace{-.48cm}

{\fontsize{8.7}{8}\selectfont
\begin{flalign}
\SubProblem  &  = \{\mathbf{x}^s,\boldsymbol{\tau}^s,\mathbf{t}^s,\boldsymbol{\delta}^s, \boldsymbol{\gamma}^s, \mathbf{zp}^s, \mathbf{zq}^s,\mathbf{u} :  \nonumber \\
%----------------------------------------------------------------------
& \textbf{On/Off AC power flow equations:} \nonumber \\
& p^s_{ij} = \tilde{x}^s_{ij} (\tilde{G}_{ij} {v^s}^{2}_{i} - \tilde{G}_{ij}v^s_i v^s_j \cos(\theta^s_{ij}) \nonumber \\  
& -\tilde{B}_{ij}v^s_i v^s_j \sin(\theta^s_{ij}))  \quad  \forall ij\in \Edges,  \label{eqn:acopf_real_ij} \\
&q^s_{ij} = \tilde{x}^s_{ij} (-\tilde{B}_{ij} {v^s}^{2}_{i} + \tilde{B}_{ij}v^s_i v^s_j \cos(\theta^s_{ij})\nonumber \\ 
& - \tilde{G}_{ij}v^s_i v^s_j \sin(\theta^s_{ij}) )  \quad  \forall ij\in \Edges,  \label{eqn:acopf_reac_ij} \\
&p^s_{ji} = \tilde{x}^s_{ij} (\tilde{G}_{ij} {v^s}^{2}_{j} - \tilde{G}_{ij}v^s_i v^s_j \cos(\theta^s_{ij}) \nonumber \\  
& +\tilde{B}_{ij}v^s_i v^s_j \sin(\theta^s_{ij}) ) \quad  \forall ij\in \Edges,  \label{eqn:acopf_real_ji} \\
&q^s_{ji} = \tilde{x}^s_{ij} (-\tilde{B}_{ij} {v^s}^{2}_{j} + \tilde{B}_{ij}v^s_i v^s_j \cos(\theta^s_{ij})\nonumber \\ 
& + \tilde{G}_{ij}v^s_i v^s_j \sin(\theta^s_{ij}) ) \quad  \forall ij\in \Edges,  \label{eqn:acopf_reac_ji} \\
& \theta^s_{ij} =  \theta^s_{i} -  \theta^s_{j} + \phi^s_{ij} \quad  \forall ij\in \Edges, \label{eqn:ph_ang}\\
& \tilde{G}_{ij} = (\bar{G}_{ij} - G_{ij})\delta^s_{ij} + G_{ij} \quad \quad \forall ij \in \Edges \label{eqn:conduct_disj}, \\
& \tilde{B}_{ij} = (\bar{B}_{ij} - B_{ij})\delta^s_{ij} + B_{ij} \quad \quad \forall ij \in \Edges \label{eqn:suscept_disj}, \\
& p^s_{ij} + p^s_{ji} = \resist_{ij} \loss_{ij} \quad  \forall ij\in \Edges,  \label{eqn:acopf_loss1} \\
& q^s_{ij} + q^s_{ji} = \react_{ij} \loss_{ij} \quad  \forall ij\in \Edges,  \label{eqn:acopf_loss2} \\
& {p^s}^{2}_{ij} + {q^s}^{2}_{ij} = \loss_{ij} {v^s}^{2}_{i} \quad  \forall ij\in \Edges,  \label{eqn:acopf_loss3} \\
%----------------------------------------------------------------------
% & \textbf{Flow Balance:} \nonumber \\
& gp^s_{i} - lp^s_{i} = \sum_{ij\in \Edges} p^s_{ij} + \sum_{ji\in \Edges} p^s_{ij} \quad \quad \forall i\in \Nodes,  \label{eqn:balance1}\\
& gq^s_{i} - lq^s_{i} = \sum_{ij\in \Edges} q^s_{ij} + \sum_{ji\in \Edges} q^s_{ij} \quad \quad \forall i\in \Nodes, \label{eqn:balance2}\\
%----------------------------------------------------------------------
& \textbf{Operational limits and topology constraints:} \nonumber \\
& {p^s_{ij}}^2+{q^s_{ij}}^2 \leq \tilde{x}^s_{ij} \capacity^2_{ij} \quad \quad \forall ij\in \Edges,  \label{eqn:thermalij} \\  
& {p^s_{ji}}^2+{q^s_{ji}}^2 \leq \tilde{x}^s_{ij} \capacity^2_{ij} \quad \quad \forall ij\in \Edges,  \label{eqn:thermalji} \\
& |\theta^s_{ij}| \leq \tilde{x}^s_{ij}\theta^u + (1-\tilde{x}^s_{ij})\theta^M \quad \quad  \forall ij \in \Edges, \label{eqn:theta_bnd}\\
& v^l_i \leq v^s_i \leq v^u \quad  \forall i\in \Nodes,  \label{eqn:v_bnd} \\
& -\phi^u \gamma^s_{ij} \leq \phi^s_{ij} \leq \phi^u \gamma^s_{ij} \quad \forall ij\in \Edges,  \label{eqn:phi_bnd} \\
%----------------------------------------------------------------------
% & \textbf{Topology Constraints:} \nonumber \\
& x^s_{ij} = t^s_{ij} \quad \quad \forall ij\in \Damages_s  \label{eqn:damage}\\
%& \tau^s_{ij} \leq x^s_{ij}  \quad \quad \forall ij\in \Edges  \label{eqn:switchable}\\
& \tilde{x}^s_{ij} = x^s_{ij} - \tau^s_{ij}\geq 0 \quad \quad \forall ij\in \Edges  \label{eqn:lineminusswitch}\\
& \delta^s_{ij} \leq \tilde{x}^s_{ij}, \ \gamma^s_{ij} \leq \tilde{x}^s_{ij} \quad \quad \forall ij \in \Edges, \\
%----------------------------------------------------------------------
& \textbf{Generation constraints}\quad  \forall i \in \Nodes: \nonumber \\
& 0 \leq zp^s_{i} \leq zp^u_{i}u_{i}, \ |zq^s_{i}| \leq zq^u_{i}u_{i}, \ zp^s_{i} \geq 2|zq^s_{i}| \label{eqn:mg_a}\\
% & -zq^u_{i}u_{i} \leq zq^s_{i} \leq zq^u_{i}u_{i}  \quad \quad \forall i\in \Nodes, \label{eqn:mg_b}\\
% & zp^s_{i} \geq 2|zq^s_{i}|, \label{eqn:zpzq_tie}\\
& 0 \leq gp^s_{i} \leq  gp^u_{i} + zp^s_{i},  \label{eqn:gp_a}\\
& (gq^l_{i} - |zq^s_{i}|) \leq gq^s_{i} \leq (gq^u_{i} + |zq^s_{i}|),  \label{eqn:gq_a}\\
%----------------------------------------------------------------------
& \textbf{Resilience constraints:} \nonumber \\
& lp^s_{i} =  yp^s_{i} dp_{i}, \ lq^s_{i} =  yq^s_{i} dq_{i} \quad \quad \forall i\in \Nodes,  \label{eqn:shed}\\
& \sum_{i \in \Load} lp^s_{i} \geq lp_{cr} \sum_{i \in \Load} dp_{i}, \label{eqn:p_cr}\\
& \sum_{i \in \Nodes \setminus \Load} lp^s_{i} \geq lp_{ncr} \sum_{i \in \Nodes \setminus \Load} dp_{i} \label{eqn:p_ncr}\\
& \sum_{i \in \Load} lq^s_{i} \geq lq_{cr} \sum_{i \in \Load} dq_{i} \label{eqn:q_cr}\\
& \sum_{i \in \Nodes \setminus \Load} lq^s_{i} \geq lq_{ncr} \sum_{i \in \Nodes \setminus \Load} dq_{i} \label{eqn:q_ncr}\\
& \mathbf{x}^s,\boldsymbol{\tau}^s,\mathbf{t}^s,\boldsymbol{\delta}^s, \boldsymbol{\gamma}^s \in \{0,1\}; \quad 0 \leq yp^s, yq^s \leq 1 \nonumber \}
% & 0 \leq yp^s, yq^s \leq 1 \} \label{eqn:cont}
\end{flalign}
}
\vspace{-1cm}

\subsection{Additional Investment Options}
\label{subsec:facts}
In this paper, we also include FACTS and PST devices and  as investment options. These devices are often useful for addressing congestion in overloaded transmission systems, a situation that can occur during major disruptions.
%and when the power system is outside normal operating conditions, which is typically the case after major disruptions. 
In addition, these devices are cost-effective  \cite{Baldick2014} and may reduce resiliency costs significantly by replacing the need for new transmission lines or hardening the existing damaged lines. 
%This paper explores this option but restricts attention to FACTS devices for series compensation and phase-shifting transformers.

\textbf{Series Compensators (FACTS)} We model series compensation with reactance reduction. As described in the nomenclature, $\delta^s_{ij}$ indicates if a compensation device is used on line $ij$ during scenario $s \in \Disasters$. To model this disjunction we introduce $(\tilde{G}_{ij}, \tilde{B}_{ij})$, a tuple representing the following conditional expression: 
\[
    (\tilde{G}_{ij},\tilde{B}_{ij})= 
\begin{cases}
    (\bar{G}_{ij},\bar{B}_{ij}),& \text{if used } (\delta^s_{ij}=1)\\       
    (G_{ij}, B_{ij}), & \text{otherwise } (\delta^s_{ij}=0)
\end{cases}
\]
This expression is modeled as constraints  \eqref{eqn:conduct_disj} and \eqref{eqn:suscept_disj}.
%where, $(\tilde{G}_{ij},\tilde{B}_{ij})$ admits reduced value ($\bar{G}_{ij},\bar{B}_{ij}$) on line $ij$ if the corresponding $\delta^s_{ij}=1.0$; otherwise, $(\tilde{G}_{ij},\tilde{B}_{ij})$ admits the original value.
Further, this disjunction is captured in the modified power flow equations \eqref{eqn:acopf_real_ij}-\eqref{eqn:acopf_reac_ji}. 
The non-convexity of this disjunction is discussed later in this section. 

\textbf{Phase-Shifting Transformers (PSTs)} PSTs are devices that adjust voltage phase angles. %They are particularly useful when lines have different susceptances. 
As described in the nomenclature, $\gamma^s_{ij}$ indicates if a phase shifter
is used on line $ij$ during scenario $s \in \Disasters$. This disjunction is described in the following conditional expression:
\[
    \theta^s_{ij}= 
\begin{cases}
    \begin{aligned}
    &\theta^s_{i}-\theta^s_{j}+\phi^s_{ij}, \\
    &-\phi^u \leq \phi^s_{ij} \leq \phi^u.
    \end{aligned} & \text{if installed } (\gamma^s_{ij}=1), \\
    \theta^s_{i}-\theta^s_{j}, &\text{otherwise } (\gamma^s_{ij}=0).
\end{cases}
\]
This expression is modeled through constraints  \eqref{eqn:ph_ang} and \eqref{eqn:phi_bnd}. 

% The phase angle difference $\theta_i^s - \theta_j^s$ can then be replaced with $\theta_i^s - \theta_j^s + \pstuse_{ij} \phi_{ij}^s$ where $\phi_{ij}^s$ represents the phase angle modification introduced by the phase shifter on line $ij$ during scenario $s \in \Disasters$. This problem is referred to as $\ORGDTdev$.


\subsection{Convex Relaxations for On/Off AC Power Flow}
\label{subsec:convex}
Numerous convex relaxations of AC power flow equations have been proposed in the literature to obtain tight lower bounds on the original nonconvex formulation. These relaxations have various trade-offs in complexity, quality, and performance. However, very few relaxation methods can be applied in the context of transmission switching problems as the addition of binary variables increases the complexity of the problems significantly. Hence, given the mixed-integer nature of ORGDT, which leads to an even more complex formulation with various resiliency options, we extend the recently developed convex quadratic relaxations of the power flow equations to ORGDT \cite{hijazi2014convex}. For completeness, we discuss the necessary convex relaxations for multilinear and trigonometric expressions which appear in nonconvex constraints \eqref{eqn:acopf_real_ij}-\eqref{eqn:acopf_reac_ji} and \eqref{eqn:acopf_loss3}. 


\label{subsec:qc}
\noindent
\textbf{Multilinear expressions} 
Given any two variables $x_i$, $x_j \in \left[\underline{x}_i \overline{x}_i\right]\times \left[\underline{x}_j\overline{x}_j\right]$, the  McCormick relaxation is used to linearize the bilinear product $x_ix_j$ by introducing a new variable $\widehat{x_{ij}} \in {\langle x_i, x_j \rangle}^{MC}$. The feasible region of this variable is defined by equations (\ref{eq:SMC}). This relaxation is exact if one variable is binary \cite{nagarajan2016tightening}. 
\vspace{-0.5cm}

{\fontsize{9}{8}\selectfont
\begin{subequations} \label{eq:SMC}
\small
\begin{align}
& \label{McCormick}\widehat{x_{ij}} \geq \underline{x}_ix_j+\underline{x}_jx_i -\underline{x}_i \hspace{2pt} \underline{x}_j \\ 
&\widehat{x_{ij}} \geq \overline{x}_ix_j+\overline{x}_jx_i - \overline{x}_i \hspace{2pt} \overline{x}_j \\ 
&\widehat{x_{ij}} \leq \underline{x}_ix_j+\overline{x}_jx_i-\underline{x}_i \hspace{2pt} \overline{x}_j \\ 
&\widehat{x_{ij}} \leq \overline{x}_ix_j+\underline{x}_jx_i-\overline{x}_i \hspace{2pt} \underline{x}_j
% & \underline{x}_i \leq x_i \leq \overline{x}_i, \ \underline{x}_j \leq x_j \leq \overline{x}_j
\end{align}
\end{subequations}%
}
\noindent
Multilinear expressions of the form $x_ix_j\ldots x_p$ are relaxed using a sequential bilinear approach as follows: $\langle \langle x_i,x_j  \rangle^{MC},\ldots, x_p\rangle^{MC}$. 

\noindent \textbf{Quadratic terms}
Given a variable $x_i \in \left[\underline{x}_i \overline{x}_i\right]$, a second-order conic relaxation is applied to $x_i^2$ by introducing a new variable $\widehat{x_{i}} \in {\langle x_i\rangle}^{MC-q}$, such that,
%
\vspace{-.25cm}
\begin{subequations} 
\label{eq:MC-q}
\small
\allowdisplaybreaks
\begin{align}
& \label{McC-q}\widehat{x_{i}} \geq x_i^2\\ 
&\widehat{x_{i}} \leq (\overline{x}_i+\underline{x}_i)x_i - \overline{x}_i\underline{x}_i
% & \underline{x}_i \leq x_i \leq \overline{x}_i
\end{align}
\end{subequations}%

By applying the above approaches and disjunctive convex relaxations for the trigonometric functions with on/off variables \cite{hijazi2014convex}, we outer-approximate  
the non-convex power flow constraints. For brevity, we focus on constraint \eqref{eqn:acopf_real_ij} but similar relaxations hold  for constraints \eqref{eqn:acopf_reac_ij}-\eqref{eqn:acopf_reac_ji}.

\vspace{-0.2cm}
{\fontsize{9}{8}\selectfont
\begin{align*}
p^s_{ij} = &(\overline{G}_{ij} - G_{ij})\left( \widehat{\delta xv}^s_{ij} - \widehat{\delta wc}^s_{ij}\right) + (\overline{B}_{ij} - B_{ij}) \widehat{\delta ws}^s_{ij} \\
&+ G_{ij} \left(\widehat{xv}^s_{ij} - \widehat{wc}^s_{ij}\right) - B_{ij} \widehat{ws}^s_{ij} \quad  \forall ij\in \Edges    
\end{align*}}
\vspace{-0.1cm}
where $\widehat{xv}^s_{ij},\widehat{wc}^s_{ij},\widehat{ws}^s_{ij},\widehat{\delta xv}^s_{ij},\widehat{\delta wc}^s_{ij}$ and $\widehat{\delta ws}^s_{ij}$ satisfy the following constraints:

\vspace{-0.1cm}
{\fontsize{9}{8}\selectfont
\begin{subequations}
\footnotesize
   \label{eq:relaxation}
\begin{align}
   &\widehat{v}^s_{i} \in \langle v^s_i\rangle^{MC-q}, \  \widehat{xv}^s_{ij} \in \langle \tilde{x}^s_{ij}, \widehat{v}^s_{i} \rangle^{MC} \\  
   &\widehat{w}^s_{ij} \in \langle v^s_i, v^s_j \rangle^{MC} \\
   &\widehat{wc}^s_{ij} \in \langle \widehat{w}^s_{ij}, \widehat{cs}^s_{ij} \rangle^{MC}, \
   \widehat{ws}^s_{ij} \in \langle \widehat{w}^s_{ij}, \widehat{sn}^s_{ij} \rangle^{MC} \\ 
   &\widehat{\delta wc}^s_{ij} \in \langle \delta^s_{ij}, \widehat{wc}^s_{ij} \rangle^{MC}, \ 
   \widehat{\delta ws}^s_{ij} \in \langle \delta^s_{ij}, \widehat{ws}^s_{ij} \rangle^{MC} \\ 
 &\widehat{cs}^s_{ij} \leq \tilde{x}^s_{ij} - \left(\frac{1- \cos(\theta^u)}{(\theta^u)^2}\right) ({\theta^s_{ij}}^2+(\tilde{x}^s_{ij}-1)(\theta^{M})^2) \label{eqn:cos_qc}\\
&\tilde{x}^s_{ij}\cos(\theta^u) \leq \widehat{cs}^s_{ij} \leq \tilde{x}^s_{ij}\\
& \widehat{sn}^s_{ij} \leq \cos(\theta^u/2)\theta^s_{ij}+\tilde{x}^s_{ij}(\sin(\theta^u/2)-\theta^u/2\cos(\theta^u/2)) \\ & \nonumber \hspace{20pt}+ (1-\tilde{x}^s_{ij})(\cos(\theta^u/2)\theta^{M} +1)\\
&\widehat{sn}^s_{ij} \geq \cos(\theta^u/2)\theta^s_{ij}-\tilde{x}^s_{ij}(\sin(\theta^u/2)-\theta^u/2\cos(\theta^u/2)) \\ & \nonumber \hspace{20pt}- (1-\tilde{x}^s_{ij})(\cos(\theta^u/2)\theta^{M} +1)\\
&\tilde{x}^s_{ij}\sin(-\theta^u) \leq \widehat{sn}^s_{ij} \leq \tilde{x}^s_{ij}\sin(\theta^u)
\end{align}
\end{subequations}}
Finally, we relax constraint \eqref{eqn:acopf_loss3} into a second-order conic constraint of the form ${p^s}^{2}_{ij} + {q^s}^{2}_{ij} \leq \loss_{ij} {\hat{v}^s}_{i}$.

\vspace{-0.18cm}
\subsection{DC approximation}
% {\color{blue} Harsha, did we allow FACTS and PST in the DC model? If so, it needs to beadded here.  If not we need to say why }
% {\color{orange}H: Pascal said we don't need to bring in FACTS devices for DC approximation. Though I had it ready in the code and had results too.}

%{\color{blue} This text may need some stronger sentences. I wrote a little bit as a place holder}
Within the expansion planning literature, the DC approximation is often used for designing power systems. For comparison purposes, 
we develop a version of the ORGDT problem based on the DC model in \eqref{eq:master_dc}. 
%Though the basic DC transmission switching problem has been studied well in the literature \cite{fisher2008optimal}, we extend it to the ORGDT problem in this paper. 
%Under the three basic assumptions such as negligible line resistances, small phase-angle differences on all lines and a constant unit voltage on all buses, mixed-integer nonlinear program in $\ORGDT$ reduces to a mixed-integer linear program as shown in formulation \eqref{eq:master_dc}. 


{\fontsize{9}{8}\selectfont
\begin{subequations}
\begin{flalign}
\SubProblemDC  &  = \{\mathbf{x}^s,\boldsymbol{\tau}^s,\mathbf{t}^s,\boldsymbol{\delta}^s, \boldsymbol{\gamma}^s, \mathbf{zp}^s, \mathbf{u} :  \nonumber \\
%----------------------------------------------------------------------
& \textbf{On/Off DC power flow equations:} \nonumber \\
& p^s_{ij} \leq -\tilde{B}_{ij} \left(\theta^s_{ij} + \theta^M(1 - \tilde{x}^s_{ij}) \right) \quad  \forall ij\in \Edges,\\
& p^s_{ij} \geq -\tilde{B}_{ij} \left(\theta^s_{ij} - \theta^M(1 - \tilde{x}^s_{ij}) \right) \quad  \forall ij\in \Edges,\\
& \tilde{G}_{ij} = (\bar{G}_{ij} - G_{ij})\delta^s_{ij} + G_{ij} \quad \quad \forall ij \in \Edges \label{eqn:dc_conduct_disj}, \\
& \tilde{B}_{ij} = (\bar{B}_{ij} - B_{ij})\delta^s_{ij} + B_{ij} \quad \quad \forall ij \in \Edges \label{eqn:dc_suscept_disj}, \\
& p^s_{ij} + p^s_{ji} = 0, \ \theta^s_{ij} =  \theta^s_{i} -  \theta^s_{j} + \phi^s_{ij} \quad  \forall ij\in \Edges, \\
% & \theta^s_{ij} =  \theta^s_{i} -  \theta^s_{j} + \phi^s_{ij} \quad  \forall ij\in \Edges, \label{eqn:dc_ph_ang}\\
%----------------------------------------------------------------------
% & \textbf{Flow Balance:} \nonumber \\
& gp^s_{i} - lp^s_{i} = \sum_{ij\in \Edges} p^s_{ij}  \quad \quad \forall i\in \Nodes,  \label{eqn:dc_balance_constraint1}\\
%----------------------------------------------------------------------
& \textbf{Operational limits and topology constraints:} \nonumber \\
& |p^s_{ij}| \leq \tilde{x}^s_{ij} \capacity_{ij} \quad \quad \forall ij\in \Edges,  \label{eqn:dc_line_constraint_ij} \\  
& |\theta^s_{ij}| \leq \tilde{x}^s_{ij}\theta^u + (1-\tilde{x}^s_{ij})\theta^M \quad \quad  \forall ij \in \Edges, \label{eqn:dc_theta_bound}\\
& -\phi^u \gamma^s_{ij} \leq \phi^s_{ij} \leq \phi^u \gamma^s_{ij} \quad \forall ij\in \Edges,  \label{eqn:dc_phi_bnd} \\
%----------------------------------------------------------------------
% & \textbf{Topology Constraints:} \nonumber \\
& x^s_{ij} = t^s_{ij} \quad \quad \forall ij\in \Damages_s  \label{eqn:dc_damage}\\
%& \tau^s_{ij} \leq x^s_{ij}  \quad \quad \forall ij\in \Edges  \label{eqn:dc_switchable}\\
& \tilde{x}^s_{ij} = x^s_{ij} - \tau^s_{ij}\geq 0 \quad \quad \forall ij\in \Edges  \label{eqn:dc_lineminusswitch}\\
%----------------------------------------------------------------------
& \textbf{Generation and demand constraints} \ \forall i \in \Nodes: \nonumber \\
& 0 \leq zp^s_{i} \leq zp^u_{i}u_{i}, \ 0 \leq gp^s_{i} \leq  gp^u_{i} + zp^s_{i}, \label{eqn:dc_microgrid_constraint_a}\\
% & 0 \leq gp^s_{i} \leq  gp^u_{i} + zp^s_{i}  \quad \quad \forall i\in \Nodes,  \label{eqn:dc_generation_constraint_a}\\
%----------------------------------------------------------------------
% & \textbf{Demand constraints:} \nonumber \\
% & lp^s_{i} =  yp^s_{i} dp_{i},  \label{eqn:dc_load_constraint_a}\\
& \sum_{i \in \Load} lp^s_{i} \geq lp_{cr} \sum_{i \in \Load} dp_{i},  \label{eqn:dc_critical_constraint1}\\
& \sum_{i \in \Nodes \setminus \Load} lp^s_{i} \geq lp_{ncr} \sum_{i \in \Nodes \setminus \Load} dp_{i} \label{eqn:dc_loadserve_constraint1}\\
& lp^s_{i} =  yp^s_{i} dp_{i}; \ 0 \leq yp^s_i \leq 1, \ \mathbf{x}^s,\boldsymbol{\tau}^s, \mathbf{t}^s \in \{0,1\}  \nonumber \}
\end{flalign}
\label{eq:master_dc}
\end{subequations}
}

\subsection{AC-feasibility formulation}
\label{subsec:ac_feas}

It is important to note that the convex relaxations introduced in section \ref{subsec:convex} 
and the DC approximation of \eqref{eq:master_dc}
%in addition to equations \eqref{eq:relaxation} 
might violate the nonconvex AC-equations described in
Eqs. \eqref{eqn:acopf_real_ij}-\eqref{eqn:acopf_reac_ji}, \eqref{eqn:acopf_loss3}. 
%Also, the DC approximation constraints in \eqref{eq:master_dc} can violate the original AC constraints due to numerous simplifying assumptions.
%In order to quantify the quality of these approximations, 
In this section, we develop a model for recovering primal feasible solutions to asses the quality of these approximations. This model fixes all the build decisions from the second stage solutions for every disaster $s \in \Disasters$ and solves for the feasibility of the nonconvex equations based on a primal-dual interior point method \cite{wright1997primal}. Formally, let $(\mathbf{\tilde{x}}^s,\mathbf{t}^s,\boldsymbol{\delta}^s, \boldsymbol{\gamma}^s, \mathbf{zp}^s, \mathbf{zq}^s,\mathbf{u})$ be the build decisions obtained in the first stage.
%by solving the MISOCP. 
Under the DC model, at bus $i$,
%In the case of DC approximation's build solutions, to accommodate for the reactive power capacity ($zq^s$), 
we assume that $zq^s_i$ is equal to half of the real power capacity, that is $|zq^s_i| = zp^s_i/2$. We then minimize a weighted sum of the violations to the resilience constraints which serves as a proxy to measure the distance to AC feasibility.
%solve the following problem to evaluate the minimum amount of load shed (critical and noncritical) required to keep the network AC feasible. 

% {\color{blue}How does $delta$ and $\gamma$ fit here?}
% {\color{orange} We again don't do feasibility analysis with FACTS devices.}
% {\color{blue}ok, but we should keep this general, but then remark in the results sections that the problems we do the analysis on don't have these devices}
\vspace{-0.14cm}
{\fontsize{8.5}{8}\selectfont
\begin{subequations}
\begin{flalign}
p_0 := \min \ \ & M\left(\lambda^p_{cr} + \lambda^q_{cr}\right) + \lambda^p_{ncr} + \lambda^q_{ncr}\nonumber \\
\text{s.t.}\quad & \mathrm{Given} \ (\mathbf{\tilde{x}}^s,\mathbf{t}^s,\boldsymbol{\delta}^s, \boldsymbol{\gamma}^s, \mathbf{zp}^s, \mathbf{zq}^s,\mathbf{u}) \nonumber\\ 
&\mathrm{Constraints} \ \eqref{eqn:acopf_real_ij}-\eqref{eqn:ph_ang}, \eqref{eqn:acopf_loss1}-\eqref{eqn:v_bnd}, \eqref{eqn:mg_a}-\eqref{eqn:shed}, \nonumber \\
& \lambda^p_{cr} \geq \left(lp_{cr} \sum_{i \in \Load} dp_{i}\right) - \sum_{i \in \Load} lp^s_{i} , \label{eqn:p_cr_slack}\\
& \lambda^p_{ncr}   \geq \left(lp_{ncr} \sum_{i \in \Nodes \setminus \Load} dp_{i}\right) - \sum_{i \in \Nodes \setminus \Load} lp^s_{i} \label{eqn:p_ncr_slack}\\
& \lambda^q_{cr}  \geq \left(lq_{cr} \sum_{i \in \Load} dq_{i}\right) - \sum_{i \in \Load} lq^s_{i} \label{eqn:q_cr_slack}\\
& \lambda^q_{ncr}  \geq \left(lq_{ncr} \sum_{i \in \Nodes \setminus \Load} dq_{i}\right) - \sum_{i \in \Nodes \setminus \Load} lq^s_{i} \label{eqn:q_ncr_slack}\\
&\lambda^p_{cr}, \lambda^q_{cr}, \lambda^p_{ncr}, \lambda^q_{ncr} \ge 0
\end{flalign}\label{eqn:primal_r}
\end{subequations}}
\vspace{-0.14cm}

\noindent
where $M$ is used to weight satisfaction of critical load and
and $\lambda^p_{cr}, \lambda^q_{cr}, \lambda^p_{ncr}, \lambda^q_{ncr}$ represent the load shed variables on critical and non-critical load constraints (up to the resilience criteria) as shown in constraints \eqref{eqn:p_cr_slack}-\eqref{eqn:q_ncr_slack}. 



% {\fontsize{9}{8}\selectfont
% \begin{subequations}
% \begin{flalign}
% p_0 := \min & \sum_{i=1}^n (lp^s_i - {lp^s_i}^{AC})_+ \nonumber \\
% \text{s.t.}\quad & ({lp^s_i}^{AC}, {lq^s_i}^{AC}) \in \text{ Eqs. } \eqref{eqn:acopf_real_ij} \text{ through }
% \eqref{eqn:voltage_bound} \nonumber \\
% & ({lp^s_i}^{AC}, {lq^s_i}^{AC}) \in \text{ Eqs. }\eqref{eqn:generation_constraint_a} \text{ through }
% \eqref{eqn:load_constraint_b} \nonumber
% \end{flalign}
% \end{subequations}}

Intro.tex
%!TEX root = ../pscc2016_orgdt.tex

%+++++++++++++++++++++++++++++;
%  Section: Introduction      ;
%+++++++++++++++++++++++++++++;
\section{Introduction}
\label{Sec:intro}

%{\color{blue} Russ note to self, restructure this to reflect focus on DC vs. AC modeling.}

\IEEEPARstart{T}{he} modern electrical system is designed for transportation of large amounts of power from sources of supply to distant points of demand. Within these systems, the underlying high-voltage transmission networks play a vital role in achieving this mission. However, when transmission networks are exposed to extreme event conditions, the ability to deliver power is degraded because of physical damage to overhead 
transmission lines and towers. One example of such events are ice storms. During an ice storm, transmission towers can fail due to leg buckling
%(Figure \ref{Fig:ice_storm}) 
and lines can fail due to the combined stress of ice accumulation and wind \cite{Brostrom2007,Eidinger,wang2010influence}. 

When such events occur on large scales, outages and impacts can be extreme. For example, in the winter of 1998, an ice storm in northeastern North America toppled over 1000 transmission towers and 30,000 wooden utility poles. Over 5 million people were without power and the economic impacts were estimated at \$2.6 billion \cite{Mills2012}. Thus, given the potential social and economic impacts of these events, it is important to consider how to upgrade the design of transmission systems to improve their performance under these conditions.


In our preliminary work \cite{7540988}, we formulated the Optimal Resilient Grid Design problem for Transmission systems (ORGDT) 
as a two-stage mixed-integer stochastic optimization problem and developed an algorithm 
to solve this problem. The first stage selects from a set of potential upgrades to the network. The
second stage evaluates the network performance benefit of the upgrades 
with a set of damage scenarios sampled from a stochastic distribution of events of concern. 
One of the important questions that was unanswered in this early work centered on the use of recently developed convex relaxations of AC power flow physics \cite{hijazi2014convex} in resiliency planning. Here, we address this concern and convincingly show that using DC approximation of power flow physics leads to severely suboptimal and infeasible solutions under certain operating conditions of the grid. These results establish the necessity of using models of AC physics in resilient design problems and analysis.  

Like \cite{7540988}, we adopt the methods discussed in \cite{wang2010influence} in order to sample 
realistic damage scenarios for transmission systems. The ORGDT upgrade options include: a) Build new lines, b) Build switches, FACTS devices and Phase-shifting transformers (PST) to provide operational flexibility, c) Harden existing lines to lower the 
probability of damage, and d) Build new generation facilities.
Minimal network (resiliency) performance is measured by satisfying a minimum fraction of critical and non-critical loads. We use the exact algorithms developed in \cite{7540988} to compare solutions found using the DC approximation and the convex quadratic relaxations of \cite{hijazi2014convex}.

%For the purposes of this study, we adopt the methods discussed in \cite{wang2010influence} in order to sample 
%realistic damage scenarios for transmission systems. The ORGDT upgrade options include: a) Build new lines, b) Build switches and FACTS devices to provide operational flexibility, c) Harden existing lines to lower the 
%probability of damage, and d) Build distributed generation facilities.
%Minimal network (resiliency) performance is measured by satisfying a minimum fraction of critical and non critical load for each scenario under the convex quadratic AC power flow constraints discussed in \cite{hijazi2014convex}. Given that this problem is a non-trivial Mixed Integer Problem (MIP), we develop exact and heuristic decomposition algorithms  exploit the block diagonal structure of the MIPs. 

\noindent \textbf{Literature Review} 
Our recent paper \cite{7540988} is the most closely related work. This paper develops the ORGDT model and algorithm used here and is based on the resilient distribution system design work developed in \cite{yamangil2015resilient}. 
%Our recent paper \cite{yamangil2015resilient} is the most closely related work, that we are aware of, in the literature.
%This paper considers the problem of 
%resilient upgrades of electrical distribution grids (radial networks). Their approach minimizes the upgrade 
%budget while satisfying a minimum standard of service by selecting from a set of 
%various potential upgrades. Like us, they pose the optimal grid design problem as 
%a two-stage, stochastic mixed-integer program with damage scenarios. They describe decomposition based optimal and heuristic algorithms that are the basis for our approach.
%We generalize their approaches in three fundamental ways. First, we extend the optimization framework to transmission networks which are non-radial by nature. Second, we add AC power flow physics, an important feature of loopy networks under stressed conditions \cite{Coffrin2012}.
%Third, we consider the use of transformers and FACTS devices to improve operations under stressed conditions and improve resilience. 

Another important area of related work is interdiction modeling and optimization. Here, the
goal is to operate or design a system to make it as resilient as possible to an adversary who can damage up to $k$ elements \cite{Chen:13,Chen2014a,Salmeron2009,Delgadillo2010}. These models are a generalization of our model when $k$ is chosen to bound a worst-case event.
However, given their min-max structure, these models are computationally challenging and are solvable only for small $k$. Instead, we exploit the
probabilistic nature of the adversary and we are able to address larger problems. In addition, existing interdiction models also do not generally include AC physics. 
Closely related to interdiction and our model is reference \cite{nezamoddini2017risk}.  They develop a model for hardening and upgrading transmission systems that is resilient to all possible physical damages bounded by a certain budget (i.e. weighted $k$). Instead of developing an attacker-defender model as in interdiction modeling, 
they include each possible damage as a scenario (like our model). Given the combinatorial number of possible scenarios, they develop heuristics for reducing the number of scenarios that are included. Unlike our approach, they focus only on the DC approximation and leave the development of decomposition algorithms for future work.
A third area of related work is stochastic transmission and generation expansion planning, where a recent survey describes some of the state-of-the-art \cite{Krishnan2015}. 

Overall, most papers in network design and expansion planning use the linearized DC model and few studies consider FACTS devices and PSTs, although they may have significant benefits. Some notable exceptions include the use of PSTs in network expansion \cite{miasaki2012transmission}, which uses
a genetic algorithm over the DC model. See also the recent work in \cite{Frolov2014} for the optimal placement of these devices to avoid congestion. 

The key contributions of this paper include:

\begin{itemize}
\item A detailed case study and extensions of the ORGDT model (from \cite{7540988}) that compares DC approximation based solutions and solutions based on state-of-the-art AC power flow relaxations. 
\item An approach for recovering AC feasible solutions from solutions found using the DC approximation and AC relaxations.
\item A detailed expansion model with appropriate generator capacity modeling, FACTS devices and PST devices, and cutting-plane-based algorithmic enhancements that was omitted from \cite{7540988}.
%\item Computationally efficient algorithms for designing resilient transmission systems.
%\item An analysis of the resiliency benefits of transmission-level control devices. We show that the stability criteria related to phase angle differences is a driving force in costs associated with resilient design and FACTS devices can be used to control these costs.
\end{itemize}


% RBENT -- drop this if we need space
The rest of this paper is organized as follows. In section \ref{Sec:formulation} we describe the ORGDT model. Section \ref{Sec:algo} describes our decomposition algorithms and Section \ref{Sec:case} our case studies. We conclude with Section \ref{Sec:conc}.

Nomenclature.tex
%!TEX root = ../pscc2016_orgdt.tex
%+++++++++++++++++++++++++++++;
%  Section: Nomenclature      ;
%+++++++++++++++++++++++++++++;

\nomenclature[AA]{$\Nodes$}{set of nodes (buses)}%
\nomenclature[AB]{$\Edges$}{set of edges (lines and transformers) %When context is clear, we refer to an edge as $i,j$, where $i$ and $j$ are the nodes it is connected to  However, more than one edge may exist between nodes
}%
\nomenclature[AC]{$\Disasters$}{set of disaster scenarios}%
\nomenclature[AD]{$\Damages_s$}{set of edges that are inoperable during $s\in S$}
\nomenclature[AE]{$\mathbf{i}$}{imaginary number constant}
\nomenclature[AEb]{$\mid.\mid$}{absolute value of the input argument}
%\nomenclature[AF]{$\HardenedDamages_s$}{set of edges that are inoperable even though they are hardened during disaster $s\in S$}
\nomenclature[AF]{$c^{x}_{ij}$}{cost to build a line ($i$,$j$) if the line does not already exist}
\nomenclature[AG]{$c^{\tau}_{ij}$}{cost to build a switch on line ($i$,$j$)}
\nomenclature[AH]{$c^t_{ij}$}{cost to harden a line ($i$,$j$)}
\nomenclature[AHb]{$c^{\delta}_{ij}$}{cost to build a FACTS device on line ($i$,$j$)}
\nomenclature[AHc]{$c^{\gamma}_{ij}$}{cost to build a PST device on line ($i$,$j$)}
\nomenclature[AI]{$c^{zp}_{i}$}{cost of real generation capacity at bus $i$}
\nomenclature[AJ]{$c^u_{i}$}{cost to build a generation facility at bus $i$}
\nomenclature[AK]{$T_{ij}$}{apparent power thermal limit on line ($i$,$j$)}
\nomenclature[AL]{$\Load$}{set of buses whose loads are critical}
\nomenclature[AM]{$G_{ij}+\mathbf{i}B_{ij}$}{admittance of line ($i$,$j$)}
\nomenclature[AN]{$R_{ij}+\mathbf{i}X_{ij}$}{impedance of line ($i$,$j$)}
\nomenclature[AM1]{$\bar{G}_{ij}+\mathbf{i}\bar{B}_{ij}$}{modified admittance of line ($i$,$j$) due to top transformer}
\nomenclature[AN1]{$R_{ij}+\mathbf{i}\bar{X}_{ij}$}{modified impedance of line ($i$,$j$) due to top transformer}
\nomenclature[AO]{$\theta^u$}{phase angle difference limit}
\nomenclature[AP]{$\theta^M$}{big-M value on phase angle difference limit}
\nomenclature[AQ]{$\phi^u$}{phase shift limit}
\nomenclature[AR]{$v^l_i,v^u_i$}{lower and upper bound on voltage at bus $i$, respectively}
\nomenclature[AS]{$dp_{i} + \mathbf{i} dq_{i}$}{AC power demand at bus $i$}
% \nomenclature[AO]{$\demandreac_{i}$}{reactive power demand at bus $i$}
\nomenclature[AT]{$gp^u_{i} + \mathbf{i} gq^u_{i}$}{max. existing AC generation capacity at bus $i$}
\nomenclature[AU]{$zp^u_{i} + \mathbf{i} zq^u_{i}$}{max. AC generation capacity that can be built at bus $i$}
\nomenclature[AW]{$lp_{ncr} + \mathbf{i}lq_{ncr}$}{fraction of non-critical AC power loads that must be served}
\nomenclature[AV]{$lp_{cr} + \mathbf{i} lq_{cr}$}{fraction of critical AC power loads that must be served}
% \nomenclature[AW]{$\criticalLoadreac$}{fraction of critical reactive power loads that must be served}



\nomenclature[BA]{$x_{ij}$}{determines if line ($i$,$j$) is built}%
\nomenclature[BB]{$\tau_{ij}$}{determines if line ($i$,$j$) has a switch built}
\nomenclature[BC]{$t_{ij}$}{determines if line ($i$,$j$) is hardened}
\nomenclature[BD]{$u_{i}$}{determines the generation capacity built at bus $i$}
\nomenclature[BE]{$\delta_{ij}$}{determines if FACTS device is built on line ($i$,$j$)}
\nomenclature[BF]{$\gamma_{ij}$}{determines if PST device is built on line ($i$,$j$)}
\nomenclature[BG]{$x^s_{ij}$}{determines if line ($i$,$j$) is used during disaster $s \in \cal{S}$}
\nomenclature[BH]{$\tau^s_{ij}$}{determines if switch ($i$,$j$) is used during disaster $s \in \cal{S}$}
\nomenclature[BI]{$t^s_{ij}$}{determines if line ($i$,$j$) is hardened during disaster $s \in \cal{S}$}
\nomenclature[BJ]{$\delta^s_{ij}$}{determines if FACTS device on line ($i$,$j$) is used during disaster $s \in \cal{S}$}
\nomenclature[BK]{$\gamma^s_{ij}$}{determines if PST device on line ($i$,$j$) is used during disaster $s \in \cal{S}$}



\nomenclature[CAa]{$\theta^s_{i}$}{phase angle at bus $i$ during disaster $s \in \cal{S}$}
\nomenclature[CAb]{$\phi^s_{ij}$}{phase angle shift introduced by PST on line $ij$ during disaster $s \in \cal{S}$}
\nomenclature[CB]{$v^s_{i}$}{voltage at bus $i$ during disaster $s \in \cal{S}$}
\nomenclature[CC]{$l^s_{ij}$}{current magnitude squared ($\lvert I_{ij}^s\rvert^2$) on line ($i$,$j$) during disaster $s \in \cal{S}$}
\nomenclature[CD]{$p^s_{ij}+\mathbf{i}q^s_{ij}$}{AC power flow on line ($i$,$j$) during disaster $s \in \cal{S}$}
\nomenclature[CE]{$zp_{i} +\mathbf{i} zq_{i}$}{determines the capacity for AC generation at bus $i$}
\nomenclature[CF]{$zp^s_{i}+\mathbf{i}zq^s_{i}  $}{determines the capacity for AC generation at bus $i$ during disaster $s \in \cal{S}$}
\nomenclature[CG]{$gp^s_{i} + \mathbf{i}gq^s_{i} $}{AC power generated at bus $i$  during disaster $s \in \cal{S}$}
\nomenclature[CH]{$lp^s_{i} + \mathbf{i} lq^s_{i}$}{AC power load delivered at bus $i$ during disaster $s \in \cal{S}$}
\nomenclature[CI]{$yp^s_{i} +\mathbf{i}yq^s_{i}$}{determines the fraction of AC power load served at bus $i$ during disaster $s \in \cal{S}$}


% \nomenclature[CD]{$\microgridVariablereac_{i}$}{determines the capacity for reactive generation at bus $i$}
% \nomenclature[CB]{$\microgriduseVariablereac_{i}$}{determines the capacity for reactive generation at bus $i$ during disaster $s \in \cal{S}$}
% \nomenclature[CE]{$\facilityuseVariable_{i}$}{indicates if the generation capacity is used at bus $i$ during disaster $s \in \cal{S}$}
% \nomenclature[CH]{$\generatorVariablereac_{i}$}{reactive power generation produced at bus $i$  during disaster $s \in \cal{S}$}
% \nomenclature[CI]{$\loadVariablereac_{i}$}{reactive power load delivered at bus $i$ during disaster $s \in \cal{S}$}
% \nomenclature[CJ]{$\loadServeVariablereac_{i}$}{determines if the reactive load at bus $i$ is served or not during disaster $s \in \cal{S}$}
% \nomenclature[CM]{$\reacflowVariable_{ij}$}{reactive flow between buses $i$ and $j$ during disaster $s \in \cal{S}$}
%\nomenclature[BP]{$\lineCycleVariable_{ij}$}{determines if at least one edge between $i$ and $j$ is used during disaster $s \in \cal{S}$}
%\nomenclature[BQ]{$\switchCycleVariable_{ij}$}{determines if at least one switch between $i$ and $j$ is used during disaster $s \in \cal{S}$}
% \nomenclature[CP]{$\lineuseVariable_{ij,0}$}{determines if there exists flow on line $(i,j)$ from $j$ to $i$, during disaster $s \in \cal{S}$}
% \nomenclature[CQ]{$\lineuseVariable_{ij,1}$}{determines if there exists flow on line $(i,j)$ from $i$ to $j$, during disaster $s \in \cal{S}$}

NR.tex
%!TEX root = ../orgdt_tps.tex

%+++++++++++++++++++++++++++;
%  Section: Case study      ;
%+++++++++++++++++++++++++++;
\section{Case Study}
\label{Sec:case}
In this section, we discuss three numerical studies for the ORGDT based on the modified single area IEEE RTS-96 system \cite{wong1999ieee}.
%demonstrate the benefits of the algorithms when solving the ORGDT to optimality and sub-optimality on the modified single area IEEE RTS-96 system \cite{wong1999ieee}. 
We first discuss the computational benefits associated with using the SBD algorithm. We second compare solutions based on the DC approximation with solutions based on the QC relaxation. Finally, we examine the benefits of FACTS and PST devices. 
%The comparison is based on the upgrade costs for resilient transmission grids and the additional cost benefit obtained due to the series compensators and phase-shift transformers in the system. We also investigate the performance of proposed algorithms with regards to computation time and scalability.
All modeling was done using JuMP \cite{dunning2015jump} on a computer with 32 threads, a 2.6GHz Intel 64 bit processor, 25.6MB L3 cache and 64GB of memory. Gurobi 6.5.2 was used for solving the MIQCPs and MILPs (optimality gap 0.1\%) and Knitro 9.1.0 was used for solving the nonlinear programs.
%All computations were implemented with KNITRO 9.1.0 as a nonlinear solver and Gurobi 6.5.2 as a MISOCP/MILP solver with optimality gap (MIPGap) set to 0.1\%. We used 32 threads on a 2.6GHz Intel 64 bit processor, 25.6MB L3 cache and 132GB of memory. All modeling was done using JuMP \cite{LubinDunningIJOC}.

%We also plan to investigate the effect of the topology and spatial distribution of the network in obtaining a trade-off between the resiliency criteria and total system damage.
%\vspace{-.1cm}
% \subsection{Test system\label{sec:testdata}}
\noindent
\textbf{Test system} We use a modified IEEE single-area RTS-96 system that has 24 buses including 17 load buses, 38 transmission lines and 32 conventional generators \cite{wong1999ieee}. The total installed capacity of the existing generators is 3405 MW. The total load in the system is 2850 MW. We labeled 1740 MW as critical loads. The network was spatially placed in an area of 4250 miles$^2$. Since this test system does not have the co-ordinate data for every bus, we approximately evaluated them based on the lengths of lines \cite{wong1999ieee}.
The admittance, impedance and apparent power thermal limit values on lines are from the standard test case. We assume that FACTS, PSTs and new generation capacities, if chosen, can be built anywhere in the network.  
The remaining parameters are outlined in Table \ref{tab:param}. 
% {\color{blue} We need a remark on where new lines are allowed, and parameters, (Maybe in Table 1), where FACTS devices are allowed, PST devices are allowed, new generation is allowed, etc.}

\begin{table}
\caption{Parameters for the test system}
\vspace{-0.2cm}
\footnotesize
\centering
% \resizebox{0.5\textipwidth}{!}{%
\begin{tabular}{rl}
\toprule
 $c^x_{ij},c^t_{ij},c^{\tau}_{ij}$ & \$1.35m/mile, \$5000/mile, \$1000 \\
 $c^{\delta}_{ij},c^{\gamma}_{ij}$ & \$50000, \$50000\\
 $c^u_i, c^{zp}_i$ & \$0.1m, \$0.817m/MW\\
 $v^l_i,v^u_i,\phi^u,\bar{X}_{ij}$ & 0.95, 1.05, $60\degree$, $0.5X_{ij}$\\
 $lp_{cr},lq_{cr},lp_{ncr},lq_{ncr}$ & 0.99, 0.99, 0.8, 0.8\\
 Total no. of disaster scenarios & 20 \\
 Critical load buses & 1,5,8,11,13,14,15,16,18,19 \\
 New lines which can be built & (1,10),(10,20),(19,24),(6,14),(20,21)\\
 \bottomrule
\end{tabular}
% }
\label{tab:param}
\vspace{-0.4cm}
\end{table}

\noindent
\textbf{Scenario generation}
The damage scenarios are based on line failure probabilities that follow a multivariate Gaussian distribution with a mean placed at the center of the network. A Bernoulli trial was applied to every line for chosen percentiles (\% damage) to generate the random scenarios.
Empirically, we observed that 20 scenarios were sufficient to represent the features of the distribution.


\subsection{Computational Performance of SBD}
\label{subsec:sbd_perf}
Table \ref{tab:times} compares the computational time of SBD 
%(in addition to cutting plane procedure described in Algorithm \ref{algo:cp})
with the computational time of solving the full model using a MIQCP (Gurobi) solver. Here, we focus on the problem where only new lines, hardened lines, and distributed generations are available as design choices.
%solving time of the full model using an MIQP solver with all the scenario constraints. 
%To illustrate the complexity of the full optimization model formulated as MISOCP, here are the statistics: 54724 constraints (78725 with FACTS), 19649 continuous variables (24649 with FACTS) amd 4144 binary variables (5184 with FACTS). For simplicity, we exclude time comparisons for the model with FACTS devices. 
On average, SBD is 17x faster when $\theta^u=15\degree$ and 26x faster when $\theta^u=45\degree$.
%The computation time speed-up is evident from the table, where SBD is 17 times faster (on an average) in the case of small phase angle differences ($\theta^u=15\degree$). SBD's superior performance in run times is even more evident in the case of larger phase angle differences ($\theta^u=45\degree$), where SBD is around 26 times faster (on an average). 
We highlight in the 80\% damage case that after 20 hours the MIQCP solver could only provide a 7\% optimality gap,
%solving the full model for 20 hours could provide only 7\% optimality gap, 
whereas SBD terminated with an optimal solution after 12 minutes. 
%This is also evident from the fact that the total number of binding scenarios for various damages thresholds were at most ``five" in both phase angle difference cases. 
It is important to note that when FACTS and PST devices are included as design options, the MIQCP solver fails to find any feasible solution after 10 hours of computation. 
%that the MISOCP solver (Gurobi) failed to find a feasible solution even after 10 hours by solving the full model with additional FACTS devices, thus necessitating for a tractable algorithm like SBD. 


\begin{table}[h]
\caption{Comparison of wall time (sec.) of SBD and MIQCP solvers on problems without FACTS and PST devices. * indicates that the best gap was 7\%. The first column indicates the level of damage. The second/fifth and fourth/seventh columns indicate the wall times of SBD and MIQCP, respectively. The third and sixth columns show the number of iterations for SBD }
\vspace{-0.2cm}
\scriptsize
\centering
% \resizebox{0.5\textipwidth}{!}{%
\begin{tabular}{crcrrcr}
\toprule
 & \multicolumn{3}{c}{$\theta^u=15\degree$} & \multicolumn{3}{c}{$\theta^u=45\degree$} \\
 \cmidrule(lr){2-4}
 \cmidrule(lr){5-7}
 Damage & SBD & Iters & Full & SBD & Iters & Full \\
% &&scenarios&&&scenarios& \\
 \midrule
 90\% & 634.8 &3 & 5008.2 & 2419.6 & 4 & 13729.2  \\
 80\% & 303.7 & 5 & 1508.7  & 730.1 &5 & 72886.3$^*$ \\
 70\% & 417.7 & 4 & 6566.9 & 485.1 & 4& 2339.0\\ 
 60\% & 188.3 & 2 & 8853.8  & 157.6 & 3 & 3537.5\\
 50\% & 90.2 & 2 & 6156.7  & 53.1 & 1 & 9372.8\\
 \bottomrule
\end{tabular}
% }
\label{tab:times}
\vspace{-0.2cm}
\end{table}


\subsection{AC vs. DC power flow models}
\label{subsec:AcDc}
In this section we scale the loads of the model between 0.8 and 1.25 at .05 intervals. We define the gap between the solutions based on the QC relaxation and solutions based on the DC approximation as

%Table \ref{tab:gaps} summarizes the $AC-DC$ solution gaps.
%evaluated over various load scaling factors ranging from 0.8x up to 1.25x at 0.5x intervals. We define $AC-DC$ gap as shown below, where opt($\mathcal{P}^{AC}_o$) and opt($\mathcal{P}^{DC}_o$) represent the optimal total upgrade costs for QC relaxation and DC approximation models, respectively. 

{\fontsize{6}{6}\selectfont
$$\zeta = \left(\frac{\mathrm{opt}(\mathcal{P}_0^{AC}) - \mathrm{opt}(\mathcal{P}_0^{DC})}{\mathrm{opt}(\mathcal{P}_0^{AC})}\right)*100$$}
\noindent where opt($\mathcal{P}^{AC}_o$) and opt($\mathcal{P}^{DC}_o$) are optimal upgrade costs for QC relaxation and DC approximation models, respectively.
As shown in table \ref{tab:gaps}, DC-approximation-based solutions are always lower bounds to solutions based on QC relaxations. In the case of small phase angle differences ($\theta^u=15$), the maximum and average $\zeta$ is at the maximum damage case (90\%). When the phase angle differences are increased, the value of $\zeta$ increases dramatically.
%However, in the case of large phase angle differences ($\theta^u=45$), the average $\zeta$ raises substantially, where the DC-based upgrade under approximates the QC-based upgrade cost up to 55\%. Moreover, the average gap between these upgrades are always greater than 22\%. 
This behavior in the solutions is not unexpected. DC approximations tend to perform poorly when there are large phase angle differences as the assumptions behind the approximation no longer hold. 
%which is contrary to one of the assumptions based on small phase angle differences. 
In Figure \ref{fig:gaps}, we expand the results of the 90\% damage row of Table \ref{tab:gaps} and see how the DC approximation consistently under estimates the design costs.
%For the case of 90\% damage in Table \ref{tab:gaps}, Figure \ref{fig:gaps} illustrates the the actual upgrade cost curves plotted at various load scales. 
\vspace{-0.3cm}
\begin{table}[h]
\caption{Each row describes the average, min, and max $\zeta$ across load scaling values between 0.8 and 1.25x for a specified damage level.
}
\vspace{-0.2cm}
\scriptsize 
\centering
% \resizebox{0.5\textipwidth}{!}{%
\begin{tabular}{ccccccc}
\toprule
 & \multicolumn{3}{c}{$\theta^u=15\degree$} & \multicolumn{3}{c}{$\theta^u=45\degree$} \\
 \cmidrule(lr){2-4}
 \cmidrule(lr){5-7}
 Dam. & avg & min & max & avg & min & max \\
 &gap (\%)&gap (\%)&gap (\%)&gap (\%)&gap (\%)&gap (\%) \\
 \midrule
 90\% & 9.7& 0 & 35.4 & 22.4& 5.5&47.3 \\
 80\% & 9.2& 0 & 25.3& 22.1& 11.7&37.6 \\
 70\% & 7.0& 0& 29.3& 24.0& 15.4&55.7 \\ 
 60\% & 5.6& 0& 29.4& 24.9& 7.4&43.1\\
 50\% & 2.1& 0& 11.3& 35.9& 0&53.1 \\
 \bottomrule
\end{tabular}
% }
\label{tab:gaps}
\end{table}
\vspace{-0.5cm}
\begin{figure}[htp]
	\centering
	\subfigure[90\% damage, $\theta^u=15^o$]{
	\includegraphics[scale=0.484]{acdc_15deg_p1}}
	%
	\subfigure[90\% damage, $\theta^u=45^o$]{
	\includegraphics[scale=0.484]{acdc_45deg_p1}}
	\vspace*{-0.18in}\\
	\caption{Comparison of total upgrade cost for the QC relaxation with the DC approximation}
	\label{fig:gaps}
\end{figure}

\noindent \textbf{AC feasibility analysis} $\zeta$ is not the only metric for comparing the solutions based on the DC and QC models.  It is also important to understand distance to AC feasibility (section \ref{subsec:ac_feas}) in each scenario.
%Based on the formulation discussed in section \ref{subsec:ac_feas}, we perform a detailed study of the quality of upgrade options obtained based on the QC relaxation and DC approximation. 
%In order to provide a compact representation, 
For a given scenario $s$, we calculate the percent apparent-load shed with: 

{\fontsize{3}{3}\selectfont
$$\mu = \left(\frac{\sqrt{(\lambda^p_{cr})^2 + (\lambda^q_{cr})^2}}{\sqrt{\left(\sum_{i\in \mathcal{L}}dp_i\right)^2+\left(\sum_{i\in \mathcal{L}}dq_i\right)^2}}, \frac{\sqrt{(\lambda^p_{ncr})^2 + (\lambda^q_{ncr})^2}}{\sqrt{\left(\sum_{i\in \mathcal{N}\setminus\mathcal{L}}dp_i\right)^2+\left(\sum_{i\in \mathcal{N}\setminus\mathcal{L}}dq_i\right)^2}}\right) \times 100$$}
\noindent for critical and noncritical loads, respectively. 

Table \ref{tab:shed} summarizes the average and maximum $\mu$ for each scenario for the solution found with the DC approximation for a load scale of 1.0.
%apparent-load sheds evaluated based on the above mentioned equations for DC approximation upgrades. Average values were evaluated over the total number of damage scenarios (20).
The QC relaxation is tight and no (additional) load shedding is required. 
%In the case of QC relaxation upgrades, the fact that this convex relaxation is very tight manifests in the form of identically ``zero percent" sheds for every damage scenario and damage probability thresholds on both critical and non-critical loads (thus we do not tabulate this data). 
%However, as expected, the DC approximation upgrades lead to both critical and non-critical load sheds as seen in table \ref{tab:shed}. 
Though the DC approximation is considered a good approximation when phase angle differences are small,
%known to perform well on smaller phase difference limits, 
we still observe load sheds of up to 3.7\% of critical loads and 6.7\% of non-critical loads in the maximum damage case. Again, the DC approximation's performance degrades further for larger phase angle differences
%at larger phase angle difference limits is very poor 
with a maximum $\mu$ of 4.5\% (critical) and 9.4\% (non-critical). 
Figure \ref{fig:shed} further expands on $\mu$ for each scenario in the 90\% damage case. 
%For a better clarity on the performance of DC approximation at every damage scenario, we plot the apparent-load shed for the 90\% damage case (row 1 of table \ref{tab:shed}) in Figure \ref{fig:shed}. 

To summarize, OGRDT solutions based on the DC approximation perform reasonable well only on systems with small phase angle differences and limited damage. When the grid is damaged heavily, DC approximation performs poorly, even at smaller $\theta^u$ values. At larger $\theta^u$ values, the DC approximation drastically under approximates the required upgrades, and results in solutions that severely violate resilience criteria based on load shedding.
%thus leading to higher load shedding. 
However, these issues can be addressed by applying AC relaxations in lue of DC approximations. 
\vspace{-0.2cm}
\begin{table}[htp]
\caption{Average and maximum $\mu$ (critical, non-critical) for DC solutions and a load scale of 1.0.}
\vspace{-0.2cm}
\scriptsize 
\centering
% \resizebox{0.5\textipwidth}{!}{%
\begin{tabular}{cl|l|l|l}
\toprule
 & \multicolumn{2}{c}{$\theta^u=15\degree$} & \multicolumn{2}{c}{$\theta^u=45\degree$} \\
 \cmidrule(lr){2-3}
 \cmidrule(lr){4-5}
 Damage & avg. $\mu$ & max. $\mu$ & avg. $\mu$ & max $\mu$ \\
% &shed (\%)&scenario (\%)&shed (\%)&scenario (\%) \\
 \midrule
 90\% & 0.3, 0.7 & 3.7, 6.7 & 0.5, 1.6 &4.5, 8.9  \\
 80\% & 0.1, 1.2& 0.3, 7.3 & 0.2, 3.0& 1.7, 9.2 \\
 70\% & 0.1, 0.1& 0.3, 2.7& 0.15, 2.5& 0.9, 9.4  \\ 
 60\% & 0.1, 0.1& 0.3, 0.2& 0.13, 2.4&0.3, 9.0 \\
 50\% & 0, 0& 0, 0.03& 0, 1.8&0, 7.5 \\
 \bottomrule
\end{tabular}
% }
\label{tab:shed}
\end{table}
\vspace{-0.3cm}
\begin{figure}[!h]
	\centering
	\subfigure[90\% damage, $\theta^u=15^o$]{
	\includegraphics[scale=0.503]{test_ang15_p1.pdf}}
	%
	\subfigure[90\% damage, $\theta^u=45^o$]{
	\includegraphics[scale=0.503]{test_ang45_p1.pdf}}
	\vspace*{-0.11in}\\
	\caption{Values of $\mu$ for the QC relaxation and DC approximation for 90\% damage scenarios when the load scale is 1.0.}
	\label{fig:shed}
\end{figure}

% $$\% \ \mathrm{savings} = \left(\frac{\mathrm{Upgrade\ without\ devices} - \mathrm{Upgrade\ with\ devices}}{\mathrm{Upgrade\ without\ devices}}\right)*100$$
\vspace{-0.5cm}
\subsection{Benefits of FACTS and PST devices}
In this section, we summarize the benefits of including FACTS and PST device options in the ORGDT. To quantify the benefits, we define the following: 

{\fontsize{6}{6}\selectfont
$$\psi = \left(\frac{\mathrm{opt}(\mathcal{P}_0^{AC}) - \mathrm{opt}(\mathcal{P}_0^{AC-dev})}{\mathrm{opt}(\mathcal{P}_0^{AC})}\right)\times 100$$}

\noindent
\textbf{FACTS devices}
Table \ref{tab:facts} summarizes the benefits of including FACTS devices to achieve advantages of operating transmission grids at larger $\theta^u$ values without the stability issues. Overall, we observe that savings is a much as 17.2\% when $\theta^u=15\degree$.
%it can be observed that the savings in the total upgrade cost is better at $\theta^u=15\degree$, ranging from 9.3\% up to 17.2\% (averaged over load scales). As indicated under max. savings in Table \ref{tab:facts}, values with parenthesis represent the load scale values at which the maximum saving occurs. 
%As expected, it may be noted that 
The benefit of FACTS devices is smaller when $\theta^u$ is larger, since the operating points are less congested and the need for such devices is reduced.
%value and at larger load scale values ($>$1.2), which is an artifact of congested networks. At $\theta^u=45\degree$, since the operating conditions are lesser congested, though not necessarily stable, the average benefit of FACTS devices is lesser and moreover, the maximum savings occur at only 1.25 load scale. 
We use Figure \ref{fig:facts} to expand on $\psi$ for the 90\% damage case.
%For better clarity in the performance of FACTS devices at various load scales, Figure \ref{fig:facts} illustrates the trend of upgrade costs and \% savings for the 90\% damage case in Table \ref{tab:facts}.

\begin{table}[htp]
\caption{Average and maximum (parenthesis indicate load scale value for the max) $\psi$ across load scales between 0.8 and 1.25 when FACTS devices are included in the OGRDT.
%in the total upgrade cost based on QC relaxations with FACTS device options. Note that the analysis in this table is over load scales from 0.8x up to 1.25x at 0.5 intervals.
}
\vspace{-0.2cm}
\scriptsize 
\centering
% \resizebox{0.5\textipwidth}{!}{%
\begin{tabular}{cl|l|l|l}
\toprule
 & \multicolumn{2}{c}{$\theta^u=15\degree$} & \multicolumn{2}{c}{$\theta^u=45\degree$} \\
 \cmidrule(lr){2-3}
 \cmidrule(lr){4-5}
 Damage & avg. $\psi$ & max. $\psi$ & avg.  $\psi$ & max. $\psi$ \\
% &\% savings & \% savings & \% savings & \% savings \\
 \midrule
90\%	&	9.3	&	31.6	(1.25)	&	5.6	&	20.7	(1.25) \\
80\%	&	9.6	&	31.3	(1.2)	&	5.5	&	20.2	(1.25) \\
70\%	&	10.5	&	28.0	(1.2)	&	4.5	&	14.6	(1.25) \\
60\%	&	10.3	&	29.7	(1.2)	&	4.6	&	28.0	(1.25) \\
50\%	&	17.2	&	33.2	(1.2)	&	11.8	&	29.4	(1.25) \\
 \bottomrule
\end{tabular}
% }
\label{tab:facts}
\end{table}

\begin{figure}[h]
   \centering
   \includegraphics[scale=0.459]{facts_15_45_p1}
   	\vspace*{-0.3in}\\
   \caption{Comparison of total upgrade costs and $\psi$ for optimal upgrade solutions with and without FACTS devices.} 
   \label{fig:facts}
\end{figure}
\vspace{-0.2cm}
In Figure \ref{fig:facts}, it is interesting to note the (sometimes) non-increasing trend in $\psi$ as load scale increases.
%another interesting trend which is noteworthy to mention is the non-monotonic behavior in the values of \% savings as the load scale increases. 
Intuitively, this indicates that upgrades required for certain load scales are also sufficient for further increases in load. More importantly, perhaps, the OGRDT solution with and without FACTS devices is sometimes identical. This phenomena indicates that at certain load scales, FACTS devices cannot eliminate the need for other upgrades (line hardening, new lines, etc.) and once those upgrades are installed, the need for FACTS devices is eliminated. Of course, as the load scale increases, FACTS devices are once again required.  

%This can be explained based on the total upgrade costs curve (QC objective without FACTS) vs. load scales. A certain consequent set of load scales do not contain any upgrade or has  minimal upgrade due to the fact that the upgrade solutions built up to the previous load scale was sufficient to maintain a feasible network. The upgrade cost increases only after the load scale value crosses a specific threshold. This is also evident in the solutions for FACTS, where \% savings drop to almost zero values because of the constant upgrade costs curve. From the grid operator's perspective, we consider this as a very encouraging result as the upgrade solutions obtained at normally operable load scales (0.95,1.0) also guarantee to be feasible and resilient (including the FACTS' benefits) even at next few discretizations of larger load scales.   


\noindent
\textbf{PST devices} Generally speaking, we only found PST devices to be of benefit when the network is congested and has larger load scales.
%Advantages of having phase shifting transformers as an additional resiliency option could be seen only in congested networks at larger load scales. 
Under the 90\% damage cases, $\theta^u=15\degree$ and load scales greater than 1.20, we found solutions with savings of around 10\%. Though the improvement in $\psi$ when PSTs are included is not as impressive as when FACTS devices are included, we believe that the savings are non trival when considering large-scale grids.





orgdt_tps.tex
\documentclass{IEEEtran}


% *** GRAPHICS RELATED PACKAGES ***
%
\ifCLASSINFOpdf
  % \usepackage[pdftex]{graphicx}
  % declare the path(s) where your graphic files are
  % \graphicspath{{../pdf/}{../jpeg/}}
  % and their extensions so you won't have to specify these with
  % every instance of \includegraphics
  % \DeclareGraphicsExtensions{.pdf,.jpeg,.png}
\else
  % or other class option (dvipsone, dvipdf, if not using dvips). graphicx
  % will default to the driver specified in the system graphics.cfg if no
  % driver is specified.
  % \usepackage[dvips]{graphicx}
  % declare the path(s) where your graphic files are
  % \graphicspath{{../eps/}}
  % and their extensions so you won't have to specify these with
  % every instance of \includegraphics
  % \DeclareGraphicsExtensions{.eps}
\fi



% *** MATH PACKAGES ***
%
\usepackage{amsmath}
\usepackage{psfrag}
\usepackage{multicol}
\usepackage{epstopdf}
\usepackage{graphicx}
% \usepackage{hyperref}
\PassOptionsToPackage{bookmarks=false}{hyperref}
\usepackage{subfigure,gensymb,url}
\usepackage{color}
%\usepackage{algorithmic}
%\usepackage{algorithm}
\usepackage{anysize, booktabs, setspace}
\marginsize{0.8in}{0.8in}{0.8in}{0.8in}
\usepackage{amsfonts, amssymb, amsthm}
\usepackage{cleveref, nomencl, blindtext, mathtools, etoolbox}
\crefname{section}{§}{§§}   
\Crefname{section}{§}{§§}
\usepackage{pgfplots,tikz}

% \hypersetup{bookmarks=false}

% *** ALGORITHMS PACKAGE ***
% \usepackage{algorithm,algorithmic}
\usepackage[ruled,vlined,linesnumbered]{algorithm2e}
% \usepackage{algorithmic}
% \usepackage[algo2e,vlined,ruled,linesnumbered]{algorithm2e}


% nomenclature
\newcommand\Nomenclature[3][X]{\nomenclature[#1#3]{#2}{#3}}
\renewcommand\nomgroup[1]{%
  \item[\bfseries
  \ifstrequal{#1}{A}{Parameters and notation}{%
  \ifstrequal{#1}{B}{Binary variables}{%
  \ifstrequal{#1}{C}{Continuous variables}{%
  \ifstrequal{#1}{D}{Parameters}{%
  \ifstrequal{#1}{E}{Random variables}{%
  \ifstrequal{#1}{F}{Scenario parameters}{%
   \ifstrequal{#1}{X}{Other Symbols}{}}}}}}}]%
  }

\immediate\write18{makeindex \jobname.nlo -s nomencl.ist -o \jobname.nls}
\makenomenclature
\setlength{\nomitemsep}{0.25em} % no spacing in the nomenclature
\setlength\nomlabelwidth{1.35cm}

\newcommand{\ORGDT}{{\mathcal{P}_0^{AC}}}
\newcommand{\ORGDTdc}{{\mathcal{P}_0^{DC}}}
% \newcommand{\ORGDT1}{\mathcal{P}_1}
\newcommand{\ORGDTdev}{\mathcal{P}_1}

\newcommand{\SubProblem}{{\cal{Q}^{AC}}(s)}
\newcommand{\SubProblemDC}{{\cal{Q}^{DC}}(s)}
\newcommand{\Graph}{{\cal{G}}}
\newcommand{\Nodes}{{\cal{N}}}
\newcommand{\Edges}{{\cal{E}}}
\newcommand{\Disasters}{{\cal{S}}}
\newcommand{\Damages}{{\cal{D}}}
\newcommand{\HardenedDamages}{{\cal{D}^\prime}}
\newcommand{\Sim}{{\cal{S}}}
\newcommand{\Var}{{\cal{X}}}
\newcommand{\From}{{\cal{F}}}
\newcommand{\To}{{\cal{T}}}
\newcommand{\VarE}{\Edges(\Var)}
\newcommand{\VarB}{\Nodes(\Var)}
\newcommand{\variable}{{y}}
\newcommand{\flow}{{\digamma}}
\newcommand{\suscept}{{B}}
\newcommand{\conduct}{{G}}
\newcommand{\resist}{{R}}
\newcommand{\react}{{X}}
\newcommand{\thetabound}{{{\theta}^u}}
\newcommand{\voltagelbnd}{{{v}^l}}
\newcommand{\voltageubnd}{{{v}^u}}


% FIRST STAGE VARIABLES
\newcommand{\lineVariable}{x}
\newcommand{\switchVariable}{\tau}
\newcommand{\hardenVariable}{t}
\newcommand{\microgridVariablereal}{zp}
\newcommand{\microgridVariablereac}{zq}
\newcommand{\facilityVariable}{u}
\newcommand{\ChanceVariable}{v}
\newcommand{\compensator}{\delta}
\newcommand{\pst}{\gamma}


% SECOND STAGE VARIABLES
\newcommand{\lineuseVariable}{x^s}
\newcommand{\switchuseVariable}{\tau^s}
\newcommand{\lineCycleVariable}{\bar{x}^s}
\newcommand{\switchCycleVariable}{\bar{\tau}^s}
\newcommand{\linehardenVariable}{t^s}
\newcommand{\microgriduseVariablereal}{{zp^s}}
\newcommand{\microgriduseVariablereac}{{zq^s}}
\newcommand{\facilityuseVariable}{u^s}
\newcommand{\generatorVariablereal}{{gp^s}}
\newcommand{\loadVariablereal}{{lp^s}}
\newcommand{\generatorVariablereac}{{gq^s}}
\newcommand{\loadVariablereac}{{lq^s}}
\newcommand{\loadServeVariablereal}{{yp^s}}
\newcommand{\loadServeVariablereac}{{yq^s}}
\newcommand{\realflowVariable}{{p^s}}
\newcommand{\reacflowVariable}{{q^s}}
\newcommand{\thetaVariable}{{\theta^s}}
\newcommand{\voltageVariable}{{v^s}}
\newcommand{\linearXv}{{\widehat{xv}^s}}
\newcommand{\linearxVVcos}{{\widehat{xwc}^s}}
\newcommand{\linearxVVsin}{{\widehat{xws}^s}}
\newcommand{\linearVVcos}{{\widehat{wc}^s}}
\newcommand{\linearVVsin}{{\widehat{ws}^s}}
\newcommand{\loss}{{l^s}}
\newcommand{\compensatoruse}{{\delta^s}}
\newcommand{\pstuse}{{\gamma^s}}


\newcommand{\lineCost}{c}
\newcommand{\switchCost}{\kappa}
\newcommand{\hardenCost}{\psi}
\newcommand{\facilityCost}{\alpha}
\newcommand{\microgridCostreal}{\zeta p}
\newcommand{\microgridCostreac}{\zeta q}
\newcommand{\scenario}{\omega}
\newcommand{\compensatorCost}{{\Delta}}
\newcommand{\pstCost}{{\Gamma}}


\newcommand{\capacity}{T}
\newcommand{\phases}{{\cal{P}}}
\newcommand{\phaseVariation}{\beta}
\newcommand{\demandreal}{dp}
\newcommand{\demandreac}{dq}
\newcommand{\availablecapacityreal}{{gp^u}}
\newcommand{\availablecapacityreac}{{gq^u}}
\newcommand{\microgridcapacityreal}{{zp^u}}
\newcommand{\microgridcapacityreac}{{zq^u}}
%\newcommand{\cycles}{\cal{C}}
\newcommand{\criticalLoadreal}{{lp_{crit}}}
\newcommand{\totalLoadreal}{{lp_{ncrit}}}
\newcommand{\criticalLoadreac}{{lq_{crit}}}
\newcommand{\totalLoadreac}{{lq_{ncrit}}}
\newcommand{\Load}{{\cal{L}}}





% % Shortcuts for formulation
% % Sets
% \newcommand{\buses}{\mathcal{B}}
% \newcommand{\lines}{\mathcal{L}}
% \newcommand{\generators}{\mathcal{G}}
% \newcommand{\gb}{\ensuremath{G_b}}
% \newcommand{\startupblocks}{\mathcal{J}}
% \newcommand{\farms}{\mathcal{W}}
% \newcommand{\hours}{\mathcal{T}}
% \newcommand{\scenarios}{\mathcal S}

% % Binary Variables
% \renewcommand{\xi}[1]{\ensuremath{x_i(#1)}}
% \newcommand{\yi}[1]{\ensuremath{y_i(#1)}}
% \newcommand{\zi}[1]{\ensuremath{z_i(#1)}}
% \newcommand{\wij}[1]{\ensuremath{w_{ij}(#1)}}


% % Continuous Variables
% \newcommand{\phase}[1]{\ensuremath{\theta_{#1}(t)}}
% \newcommand{\droop}[1]{\ensuremath{\alpha_{#1}(t)}}
% \newcommand{\flow}{\ensuremath{f_{\ell}(t)}}

% % Parameters
% \newcommand{\sus}{\ensuremath{\beta_{\ell}}}
% \newcommand{\pmin}{\ensuremath{p_i^{\min}}}
% \newcommand{\pmax}{\ensuremath{p_i^{\max}}}
% \newcommand{\flowmax}{\ensuremath{f_{\ell}^{\max}}}
% \newcommand{\Lup}{\ensuremath{\overline L_i}}
% \newcommand{\Ldown}{\ensuremath{\underline L_i}}
% \newcommand{\uptime}{UT_i}
% \newcommand{\downtime}{DT_i}
% \newcommand{\onoff}{p_i^{\operatorname{on-off}}}
% \newcommand{\upinit}{p_i^{\operatorname{up,init}}}
% \newcommand{\downinit}{p_i^{\operatorname{down,init}}}
% \newcommand{\opinit}{p_i(0)}
% \newcommand{\rampup}{RU_i}
% \newcommand{\rampdown}{RD_i}
% \newcommand{\sccost}{\ensuremath{c_{ij}^{sc}}}

% % Random variables
% \newcommand{\rvw}{\ensuremath{\boldsymbol\omega_b(t)}}
% \newcommand{\rvwvector}{\ensuremath{\boldsymbol\omega(t)}}
% \newcommand{\rvomega}{\ensuremath{\boldsymbol\Omega(t)}}
% \newcommand{\rvp}{\ensuremath{\boldsymbol p_i(t)}}
% \newcommand{\rvtheta}[1]{\ensuremath{\boldsymbol \theta_{#1}(t)}}
% \newcommand{\rvthetavector}{\ensuremath{\boldsymbol \theta(t)}}
% \newcommand{\rvf}{\ensuremath{\boldsymbol f_{\ell}(t)}}
% \newcommand{\rvfs}{\ensuremath{\boldsymbol f^k_{\ell}(t)}}

% Line scenario variables

% \renewcommand{\exp}{\mathbb{E}}
% \newcommand{\defeq}{\vcentcolon=}
% \newcommand{\var}{\operatorname{var}}
% \newcommand{\prob}{\operatorname{Pr}}


\allowdisplaybreaks
\newcommand{\deriv}[3][]{% \deriv[<order>]{<func>}{<var>}
  \ensuremath{\frac{\partial^{#1} {#2}}{\partial {#3}^{#1}}}}
  
  
\newtheorem{thm}{Theorem}[section]
\newtheorem{cor}[thm]{Corollary}
\newtheorem{lem}[thm]{Lemma}
\newtheorem{conj}[thm]{Conjecture}
\newtheorem{clm}[thm]{Claim}
\newtheorem{rmk}[thm]{Remark}
\newtheorem{obs}[thm]{Observation}

\begin{document}
%\title{Optimal Resilient Grid Design for Transmission Systems}
% \title{Resilient Transmission Grid Design: AC vs. DC power flow models}
\title{Resilient Transmission Grid Design: AC Relaxation vs. DC approximation}


%% To specify the authors when (number of affiliations <= 2)
\author{
   %\footnotesize
   %\small
%   ``The gang"
   Harsha Nagarajan$^{\dag}$, Russell Bent$^{\dag}$, Pascal Van Hentenryck$^\ddag$, Scott Backhaus$^{\dag}$, Emre Yamangil$^\dag$ \\
$^\dag$ Center for Nonlinear Studies, Los Alamos National Laboratory, NM, United States. \\
$^\ddag$ Department of Industrial Operations Engineering, University of Michigan, Ann Arbor, United States.\\
Contact: harsha@lanl.gov
}


% make the title area
\maketitle

\begin{abstract}

As illustrated in recent years (Superstorm Sandy, the Northeast Ice Storm of 1998, etc.), extreme weather events pose an enormous threat to the electric power transmission systems and the associated socio-economic systems that depend on reliable delivery of electric power. Besides inevitable malfunction of power grid components, deliberate malicious attacks can cause high risks to the service. These threats motivate the need for approaches and methods that improve the resilience of power systems. In this paper, we develop a model and tractable methods for optimizing the upgrade of transmission systems through a combination of hardening existing components, adding redundant lines, switches, generators, and FACTS and phase-shifting devices. While many of these controllable components are included in traditional design (expansion planning) problems, we uniquely assess their benefits from a resiliency point of view. 
More importantly, perhaps, we evaluate the suitability of using state-of-the-art AC power flow relaxations versus the common DC approximation in resilience improvement studies. The resiliency model and algorithms are tested on a modified version of the RTS-96 (single area) system. 


\end{abstract}


 \input{Intro}

\input{Formulation}

 \input{Algo}

 \input{NR}

 \input{Conclusions}




%==============================;
%  Include all the references  ;
%==============================;
\bibliographystyle{IEEEtran}
\bibliography{references.bib}


% that's all folks
\end{document}



results.tex
\documentclass{article}
\usepackage[cmex10]{amsmath}
\usepackage{psfrag}
\usepackage{multicol}
\usepackage{epstopdf}
\usepackage{graphicx,subfigure}
\usepackage{hyperref}
\usepackage{subfigure,gensymb,url}
\usepackage{color}
%\usepackage{algorithmic}
%\usepackage{algorithm}
\usepackage{anysize, booktabs, setspace}
\marginsize{1in}{1in}{0.5in}{0.5in}
\usepackage{amsfonts, amssymb, amsthm}
\usepackage{cleveref, nomencl, blindtext, mathtools, etoolbox}
\crefname{section}{§}{§§}   
\Crefname{section}{§}{§§}
\usepackage{pgfplots,tikz}


% *** ALGORITHMS PACKAGE ***
\usepackage[ruled,vlined,linesnumbered]{algorithm2e}

\title{Updated results for ORGDT}
\date{\today}

\begin{document}

\maketitle

\section{Upgrade cost analysis - without FACTS devices}
\textbf{Note}: The base network is feasible only at load factor = 1.0

%---------------------------------------------------------------;
%   Figure : Comparison of upgrade costs up to loadfactor=1.25  ;
%---------------------------------------------------------------;
\begin{figure}[htp]
	\centering
	\subfigure[$p=0.1$, $\theta^u=15^o$]{
	\includegraphics[scale=0.6]{Figures/test_p1_ang15_n.pdf}}
	%
	\subfigure[$p=0.1$, $\theta^u=45^o$]{
	\includegraphics[scale=0.6]{Figures/test_p1_ang45_n.pdf}}
	%
	\subfigure[$p=0.1$, Load factor = 1.30]{
	\begin{tikzpicture}
    \begin{axis}[
    scale=0.57,
    ybar,
    enlarge x limits=0.25,
    legend style={at={(0.6,1.4)},draw=none,
      anchor=north},
    % ylabel={Upgrade cost (\$)},
    symbolic x coords={$\theta^u=15^o$,$\theta^u=45^o$},
    xtick=data,
    ymajorgrids,
    %nodes near coords,
    %nodes near coords align={vertical},
    ]
    \addplot coordinates {($\theta^u=15^o$,3.8087E7) ($\theta^u=45^o$,3.7764E7)};
    \addplot coordinates {($\theta^u=15^o$,1.075E7) ($\theta^u=45^o$,1.031E7)};
    \legend{QC relax.,DC approx.}
    \end{axis}
    \end{tikzpicture}
	}
	%
    \subfigure[$p=0.2$, $\theta^u=15^o$]{
	\includegraphics[scale=0.6]{Figures/test_p2_ang15_n.pdf}}
	%
	\subfigure[$p=0.2$, $\theta^u=45^o$]{
	\includegraphics[scale=0.6]{Figures/test_p2_ang45_n.pdf}}
	%
	%
	\subfigure[$p=0.2$, Load factor = 1.30]{
	\begin{tikzpicture}
    \begin{axis}[
    scale=0.57,
    ybar,
    enlarge x limits=0.25,
    legend style={at={(0.6,1.4)},draw=none,
      anchor=north},
    % ylabel={Upgrade cost (\$)},
    symbolic x coords={$\theta^u=15^o$,$\theta^u=45^o$},
    xtick=data,
    ymajorgrids,
    %nodes near coords,
    %nodes near coords align={vertical},
    ]
    \addplot coordinates {($\theta^u=15^o$,3.7749E7) ($\theta^u=45^o$,3.7447E7)};
    \addplot coordinates {($\theta^u=15^o$,1.067E7) ($\theta^u=45^o$,1.038E7)};
    \legend{QC relax.,DC approx.}
    \end{axis}
    \end{tikzpicture}
	}
	%
		%
	\subfigure[$p=0.3$, $\theta^u=15^o$]{
	\includegraphics[scale=0.6]{Figures/test_p3_ang15_n.pdf}}
	%
	\subfigure[$p=0.3$, $\theta^u=45^o$]{
	\includegraphics[scale=0.6]{Figures/test_p3_ang45_n.pdf}}
	%
	%
	\subfigure[$p=0.3$, Load factor = 1.30]{
	\begin{tikzpicture}
    \begin{axis}[
    scale=0.57,
    ybar,
    enlarge x limits=0.25,
    legend style={at={(0.6,1.4)},draw=none,
      anchor=north},
    % ylabel={Upgrade cost (\$)},
    symbolic x coords={$\theta^u=15^o$,$\theta^u=45^o$},
    xtick=data,
    ymajorgrids,
    %nodes near coords,
    %nodes near coords align={vertical},
    ]
    \addplot coordinates {($\theta^u=15^o$,3.74287E7) ($\theta^u=45^o$,3.7135E7)};
    \addplot coordinates {($\theta^u=15^o$,1.0453E7) ($\theta^u=45^o$,1.0169E7)};
    \legend{QC relax.,DC approx.}
    \end{axis}
    \end{tikzpicture}
	}
	%
% 	\caption{}
	\label{Fig:fig1}
\vspace*{-0.1in}
\end{figure}

\begin{figure}[htp]
	\centering
	%
	\subfigure[$p=0.4$, $\theta^u=15^o$]{
	\includegraphics[scale=0.6]{Figures/test_p4_ang15_n.pdf}}
	%
	\subfigure[$p=0.4$, $\theta^u=45^o$]{
	\includegraphics[scale=0.6]{Figures/test_p4_ang45_n.pdf}}
	%
	%
	\subfigure[$p=0.4$, Load factor = 1.30]{
	\begin{tikzpicture}
    \begin{axis}[
    scale=0.57,
    ybar,
    enlarge x limits=0.25,
    legend style={at={(0.6,1.4)},draw=none,
      anchor=north},
    % ylabel={Upgrade cost (\$)},
    symbolic x coords={$\theta^u=15^o$,$\theta^u=45^o$},
    xtick=data,
    ymajorgrids,
    %nodes near coords,
    %nodes near coords align={vertical},
    ]
    \addplot coordinates {($\theta^u=15^o$,3.7311E7) ($\theta^u=45^o$,3.7023E7)};
    \addplot coordinates {($\theta^u=15^o$,1.044E7) ($\theta^u=45^o$,1.0169E7)};
    \legend{QC relax.,DC approx.}
    \end{axis}
    \end{tikzpicture}
	}
	%
	\subfigure[$p=0.5$, $\theta^u=15^o$]{
	\includegraphics[scale=0.6]{Figures/test_p5_ang15_n.pdf}}
	%
	\subfigure[$p=0.5$, $\theta^u=45^o$]{
	\includegraphics[scale=0.6]{Figures/test_p5_ang45_n.pdf}}
	%
	%
	\subfigure[$p=0.5$, Load factor = 1.30]{
	\begin{tikzpicture}
    \begin{axis}[
    scale=0.57,
    ybar,
    enlarge x limits=0.25,
    legend style={at={(0.6,1.4)},draw=none,
      anchor=north},
    % ylabel={Upgrade cost (\$)},
    symbolic x coords={$\theta^u=15^o$,$\theta^u=45^o$},
    xtick=data,
    ymajorgrids,
    %nodes near coords,
    %nodes near coords align={vertical},
    ]
    \addplot coordinates {($\theta^u=15^o$,3.508E7) ($\theta^u=45^o$,3.6798E7)};
    \addplot coordinates {($\theta^u=15^o$,1.0252E7) ($\theta^u=45^o$,9.9419E6)};
    \legend{QC relax.,DC approx.}
    \end{axis}
    \end{tikzpicture}
	}
	%
% 	\caption{}
	\label{Fig:fig1}
\vspace*{-0.1in}
\end{figure}

% Table ----- 
% \begin{table}[h]
%     \centering
%     \scriptsize
%     \caption{Total real and reactive power installed at damage probability = 0.1 (maximum damage case).}
%     \begin{tabular}{cllll}
%          \toprule
%          Load factor & $\theta^u=15^o$ & & $\theta^u=45^o$ &\\
%          &mg-real&mg-reac. & mg-real & mg-reac. \\
%          \cmidrule{1-5}
%          1.00 &0 &310 &0 &310\\
%          1.05 &0  & 0 & 0 & 310\\ 
%          1.10 & 0 & 0 & 0 & 0\\
%          1.15 & 0 & 310 & 0 & 310\\
%          1.20 & 0  & 0 & 0 & 310\\
%          1.25 & 0 & 930 & 0 & 310\\
%          1.30 & 40.8 & 15300 & 40.8 & 12400\\
%          \bottomrule
%     \end{tabular}
%     \label{tab:my_label}
% \end{table}


\section{AC-feasibility analysis}
For given QC relaxations, we fix $x^*, \tau^*, t^*, u^*, zp^*,zq^*$. In case of DC approximations, we fix $x^*, \tau^*, t^*, u^*, zp^*$ and based on $zp^*$,we fix the bounds $-zp^*/2 <= zq <= zp^*/2$. Based on these build variables, we solve the following problem which quantifies the minimum amount of critical and non-critical loads which have to be shed to ensure that the AC power flow equations are feasible:
\begin{equation}
\begin{aligned}
& \text{minimize} & & \left(
\lambda^p_{cr} + \lambda^q_{cr} + \lambda^p_{ncr} + \lambda^q_{ncr}\right) \\
&\text{subject to} & & \mathrm{AC \ power \ flow \ (nonconvex) \ constraints}, \\
& & & \sum_{i\in \mathcal{L}} lp_i^s \geq lp_{crit} \left(\sum_{i \in \mathcal{L}} dp_i\right) - \lambda^p_{cr}, \\
& & & \sum_{i\in \mathcal{L}} lq_i^s \geq lq_{crit} \left(\sum_{i \in \mathcal{L}} dq_i\right) - \lambda^q_{cr}, \\
& & & \sum_{i\in \mathcal{N}\setminus \mathcal{L}} lp_i^s \geq lp_{ncrit} \left(\sum_{i \in \mathcal{N}\setminus\mathcal{L}} dp_i\right) - \lambda^p_{ncr}, \\
& & & \sum_{i\in \mathcal{N}\setminus \mathcal{L}} lq_i^s \geq lq_{ncrit} \left(\sum_{i \in \mathcal{N}\setminus\mathcal{L}} dq_i\right) - \lambda^q_{ncr}, \\
& & & \lambda^p_{cr}, \lambda^q_{cr}, \lambda^p_{ncr},  \lambda^q_{ncr} \geq 0.
\end{aligned}
\end{equation}
In the following figures, for a given scenario $s$, percent apparent-load shed is given by 
$$\frac{\sqrt{(\lambda^p_{cr})^2 + (\lambda^q_{cr})^2}}{\sqrt{\left(\sum_{i\in \mathcal{L}}dp_i\right)^2+\left(\sum_{i\in \mathcal{L}}dq_i\right)^2}}*100, \ \frac{\sqrt{(\lambda^p_{ncr})^2 + (\lambda^q_{ncr})^2}}{\sqrt{\left(\sum_{i\in \mathcal{N}\setminus\mathcal{L}}dp_i\right)^2+\left(\sum_{i\in \mathcal{N}\setminus\mathcal{L}}dq_i\right)^2}}*100$$

\subsection{AC-feasibility at load scale = 1.0}
\vspace{-0.5cm}
\begin{figure}[htp]
	\centering
	\subfigure[$p=0.1$, $\theta^u=15^o$]{
	\includegraphics[scale=0.8]{Figures/test_ang15_p1.pdf}}
	%
	\subfigure[$p=0.1$, $\theta^u=45^o$]{
	\includegraphics[scale=0.8]{Figures/test_ang45_p1.pdf}}
	%
	\subfigure[$p=0.2$, $\theta^u=15^o$]{
	\includegraphics[scale=0.8]{Figures/test_ang15_p2.pdf}}
	%
	\subfigure[$p=0.2$, $\theta^u=45^o$]{
	\includegraphics[scale=0.8]{Figures/test_ang45_p2.pdf}}
	%
	\subfigure[$p=0.3$, $\theta^u=15^o$]{
	\includegraphics[scale=0.8]{Figures/test_ang15_p3.pdf}}
	%
	\subfigure[$p=0.3$, $\theta^u=45^o$]{
	\includegraphics[scale=0.8]{Figures/test_ang45_p3.pdf}}
	%
	\subfigure[$p=0.4$, $\theta^u=15^o$]{
	\includegraphics[scale=0.8]{Figures/test_ang15_p4.pdf}}
	%
	\subfigure[$p=0.4$, $\theta^u=45^o$]{
	\includegraphics[scale=0.8]{Figures/test_ang45_p4.pdf}}
\end{figure}

\begin{figure}[htp]
	\centering
	\subfigure[$p=0.5$, $\theta^u=15^o$]{
	\includegraphics[scale=0.8]{Figures/test_ang15_p5.pdf}}
	%
	\subfigure[$p=0.5$, $\theta^u=45^o$]{
	\includegraphics[scale=0.8]{Figures/test_ang45_p5.pdf}}
	\label{Fig:fig_feas1}
\end{figure}


% \subsection{AC-feasibility at load scale = 1.3}






% \section{Base network feasibility analysis}
% \begin{table}[!h]
%     \centering
%     \scriptsize
%     \caption{Total upgrade cost (\$) to achieve feasibility of base network at damage probability = 0.1 (maximum damage case).}
%     \begin{tabular}{cll}
%          \toprule
%          Load factor & $\theta^u=15^o$ & $\theta^u=15^o$ \\
%           & (all loads=100\%) & (critical=100\%, non-critical=80\%)\\
%          \cmidrule{1-3}
%          1.00 &0 &0 \\
%          1.05 &3000  & 0\\  
%           & {\color{gray}(3 sw)} & \\
%          1.10 &100000 & 0\\
%          & {\color{gray}(1 mg - reac. support)} & \\
%          1.15 &1.217E7 & 0\\
%           &{\color{gray}(2 sw, 8 mg - huge reac. support)} & \\
%          1.20 & 6.497E7 & 0\\
%           &{\color{gray}(5 sw, 7 mg - huge reac. support)}& \\
%          1.25 & 1.55E8 & 0 \\
%           & {\color{gray}(2 sw, 6 mg - real+reac. support)}& \\
%          1.30 & 2.69E8 & 3.37E7\\
%          &{\color{gray}(8 mg - real+huge reac. support)} & {\color{gray}(1 sw, 4 mg- real+huge reac. support)} \\
%          \bottomrule
%     \end{tabular}
%     \label{tab:my_label}
% \end{table}




\newpage
\section{Analysis with FACTS device (var compensator)}
Following results are at $p=0.1$ and $\theta^u=15^o$.  \\
$$\% \ \mathrm{savings} = \left(\frac{\mathrm{Upgrade\ without\ devices} - \mathrm{Upgrade\ with\ devices}}{\mathrm{Upgrade\ without\ devices}}\right)*100$$

\begin{figure}[htp]
\centering
%
\subfigure[Reactance reduction factor = 0.5]{
\begin{tikzpicture}

    \begin{axis}[
    scale=0.88,
    ybar,
    axis y line*=right,
    % enlarge x limits=0.25,
    % legend style={at={(1.2,0.6)},draw=none,anchor=north},
    % ylabel={\% savings},
    xlabel={load scale},
    symbolic x coords={1.0, 1.05, 1.1, 1.15, 1.2, 1.25},
    xtick=data,
    ymajorgrids,
    ]
    \addplot coordinates {(1.0,6.07)(1.05,0) (1.1,1.19) (1.15,19.89) (1.2,18.58) (1.25,31.59)};
    \addplot coordinates {(1.0,0)(1.05,0) (1.1,1.02) (1.15,20.96) (1.2,16.64) (1.25,22.78)};
    % \legend{QC relax.,DC approx.}
    \end{axis}
    %
    \begin{axis}[
    scale=0.88,
    axis y line*=left,
    ylabel=Total upgrade cost,
    axis x line=none,
    ymin=0,
    xtick=data,
    % ymajorgrids,
    legend style={at={(0.8,1.42)},draw=none,anchor=north},
    ]
    \addplot[mark=square*,color=blue] coordinates {(1.0,743098)(1.05,743098) (1.1,802627) (1.15,990048) (1.2,1.14E+06) (1.25,1.42E+06)};
    \addplot[mark=*,color=red] coordinates {(1.0,480241)(1.05,726112) (1.1,802627) (1.15,981325) (1.2,9.91E+05) (1.25,1.26E+06)};
    \addplot[mark=square*,dashed,color=blue] coordinates {(1.0,698028)(1.05,743098) (1.1,793098) (1.15,793098) (1.2,931977) (1.25,972769)};
    \addplot[mark=*,dashed,color=red] coordinates {(1.0,480241)(1.05,726112) (1.1,775112) (1.15,775612) (1.2,826112) (1.25,971769)};
    \legend{QC (w/o FACTS),DC (w/o FACTS), QC (with FACTS),DC (with FACTS) }
    \end{axis}
\end{tikzpicture}}
%
\subfigure[Reactance reduction factor = 0.8]{
\begin{tikzpicture}
    \begin{axis}[
    scale=0.88,
    axis y line*=right,
    ybar,
    % enlarge x limits=0.25,
    legend style={at={(0.2,1.29)},draw=none,anchor=north},
    % ylabel={\% savings},
    xlabel={load scale},
    symbolic x coords={1.0, 1.05, 1.1, 1.15, 1.2, 1.25},
    xtick=data,
    ymajorgrids,
    ]
    \addplot coordinates {(1.0,0)(1.05,0) (1.1,0) (1.15,1.85) (1.2,10.7) (1.25,9.67)};
    \addplot coordinates {(1.0,0)(1.05,0) (1.1,0.9) (1.15,2.35) (1.2,0.05) (1.25,3.09)};
    \legend{QC (\% savings),DC (\% savings)}
    \end{axis}
    %
    \begin{axis}[
    scale=0.88,
    axis y line*=left,
    % ylabel=Total upgrade cost,
    axis x line=none,
    ymin=0,
    xtick=data,
    % ymajorgrids,
    legend style={at={(0.8,1.4)},draw=none,anchor=north},
    ]
    \addplot[mark=square*,color=blue] coordinates {(1.0,743098)(1.05,743098) (1.1,802627) (1.15,990048) (1.2,1.14E+06) (1.25,1.42E+06)};
    \addplot[mark=*,color=red] coordinates {(1.0,480241)(1.05,726112) (1.1,802627) (1.15,981325) (1.2,9.91E+05) (1.25,1.26E+06)};
    \addplot[mark=square*,dashed,color=blue] coordinates {(1.0,743098)(1.05,743098) (1.1,802627) (1.15,971769) (1.2,1.02E+06) (1.25,1.28E+06)};
    \addplot[mark=*,dashed,color=red] coordinates {(1.0,480241)(1.05,726112) (1.1,776112) (1.15,958312) (1.2,990548) (1.25,1.22E+06)};
    \end{axis}
\end{tikzpicture}}
% \caption{test}
\end{figure}





\end{document}

orgdt_tps.bbl
% Generated by IEEEtran.bst, version: 1.13 (2008/09/30)
\begin{thebibliography}{10}
\providecommand{\url}[1]{#1}
\csname url@samestyle\endcsname
\providecommand{\newblock}{\relax}
\providecommand{\bibinfo}[2]{#2}
\providecommand{\BIBentrySTDinterwordspacing}{\spaceskip=0pt\relax}
\providecommand{\BIBentryALTinterwordstretchfactor}{4}
\providecommand{\BIBentryALTinterwordspacing}{\spaceskip=\fontdimen2\font plus
\BIBentryALTinterwordstretchfactor\fontdimen3\font minus
  \fontdimen4\font\relax}
\providecommand{\BIBforeignlanguage}[2]{{%
\expandafter\ifx\csname l@#1\endcsname\relax
\typeout{** WARNING: IEEEtran.bst: No hyphenation pattern has been}%
\typeout{** loaded for the language `#1'. Using the pattern for}%
\typeout{** the default language instead.}%
\else
\language=\csname l@#1\endcsname
\fi
#2}}
\providecommand{\BIBdecl}{\relax}
\BIBdecl

\bibitem{Brostrom2007}
E.~Brostrom, ``{Ice Storm Modelling in Transmission System Reliability
  Calculations},'' Ph.D. dissertation, Royal Institute of Technology, 2007.

\bibitem{Eidinger}
J.~Eidinger and L.~Kempner, ``{Reliability of Transmission Towers under Extreme
  Wind and Ice Loading},'' \emph{G\&E Engineering Systems Inc. and Bonneville
  Power Administration}, 2013.

\bibitem{wang2010influence}
S.~Wang, ``Influence of ice accretions on mechanical and electrical performance
  of overhead transmission lines,'' \emph{Journal of Electrical and Mechanical
  Engineering}, vol.~1, no.~1, 2010.

\bibitem{Mills2012}
E.~Mills, ``{Electric Grid Disruptions and Extreme Weather},'' in \emph{US
  Disaster Reanalysis Workshop National Climatic Data Center}, Asheville, NC,
  2012.

\bibitem{7540988}
H.~Nagarajan, E.~Yamangil, R.~Bent, P.~V. Hentenryck, and S.~Backhaus,
  ``Optimal resilient transmission grid design,'' in \emph{2016 Power Systems
  Computation Conference (PSCC)}, June 2016, pp. 1--7.

\bibitem{hijazi2014convex}
H.~Hijazi, C.~Coffrin, and P.~Van~Hentenryck, ``Convex quadratic relaxations
  for mixed-integer nonlinear programs in power systems,'' \emph{Mathematical
  Programming Computation}, pp. 1--47, 2014.

\bibitem{yamangil2015resilient}
E.~Yamangil, R.~Bent, and S.~Backhaus, ``Resilient upgrade of electrical
  distribution grids,'' in \emph{Twenty-Ninth AAAI Conference on Artificial
  Intelligence}, 2015.

\bibitem{Chen:13}
R.~L.-Y. Chen and C.~A. Phillips, ``k-edge failure resilient network design,''
  \emph{Electronic Notes in Discrete Mathematics}, vol.~41, no.~0, pp. 375 --
  382, 2013.

\bibitem{Chen2014a}
R.~L.-Y. Chen, A.~Cohn, N.~Fan, and A.~Pinar, ``{Contingency-Risk Informed
  Power System Design},'' \emph{IEEE Transactions on Power Systems}, vol.~29,
  no.~5, pp. 2087--2096, Sep. 2014.

\bibitem{Salmeron2009}
J.~Salmeron, K.~Wood, and R.~Baldick, ``{Worst-case interdiction analysis of
  large-scale electric power grids},'' \emph{IEEE Transactions on Power
  Systems}, vol.~24, no.~1, pp. 96--104, 2009.

\bibitem{Delgadillo2010}
A.~Delgadillo, J.~Arroyo, and N.~Alguacil, ``{Analysis of Electric Grid
  Interdiction with Line Switching},'' \emph{IEEE Transactions on Power
  Systems}, vol.~25, no.~2, pp. 633--641, 2010.

\bibitem{nezamoddini2017risk}
N.~Nezamoddini, S.~Mousavian, and M.~Erol-Kantarci, ``A risk optimization model
  for enhanced power grid resilience against physical attacks,'' \emph{Electric
  Power Systems Research}, vol. 143, pp. 329--338, 2017.

\bibitem{Krishnan2015}
V.~Krishnan, J.~Ho, B.~Hobbs, A.~Liu, J.~McCalley, M.~Shahidehpour, and
  Q.~Zheng, ``{Co-Optimization of Electricity Transmission and Generation
  Resources for Planning and Policy Analysis: Review of Concepts and Modeling
  Approaches},'' \emph{Energy Systems}, 2015.

\bibitem{miasaki2012transmission}
C.~T. Miasaki, E.~Franco, and R.~A. Romero, ``Transmission network expansion
  planning considering phase-shifter transformers,'' \emph{Journal of
  Electrical and Computer Engineering}, vol. 2012, 2012.

\bibitem{Frolov2014}
V.~Frolov, S.~Backhaus, and M.~Chertkov, ``Efficient algorithm for locating and
  sizing series compensation devices in large power transmission grids: Ii.
  solutions and applications,'' \emph{New Journal of Physics}, vol.~16, no.~10,
  p. 105016, 2014.

\bibitem{baran1989optimal}
M.~Baran and F.~F. Wu, ``Optimal sizing of capacitors placed on a radial
  distribution system,'' \emph{IEEE Trans. on power Delivery}, vol.~4, no.~1,
  pp. 735--743, 1989.

\bibitem{Baldick2014}
R.~Baldick and R.~O'Neill, ``Estimates of comparative costs for uprating
  transmission capacity,'' \emph{Power Delivery, IEEE Trans. on}, vol.~24,
  no.~2, pp. 961--969, April 2009.

\bibitem{nagarajan2016tightening}
H.~Nagarajan, M.~Lu, E.~Yamangil, and R.~Bent, ``Tightening mccormick
  relaxations for nonlinear programs via dynamic multivariate partitioning,''
  in \emph{International Conference on Principles and Practice of Constraint
  Programming}.\hskip 1em plus 0.5em minus 0.4em\relax Springer, 2016, pp.
  369--387.

\bibitem{wright1997primal}
S.~J. Wright, \emph{Primal-dual interior-point methods}.\hskip 1em plus 0.5em
  minus 0.4em\relax SIAM, 1997.

\bibitem{wong1999ieee}
P.~Wong, P.~Albrecht, R.~Allan, R.~Billinton, Q.~Chen, C.~Fong, S.~Haddad,
  W.~Li, R.~Mukerji, D.~Patton \emph{et~al.}, ``The ieee reliability test
  system-1996. a report prepared by the reliability test system task force of
  the application of probability methods subcommittee,'' \emph{Power Systems,
  IEEE Trans. on}, vol.~14, pp. 1010--1020, 1999.

\bibitem{dunning2015jump}
I.~Dunning, J.~Huchette, and M.~Lubin, ``Jump: A modeling language for
  mathematical optimization,'' \emph{arXiv:1508.01982}, 2015.

\bibitem{Mak2014}
T.~Mak, C.~Coffrin, P.~Van~Hentenryck, I.~Hiskens, and D.~Hill, ``{Power System
  Restoration Planning with Standing Phase Angle and Voltage Difference
  Constraints},'' in \emph{Power Systems Computation Conference}, Wroclaw,
  Poland, 2014, pp. 1--8.

\end{thebibliography}
