Simple Certificate of Solvability of Power Flow Equations for Distribution Systems
\IEEEauthorblockN{Suhyoun Yu
Hung D. Nguyen
and Konstantin S. Turitsyn
%\boldmath Power flow solvable boundary plays an important role in contingency analysis, security assessment, and planning processes. However, to construct the real solvable boundary in multidimensional parameter space is burdensome and time consuming. In this paper, we develop a new technique to approximate the solvable boundary of distribution systems based on Banach fixed point theorem. Not only the new technique is fast and non-iterative, but also the approximated boundary is more valuable to system operators in the sense that it is closer to the feasible region. Moreover, a simple solvable criterion is also introduced that can serve as a security constraint in various planning and operational problems.
main.tex
\documentclass[conference]{IEEEtran}

\hyphenation{Power System Engineering}

\usepackage{pgfplots,pgfplotstable,booktabs}

\pgfplotsset{compat=1.5}

\usepackage{amsmath, bbold}

%\usepackage{subfig}

%\graphicspath{{images/}}

\usepackage{etoolbox}

\usepackage{graphicx, epstopdf}
%\usepackage{graphicx}
\usepackage{caption}
\usepackage{subfigure}

%\usepackage{subcaption}

%\epstopdfsetup{update}

\DeclareGraphicsExtensions{.eps}

\apptocmd{\sloppy}{\hbadness 10000\relax}{}{}

\usepackage{cite}

\newtheorem{theorem}{Theorem}
\newtheorem{conjecture}{Conjecture}
\newtheorem{proof}{Proof}

\usepackage{epsf}

\begin{document}
%

\title{Simple Certificate of Solvability of Power Flow Equations for Distribution Systems}

\author{
\IEEEauthorblockN{Suhyoun Yu, Hung D. Nguyen, and Konstantin S. Turitsyn}
\IEEEauthorblockA{Department of Mechanical Engineering\\
Massachusetts Institute of Technology\\
Cambridge, MA 02139\\
Email: syu2@mit.edu, hunghtd@mit.edu, and turitsyn@mit.edu}}

% make the title area
\maketitle

\begin{abstract}
%\boldmath
Power flow solvable boundary plays an important role in contingency analysis, security assessment, and planning processes. However, to construct the real solvable boundary in multidimensional parameter space is burdensome and time consuming. In this paper, we develop a new technique to approximate the solvable boundary of distribution systems based on Banach fixed point theorem. Not only the new technique is fast and non-iterative, but also the approximated boundary is more valuable to system operators in the sense that it is closer to the feasible region. Moreover, a simple solvable criterion is also introduced that can serve as a security constraint in various planning and operational problems.

\end{abstract}

\begin{IEEEkeywords}
Distribution systems, power flow, solvability, security constraint, voltage stability.
\end{IEEEkeywords}

\section{Introduction}

% {\color{blue}
% % @Hung (from essie) Line 86: are there supposed to be some equations? No. I think Sharelarex made the gap itself so that section II begins at the second page.\\
% @Prof. Turitsyn: From \eqref{eq:solcri}, if we introduce a new variable: $\sigma = \sum_k |s_k/s_{k,max}|$. Could $1-\sigma$ be used as a security indicator?\\
% % @Essie: Does $s_{k,max}$ in \eqref{eq:solcri} depend on the slack bus voltage $V_0$? or you assume that $V_0 = 1$.
% }

The penetration of renewable energy sources is expected to increase substantially in the near future. Besides numerous advantages of such more environmentally friendly resources, additional risks may make the distribution system more vulnerable to instability \cite{KundurDSA}. Voltage security is one of main existing issues that we need to address carefully and thoroughly in order to operate the systems in a reliable and secure fashion. Hence, the system operators also need to know the system limits which require faster and more reliable security assessment tools. In fact, static voltage stability criteria is widely used by system operators \cite{xie2007novel,LeXieQSVS}; moreover, it has been argued that static analysis is preferred over dynamic approach \cite{morison1993voltage}. Most of static analysis techniques rely on the solution of quasi-steady state analysis or power flow problem. The concept of quasi-steady state is based on the assumption of slow dynamics so that the power system can be considered in steady state.

% Power flow in a network is determined by the voltage at each bus of the network and the impedance of the lines between buses. The load flow problem consists of finding the set of voltages: magnitude and angle, which, together with the network impedance, produces the load flows that are known to be correct at the system terminals. Power flow study is a steady state analysis used to solve such unknown values in order to plan ahead and account for various hypothetical situations of networks \cite{}.

These power flow equations have been posing much difficulty to the engineers due to their inherent nonlinearity and high computational cost of finding solutions with conventional approaches \cite{HungTuritsyn}. Iterative methods such as Newton-Raphson or Gauss-Seidel methods are the most widely used by power flow solvers for the following reasons. First, the solution can be found after several iterations in most of the cases, and the iterative methods can be easily scaled to bulk systems. Second, the solutions obtained can be used for control purposes. However, these iterative methods may be subject to divergence due to either bad initial guesses, proximity to solvable boundaries, or the nonexistence of solutions \cite{JEET}.

%Two traditional iterative algorithms for solving these non-linear power flow equations are the Newton-Raphson method and Fast Decoupled methods; however, despite their popular usage, the former method requires an expensive LU factorization of power flow Jacobian at each Newton-iteration, and the latter may subject to divergence because of either bad initial guess, proximity to solvable boundaries, or the nonexistence of solutions.
%unless the power system at hand is well-behaved. 
In order to overcome such handicaps, numerous attempts have been made to design more efficient and accurate ways to solve or certify existence of solutions to those equations \cite{coffrin2014linear,bergen, Cvijic2012}. Among them, in \cite{bolognani2014existence}, Bolognani proposed a method of certifying the solvability of power flow equations using the Banach Fixed Point Theorem. The resulting criterion has a simple algebraic representation and dramatically reduces the computational of assessing system security. The downside of the approach is the conservativeness of the criterion and its inability to certify the solvability of the whole region. The motivation for this paper is to extend Bolognani’s criterion and alleviate the conservativeness observed in the original definition.

The key contributions of this paper are as follows. We analyze an extension of Bolognani's criterion and develop a new mathematical criterion for solvability certification. Then, we derive a simple algebraic representation and propose a new technique for fast and non-iterative characterization and visualization of the solvability boundaries in a multidimensional parameter space. Finally, two important applications for solvable boundary approximation and security constraint optimization are presented.

%such that the solution boundary may be manipulated as desired;
%; one, revisit Bolognani’s method of establishing a practical method for nonlinear power flow equations and propose a modification of its criterion such that the solution boundary may be manipulated as desired; two, characterize the result from this proposed modification made.
In section \ref{sec:model}, we describe the nonlinear power flow equations and the notations that be used throughout this paper, followed by a brief explanation of the criterion established by Bolognani which is the motivation for this paper. The method for extending Bolognani’s criterion, along with its proof, be presented in section \ref{sec:modify}. Then, section \ref{sec:simulation} displays the $2$-bus simulation results and related empirical observation that discuss the relationship between the resulting ellipsoids and rhombuses. This section also discusses the algorithm, or the convex hull method for constructing outermost rhombuses that can be observed from simulations. Section \ref{sec:outermost} provides with more simulations of $3$-bus and $13$-bus radial systems, which is followed by section \ref{sec:discussion} where main results are discussed. In section \ref{sec:apps}, we introduce two possible application of this modification.

\section{Power flow problem and solvable boundary approximation} \label{sec:model}

\subsection{Power Flow Equations}

In this work we use the notation introduced originally by Bolognani in \cite{bolognani2014existence}. %The considered networks are radial ones. 
Let $\mathcal{L}:=\{1,...,n\}$ be the index set denoting all load buses, and let index $0$ correspond to the slack bus. For $h\in\mathcal{L}$, let $v_h\in\mathbb{C}$ and $i_h\in\mathbb{C}$ be the voltage and current, respectively, that define each load bus, and let $v,i\in\mathbb{C}^{n+1}$ be complex vectors with entries of $v_h$ and $i_h$, respectively. Also, let $s_h\in\mathbb{C}$ be the imposed power at each load bus with $s\in\mathbb{C}^{n+1}$ as the vector with $s_h$ as entries. We assume that the voltage level at slack bus is fixed at $v_0 = V_0$.

%For the sake of simplicity, assume that the system is in steady state and that all voltages and currents are sinusoidal signals at the same frequencies. Then each signal can be characterized as a complex number of the form $y=|y|\exp(j\angle y)$; hence 
% For the slack bus, we have that:

% \begin{equation*}
% v_0=V_0\exp(j\theta_0)
% \end{equation*}
% where $V_0, \theta_0 \in \mathbb{R}$; normally, $V_0=1$ and $\theta_0=0$.

% Then the objective of power flow equations is to solve for the voltages and currents as functions of the slack bus voltage and power at each bus. In short,
% \begin{align*}
% v_h&=v_h(V_0,\theta_0,s_1,\dots,s_n)\\
% i_h&=i_h(V_0,\theta_0,s_1,\dots,s_n)
% \end{align*}
% The construction is described below.

The power flow equations establish the relation between the power flows injected and consumed at each bus and the corresponding voltage and current levels:
%\indent For $PQ$ buses, the following holds:
\begin{equation*}
v_h\bar{i}_h=s_h\qquad \forall{h}\in\mathcal{L}:={1,\dots,n}
\end{equation*}
where the bar denotes the complex conjugate of the given vector. The power flow equations need to be supplemented with Kirchhoff laws that can be represented in the following form:
\begin{equation*}
v_L=v_0\mathbb{1}+Zi_\mathcal{L}\text{, where } Z:=Y_{\mathcal{L}\mathcal{L}}^{-1}
\end{equation*}
where $\mathbb{1}$ is a vector with all components equal to $1$, i.e. $\mathbb{1}_h = 1$.

In this work, we only consider $PQ$ buses, so the values of $s$'s are assumed to be given. The extension for $PV$ buses is rather straightforward and is also described in \cite{bolognani2014existence}.

%Following Bolognani's setup, the power flow equations are reduced to the following:

\subsection{Bolognani's certificate of solvability}

A sufficient condition for existence of a particular solution to the nonlinear power flow equations was derived in \cite{bolognani2014existence} and forms the foundation of the algorithm proposed in this work. The criterion is based on reformulation of power flow equation as fixed point of nonlinear map, i.e. equation of the form $f = G(f)$, where the map is defined via the following function:
\begin{equation}\label{map}
G(f):=-\frac{1}{V^2_0}\text{diag}(f+s_\mathcal{L})Z(\bar{f}+\bar{s}_{\mathcal{L}})
\end{equation}
and the quantity $f$ is related to bus voltage and current levels via the following relation:
\begin{equation*}
f:=v_0\bar{i}_{\mathcal{L}}-s_{\mathcal{L}}=V_0\bar{i}_{\mathcal{L}}-s_\mathcal{L}
\end{equation*}
% {\color{red} What is the definition of f ? Please replace the previous equation with the definition of f in terms of voltage and currents}

The existence of solution is certified by proving that the function $G(f)$ is contractive, and thus the solution exists. The final criterion can be compactly expressed as the following condition on the power flow levels:
\begin{equation}\label{saverio}
4\|Z\|^*_2\cdot\|s_\mathcal{L}\|_2 \leq V_0^2,
\end{equation}
with the ``nuclear'' matrix norm defined as 
\begin{align}\label{two_norm}
 \|A\|_2^*& = \max_h\sqrt{\sum_j{A_{hj}^2}}
\end{align}
Similar relation can be derived also for $1$-norms. In this case one has 
\begin{align}\label{1normcrit}
4\|Z\|^*_\infty\cdot\|s_\mathcal{L}\|_1 &\leq V_0^2 \\
 \|A\|_\infty^*& = \max_{h,k}|A_{hk}|\label{one_norm}
\end{align}
Whenever the relation \eqref{saverio} is satisfied for some power flow injection vector $s$, the power flow equations a proven to have a solution. Some additional properties of this solution were proven in the original work \cite{bolognani2014existence}.


\section{Extension of the solvability certificate} \label{sec:modify}
The solvability certificate defined by \eqref{saverio} is not unique. In fact, it was noted in the original work \cite{bolognani2014existence}, that it can be generalized by rescaling the power and voltage vectors in the original power flow equations. In this work we explore this freedom and show that these rescalings can be used to extend the criterion, and make it both less conservative and more computationally tractable. 


To extend the criterion we replace $f$ with $\Lambda\varphi$, where $\Lambda$ is some diagonal matrix in $\mathbb{R}^{n \times n}$ and $\varphi$ some vertical vector in $\mathbb{R}^n$. Then we have:

\begin{equation*}
f =\Lambda\varphi =  \begin{bmatrix}\lambda_1\varphi_1\\ \lambda_2\varphi_2\\ \vdots\\ \lambda_n\varphi_n\end{bmatrix}
\end{equation*}
It follows from \eqref{map} that
\begin{equation}
\varphi_i=-\frac{1}{V_0^2}(\varphi_i+\lambda_i^{-1}s_i)\sum\limits_{j=1}^n\left(Z_{ij}\lambda_j\right)(\bar{\varphi}_j+\lambda_j^{-1}\bar{s}_j)
\end{equation}
Notice that this equation has the same structure as \eqref{map}, so the same techniques that were used to prove the solvability based on \eqref{saverio} can be applied to this equation. Then the revised criterion can be rewritten as 
\begin{equation} \label{rescaled}
4\|Z\Lambda\|^*\cdot \|\Lambda^{-1} s_\mathcal{L}\| \leq V_0^2
\end{equation}
The key idea behind the approach explored in this work is to characterize the whole union of certificates defined by the inequalities \eqref{rescaled} for all possible diagonal matrices $\Lambda$. More precisely we are interested in characterization of the set defined formally as 
\begin{align}\label{pset}
 \mathcal{P} = \left\{s_\mathcal{L}: \,\, \exists \Lambda\textrm{ such that } 4\|Z\Lambda\|^*\cdot \|\Lambda^{-1} s_\mathcal{L}\| \leq V_0^2\right\}
\end{align}
In the following sections we present the results of our analysis. In the section \ref{sec:simulation} we report the results of numerical analysis on various models, whereas in  section \ref{sec:outermost} we provide preliminary results from theoretical characterization of the set $\mathcal{P}$.

\section{Simulation results} \label{sec:simulation}
We start our analysis by considering simplest possible radial system composed of $2$ and $3$ buses. To visually characterize the union of certificates generated by different matrices $\Lambda$ in \eqref{rescaled} we can simply plot the boundaries corresponding to different $\Lambda$ on top of each other. 




In Figure \ref{twobus} we present results of analysis of a simple two bus system with one slack bus and one $PQ$ bus characterized by levels of power consumption defined by $P$ and $Q$. In this model we assume, without loss of generality, that the line is purely resistive. More realistic situations can be recovered by rotating the vector $s$ in complex plane by an angle defined by the impedance of the lines. We overlay the plots of the actual solution boundary and the certificates constructed with the original criterion \eqref{saverio} as well as the modified ones \eqref{rescaled}. Several observations can be made from Figure \ref{twobus}. First, it can be clearly seen that certificates defined by \eqref{saverio} and \eqref{rescaled} do not cover the whole region where the solution to power flow equation exists. Second, one can see that rescaling does affect the region certified by the criterion. Finally, one can see that the original criterion \eqref{saverio} is the least conservative boundary, so rescaling does not increase the size of the region certified to have a power flow solution. As we will see, this last observation does not apply to more general systems discussed below.

\begin{figure}[t!]
\centering
\includegraphics[width=0.9 \columnwidth]{twobus.eps}
\caption{The dotted curves illustrate the solvable boundary from varying rescalings for a $2$-bus network. The rescaled regions cover only fractions of the solution region certified by Bolognani's criterion, hence represents more conservative solution boundaries.}
\label{twobus}
\end{figure}
Next, we consider a more complicated case of a $3$-bus radial system with two $PQ$ buses. In the $3$-bus system, both of the power lines between buses $1-2$ and $2-3$ have the same values of resistance and inductance: $r=0.0734\, p.u.$ and $x=0.2581\, p.u.$. We consider the case in which $P_1=P_2$ and $Q_1=Q_2$. To understand the impact of the rescaling we consider diagonal matrices $\Lambda$ with elements in the range $\Lambda_{ii} \in [0.5, 25]$. The impedance matrix $Z$ in this problem is given by
\begin{equation*}
Z=\begin{bmatrix}-\frac{1}{7}&-\frac{1}{7}\\-\frac{1}{7}&-\frac{64}{203} + j\frac{2}{29}\end{bmatrix}
\end{equation*}
Figures \ref{64ellipse} and \ref{64rhombus} illustrate unions of $P$ \textit{v.s.} $Q$ plots with regions corresponding to certificates produced by different diagonal matrices $\Lambda$. Figure \ref{64ellipse} represents the union of 64-many $PQ$ plots corresponding to $||\cdot||_2^*$ norm as  and Figure \ref{64rhombus} represents the union of plots from $||\cdot||_\infty^*$ norm. The key observation that one can draw from these plots is that the outer boundary of the union of all certificates has a simple polytope, that we refer to as a ``Rhombus'' in our figures due to the shape of its projection. Second, as one can see the outer boundary improved the certificate provided by each individual $\Lambda$ and provides a less conservative way of checking the system solvability. Finally, it is important to note, that the outer boundary of both the $||\cdot||_\infty^*$-norm union and $||\cdot||_2^*$-norm union is the same. This observation suggests that this boundary can be characterized by a simple algebraic expression. 

\begin{figure}[t!]
\centering
\includegraphics[width=0.9 \columnwidth]{64ellipse.eps}
\caption{64 elliptical solvable regions for the given 3-bus network. Note that the contour of the union of the ellipses forms a geometry identical to that of the union of rhombuses in Figure \ref{64rhombus}.}
\label{64ellipse}
\end{figure}

\begin{figure}[t!]
\centering
\includegraphics[width=0.9 \columnwidth]{64rhombus.eps}
\caption{64 rhombus-shaped solvable regions for the 3-bus network.}
\label{64rhombus}
\end{figure}



%Figure \ref{3and10},
On figure  \ref{1and2} we provide yet another empirical confirmation of our observation by looking at the projections plots of the 13-bus network, assuming $P_1=P_2=P_3=P_4=P_5=\frac{1}{2}P_6=\frac{1}{2}P_7=\frac{1}{2}P_8=\frac{1}{2}P_9=\frac{1}{2}P_{10}=\frac{1}{4}P_{11}=\frac{1}{4}P_{12}$ and $Q_1=Q_2=Q_3=Q_4=Q_5=\frac{1}{2}Q_6=\frac{1}{2}Q_7=\frac{1}{2}Q_8=\frac{1}{2}Q_9=\frac{1}{2}Q_{10}=\frac{1}{4}Q_{11}=\frac{1}{4}Q_{12}$. The $13$-bus network configuration and data is described in \cite{nguyen2014voltage}. Moreover, fixed shunt capacitors are installed at buses $3$, $4$, $6$, $9$, $10$, $13$ with the capacities of $0.2$, $0.1$, $0.05$, $0.05$, $0.05$, $0.1$, respectively. The unit of the capacities are in $p.u.$. The impedance matrix $Z\in\mathbb{R}^{12\times12}$ is then computed, and the plot for $P$ \textit{v.s.} $Q$ is plotted using the modified criterion, using both norms $||\cdot||^*_2$ and $||\cdot||^*_\infty$. Again, we see that the union is always a simple geometric shape that looks like a rhombus on the projection plane. %{\color{red} Is it really true ? Please check convergence in figure 5. Or let's remove it in this paper and return to the question later. I'm also not sure why the maximal rhombus is exactly the same on Fig. 4 and 5. There is no reason why the polytope should have the same dimension in these two directions.} %Notice that in this particular case, the projection of $1^{st}$ and $2^{nd}$ variables' plane yields the largest outermost rhombus, hence in any dimension, it is sufficient to find the largest outermost rhombus of any projection.

%, the base power and base voltage are $S_{base} = 100\,MVA$, $V_{base} = 4.16\, kV$.

\begin{figure}[t!]
   \centering
  \subfigure{\includegraphics[width=4.3cm]{r1.eps}}
 \subfigure{\includegraphics[width=4.3cm]{e1.eps}}
  \caption{Projection of $1^{st}$ and $2^{nd}$ variables' plane for 13-bus network. On the left is the rescaled regions from the definition from (\ref{one_norm}), and on the right is from (\ref{two_norm}).}%
  \label{1and2}
\end{figure}

% \begin{figure}[!ht]
%     \centering
%     \subfigure{{\includegraphics[width=5cm]{r3.eps}}}%
%     \subfigure{{\includegraphics[width=5cm]{e3.eps} }}%
%     \caption{Projection of $6^{th}$ and $12^{th}$ variables' plane for 13-bus network}%
%     \label{6and12}
% \end{figure}

In the following section, based on our empirical observations we introduce a conjecture that defines the solvability certificate in terms of of the simple piecewise linear function defining the convex hull of the critical points. 

\section{Constructing the outermost rhombus} \label{sec:outermost}

In this section, we propose a method for constructing an algebraic representation of the ``Maximal rhombus'' shape that we observed in our numerical experiments. This shape is most easily constructed using the representation \eqref{1normcrit}. Consider the rescaling defined by $\lambda_k^{-1} = \max_h |Z_{hk}|$. In this case the nuclear norm of the matrix $\|Z\Lambda\|_\infty^* = 1$ that can be easily checked from the definition \eqref{one_norm}. The certificate is then reduced to 
\begin{equation} \label{convexhull}
 \sum_k \left|\frac{s_k}{s_k^{\max}}\right| \leq 1
\end{equation}
where the value of $s_k^{\max}$ is given by the following expression:
\begin{equation}
 s_k^{\max}  = \frac{V_0^2}{4 \max_h |Z_{hk}|}
\end{equation}
Note that the expression \eqref{convexhull} defines the convex hull of critical points $\mathbf{s}^{(k)} = s_k^{\max} \mathbf{e}_k$, where the unit vectors $\mathbf{e}_k$ define the rectangular coordinate system in the $s$-space. One can easily show that the polytope defined by \eqref{convexhull} is the envelope of union of polytopes defined by the rescaled version of \eqref{1normcrit}. It is also easy to check that it encompasses all the critical ellipsoids defined by the rescaled version of \eqref{saverio} with $\Lambda_k \to \infty$. At this stage we don't know the conditions under which the union of all the rescaled ellipsoids coincides with the convex hull defined by \eqref{convexhull}, however our numerical experiments show that it is often the case. 


\section{Comparison with real boundary} \label{sec:discussion}

% Figure \ref{convexhull} illustrates the procedure described above for the $3$-bus case. Points A and B indicate the two endpoints of the extreme ellipse as $\lambda_2 \rightarrow 0$, and points C and D indicate the two endpoints of the extreme ellipse as $\lambda_1 \rightarrow 0$. Taking the convex hull of point A, B, C, D produces the rhombus identical to the outermost rhombus from Figure \ref{64rhombus}.

% \begin{figure}[!ht]
% \centering
% \includegraphics[width=0.7 \columnwidth]{convexhull.eps}
% \caption{Convex hull of the four endpoints}
% \label{convexhull}
% \end{figure}

As the certificate can be naturally used to guarantee the existence of power flow equation solutions it is reasonable to assess its conservativeness by comparing the certified region to the actual solution boundary. On Figure \ref{realboundary} we perform this comparison for the 3-bus example system. Here the outermost curve indicates the actual solution boundary found via Groebner basis approach explained in \cite{HungTuritsyn}. The polytope shown represents the first quadrant part of the convex hull as defined in \eqref{convexhull}. As expected, the polytope lies inside the solvability region, and provides a remarkably tight approximation of this boundary. %{\color{red} I suggest to remove the references to Vertical and Horizontal ellipse references from the figure, as we don't use them anymore in our construction }{\color{green} From Essie: I got rid of the two lines (extreme ellipses). Let me know if this was not the correct modification.}% Lastly, the proposed convexhull method was implemented in order to find the extreme points in the first quadrant. One can see that not only is this new linear boundary contained inside the actual solution boundary, but also that a new solution boundary can be found around a desired solution. 

% {\color{red} Hung, insert the right citation please, also correct the bib file to refer to our conference paper}

\begin{figure}[t!]
\centering
\includegraphics[width=0.9 \columnwidth]{realbound2.eps}
\caption{The plot displays two different solvable boundaries of $3$-bus system; one, the actual solution boundary for a 3-bus network, and two, the contour of the union of elliptical rescaled solvable boundaries. The contour was plotted using the method for constructing the "Outermost Rhombus".}
\label{realboundary}
\end{figure}

\begin{figure}[ht]
\centering
\includegraphics[width=0.85
\columnwidth]{PVcurve.eps}
\caption{The $PV$ curves at load bus $2$ with different levels of reactive power in the $3$-bus system.}
\label{PVcurve}
\end{figure}

As a supplement to Figure \ref{realboundary}, Figure \ref{PVcurve} presents $3$ $PV$ curves corresponding to different fixed values of reactive power $Q_2=Q_3$. $P_A$ and $P_E$ indicate the actual value and estimated value of maximum active power, respectively. For $Q_2=0$, the two values match as shown at point $B$ in Figure \ref{realboundary}. Interestingly enough, the certificate is tight for purely active power flows or DC power flow, but is more conservative as the amount of reactive power consumption is increased. %As one can observe, at $P_A$ the system encounters saddle-node bifurcation.

\section{Practical applications} \label{sec:apps}
In this section, we propose two important applications of the solvability certificate derived above.

\subsection{Non-iterative approximation of solvable boundaries}
In steady state analysis or power flow problem, the solvable boundary consists points corresponding to saddle-node bifurcation where the Jacobian matrix becomes singular. Knowing the solvable boundary of the system is essential to system operators because it indicates the limit of the system and the threats to viability of the system operation \cite{overbye1994power}. However, constructing the solvable boundary in multi-dimensional parameter space is difficult due to the complexity of the system \cite{venkatasubramanian1995local}. It has been drawing much attention for decades, and most of the proposed techniques are computationally prohibitive iterative methods \cite{overbye1994power,overbye1995computation,makarov2000computation, Chiang2008continuation}, which may suffer from the divergence problem. On the other hand, the simple linear criterion derive in this work a reasonable alternative to characterize non-iteratively the solvable boundary. Besides computational efficiency provided, a prompt technique for visualization of the solvable boundary is extremely helpful not only in contingency analysis and planning processes, but also in emergency controls where faster decision-making is required. Moreover, by construction, the approximated region bounded by the outermost rhombus is  inscribed in the actual solvability boundary; hence, the approximated boundary is closer to the feasible region in which the operational constraints are satisfied. %We envision that the approximated region is the intermediate one between the actual and the feasible ones. In this sense, the approximated boundary is more meaning to the system operators than the actual solvable boundary.

%Despite of unavoidable errors due to approximation which, however, in our test cases, are rather small;  

\subsection{Solvability criterion for security constraints}
Another important application can be developed from simple the solvability criterion described in \eqref{convexhull}. The simple criterion can be incorporated as a security constraint in planning and operational problems such as OPF. Whenever the solvability criterion is satisfied, the system is guaranteed to possess a solution to the corresponding power flow problem. The simple criterion \eqref{convexhull} is also useful for corrective control. For example, one can use this to find optimal load shedding.

Moreover, note that the proposed simple criterion can be expressed in terms of the slack bus voltage, $V_0$. Other forms of solvable criteria are described in \cite{thorp1986reactive, molzahn2012sufficient}. In the future work, we plan to extend the proposed criterion to take in account the feasibility of the solutions, including voltage and current criteria.  We expect that the new criterion could be applied directly to assess quasi-static voltage stability and would help to avoid voltage collapse in static scenarios.

\section{Conclusions}

In this work, we extended the solvability criterion derived in \cite{bolognani2014existence} by exploring the freedom associated with rescaling matrices $\Lambda$. We observed, that the union of the regions described by these certificate has often a simple polynomial shape and can be represented as a convex hull of relatively small number of points represented by the expression \eqref{convexhull}. This simple solvability certificate  depends only on the system configuration and can be naturally used as voltage security constraint.

In future works, we plan to develop a new method to identify feasible load change and design secondary voltage control based on the definition of the polytopic solvability certificate. Moreover, we plan to explore the possibility of using this criterion as static voltage stability indicator.

% \section{Meeting from 10/14/2014 Paper Outline}

% \noindent Introduction

% \begin{itemize}
%     \renewcommand{\labelitemi}{$\cdot$}
% \item Penetration of renewables
% \item Voltage stability problem in presence of renewables
% \item Previous studies: Overbye, Conejo, Ilic, Bolognani
% \item Our contribution: computationally tractable extended solvability certificates
% \end{itemize}

% \noindent Model
% \begin{itemize}
%     \renewcommand{\labelitemi}{$\cdot$}
% \item Power flow equations + motivations
% \item Bolognani's bound: review
% \item Extension ($\Lambda$) - set of certificates
% \item Union over many $\Lambda$'s $\rightarrow$ envelopes
% \end{itemize}

% \noindent Simulations for 2-bus system
% \begin{itemize}
%     \renewcommand{\labelitemi}{$\cdot$}
% \item Ellipsoids, Rhombus 
% \end{itemize}
% \noindent Algorithm for constructing outermost Rhombus
% \begin{itemize}
%     \renewcommand{\labelitemi}{$\cdot$}
% \item Derivation
% \item Prove that every ellipsoid is inside convex hull
% \end{itemize}
% \noindent Simulations
% \noindent Conclude
% \begin{itemize}
%     \renewcommand{\labelitemi}{$\cdot$}
% \item New technique for voltage stability assessment 
% \end{itemize}


% \section{Power flow approximation}
% The summation of Bolognani's paper goes here.

% \section{Solvability boundary approximation}

% Put your $\Lambda$ modified criterion there.

% \section{Simulations}

% \subsection{Solvability boundary approximation}
% \subsection{Feasible load change and secondary voltage control}
% \subsection{Probability of losing solvability}

% \section{Discussions}

\section{Acknowledgement}

We thank the NSF, MIT/Skoltech, and Masdar Initiative, as well as the Vietnam Education Foundation for their supports. 

% \appendix \label{app}

% \begin{figure}[!ht]
%     \centering
%     \includegraphics[width=0.5 \columnwidth]{figures//13bus.eps}
% 	\caption{The thirteen-bus feeder test case}
%     \label{fig:13busfeeder}
% \end{figure}

% \begin{table}[ht]
%     \caption{The thirteen-bus test feeder branch data}
%     \centering
%     \label{branchdata}
%     \begin{tabular}{|c |c |c||c |c |c|}
%     \hline
%     \textbf{From} & \textbf{To} & \textbf{Line impedance}& \textbf{From} & \textbf{To} & \textbf{Line impedance} \\
%     \textbf{bus} & \textbf{bus } & \textbf{$(p.u.)$}&  \textbf{bus } & \textbf{bus} & \textbf{$(p.u.)$}\\
%               \hline
%     $1$ & $4$ & $0.041+j0.131$&$8$ & $9$ & $0.027+j0.070$ \\
%     $2$ & $3$ & $0.027+j0.070$ &$4$ & $9$ & $0.041+j0.131$  \\
%     $3$ & $4$ & $0.045+j0.100$ &$9$ & $10$ & $0.006+j0.020$\\
%     $4$ & $5$ & $0.032+j0.070$ &$10$ & $11$ & $0.027+j0.070$\\
%     $5$ & $6$ & $0.100+j0.240$ &$8$ & $12$ & $0.039+j0.090$\\
%     $7$ & $8$ & $0.015+j0.050$ &$9$ & $13$ & $0.020+j0.065$\\
    
%     \hline
%     \end{tabular}
% \end{table}
% The $13$-bus radial system is constructed based on the IEEE $13$-bus distribution feeder test case described in \cite{IEEEtestcase}. The network configuration is shown in Figure \ref{fig:13busfeeder} and the branch data can be found in Table \ref{branchdata}. Moreover, fixed shunt capacitors are installed at buses $3$, $4$, $6$, $9$, $10$, $13$ with the capacities of $0.2$, $0.1$, $0.05$, $0.05$, $0.05$, $0.1$, respectively. The unit of the capacities are in $p.u.$, the base power and base voltage are $S_{base} = 100\,MVA$, $V_{base} = 4.16\, kV$.

\bibliographystyle{IEEEtran}
\bibliography{main.bbl}
\end{document}

\section{Meeting from 8/27/2014 Paper Outline}

\renewcommand\labelitemi{$-$}
\renewcommand\labelitemii{$\hookrightarrow$}
\begin{itemize}
\item Intro
    \begin{itemize}
    \item Voltage Stability
    \item Security assessment
    \item Renewable $\rightarrow$ Probabilities
   % \item Secondary voltage
    \end{itemize}
\item Certificate
    \begin{itemize}
    \item Computationally fast, approximate
    \item Bolognani
    \item Extension
        \begin{itemize}
        \item $\Lambda$
        \item Different operating point
        \end{itemize}
    \item Envelope, Analysis
    \end{itemize}
\item Algorithm for certificate construction:
    \begin{itemize}
    \item Closest ellipsoid (understand zero method)
    \item Convex hull approach
    \end{itemize}
\item Applications
    \begin{itemize}
    \item Solvability boundary
    \item Lower bound on loading
    \item Probabilistic guarantees
   % \item Voltage regulation
    \end{itemize}
\item Simulations:
    \begin{itemize}
    \item 2 bus
    \item 3 bus + \underline{real boundary}
    \item 39 bus: comparison with NR/Continuation, time and margin
    \end{itemize}
\end{itemize}










%polynomial that we named as the "outermost rhombus". Empirical observation of the plots indicate that there exists the outermost rhombus that encompasses all other rhombuses and ellipsoids produced from varying $\Lambda$ matrices, directing us to the following conjecture:

% \vspace{2mm}

% \begin{conjecture}
% The convex hull of the four endpoints of two extreme ellipses encompasses the outermost solution area that may be defined by varying lambda.
% \end{conjecture}

% % \vspace{4mm}
% % \textbf{\textit{Conjecture}} \textit{The convex hull of the four endpoints of two extreme ellipses encompasses the outermost solution area that may be defined by varying lambda.}

% % \vspace{4mm}

% Let us be given a diagonal matrix:

% \begin{equation*}
% \Lambda=\begin{bmatrix}\lambda_1&0\\0&\lambda_2\end{bmatrix}
% \end{equation*}

% We claim that $||Z\Lambda||_2^*||\Lambda^{-1}s||_2=c$ produces an ellipsoid, and that the union of such ellipsoids is a rhombus. Moreover, this rhombus is the convex hull of the end-points of the two extreme ellipsoids in the set of ellipsoids defined above. First consider the case for dimensions in $\mathbb{R}^2$ and $\mathbb{R}^{2\times 2}$. 

% \begin{align}\label{ellipse 1}
% c&=||Z\Lambda||_2^*||\Lambda^{-1}s||\\
% &=\max_h{\sqrt{Z_{h1}^2\lambda_1^2+Z_{h2}^2\lambda_2^2}}\sqrt{\frac{P^2}{\lambda_1^2}+\frac{Q^2}{\lambda_2^2}}
% \end{align}

% Since $\max_h{\sqrt{Z_{h1}^2\lambda_1^2+Z_{h2}^2\lambda_2^2}}$ is just some constant, equation (\ref{ellipse 1}) can be rewritten as follows:

% \begin{equation}\label{ellipse 2}
% \frac{P^2}{\lambda_1^2}+\frac{Q^2}{\lambda_2^2}=d
% \end{equation}

% where $d=c/\max_h{\sqrt{Z_{h1}^2\lambda_1^2+Z_{h2}^2\lambda_2^2}}$ is some constant. Then (\ref{ellipse 2}) represents a non-tilted 2D ellipse with axis of lengths $\lambda_1$ and $\lambda_2$. 

% WLOG, fix $\lambda_1$ and let $\lambda_2\rightarrow 0$. By doing so, we can eliminate one of the two radii so that the ellipse practically becomes a line between two points, and obtain the maximum value for $\lambda_1$. Repeat the same procedure for $\lambda_2$. Since the rhombus must contain all ellipse of the form \eqref{ellipse 2}, it must be the case that the rhombus also encompasses the two extreme ellipses as well; hence the rhombus is exactly the convex hull of the end points of the two extreme ellipses.

% We can apply the idea above into a more general dimension and use the definition from (\ref{one_norm}) in order to establish a simpler criterion for the multi-dimensional convex hull. Given

% \begin{equation*}
% \begin{aligned}
% s&=\{s_1,\dots,s_m\}\\
% \Lambda&=\text{diag}(\lambda_1,\dots,\lambda_m)
% \end{aligned}
% \end{equation*}

% the modified criterion can be rewritten as:

% \begin{equation}\label{general case}
% \sqrt{\sum \frac{s_k^2}{\lambda_k^2}}=\frac{c^*}{\max_h \sqrt{\sum_k Z_{hk}^2}\lambda_k^2}
% \end{equation}

% where 

% \begin{equation*}
% \begin{aligned}
% c^*&=||Z\Lambda||_2^*||\Lambda^{-1}s||_2^*\quad\text{and}\\
% ||Z\Lambda||^*_2&=\max_h \sqrt{\sum_k Z_{hk}^2\lambda_k^2}
% \end{aligned}
% \end{equation*}

% In order to find the outermost rhombus, extreme ellipses need to be found by setting $\lambda_{k\neq l}\rightarrow0$. Then (\ref{general case}) can be rewritten as

% \begin{equation} \label{eq:smax}
% s_{k,max}=\frac{c^*}{\max_h|Z_hk|}
% \end{equation}

% Let $\mathcal{R}$ be the convex hull of points $s_{i,max}$ ($i=1,\dots,m$). Then by the definition of a convex hull, we have

% \begin{equation*}
% \begin{aligned}
% \mathcal{R}=\{&\mathcal{S}|\mathcal{S}=s_{k,max}e_k(x_k-y_k),\\
% &\text{ }x_k,y_k\geq0 \text{ and } \sum(x_k-y_k)\leq 1\}
% \end{aligned}
% \end{equation*}

% By combining the expression from (\ref{eq:solcri}) and the definition of a matrix norm in (\ref{one_norm}), we can define the set

% \begin{equation*}
% \Big\{s|\text{ }||\Lambda^{-1}s||_1<\frac{c^*}{||Z\Lambda||_1}\Big\}
% \end{equation*}

% where

% \begin{equation*}
% x_k-y_k=\frac{e_ks}{s_{k,max}}
% \end{equation*}

% We can then condense to a simpler expression and summarize into the following conjecture:

% \begin{conjecture} The outermost rhombus $\mathcal{R}$ can be constructed by taking the convex hull of $s_k$:
% \begin{equation} \label{eq:solcri}
% \mathcal{R}=\Big\{s_k|\sum|s_k/s_{k,max}|\leq1\Big\}
% \end{equation}
% \end{conjecture}

% {\color{blue}
% \section{Solvability certificate analysis}
% Consider the power flow equation solvability certificate defined in the work of Bolognani. 

% \begin{equation}
%  \|Z \Lambda\|^*_2 \cdot\|\Lambda^{-1} s\|_2 \leq c
% \end{equation}
% Each choice of diagonal matrix $\Lambda$ defines an interior of an ellipsoid with primary axes aligned with Cartesian unit vectors $\mathbf{e}_k$. The union of the ellipsoid interiors can be formally defined as a the following set:
% \begin{equation}
%  \mathcal{P} = \left\{s: \exists \Lambda \textrm{ such that }  \|Z \Lambda\|^*_2 \cdot\|\Lambda^{-1} s\| \leq c\right\}
% \end{equation}
% Our goal is to characterize the geometry of the set $\mathcal{P}$ and prove that it can be represented as a convex hull of some simple set of points. To characterize the outer boundary of the set $\mathcal{P}$ we will use the method of supporting hyperplanes. Within this approach we will characterize the set of hyperplanes defined by $\mu^\top \left[s-s^*(\mu)\right] = 0$ that touch the boundary of the region $\mathcal{P}$ at points $s^*(mu)$. These hyperplanes are characterized by the fact that there are no points in $\mathcal{P}$ for $\mu^\top s > \mu^\top s^*$. In this case the $\mu$ we can define the following optimization problem to find the values of $s^*$:
% \begin{align}
% \mu^\top s^*(\mu) = & \max_{s,\Lambda} \mu^\top s \\
% \textrm{subject to: } &\|Z \Lambda\|^*_2 \cdot\|\Lambda^{-1} s\|_2 \leq c
% \end{align}
% Note that the actual point of the boundary $s^*(\mu)$ solves the optimization problem and can be recovered via analysis of its solution. A special case of interest in the context of power system corresponds to situation when the optimum is achieved on some values of $s$ that do not depend on the hyperplane normal vector direction $\mu$ at least in some region of $\mu$. In this case the hyperplanes touch the set $\mathcal{P}$ in some countable set of points $s^*$ and the set itself can be characterized as convex hull of these points yielding a simple polytope geometry.

% In the following we will transform the problem into simple linear optimization and discuss the conditions under which the set $\mathcal{P}$ is a simple polytope. We will solve this problem in two steps first optimizing with respect to $s$ and then performing the optimization with respect to $\Lambda$. Optimization with respect to $s$ can be easily performed using the standard Lagrangian techniques. Using the Lagrangian
% \begin{equation}
%  \mathcal{L}(s,\kappa) = \mu^\top s + \kappa \|\Lambda^{-1} s\|_2^2 =\mu^\top s + \kappa s^\top \Lambda^{-2} s
% \end{equation}
% we get the following expression for the optimal $s^*(\Lambda)$:
% \begin{equation}
%  s^*(\Lambda,\mu) = \frac{c}{\|Z \Lambda\|^*_2} \frac{\Lambda^2 \mu}{\sqrt{\mu^\top \Lambda^2 \mu}}
% \end{equation}
% The supporting hyperplane can be now found by optimizing the value of $\mu^\top s^*(\Lambda,\mu)$ with respect to $\Lambda$. This yields the following optimization problem:
% \begin{align}
% \mu^\top s^*(\mu) = \max_\Lambda c \frac{\sqrt{\mu^\top \Lambda^2 \mu}}{\|Z \Lambda\|^*_2}
% \end{align}
% Note, that the nuclear norm of the matrix can be itself represented as $\|Z \Lambda\|^*_2 = \max_h \sqrt{Z_h^\top \Lambda^2 Z_h}$ which leads us to the following problem
% \begin{align} \label{lambdaproblem}
% \mu^\top s^*(\mu) =  \max_\Lambda\min_h c \sqrt{\frac{\mu^\top \Lambda^2 \mu}{Z_h^\top \Lambda^2 Z_h}} = 
%  \max_\Lambda\min_h 
%  c \sqrt{\frac{\sum_k \mu_k^2 \Lambda_k^2}{\sum_k Z_{hk}^2 \Lambda_k^2}}
% \end{align}
% Let us introduce a new variable 
% \begin{equation}
%  x_k = \frac{Z_{hk}^2\Lambda_k^2}{\sum_l Z_{hl}^2 \Lambda_l^2}
% \end{equation}
% Then, the set of all possible diagonal matrices $\Lambda$ can be represented as a simplex $x_k \geq 0$ and $\sum_k x_k = \mathbf{1}^\top x = 1$. The problem \eqref{lambdaproblem} can be then rewritten as
% \begin{align}
% \mu^\top s^*(\mu) =  & \max_x\min_h c \sqrt{\sum_k Q_{hk} x_k}\\
%  \textrm{subject to: } 
% & \mathbf{1}^\top x  = 1, \\
% & x\geq 0
% \end{align}
% where $Q_{hk} = \left(\mu_k/Z_{hk}\right)^2$. The inner program $\min_h \left[\dots\right]$ can be represented as the solution of the following maximization problem:
% \begin{align}
%  \min_h  c \sqrt{\sum_k Q_{hk} x_k}  = & c \sqrt{t^*} \\
%  t^* = &\max_{t} t \\
%  \textrm{subject to: } &
%  Q x \geq t \mathbf{1}, 
% \end{align}
% Hence the overall minimization problem reduces to
% \begin{align}\label{linear}
% t^* = & \max_{x,t} t \\
% \textrm{subject to: } & \label{Qx}
%   Q x \geq t \mathbf{1}, \\
% & \mathbf{1}^\top x  = 1, \label{sumx} \\
% & x\geq 0 \label{xpos} 
% \end{align}
% This is a linear programming problem in $\mathbb{R}^{n+1}$, therefore the solution lies on intersection of $n+1$ hyperplanes. There are a total of $2n+1$ hyperplanes defined in \eqref{Qx},\eqref{sumx},\eqref{xpos} and the intersection point depends on the vector $\mu$. 

% Of special interest is the case when the solution lies on the intersection of hyperplanes forming the original simplex in $x$ space, i.e. critical points of  \eqref{sumx} and \eqref{xpos}). This set of points simply coincides with unit vectors $\mathbf{e}_k$, and results in the the following expression for $t^*$:
% \begin{align}
%  t^* = \min_{h,k} Q_{hk} = \min_{h,k} \frac{\mu_k^2}{Z_{hk}^2}. = \frac{\mu_{k^*}^2}{Z_{h^*k^*}^2}
% \end{align}
% The original hyperplanes in this case are defined by 
% \begin{align}
%  \mu^\top s  = \mu_{k^*} \frac{c}{Z_{h^*k^*}} = \mu^\top \frac{c\mathbf{e}_{k^*}}{Z_{h^*k^*}}
% \end{align}
% The set $\mathcal{P}$ is then formed as a convex hull of critical points 
% \begin{equation}
%  s^* = \frac{c}{Z_{h^*k^*}}\mathbf{e}_{k^*}
% \end{equation}
% that are the same as constructed in empirical studies of $\mathcal{P}$. 

% To conclude, we have provided that the set $\mathcal{P}$ has simple ``convex hull'' structure whenever the linear program \eqref{linear} is solved by vertices of the polytope \eqref{sumx},\eqref{xpos}. It is still not clear at this stage what are the necessary and sufficient conditions responsible for this property. 
% }
main.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{KundurDSA}
K.~Morison, L.~Wang, and P.~Kundur, ``Power system security assessment,''
  \emph{Power and Energy Magazine, IEEE}, vol.~2, no.~5, pp. 30--39, Sept 2004.

\bibitem{xie2007novel}
L.~Xie, J.~Ilic, and M.~Ilic, ``Novel performance index and multi-layered
  information structure for monitoring quasi-static voltage problems,'' in
  \emph{PESGM, 2007. IEEE}.\hskip 1em plus 0.5em minus 0.4em\relax IEEE, 2007,
  pp. 1--7.

\bibitem{LeXieQSVS}
L.~Xie, Y.~Chen, and H.~Liao, ``Distributed online monitoring of quasi-static
  voltage collapse in multi-area power systems,'' \emph{Power Systems, IEEE
  Transactions on}, vol.~27, no.~4, pp. 2271--2279, Nov 2012.

\bibitem{morison1993voltage}
G.~Morison, B.~Gao, and P.~Kundur, ``Voltage stability analysis using static
  and dynamic approaches,'' \emph{Power Systems, IEEE Transactions on}, vol.~8,
  no.~3, pp. 1159--1171, 1993.

\bibitem{HungTuritsyn}
{Hung D. Nguyen and Konstantin S. Turitsyn}, ``Appearance of multiple stable
  load flow solutions under power flow reversal conditions,'' in \emph{Power \&
  Energy Society General Meeting}.\hskip 1em plus 0.5em minus 0.4em\relax IEEE,
  2014.

\bibitem{JEET}
{Hung Nguyen Dinh, Yong Tae Yoon, Minh Y Nguyen},
  ``\BIBforeignlanguage{English}{A new approach for corrective and preventive
  control to unsolvable case in power networks having ders},''
  \emph{\BIBforeignlanguage{English}{JEET}}, vol.~8, no.~3, pp. 411--420, May
  2013.

\bibitem{coffrin2014linear}
C.~Coffrin and P.~Van~Hentenryck, ``A linear-programming approximation of ac
  power flows,'' \emph{INFORMS Journal on Computing}, 2014.

\bibitem{bergen}
{Arthur R. Bergen, Vijay Vittal}, ``Power systems analysis,'' 2000.

\bibitem{Cvijic2012}
S.~Cvijic, P.~Feldmann, and M.~Hie, ``Applications of homotopy for solving ac
  power flow and ac optimal power flow,'' in \emph{Power and Energy Society
  General Meeting, 2012 IEEE}, July 2012, pp. 1--8.

\bibitem{bolognani2014existence}
S.~Bolognani and S.~Zampieri, ``On the existence and linear approximation of
  the power flow solution in power distribution networks,'' \emph{arXiv
  preprint arXiv:1403.5031}, 2014.

\bibitem{nguyen2014voltage}
H.~Nguyen, D. and K.~Turitsyn, ``Voltage multistability and pulse emergency
  control for distribution system with power flow reversal,'' \emph{arXiv
  preprint arXiv:1407.1355}, 2014.

\bibitem{overbye1994power}
T.~J. Overbye, ``A power flow measure for unsolvable cases,'' \emph{Power
  Systems, IEEE Transactions on}, vol.~9, no.~3, pp. 1359--1365, 1994.

\bibitem{venkatasubramanian1995local}
V.~Venkatasubramanian, H.~Schattler, and J.~Zaborszky, ``Local bifurcations and
  feasibility regions in differential-algebraic systems,'' \emph{Automatic
  Control, IEEE Transactions on}, vol.~40, no.~12, pp. 1992--2013, 1995.

\bibitem{overbye1995computation}
{T. J. Overbye}, ``Computation of a practical method to restore power flow
  solvability,'' \emph{Power Systems, IEEE Transactions on}, vol.~10, no.~1,
  pp. 280--287, 1995.

\bibitem{makarov2000computation}
Y.~V. Makarov, D.~J. Hill, and Z.-Y. Dong, ``{Computation of bifurcation
  boundaries for power systems: a new $\Delta$-plane method},'' \emph{Circuits
  and Systems I: Fundamental Theory and Applications, IEEE Transactions on},
  vol.~47, no.~4, pp. 536--544, 2000.

\bibitem{Chiang2008continuation}
S.-H. Li and H.-D. Chiang, ``Continuation power flow with nonlinear power
  injection variations: a piecewise linear approximation,'' \emph{Power
  Systems, IEEE Transactions on}, vol.~23, no.~4, pp. 1637--1643, 2008.

\bibitem{thorp1986reactive}
J.~Thorp, D.~Schulz, and M.~Ili{\'c}-Spong, ``Reactive power-voltage problem:
  conditions for the existence of solution and localized disturbance
  propagation,'' \emph{IJEPES}, vol.~8, no.~2, pp. 66--74, 1986.

\bibitem{molzahn2012sufficient}
D.~K. Molzahn, B.~C. Lesieutre, and C.~L. DeMarco, ``A sufficient condition for
  power flow insolvability with applications to voltage stability margins,''
  \emph{arXiv preprint arXiv:1204.6285}, 2012.

\end{thebibliography}
