\section{Introduction}

In this note, we only consider finite simple graphs. Let $G$ be a
graph with vertex set $V(G)$ and edge set $E(G)$. For any vertex
$v\in V(G)$, we use $N_G(v)$ and $\deg_G(v)$ to denote the set of neighbors of $v$ and the
degree of $v$ in $G$, respectively. For any $X\subseteq V(G)$, we denote by
$|X|$ the cardinality of $X$. We define $G-uv$ to be the graph obtained from $G$ by
deleting the edge $uv\in E(G)$, and $G+uv$ to be the graph obtained
from $G$ by adding an edge $uv$ between two non-adjacent vertices
$u$ and $v$ of $G$. For two vertices $u$ and $v$ of $G$, the distance between $u$ and $v$ in $G$ is denoted
by $d_{G}(u, v)$.\\
For an integer $m\geqslant 2,$ let $\alpha^{m}(G)$ denote the number defined by
$$\alpha^{m}(G)=\max\{ |S|:S\subset V(G),d_{G}(x,y)\geqslant m\,\quad\text{for all vertices }\,x,y\in S\}.$$
For an integer $p\geqslant 2$, we define
$$\sigma_p^{m}(G)=\min\left\lbrace \sum_{a\in S}\deg_{G}(a):S\subset V(G),|S|=p, d_{G}(x,y)\geqslant m\quad\text{for}\;x,y\in S\right\rbrace.$$
For convenience, we define $\sigma^{m}_{p}(G)=+\infty$ if $\alpha^{m}(G)<p$. We note that,  $\alpha^{2}(G)$ is often written $\alpha(G)$, which  is the independent number of $G,$ and $\sigma_p^{2}(G)$ is often written $\sigma_{p}(G)$, which is the minimum degree sum of $p$ independent vertices.

Let $T$ be a tree, a vertex of degree one is called a leaf. The set of leaves of $T$
is denoted by $Leaf(T)$. The subtree $T-Leaf(T)$ of $T$ is called the $stem$ of $T$ and is denoted by $Stem(T)$. A tree having at most $l$ leaves is called $l-$ended tree and a stem having at most $l$ leaves is  called $l-$ended stem.
There are several well-known conditions (such as the independence
number conditions and the degree sum conditions) ensuring that a
graph $G$ contains a spanning tree with a bounded number of leaves
or branch vertices (see the survey paper~\cite{OY} and the
references cited therein for details). Win~\cite{Wi} obtained a
sufficient condition related to the independence number for
$k$-connected graphs, which confirms a conjecture of Las
Vergnas~\cite{LV}. Broersma and Tuinstra~\cite{BT} gave a degree
sum condition for a connected graph to contain a spanning tree with
at most $l$ leaves.
\begin{theorem}\label{theo1.1}{\rm (Win~\cite{Wi})}
	Let $G$ be a $k$-connected graph and let $l\geq 2$.  If
	$\alpha(G)\leq k+l-1$, then $G$ has a spanning tree with at most $l$
	leaves.
\end{theorem}
\begin{theorem}\label{theo1.2}{\rm (Broerma and Tuinstra~\cite{BT})}
	Let $G$ be a connected graph with $n$ vertices and let $l\geq 2$. If
	$\sigma_2(G)\geq n-l+1$, then $G$ has a spanning tree with at most
	$l$ leaves.
\end{theorem}

Recently, many researches are studied on spanning trees in connected graphs whose stems having a bounded number of leaves or branch vertices (see \cite{TZ}, \cite{KY}, \cite{KY15} and \cite{Yan} for more details). We introduce here some results on spanning trees whose stems have a few leaves.
\begin{theorem}{\rm (Tsugaki and Zhang~\cite{TZ})}
	Let $G$ be a connected graph and let $l\geqslant 2$ be an integer. If $\sigma_{3}(G)\geqslant|G|-2l+1$, then $G$ have a spanning tree whose stem has at most $l$ leaves.
\end{theorem}
\begin{theorem}{\rm (Kano and Yan~\cite{KY})}\label{KY1}
	Let $G$ be a connected graph and let $l\geqslant 2$ be an integer. If $\sigma_{l+1}(G)\geqslant|G|-l-1$, then G has a spanning tree with $l$-ended stem.
\end{theorem}
\begin{theorem}{\rm (Kano and Yan~\cite{KY})}\label{KY2}
	Let $G$ be a connected claw-free graph and let $l\geqslant 2$ be an integer. If $\sigma_{l+1}(G)\geqslant|G|-2l-1$, then G have a spanning tree with $l$-ended stem.
\end{theorem}
On the other hand, for a positive integer $t \geq 3,$ a graph $G$ is said to be  $K_{1,t}-$ free graph if it contains no $K_{1,t}$ as an induced subgraph. If $t=3,$ the $K_{1,3}-$ free graph is also called the claw-free graph. Moreover, if the maximum degree of a graph $G$ is
denoted by $\Delta(G)$ then $G$ is nothing but a $K_{1,t}-$free graph for all $t\geq \Delta(G)+1.$
In this note, we would like to introduce a generalization of above theorems. We study on the spanning tree of a $K_{1,t}-$free graph whose stem has a bounded number of leaves. In particular, we prove the following result.
\begin{theorem}\label{thm-main}
	For a positive integer $t \geq 3,$ let $G$ be a connected $K_{1,t}$-free graph and let $l(\not= t-2)$ be an integer. If
\begin{equation}\label{eq1}
\sigma^{4}_{l+1}(G)\geqslant |G|-\lfloor \dfrac{l(t-1)}{t-2}\rfloor-1
\end{equation}
	then $G$ has a spanning tree with $l$-ended stem. Here, the notation $\lfloor r\rfloor$ stands for the biggest integer not exceed the real number $r.$
\end{theorem}
We also note that the reason why we consider $\sigma^{4}_{l+1}(G)$ is based on the following theorem of Kano and Yan.
\begin{theorem}{\rm (Kano and Yan~\cite{KY})}\label{KY}
	Let $G$ be a connected graph and let $l\geqslant 2$ be an integer. If $\alpha^{4}(G)\leq l$, then G has a spanning tree with $l$-ended stem.
\end{theorem}
By using Theorem \ref{thm-main} when $t=3,$ we have Theorem \ref{KY2}. Moreover, Theorem \ref{thm-main} is an improvement of Theorem \ref{KY1} when we consider the positive integer $t$ big enough.

 We end this section by constructing two examples to show that
the conditions of Theorem \ref{thm-main}
are sharp. Let $t, k, m$ be integers such that $t\geq 3, k\geq 2, m \geq 1$ and let $l=k(t-2)$. Let $D$ be a complete graph with $k+1$ vertices $u_1, u_2, ..., u_{k+1}$. Let $D_1, D_2, ..., D_{k(t-2)+1}$ be copies of the graph $K_m$. Let $v_1, v_2, ..., v_{k(t-2)+1}$ be vertices which are not in $D\cup D_1\cup D_2\cup \cdots \cup D_{k(t-2)+1}$. For each $i\in\{1,2,...,k\},$ join $u_i$ to all vertices of the graphs $D_{(i-1)(t-2)+1}, D_{(i-1)(t-2)+2}, ..., D_{i(t-2)} $ and join $u_{k+1}$ to all vertices of the graph $D_{k(t-2)+1}$. Join $v_j$ to all vertices of $D_j$ for all $j\in \{1;2;\cdots; k(t-2)+1\}$. Then the resulting graph $G$ is a $K_{1,t}-$free graph. Moreover, we have $|G|=k+1+(k(t-2)+1)(m+1)$ and \begin{align*}
	\sigma^4_{l+1}(G)&=\sigma^4_{k(t-2)+1}(G)=\sum\limits_{i=1}^{k(t-2)+1}\deg_G(v_i)\\
	&=(k(t-2)+1).m=|G|-k(t-1)-2=|G|-\lfloor \dfrac{l(t-1)}{t-2}\rfloor-2.
\end{align*}
But $G$ has no spanning tree with $l$-ended stem. Hence the condition (\ref{eq1}) is sharp.\\
On the other hand, when $l=t-2,$ let $D_1, D_2, ..., D_{l+1}$ be copies of the graph $K_m$. Let $w, v_1, v_2, ..., v_{l+1}$ be distinct vertices which are not in $ D_1\cup D_2\cup \cdots \cup D_{l+1}$. For each $i\in\{1,2,...,l+1\}$ join $v_i$ to all vertices of the graph $D_{i} $ and join $w$ to all vertices of the graphs $D_1, D_2, \cdots D_{l+1}$. The resulting graph is denoted by $H.$ Then $H$ is a $K_{1,t}-$free graph. Moreover, we may obtain that $|H|=1+(l+1)(m+1)$ and \begin{align*}
	\sigma^4_{l+1}(H)=\sum\limits_{i=1}^{l+1}\deg_{H}(v_i)=(l+1).m=|H|-l-2=|H|-\lfloor \dfrac{l(t-1)}{t-2}\rfloor-1.
\end{align*}
But $H$ has no spanning tree with $l$-ended stem. This implies that the condition $l\not=t-2$ is necessary.

\section{Proof of Theorem \ref{thm-main}}
In this section, we extend the idea of Kano and Yan in~\cite{KY} to prove
Theorem~\ref{thm-main}.

 We prove the theorem
by contradiction. Suppose to the contrary that $G$ contains no spanning tree with $l-$ended stem. We choose a maximal tree $T$ with $l-$ended stem in $G$ so that

(C) $|Leaf(T)|$ is as large as possible.\\
By the maximality of $T$, we have the following claim.
\begin{claim}\label{claim1}
	For every $v\in V(G)-V(T), N_{G}(v)\subseteq Leaf(T)\cup(V(G)-V(T))$.
\end{claim}
 Because $G$ is connected and $T$ is not a spanning tree of $G,$ there exist two vertices $v_1\in V(G)-V(T)$ and $v_2\in Leaf(T)$ such that $v_1v_2\in E(G)$.
We may obtain that $Stem(T)$ has exactly $l$ leaves. Indeed, otherwise we consider the tree $T'=T+v_1v_2.$ Then $T'$ has $l$-ended stem and $|T'|>|T|$, this implies a contradiction with the maximality of $T$. Let $\{x_1, x_2, ..., x_l\}$ be the leaf set of $Stem(T)$.
\begin{claim}\label{claim2}
	For every
	$x_{i}(1 \leqslant i \leqslant l)$, there exists a vertex $y_{i}\in Leaf(T)$ such that $y_{i}$ is adjacent to $x_{i}$
	and $N_{G}(y_{i})\subset Leaf(T)\cup\left\{x_{i}\right\}$.
\end{claim}
\noindent {\it Proof.}
	 By the maximality of $T,$ it is easy to see that for each leaf $x\in Leaf(Stem(T)),$ there exists at least a vertex $y$ in $Leaf(T)$ such that $y$ is adjacent to $x.$ Suppose that for some $1\leqslant i\leqslant l$, each leave $y_{i_j}$ of $T$ adjacent to $x_{i},$ is also adjacent to a vertex $z_{i_j}\in (Stem(T)-\left\{x_{i}\right\})$. Then we consider $T'$ to be the tree obtained from $T$ by removing the edge $y_{i_j}x_{i}$ and adding the edge $y_{i_{j}}z_{i_{j}}$. Hence $T'$ is a tree with $l-$end stem such that $|T'|=|T|$ and $Leaf(T')=Leaf(T)+\{x_i\}$, which contradicts the condition (C). Therefore, for each $x_i$, there exists a leaf $y_i\in N_G(x_i)$ such that $N_G(y_i)\cap (Stem(T)-\{x_i\})=\emptyset$. By the maximality of $T$ we also see that  $N_{G}(y_i)\cap (V(G)-V(T))=\emptyset.$ The claim holds.
\hfill \rule{4pt}{7pt}
\begin{claim}\label{claim3}
	For any two dictinct vertices $y,z\in \left\{v_1,y_{1},y_{2},\dots,y_{l}\right\},d_{G}\left(y,z\right)\geqslant 4$.
\end{claim}
\noindent {\it Proof.}
	First, we show that $d_{G}\left(v_1,y_{i}\right)\geqslant 4$ for every $1\leqslant i\leqslant l$. Let $P_i$ be the shortest path connecting $v_1$ and $y_{i}$ in $G$.  If all the vertices of $P_i$ between $v_1$ and $y_{i}$ are contained in $Leaf(T)\cup (V(G)-V(T))\cup\left\{x_{i}\right\},$ add $P_{i}$ to $T$ (if $P_{i}$ passes through $x_{i}$, we just add the segment of $P_{i}$ between $v_1$ and $x_{i})$ and remove the edges of $T$ joining $V(P_{i}\cap Leaf(T))$ to $V(Stem(T))$ except the edge $y_{i}x_{i}.$ The resulting tree is denoted by $T'.$ Then $T'$ is a tree in $G$ with $l$-ended stem  and $|T'|>|T|$, which contradicts to the maximality of $T.$ Hence we may choose the vertex $s$ in $Stem(T)\cap P_i$ such that it is nearest to $v_1$ in $P_i$. If $s=x_j$ for some $1\leq j\leq l$, then we add the segment of $P_{i}$ between $v_1$ and $x_j$ (which is denoted by $Q$) to $T$ and remove the edges of $T$ joining $V(Q)\cap Leaf(T)$ to $V(Stem(T))$ except $x_jy_j$. Hence the resulting tree has $l$-ended stem and its order is greater then $|T|$, contradicting the maximality of $T$. Thus $s\in Stem(T)-\{x_1,...,x_l\}$. By Claims \ref{claim1} and \ref{claim2}, we have $d_G(v_1,s)\geq 2$, $d_G(s,y_i)\geq 2.$ Therefore we conclude that $d_G(v_1,y_i)=d_G(v_1,s)+d_G(s,y_i)\geq 4$.

	\indent Next, we show that $d_{G}(y_{i},y_{j})\geqslant 4$ for all $1\leqslant i<j\leqslant l$. Let $P_{ij}$ be the shortest path connecting $y_{i}$ and $y_{j}$ in $G$. We note that if $P_{ij}$ passes through $x_{i}$ (or $x_{j}$), then $y_{i}x_{i}\in E(P_{ij})$ (or $y_{j}x_{j}\in E(P_{ij})$), respectively. We consider following two cases.\\
	{\it Case 1.} All vertices of $P_{ij}$ between $y_{i}$ and $y_{j}$ are contained in $Leaf(T)\cup(V(G)-V(T))\cup\left\{x_{i},x_{j}\right\}$. Then we add $P_{ij}$ to $T$ and remove the edges of $T$ joining $V(P_{ij }\cap Leaf(T))$ to $V(Stem(T))$ except the edges $y_{i}x_{i}$ and $y_{j}x_{j}$. Hence the resulting graph has exactly a cycle, which contains an edge $e$ of $Stem(T)$ incident with a branch vertex in $stem(T)$. By removing the edge $e$ and by adding an edge $v_1v_2$, we have a resulting tree $T'$ with $l$-ended stem of order greater than $|T|$, which contradicts the maximality of $T$. \\
	{\it Case 2.} There exists a vertex $s \in P_{ij}\cap (Stem(T)-\{x_i,x_j\}).$ Then $d_G(y_i,s)\geq 2, d_G(s,y_j)\geq 2$ by Claim \ref{claim2}. This concludes that $d_G(y_i,y_j)=d_G(y_i,s)+d_G(s,y_j)\geq 4$. \\
	So the assertion of the claim holds.
\hfill \rule{4pt}{7pt}

Denote $Y=\left\{y_{1},y_2,\dots,y_{l}\right\}$. Since Claims \ref{claim1}-\ref{claim3}, we have
\begin{center}
	$N_{G}(v_1)\subseteq (V(G)-V(T)-\left\{v_1\right\})\cup (N_{G}(v_1)\cap(Leaf(T)-Y)),$\\
	$\displaystyle\bigcup_{i=1}^{l}N_{G}(y_{i})\subseteq (Leaf(T)-Y-N_{G}(v_1))\cup\left\{x_{1},\dots,x_{l}\right\}.$
\end{center}
Hence by setting $q=|N_{G}(v_1)\cap(Leaf(T)-Y)|$, we obtain
\begin{equation*}
	\begin{split}
		\deg_{G}(v_1)+\displaystyle\sum\limits_{i=1}^{l}\deg_{G}(y_{i})&\leqslant (|G|-|T|-1+q)+(|Leaf(T)|-l-q)+l\\
		&=|G|-|Stem(T)|-1.
	\end{split}
\end{equation*}
On the other hand, by the assumption of Theorem \ref{thm-main} we have $$|G|-\lfloor \dfrac{l(t-1)}{t-2}\rfloor-1\leq \sigma^4_{l+1}(G)\leq \deg_{G}(v_1)+\displaystyle\sum\limits_{i=1}^{l}\deg_{G}(y_{i}).$$
Therefore we obtain $|Stem(T)|\leq \lfloor \dfrac{l(t-1)}{t-2}\rfloor$. By combining with  $|Leaf(Stem(T))|=l$ we conclude that
\begin{align} \label{1}
	|Stem(Stem(T))|\leq \lfloor \dfrac{l}{t-2}\rfloor.
\end{align}
\begin{claim}\label{claim5}
	$N_G(v_2)\cap \{x_1,x_2,...,x_l\}=\emptyset$.
\end{claim}
\noindent {\it Proof.}
	Suppose the assertion of the claim is false. Then there exists a vertex $z\in N_G(v_2)\cap \{x_1,x_2,...,x_l\}$. Remove the edge of $T$ joining $v_2$ to $Stem(T)$ except $zv_2$ and add the edges $zv_2, v_1v_2.$ Hence the resulting tree has $l$-ended stem and $|T'|>|T|$, this contradicts to the condition (T1). Claim \ref{claim5} is proved.
\hfill \rule{4pt}{7pt}

Now, we complete the proof of Theorem \ref{thm-main} by considering following two steps.

{\bf Step 1}. $|Stem(Stem(T))|=1$.

We assume that $Stem(Stem(T))=\{u\}$. It follows from $t\geq 3$ and (\ref{1}) that $l\geq t-2.$ But $l \not= t-2$ we have $l \geq t- 1.$ By combining with Claims \ref{claim3}- \ref{claim5}, $G$ induced a $K_{1,t}$ subgraph with the vertex set $\{u, x_1, x_2,...,x_{t-1},v_2\}$, this gives a contradiction.

{\bf Step 2}. $|Stem(Stem(T))|\geq 2$.

By Claim \ref{claim5}, there exists a vertex $v_3\in N_G(v_2)\cap Stem(Stem(T))$.\\
Now, we conclude that $|N_T(v_3)\cap\{x_1,x_2,...,x_l\}|<t-2$. Indeed, otherwise, without loss of generality we may assume $x_1,x_2,...,x_{t-2}\in N_G(v_3)$. Setting $s\in Stem(Stem(T))\cap N_T(v_3)$. We consider the subgraph with the vertex set $\{v_3, v_2, s, x_1, x_2,...,x_{t-2}\}$ in $G.$ By combining with the fact that $G$ is $K_{1,t}-$free we have following two cases.\\
{\it Case 1.}  $sv_2\in E(G)$. This implies that the tree  $T'=T+sv_2+v_2v_1-sv_3$ has $l$-ended stem and $|T'|>|T|$, this contradicts to the maximality of $G$.\\
{\it Case 2.} $x_js \in E(G)$ for some $j \in \{1; ...; t-2\}.$ Then we consider the tree  $T'=T+x_js+x_jv_3+v_2v_1-sv_3.$ Hence $T'$ has $l$-ended stem and $|T'|>|T|$, this also contradicts to the maximality of $G$.\\
Therefore $|N_T(v_3)\cap\{x_1;x_2;...;x_l\}|<t-2$. \\
On the other hand, since $|Leaf(Stem(T))|=l$ and $|Stem(Stem(T))|\leq \lfloor \dfrac{l}{t-2}\rfloor$ there exists a vertex $u\in Stem(Stem(T))$ such that $|N_T(u)\cap\{x_1,x_2,...,x_l\}|\geq t-1$. Without loss of generality, we may assume $x_1,x_2,...,x_{t-1}\in N_T(u)$. Set $s\in Stem(Stem(T))\cap N_T(u)$. By repeating the same arguments as in {\it Case 2} we may conclude that $x_js \not\in E(G)$ for all $j \in \{1; ...; t-1\}.$ Then, $G$ induces a $K_{1,t}$ subgraph with vertex set $\{u, s, x_1,x_2,...,x_{t-1}\}$. This gives a contradiction with the assumption of Theorem \ref{thm-main}.\\
Therefore we complete the proof of Theorem \ref{thm-main}.

\bigskip

{\bf Acknowledgements.} The research is supported by the NAFOSTED Grant of Vietnam (No. 101.04-2018.03).


\section{Introduction}

\IEEEPARstart{A}{vailable} transfer capability (ATC) is a crucial index designed to describe how much more electric power (MW) is available to buy or sell in a specified period in a competitive electric power market. It is calculated based on a set of assumed operating conditions well before the system approaches that operational state. The U.S. Federal Energy Regulatory Commission (FERC) requires the ATC to be periodically calculated and posted on the Open Access Same-time Information System (OASIS) for public access \cite{FERC96}.

In 1996, the North American Electric Reliability Council (NERC) proposed a uniform definition for ATC and the associated terminologies \cite{NERC96a}.
By definition, ATC is a measure of the transfer capability remaining in the physical transmission network for further commercial activity over and above already committed uses. Mathematically, ATC is defined as:
\begin{equation}
\mbox{ATC = TTC -- TRM -- ETC -- CBM}
\label{eq:ATCdefinition}
\end{equation}
where TTC denotes the total transfer capability, TRM denotes the transmission reliability margin, ETC denotes the existing transmission commitments (base case) which include retail customer service, and CBM denotes the capacity benefit margin (CBM).
More detailed explanations for the terminologies are presented below.
\begin{itemize}
\item TTC is the amount of electric power that can be transferred over the interconnected transmission network in a reliable manner while meeting all of a specific set of defined pre- and post-contingency system conditions (i.e., thermal, voltage, and stability limits).
The TTC is the same as the first contingency total transfer capability (FCTTC) defined in \cite{NERC95a}.
\item TRM is the amount of transmission transfer capability necessary to ensure that the interconnected transmission network is secure under a reasonable range of \textbf{uncertainties} in system conditions. The TRM can be considered as the difference of TTC between cases with and without consideration of uncertainties brought about by, for instance, forecasted load demand and the integration of renewable energy resources (RES).
TRM can be determined by assuming a fixed percentage reduction in TTC (typically too conservative).
Reference \cite{PWSauer98a} summarized four methods to determine the TRM, of which the probabilistic approaches is preferable considering the time-variant system conditions.
\item CBM is the amount of transmission transfer capability reserved by load-serving entities to ensure access to generation from interconnected systems to meet generation reliability requirements.
Typical CBM could be a multiple of the largest generation unit within the transmission system \cite{NERC99a}.
\item ETC is the sum of existing transfer capability between the source and the sink in the base case.
\end{itemize}

Since we are interested in the probabilistic characteristics of the additional transfer capability on top of the base case (i.e., TTC -- ETC). Without loss of generality, we assume $\mbox{CBM}=0$ and use the term probabilistic ATC (PATC) to represent ATC plus TRM (i.e., TTC -- ETC) in the rest of the paper.

From the above definitions, it can be seen that the nucleus of assessing ATC is to calculate the PATC and its cumulative distribution function (CDF), based on which TRM can be estimated. As a result, the ATC can be computed directly from
\begin{equation}
\label{eq:ATCcompactdefinition}
\text{ATC} = \mathbb{E}[\text{PATC}] - \text{TRM}
\end{equation}
Indeed, the main contribution of the paper is to propose a mathematical formulation for PATC and to develop a computationally efficient yet accurate algorithm for the PATC calculation.

In the previous work, traditional deterministic ATC calculation methods can be categorized into four classes: i) Linear approximation methods (e.g., \cite{GLLandgren72a}); ii) Repeated power flow (RPF) methods (e.g., \cite{YOu03a}); iii) Continuation power flow (CPF) methods (e.g., \cite{GCEjebe98a}, \cite{HDChiang05a}); and iv) Optimal power flow (OPF) methods (e.g., \cite{GDIrisarri97a}). %
However, deterministic ATC calculation does not consider the uncertainties of the renewable energy generation and their intrinsic dependencies.

To overcome the limitation of deterministic ATC methods, various probabilistic ATC assessment methods have been discussed in the previous literature, which can be categorized into 1) Monte-Carlo simulation; 2) Bootstrap method; 3) point estimation method; 4) stochastic programming.
Monte-Carlo simulation method is the most widely used method due to its simplicity. It has been used in combination with all deterministic ATC solvers mentioned above. However, it suffers from a high computational effort that may hamper its practical applications even with efficient sampling methods like the Latin Hypercube sampling method \cite{HYu09a}.
Clustering methods \cite{MRamezani09a} were applied to speed up Monte-Carlo simulation yet at the expense of accuracy.
The Bootstrap method is powerful in constructing confidence region of ATC but is time-consuming due to the re-sampling procedure \cite{RRChang02a}.
Attempting to release the computational burden, some analytical methods have been developed utilizing mathematical approximations. In \cite{WYLi13a}, the point estimation method was proposed to estimate the variance of the TTC to determine the TRM and thus the ATC.

Apart from the methods mentioned above, another class of methods for uncertainty quantification is the meta-modeling, which aims to build a statistically-equivalent functional representation for the desired response (e.g., the PATC in this study) using a small number of model evaluations. A representative method is the polynomial chaos expansion (PCE) which has been applied in the context of power systems to study the probabilistic power flow \cite{FNi17a}, the load margin problem \cite{EHaesen09a}, and the available delivery capability problem \cite{HSheng18a}.
Another emerging method, an alternative to PCE, is the canonical low-rank approximations (LRA), which employs the canonical decomposition to express the desired response as a sum of rank-one functions \cite{KKonakli16a}. The original idea of canonical decomposition dates back to 1927 \cite{FHitchcock27a}, and the method has recently become attractive for uncertainty quantification in structural vibration problems \cite{MChevreuil15a}.

In this paper, we propose a mathematical formulation of the PATC problem considering the uncertainties from the wind power, the solar PV and the loads. More importantly, we apply a novel LRA approach to solve the important PATC problem.
The main contribution of this paper are summarized as follows:
\begin{itemize}
\item A mathematical formulation of PATC problem is developed for transmission system integrating renewable energy resources and loads that may follow various marginal distributions. This formulation considers both the pre-contingency and post-contingency cases.
\item A novel efficient yet accurate algorithm---LRA is proposed to evaluate the PATC.
Particularly, random variables with diverse marginal distributions and correlation can be accommodated using proper polynomial basis and Nataf transformation.
\item Accurate probabilistic characteristics (probabilistic density function (PDF), cumulative distribution function (CDF), mean, standard deviation, etc.) of the PATC can be achieved with much less computational efforts compared to the LHS-based MCS. Such probabilistic information can be further used to compute a more realistic TRM, leading to a more fully utilized transmission assets.
\end{itemize}

The rest of the paper is organized as follows. Section II introduces the probabilistic models of the randomness in the PATC problem. Section III proposes the mathematical formulation of the PATC problem. Section IV elaborates the low-rank approximation method and its implementation in PATC calculation. The detailed algorithm to assess the PATC is summarized in Section V. The simulation results are presented in Section VI. Conclusions and perspectives are given in Section VII.

\section{Modeling Uncertainties in the Probabilistic ATC Problem}
\normalsize
The ATC calculation is typically applied to estimate the transfer capability for a pre-specified future period before the system approaches that operational state. To obtain a reasonable and dependable ATC, the base case condition, the target transactions, a credible contingency list, the physical and operational limits, and the network response should be well defined and estimated.
The variability of diverse renewable energy sources and loads result in uncertainties in the base case, which therefore requires careful modeling.

\subsection{Projection of the Wind Generation in Base Case}
\normalsize
In this study, the wind farms are modeled as an aggregated wind turbine with equivalent parameters. The Weibull distribution provides a great fit for the wind speed in many locations around the world \cite{SHJangamshetti99a}, the probability density function of which is
\small
\begin{equation}
\label{eq:wind_pdf_weibull}
f_{V}(v)=\frac{k}{c}{{\left(\frac{v}{c}\right)}^{k-1}} \exp\left[-{{\left(\frac{v}{c}\right)}^{k}}\right]
\end{equation}
\normalsize
where ${v}$ is the wind speed, ${k}$ and ${c}$ are the equivalent shape and scale parameters, respectively.
As a result, the active power output ${{P}_{w}}$ can be calculated by the piece-wise wind speed-power output relation \cite{MAien15a}
\small
\begin{equation}
\label{eq:wind_p_val}
{{P}_{w}}(v)=\left\{ \begin{array}{*{35}{l}}
0 & v\le {{v}_{in}} \mbox{ or } v>{{v}_{out}} \\
\displaystyle \frac{v-{{v}_{in}}}{{{v}_{rated}}-{{v}_{in}}}{{P}_{r}} & {{v}_{in}}<v\le {{v}_{rated}} \\
{{P}_{r}} & {{v}_{rated}}<v\le {{v}_{out}} \\
\end{array} \right.
\end{equation}
\normalsize
where ${{v}_{in}}$, ${{v}_{out}}$ and ${{v}_{rated}}$ are the cut-in, cut-out, and rated wind speed ($m/s$), ${{P}_{r}}$ is the rated wind power ($kW$). The reactive power can be determined according to the speed control type of the wind turbine \cite{EHCamm09b} since the wind turbine can be modelled as either a constant P-Q bus or a constant P-V bus with given Q-limits.

\subsection{Projection of the Solar Generation in Base Case}
Similar to the wind farms, large solar PV power plants can be modeled as an aggregated solar PV with equivalent parameters. Typically, the solar radiation is represented by the Beta distribution \cite{ZMSalameh95a}
\small
\begin{equation}
\label{eq:solar_pdf_beta}
f_{R}(r)=\frac{\Gamma(\alpha+\beta)}{\Gamma(\alpha)\Gamma(\beta)}{{\left(\frac{r}{{{r}_{\max}}}\right)}^{\alpha-1}}{{\left(1-\frac{r}{{{r}_{\max}}}\right)}^{\beta-1}}
\end{equation}
\normalsize
where ${\alpha}$ and ${\beta}$ are the shape parameters of the distribution, ${\Gamma}$ denotes the Gamma function. The parameters are typically obtained from fitting historical solar radiation data.
$r$ and ${{r}_{\max}}$ (${W/m}^{2}$) are the respective actual and maximum solar radiations. The active power ${{P}_{pv}}$ corresponding to the solar radiation ${r}$ is determined by the piece-wise function \cite{MAien15a}
\small
\begin{equation}
\label{eq:solar_p_val}
{{P}_{pv}}(r)=\left\{ \begin{array}{*{35}{l}}
\displaystyle \frac{{{r}^{2}}}{{{r}_{c}}{{r}_{std}}}{{P}_{r}} & 0\le r<{{r}_{c}} \\
\displaystyle \frac{r}{{{r}_{std}}}{{P}_{r}} & {{r}_{c}}\le r \le {{r}_{std}} \\
{{P}_{r}} & r>{{r}_{std}} \\
\end{array} \right.
\end{equation}
\normalsize
where ${{r}_{c}}$ is a certain radiation point typically set as 150 $W/m^2$, ${{r}_{std}}$ is the solar radiation in the standard environment, ${{P}_{r}}$ is the rated power of the solar PV. Solar generation is required to be injected into the power grid at unity power factor \cite{WECC10}, and hence ${{Q}_{pv}}$ is assumed to be zero in this study.

\subsection{Projection of the Forecasted Load Demand in Base Case}

By nature, load demand is uncertain in power systems.
It is a common practice to model the load uncertainty by Gaussian distribution due to the forecasting error \cite{RBillinton08a}
\small
\begin{equation}
\label{eq:load_pdf}
f({{P}_{L}})=\frac{1}{\sqrt{2\pi}{{\sigma}_{P}}} \exp \left( -\frac{{{({{P}_{L}}-{{\mu}_{P}})}^{2}}}{2\sigma_{P}^{2}} \right)
\end{equation}
\normalsize
where the mean value ${{\mu}_{P}}$ of ${{P}_{L}}$ is provided by load forecaster, and ${{\sigma}_{P}}$ denotes the forecasting error. Generally, only the active power is predicted, whereas the reactive power is calculated under the assumption of constant power factor \cite{WYLi11a}.

\section{Mathematical Formulation of Probabilistic ATC Problem} \label{ATC_Formulation}

In this section, we present a CPF-based mathematical formulation of the probabilistic ATC problem. The near-future base case conditions are formulated as a probabilistic power flow problem in which uncertainties of the RES and loads are incorporated. Besides, target inter-area transactions are modeled as the load-generation variation vector.

The deterministic power flow equations of a $N$-bus transmission system can be represented as
\small
\begin{equation}
\label{eq:power_flow}
{\bm{f}}(\bm{x})=\begin{bmatrix}
{P_{Gi}-P_{Li}-P_{i}(\bm{x})} \\
{Q_{Gi}-Q_{Li}-Q_{i}(\bm{x})} \end{bmatrix}=0
\end{equation}
\normalsize
\noindent with
\small
\begin{equation}
\label{eq:bus_flow_out}
\begin{gathered}
P_{i}(\bm{x})=V_{i}\sum\limits_{j=1}^{N}{{V_{j}(G_{ij}\cos \theta_{ij}+B_{ij}\sin \theta_{ij})}} \\
Q_{i}(\bm{x})=V_{i}\sum\limits_{j=1}^{N}{{V_{j}(G_{ij}\sin \theta_{ij}-B_{ij}\cos \theta_{ij})}}
\end{gathered}
\end{equation}
\normalsize
where ${\bm{x}}={{[{\theta},{V}]}^{T}}$, e.g., voltage angles and magnitudes for all buses; ${P_{Gi}}$ and ${Q_{Gi}}$ are the total active and reactive generation power at bus $i$; ${P_{Li}}$ and ${Q_{Li}}$ are the total active and reactive load power at bus $i$; ${G_{ij}}$ and ${B_{ij}}$ are the real and imaginary part of the entry ${Y_{ij}}$ in the bus admittance matrix.

Let ${\bm{v}}$, ${\bm{r}}$ and ${\bm{P}_{L}}$ be the random vectors representing wind speeds, solar radiations and load variations in the pre-specified period, respectively. The resulting base case condition of a ${N}$-bus system can be described as a set of probabilistic power flow (PPF) equations. Specifically, for P-Q type buses, the PPF equations are:
\small
\begin{equation}
\label{eq:pq_bus_pq}
\begin{gathered}
P_{Gi}+P_{wi}({{v}_{i}})+P_{pvi}({{r}_{i}})-P_{Li}(P_{Li})-P_{i}(\bm{x})=0 \\
Q_{Gi}+Q_{wi}({{v}_{i}})-Q_{Li}(P_{Li})-Q_{i}(\bm{x})=0
\end{gathered}
\end{equation}
\normalsize
\noindent For P-V type buses, the corresponding PPF equations are:
\small
\begin{equation}
\label{eq:pv_bus_pvq}
\begin{gathered}
P_{Gi}+P_{wi}({{v}_{i}})+P_{pvi}({{r}_{i}})-P_{Li}(P_{Li})-P_{i}(\bm{x})=0 \\
V_{i}={{V}_{i0}} \\
Q_{Gi}=-Q_{wi}({{v}_{i}})+Q_{Li}(P_{Li})+Q_{i}(\bm{x}) \\
{{Q}_{min,i}}\le Q_{Gi}\le {{Q}_{max,i}}
\end{gathered}
\end{equation}
\normalsize
\noindent where ${P_{wi}({{v}_{i}})}$, ${P_{pvi}({{r}_{i}})}$, ${P_{Li}}$ and ${P_{Gi}}$ are the real power injection from the wind farm, the solar PV power plant, the load, and the conventional generator at bus ${i}$; ${Q_{wi}({{v}_{i}})}$, ${Q_{Li}}$ and ${Q_{Gi}}$ are the corresponding reactive power injections.
If ${Q_{Gi}}$ exceeds its limits, then the terminal bus switches from P-V to P-Q with ${Q_{Gi}}$ fixed at the violated limit.

Besides, the target transactions under the study can be described as a load-generation variation vector ${\bm{b}}$ of the system in the form of
\small
\begin{equation}
\label{eq:cpf_variation}
{\bm{b}}=\begin{bmatrix}
\Delta{\bm{P}}_{G}-\Delta{\bm{P}}_{L} \\
-\Delta{\bm{Q}}_{L}
\end{bmatrix}
\end{equation}
\normalsize

In order to study how much power can be transferred along the direction specified by vector $b$, the set of equations (\ref{eq:pq_bus_pq})-(\ref{eq:cpf_variation}) can be formulated as a set of probabilistic CPF equations in the following compact form:
\small
\begin{equation}
\label{eq:cpf_equation}
{\bm{f}}({\bm{x},\bm{\mu},\bm{\lambda},\bm{U}})={\bm{f}}({\bm{x},\bm{\mu}},{\bm{U}})-\lambda{\bm{b}}=0
\end{equation}
\normalsize
\noindent where ${\bm{x}}$ is the vector of state variables, ${\bm{\mu}}$ is the vector of control parameters such as the tap ratio of transformers, ${\bm{U}=[\bm{v},\bm{r},\bm{{P}_{L}}]}$ is the random vector describing the wind speed, the solar radiation, and the load active power. It is obvious that the set of the parameterized power flow equations become the base-case power flow equation if ${\lambda=0}$.

For reliable operation of a power system, ATC calculation is required to account for both normal operating state and the state when a contingency occurs. Typically, only the ${N-1}$ contingencies are of interest in ATC calculation. However, more complex multiple contingencies may be required by the reliability criteria in a certain individual system. Enumeration of all contingencies is unnecessary and usually leads to a too conservative ATC value, hence it is a common practice to obtain a credible contingency list from the Security Analysis (SA) module in the Energy Management System (EMS). Therefore, the probabilistic ATC formulation considering a credible contingency list can be formulated as below:
\small
\begin{equation}
\label{eq:atc_final}
\lambda_{ATC} = min\{\lambda^{(0)},\lambda^{(1)},...,\lambda^{(N_{C})}\}
\end{equation}
\normalsize
\noindent in which the transfer capability $\lambda^{(v)}$ for the $v$-th case is determined by
\small
\begin{equation}
\label{eq:atc_equation}
\begin{aligned}
& \text{maximize} & & \lambda^{(v)} \\
& \text{subject to} & & {\bm{f}^{(v)}}({\bm{x}},{\bm{\mu}},{\bm{U}})-\lambda^{(v)}{\bm{b}}=0, & (a) \\
& & & {{V}_{min,i}^{(v)}}\le V_{i}({\bm{x},\bm{\mu},\lambda^{(v)},\bm{U}})\le {{V}_{max,i}^{(v)}}, & (b) \\
& & & S_{ij}^{(v)}({\bm{x},\bm{\mu},\lambda^{(v)},\bm{U}})\le {{S}_{ij,max}^{(v)}}, & (c) \\
& & & {{P}_{min,i}}\le P_{Gi}({\bm{x},\bm{\mu},\lambda^{(v)},\bm{U}})\le {{P}_{max,i}}, & (d) \\
& & & {{Q}_{min,i}}\le Q_{Gi}({\bm{x},\bm{\mu},\lambda^{(v)},\bm{U}})\le {{Q}_{max,i}}, & (e) \\
\end{aligned}
\end{equation}
\normalsize
where $v=0$ represents the base case, and $v=1,...,{{N}_{C}}$ represent the contingency cases. Specifically, ${\bm{f}^{(0)}}$ corresponds to the pre-contingency (base case) network configuration, while ${\bm{f}^{(v)}},v=1,...,{{N}_{C}}$ corresponds to the $v$-th post-contingency network configuration.
Similarly, $[{{V}_{min,i}^{(0)}},{{V}_{max,i}^{(0)}}]$ and ${{S}_{ij,max}^{(0)}}$ are the respective normal voltage and thermal limits applied to the base case; $[{{V}_{min,i}^{(v)}},{{V}_{max,i}^{(v)}}]$ and ${{S}_{ij,max}^{(v)}} (v=1,...,{{N}_{C}})$ are the emergency voltage and thermal limits applied to the contingency cases.
${\lambda}$ is the normalized load margin under the given load-generation variation vector.
The maximum value of $\lambda$ that could be achieved without the violation of (\ref{eq:atc_equation}) gives the ATC.
Note that ${\lambda}$ is a random variable due to the random input ${\bm{U}}$. Equation ($a$) specifies that the solution must satisfy the parameterized power flow equations (\ref{eq:cpf_equation}); Equations ($b$)-($e$) imply that the solution has to satisfy typical operational and electrical constraints.

For each realization of ${\bm{U}}$, there are $N_{C}+1$ ($N_{C}$ is the size of the contingency list) deterministic cases of (\ref{eq:atc_equation}) to be solved. As a result, it is extremely time-consuming to run Monte-Carlo simulations to estimate PATC using a large number of samples. It is imperative and essential to develop an efficient and accurate method to estimate the PATC.

\section{Canonical Low-rank Approximation using Polynomial Basis} \label{LRA_Section}
This section presents a general framework of the low-rank approximation of a multivariate stochastic response function. For simplicity, we first consider a scalar response function of independent inputs. The case of dependent inputs will be addressed in Section \ref{LRA_Section}--E.

\subsection{Low-rank Approximation with Polynomial Basis}

Consider a random vector ${\bm{\xi}}$ = (${{\xi}_{1},{\xi}_{2},...,{\xi}_{n},}$) with joint probability density function (PDF) ${f_{\bm{\xi}}}$ and marginal distribution functions ${f_{\xi_{i}},i=1,...,n}$ (${\xi_{i}}$ is related with the random variables ${U_{i}}$ in \eqref{eq:atc_equation}, see Section \ref{ATC_Formulation}), then the canonical rank-$r$ approximation \cite{KKonakli16a} of the target stochastic response (e.g., \textbf{PATC in this study}) ${Y = g(\bm{\xi})}$ can be represented by:
\small
\begin{equation}
\label{eq:lra_rank_R}
{Y}\approx{\hat{Y}}=\hat{g}(\bm{\xi})=\sum\limits_{l=1}^{r}{{{b}_{l}}{{\omega}_{l}}(\bm{\xi})}
\end{equation}
\normalsize
\noindent in which $b_{l},l=1,...,r$ are normalizing weighting factors, and ${{\omega}_{l}}$ is a rank-one function of ${\xi}$ in the form of
\small
\begin{equation}
\label{eq:lra_rank_one}
\omega_{l}(\bm{\xi})=\prod\limits_{i=1}^{n}{{{v}_{l}^{(i)}}({{\xi}_{i}})}
\end{equation}
\normalsize
\noindent where ${v_{l}^{(i)}}$ denotes the ${i}$-th dimensional univariate function in the ${l}$-th rank-one function.
For most applications, the number ${r}$ of rank-one terms is usually small (under 5), hence \eqref{eq:lra_rank_R} and \eqref{eq:lra_rank_one} represent a canonical low-rank approximation.

In order to obtain the rank-$r$ approximation, a natural choice is expanding ${v_{l}^{(i)}}$ on a polynomial basis ${\{\phi_{k}^{(i)}, k \in N\}}$ that is orthogonal to ${f_{X_{i}}}$, the resulting rank-$r$ approximation takes the form:
\small
\begin{equation}
\label{eq:lra_rank_r_pce}
{\hat{Y}}=\hat{g}(\bm{\xi})=\sum\limits_{l=1}^{r}{{{b}_{l}} \left[\prod\limits_{i=1}^{n}{\left(\sum\limits_{k=0}^{{{p}_{i}}}{z_{k,l}^{(i)}\phi_{k}^{(i)}({{\xi}_{i}})} \right)} \right]}
\end{equation}
\normalsize
where ${\phi}_{k}^{(i)}$ denotes the $k$-th degree univariate polynomial in the $i$-th random input, $p_{i}$ is the maximum degree of ${\phi}^{(i)}$ and ${z}_{k,l}^{(i)}$ is the coefficient of ${\phi}_{k}^{(i)}$ in the $l$-th rank-one function.

Building the low-rank approximation for desired response in \eqref{eq:lra_rank_r_pce} requires: (i) choose an appropriate univariate polynomial for each random input; (ii) solve the polynomial coefficients ${z_{k,l}^{(i)}}$ as well as the weighing factors ${b_{l}}$.
This process relies on a set of samples and their corresponding accurate response which are usually termed as experimental design (ED).

\subsection{Selection of the Univariate Polynomial Basis}

It is crucial to choose a proper polynomial ${\phi}_{i}$ for the $i$th random input ${\xi}_{i},i=1,2,...,n$ to avoid low convergence rate and/or higher degree of expansion \cite{DBXiu02a}, which will hamper the capability of LRA in dealing with high-dimensional problems.
Table \ref{tab:gpce_mapping} shows a set of typical continuous distributions and the respective optimal Wiener-Askey polynomial basis,
which can ensure the exponential convergence rate. In case that $f_{\xi_{i}}$ is out of the list in Table \ref{tab:gpce_mapping}, the discretized Stieltjes procedure is adopted to numerically construct a set of univariate orthogonal polynomial basis
\cite{UQLabPCE17}.
\begin{table}[t]
\renewcommand{\arraystretch}{1.3}
\caption{Standard Forms of Classical Continuous Distributions and their Corresponding Orthogonal Polynomials \cite{DBXiu02a}}
\label{tab:gpce_mapping}
\begin{center}
\begin{tabular}{c|c|>{\centering}p{1.5cm}|c}
\hline
\bfseries Distribution & \bfseries Density Function & \bfseries Polynomial & \bfseries Support \\
\hline
Normal & ${\frac{1}{\sqrt{2\pi}}{{e}^{-{{x}^{2}}/2}}}$ & Hermite & (-$\infty$,$\infty$) \\
\hline
Uniform & ${\frac{1}{2}}$ & Legendre & [-1,1] \\
\hline
Beta & ${\frac{{{(1-x)}^{\alpha}}{{(1+x)}^{\beta}}}{{{2}^{\alpha+\beta+1}}{B(\alpha+1,\beta+1)}}}$ & Jacobi & [-1,1] \\
\hline
Exponential & ${{e}^{-x}}$ & Laguerre & (0,$\infty$) \\
\hline
Gamma & ${\frac{{{x}^{\alpha}}{{e}^{-x}}}{\Gamma(\alpha+1)}}$ & Generalized Laguerre & [0,$\infty$) \\
\hline
\end{tabular}
\end{center}
\begin{tablenotes}
\item * The Beta function is defined as $B(p,q)=\frac{\Gamma(p)\Gamma(q)}{\Gamma(p+q)}$.
\end{tablenotes}
\end{table}

\subsection{Calculation of the Coefficients and Weighing Factors}\label{section_solve_lra}

Based on an experiment design of size $M_{C}$, i.e., a set of samples of $\bm{\xi}_{C}=\{\xi^{(1)},\xi^{(2)},...,\xi^{(M_{C})}\}$ and the corresponding response $\bm{y}_{C}=\{y^{(1)},y^{(2)},...,y^{(M_{C})}\}$ evaluated by deterministic tools, different algorithms have been proposed in the literature for solving the LRA coefficients and the weighting factors in a non-intrusive manner \cite{ADoostan13a}, \cite{PRai14a}.
The sequential correction-updating scheme presented in \cite{MChevreuil15a} is employed in this study due to its efficiency and capability of constructing low-rank approximation using less sample evaluations.
In the ${r}$-th correction step, a new rank-one function ${\omega_{r}}$ is built, while in the ${r}$-th updating step, the set of weighing factors ${\{b_{1},...,b_{r}\}}$ is determined. This process continues until the applied error index (the relative empirical error in this study) stop decreasing \cite{MChevreuil15a}.

\textbf{Correction step}:
the ${r}$-th correction step aims to find a new rank-one tensor $\omega_{r}$, which can be obtained by solving the following minimization problem:
\small
\begin{equation}
\label{eq:lra_rank_one_min}
\begin{aligned}
{{\omega}_{r}}(\bm{\xi}) &= \arg \underset{\omega \in W}{\mathop{\min}}\,\left\| {{e}_{r-1}}-\omega \right\|_{\bm{\xi}_{C}}^{2} \\
& = \arg \underset{\omega \in W}{\mathop{\min}}\,\sum_{m=1}^{M_{C}}{\left[y^{(m)}-\hat{g}_{r-1}(\bm{\xi}^{(m)})-\omega(\bm{\xi}^{(m)})\right]^{2}}
\end{aligned}
\end{equation}
\normalsize
\noindent where ${W}$ represents the space of rank-one tensors, ${{e}_{r-1}}=(g-{{\hat{g}}_{r-1}})$ is the approximation error of the response $Y$ at the ${(r-1)}$-th step, $\|.\|^{2}$ represents the norm 2 of the residual after the new rank-one tensor $w$ is applied, and the subscript ${\bm{\xi}_{C}}$ indicates that the minimization is carried over the whole set of samples in the experiment design $(\bm{\xi}_{C},\bm{y}_{C})$.

By exploiting the retained tensor-product form of the univariate polynomial basis, as shown in \eqref{eq:lra_rank_r_pce}, typical scheme for solving equation \eqref{eq:lra_rank_one_min} is the alternated least-square (ALS) minimization, which involves sequential minimization along each dimension ${i=1,...,n}$ to solve the corresponding polynomial coefficients ${\bm{z}_{r}^{(i)}=(z_{0,r}^{(i)},...,z_{p_{i},r}^{(i)})}$.
The total number of coefficients to be solved in each correction step is $\sum\nolimits_{i=1}^{n}{\left({{p}_{i}}+1\right)}$, which grows linearly as the number of random inputs $n$ increases.
Since ${{\omega }_{r}}$ is the product of $v_{l}^{(i)}(\xi)$ as shown in \eqref{eq:lra_rank_one}, $v_{l}^{(i)}({\xi}_{i})$ can be initialized by setting to $1.0$.

In the minimization along the ${i}$-th dimension, the polynomial coefficients corresponding to all other dimensions are "frozen" at their current values and the coefficients ${\bm{z}_{r}^{(i)}=(z_{0,r}^{(i)},...,z_{p_{i},r}^{(i)})}$ can be determined by:
\small
\begin{equation}
\label{eq:lra_coeff}
\bm{z}_{r}^{(i)}=\arg \underset{\bm{\zeta} \in {{R}^{({{p}_{i}}+1)}}}{\mathop{\min}}\,\left\|{{e}_{r-1}}-C_{i}\left( \sum\limits_{k=0}^{{{p}_{i}}}{{{\zeta}_{k}}\phi_{k}^{(i)}}\right)\right\|_{\bm{\xi}_{C}}^{2}
\end{equation}
\normalsize
\noindent where $C_{i}$ is a scalar
\small
\begin{equation}
\label{eq:lra_coeff_c}
C_{i}=\prod\limits_{j\ne i}{v_{r}^{(j)}(\xi_{j})}=\prod\limits_{j\ne i}{\left( \sum\limits_{k=0}^{{{p}_{j}}}{z_{k,r}^{\left( j \right)}\phi _{k}^{\left( j \right)}\left( {{\xi}_{j}} \right)} \right)}
\end{equation}
\normalsize

\textbf{Updating step}: After the ${r}$-th correction step is completed, the algorithm proceed to the ${r}$-th updating step to determine the weighing factor $b_{r}$ of the newly solved rank-one function ${{\omega }_{r}(\xi)}$, meanwhile, the set of existing weighing factors ${\bm{b}=(b_{1},...,b_{r-1})}$ are updated too. The updating step can be achieved by solving the following minimization problem
\small
\begin{equation}
\label{eq:lra_weight}
\bm{b}=\arg \underset{\bm{\beta} \in {{R}^{r}}}{\mathop{\min}}\,\left\|g-\sum\limits_{l=1}^{r}{{{\beta}_{l}}{{\omega}_{l}}}\right\|_{\bm{\xi}_{C}}^{2}
\end{equation}
\normalsize

\textbf{Stop criteria}: The correction-updating scheme successively adds new rank-one function to improve the accuracy of the approximation \eqref{eq:lra_rank_r_pce}. Hence, error reduction in two successive iterations becomes a natural stop criteria for this process. In this study, the relative empirical error is employed which is given by:
\small
\begin{equation}
\label{eq:lra_error}
{{\hat{e}}_{r}}=\frac{\left\|{{e}_{r-1}}-{{\omega}_{r}}\right\|_{\bm{\xi}_{C}}^{2}}{\mathbb{V}(\bm{y}_{C})}
\end{equation}
\normalsize
where $\mathbb{V}(\bm{y}_{C})$ denotes the empirical variance of the desired response over the experimental design.

Remark: (i) It is worth pointing out that the minimization problems in (\ref{eq:lra_coeff}) and (\ref{eq:lra_weight}) can be efficiently solved with the ordinary least-squares (OLS) technique because the dimension of unknowns are small.
(ii) When the LRA \eqref{eq:lra_rank_r_pce} for the desired response is built up, the response of any new samples can be evaluated efficiently by directly substituting to \eqref{eq:lra_rank_r_pce} instead of solving the original complex problem (e.g., the PATC problem \eqref{eq:atc_equation}).

\subsection{Selection of optimal rank and polynomial degree}
Currently, there is no systematic way to identify the optimal rank $r$ and the polynomial degree $p_i$ for individual random input beforehand. In this study, we specify a candidate set ${\{1,2,3,4,5\}}$ for the rank $r$, and another one ${\{2,3,4,5\}}$ for all the univariate polynomials. The rank selection is performed by progressively increasing the rank and applying the corrected error \eqref{eq:lra_error} to select the best one. The selection of optimal degree can be implemented similarly.

\subsection{Integration of Dependent Random Inputs}

So far, the random inputs are assumed to be mutually independent as requested by the LRA method.
To accommodate dependent random inputs with correlation matrix $\bm{\rho}$, the Nataf transformation \cite{RLebrun09a}, \cite{RLebrun09b} and the isoprobabilistic transformation can be employed to build up a mapping between $\bm{U}$ and the independent standard random variables $\bm{\xi}$: $\bm{u}=T(\bm{\xi})$, where $T$ is invertible. Therefore, the set of samples of ${\bm{\xi}}$ can be transformed back into samples of ${\bm{U}}$ to evaluate the corresponding responses ${\bm{Y}}$, after which the desired response $\bm{y}=g(T^{-1}(\bm{u}))$ can be expanded onto the polynomial basis with $\bm{\xi}$ using the aforementioned method \cite{KKonakli16a}.

\subsection{Moments of a Low-Rank Approximation}

Due to the orthogonality of the univariate polynomials that form the LRA basis (see (\ref{eq:lra_rank_r_pce})), the mean and the variance of the meta-model can be obtained \textit{analytically} in terms of the polynomial coefficients and the weighting factors. In particular, the mean and variance of the LRA response are given respectively by \cite{KKonakli16a}:
\small
\begin{equation}
\label{eq:lra_mean_var}
\begin{gathered}
{{\mu}_{y}}=E\left[{\hat{g}}(\bm{\xi})\right]=\sum\limits_{l=1}^{r}{{{b}_{l}}\left(\prod\limits_{i=1}^{n}{z_{0,l}^{(i)}}\right)} \\
\sigma_{y}^{2}=\sum\limits_{l=1}^{r}{\sum\limits_{m=1}^{r}{{{b}_{l}}{{b}_{m}}\prod\limits_{i=1}^{n}{\left[\left(\sum\limits_{k=0}^{{{p}_{i}}}{z_{k,l}^{(i)}z_{k,m}^{(i)}}\right)-z_{0,l}^{(i)}z_{0,m}^{(i)} \right]}}}
\end{gathered}
\end{equation}
\normalsize
Hence, if only mean and variance are of interests, \eqref{eq:lra_mean_var} can be applied directly without evaluating a large size of samples, which in contrast is required by most of the simulation-based methods.

\section{Computation of Probabilistic Available Transfer Capability}

In this section, a step-by-step description of the LRA method for PATC calculation is summarized below:

\noindent \textbf{Step 1}: Input the network data, the transactions to be studied, the credible contingency list, the probability distribution and the parameters of the random inputs $\bm{U}$, i.e., the wind speed, the solar radiation, the active load power, and their correlation matrix ${\bm{\rho}}$.

\noindent \textbf{Step 2}: Build the load-generation pattern (i.e., vector $b$ in \eqref{eq:atc_equation}-(a)) according to the transactions under study.

\noindent \textbf{Step 3}: Choose the independent standard variable ${\xi_{i}}$ and the corresponding univariate polynomial ${\phi_{i}}$ for each random input $U_{i}$.

\noindent \textbf{Step 4}: Generate an experimental design of size ${{M}_{C}}$:
\begin{itemize}
\item i) Generate ${M_{C}}$ samples ${\bm{\xi}_{C}=(\xi^{(1)},\xi^{(2)},...,\xi^{({M_{C}})})}$ in the standard space by the LHS method.
\item ii) Transform ${\bm{\xi}_{C}}$ into the physical space by the inverse Nataf transformation ${\bm{u}_C}=T^{-1}(\bm{\xi}_C)$.
\item iii) Apply the deterministic ATC solver to evaluate the accurate response ${\bm{y}_{C}}$ (i.e., PATC) of ${\bm{u}_{C}}$.
For each realization, the overall ATC is computed by applying the parallel scheme \cite{GCEjebe98a} which enables an early stop in CPF calculation when the first violation is encountered in either normal or contingency cases.
Pass the sample-response pairs $({\bm{\xi}_{C}},{\bm{y}_{C}})$ to Step 5.
\end{itemize}

\noindent \textbf{Step 5}: Apply the algorithm in Section \ref{section_solve_lra} to build the low-rank approximation (\ref{eq:lra_rank_r_pce}) for the PATC.
If the LRA for PATC has reached the prescribed accuracy, go to Step 7; otherwise, go to Step 6.

\noindent \textbf{Step 6}: Generate additional ${\Delta M_{C}}$ new samples and evaluate them, then go back to Step 5 using the enriched experiment design $({\bm{\xi}_{C}}+{\Delta \bm{\xi}_{C}},{\bm{y}_{C}}+{\Delta \bm{y}_{C}})$.

\noindent \textbf{Step 7}: Calculate the mean and standard deviation of the response through \eqref{eq:lra_mean_var}.

\noindent \textbf{Step 8}: Sample $\bm{\xi}$ extensively, e.g., ${{M}_{S}}$ samples, and apply the solved functional approximation (\ref{eq:lra_rank_r_pce}) to evaluate the corresponding response $\bm{y}_{S}$ for all samples. Then compute the statistics of interest (PDF/CDF of PATC in this study).

\noindent \textbf{Step 9}: Compute the TRM value and the resulting ATC ($ATC=\mu_{PATC} - TRM$) for the given confidence level $p_{cl}\%$.

\noindent \textbf{Step 10}. Generate the result report.

Remark: the number of samples ${{M}_{C}}$ in Step 4 is usually much smaller than ${{M}_{S}}$ in Step 8. Unlike MCS, LRA does not solve power flow equations for all ${{M}_{S}}$ samples in Step 8, and hence it is more efficient. The main computational effort of LRA lies in Step 4.

\section{Numerical Studies}

In this section, we apply the proposed LRA method to investigate the probabilistic ATC of the modified IEEE 24-bus reliability test system (RTS)
and IEEE 118-bus systems \cite{MATPOWER}. The LHS-based Monte Carlos simulation serves as a benchmark for validating the accuracy and the performance of the proposed method.
In addition, a comparison between the LRA and the sparse PCE \cite{HSheng18a} is also presented.

In this study, we assume that the probability distributions and the associated parameters of all random inputs are available from up-front modeling. Particularly, the wind speed follows Weibull distribution; the solar radiation follows Beta distribution; the load power follows Normal distribution. For each individual load, the mean value is set to be its base case value and the variance is equal to 5\% of its mean value.
The univariate polynomial basis used in \eqref{eq:lra_rank_r_pce} for Weibull, Beta and Normal distributions are chosen appropriately
according to Table \ref{tab:gpce_mapping}.
For simplicity, the linear correlation coefficient ${\rho_{ij}}$ between component $i$ and $j$ of wind speed $\bm{v}$, solar radiations $\bm{r}$ and load power $\bm{P_{L}}$ are 0.8040, 0.5053 and 0.4000, respectively.

\subsection{The Modified IEEE 24-Bus RTS}

The IEEE 24-bus RTS is composed of 4 areas, containing 33 generators, 32 branches and 17 loads.
In this paper, we modify this test system by adding 4 wind farms of 80 MW at bus \{15, 18, 21, 23\} and 4 solar PV power plants of 60 MW at bus \{1, 2, 7, 16\}, respectively. Their corresponding parameters $\{W1, W2, W3, W4\}$ and $\{S1, S2, S3, S4\}$ are shown in Table \ref{tab:wind_param}-\ref{tab:solar_param}. Together with 17 stochastic loads, there are totally 25 random inputs.
The existing transmission commitments (base case) is the first-day peak load as defined in \cite{PFAlbrecht79a}, and the new transaction under study is to transfer 75 MW from generator bus \{7\} in area 2 (source) to load buses \{3,4,9\} in area 1 (sink). The contingency list under study contains four $N-1$ outages $\{G1\#1, L2-4, L3-24, L9-11\}$.

\begin{table}[]
\renewcommand{\arraystretch}{1.3}
\caption{Wind speed and wind turbine parameters \cite{SHJangamshetti99a}}
\label{tab:wind_param}
\centering
\begin{tabular}{c|c|c|c|c|c}
\hline
$No.$ & ${c}$ & ${k}$ & ${V_{r}}$ & ${V_{in}}$ & ${V_{out}}$ \\
\hline
${W1}$ &  8.0063 & 2.1182 & 13.50 & 3.50 & 25.00 \\
\hline
${W2}$ & 11.5762 & 2.7022 & 13.80 & 3.50 & 25.00 \\
\hline
${W3}$ & 11.2441 & 3.6322 & 13.00 & 5.00 & 25.00 \\
\hline
${W4}$ & 12.4813 & 3.2465 & 12.90 & 5.00 & 24.00 \\
\hline
${W5}$ & 11.1533 & 3.2895 & 12.00 & 5.50 & 24.00 \\
\hline
${W6}$ &  8.8261 & 2.6511 & 10.00 & 3.50 & 20.00 \\
\hline
\end{tabular
\end{table}
\begin{table}[]
\renewcommand{\arraystretch}{1.3}
\caption{Solar radiation and solar PV parameters \cite{FYEttoumi02a}}
\label{tab:solar_param}
\centering
\begin{tabular}{c|c|c|c|c|c|c}
\hline
$No.$ & ${\alpha}$ & ${\beta}$ & $r_{min}$ & $r_{max}$ & ${R_{c}}$ & ${R_{std}}$ \\
\hline
${S1}$ & 1.110 & 0.730 & 0.0 & 1000.0 & 150.0 & 1000.0 \\
\hline
${S2}$ & 1.320 & 0.690 & 0.0 & 1000.0 & 150.0 & 1000.0 \\
\hline
${S3}$ & 1.700 & 0.740 & 0.0 & 1000.0 & 150.0 & 1000.0 \\
\hline
${S4}$ & 2.970 & 0.940 & 0.0 & 1000.0 & 150.0 & 1000.0 \\
\hline
${S5}$ & 2.540 & 0.780 & 0.0 & 1000.0 & 150.0 & 1000.0 \\
\hline
${S6}$ & 1.780 & 0.850 & 0.0 & 1000.0 & 150.0 & 1000.0 \\
\hline
\end{tabular
\end{table}

We first apply the CPF-based ATC solver to the deterministic system (i.e., without uncertainty). As shown in Table \ref{tab:deter_atc}, the binding limit is the thermal limit at branch 7-8 in the base case, leading to an overall ATC 83.0207 MW. The contingencies do not deteriorate the transfer capability potentially because the given limits of branches are larger in emergency compared to the normal operating condition.
\begin{table}[]
\renewcommand{\arraystretch}{1.3}
\caption{The determination of ATC considering normal and contingency cases}
\label{tab:deter_atc}
\centering
\begin{tabular}{c|c|>{\centering\arraybackslash}m{1.2cm}|>{\centering\arraybackslash}m{1.2cm}|>{\centering\arraybackslash}m{1.2cm}|>{\centering\arraybackslash}m{0.8cm}}
\hline
\multirow{2}{*}{Case} & \multirow{2}{*}{Outage} & \multicolumn{3}{c|}{ATC w.r.t limits (MW)} & \multirow{2}{*}{Overall} \\
\cline{3-5}
No. & Facility & Voltage limit & Thermal limit & Voltage collapse & ATC\\
\hline
${0}$ & Base case & 280.3934 & 83.0207 & 511.9046 & \multirow{5}{*}{83.0207} \\
\cline{1-5}
${1}$ & G1-1 & 469.0375 & 128.0214 & 508.8632 & \\
\cline{1-5}
${2}$ & L2-4 & 201.7742 & 127.8677 & 451.6651 & \\
\cline{1-5}
${3}$ & L3-24 & 183.7445 & 127.5457 & 355.2315 & \\
\cline{1-5}
${4}$ & L9-11 & 395.0070 & 127.8440 & 460.6714 & \\
\hline
\end{tabular
\end{table}

Next, we exploit the proposed LRA method to assess the probabilistic characteristics (e.g., mean, variance, PDF and CDF) of the PATC and compare the results with those of the LHS-based MCS and with those of the PCE method.
Applying the proposed algorithm, 125 simulations are required in Step 4-5 to build up the LRAs \eqref{eq:lra_rank_r_pce} of the PATC, which consist of 1 rank-one function with the highest polynomial degrees $p_{i}=2$. The total number of coefficients plus weighting factors is 76 (i.e., $\sum\nolimits_{l=1}^{r}{\left(\sum\nolimits_{i=1}^{n}(p_{i}+1)+1\right)}=\sum\nolimits_{i=1}^{25}(2+1)+1$). Once the coefficients and the weighting factors are computed, the mean and the standard deviation of the responses are computed in Step 7 and are compared with those of the LHS-based MCS and with those of the sparse PCE method as shown in Table \ref{tab:comp_atc_24}. Furthermore, 2000 samples are generated in Step 8 to assess the PDF and the CDF of the response. Fig. \ref{fig:total_atc_dist} shows the PDF and CDF of the PATC computed by the LHS-based MCS, the sparse PCE, and the solved LRA, respectively. These results clearly demonstrate that the LRA-PATC method can provide accurate estimation for the probabilistic characteristics of the PATC.
\begin{figure}[h]
\centering
\includegraphics[width=0.45\textwidth,height=0.18\textheight]{figure_01.eps}
\caption{The distribution of the PATC computed by the MCS, the PCE and the LRA. They are almost overlapped. The TRM for 95\% confidence level is 21.1497 MW and the resulting ATC is 61.8815 MW.}
\label{fig:total_atc_dist}
\end{figure}
\begin{table}[]
\renewcommand{\arraystretch}{1.3}
\caption{Comparison of the estimated statistics of the overall ATC by the MCS, PCE and LRA methods}
\label{tab:comp_atc_24}
\centering
\begin{tabular}{c|c|c|c|c|c}
\hline
Indices & MCS & PCE & LRA & ${\frac{\Delta PCE}{MCS}} \%$ & ${\frac{\Delta LRA}{MCS}} \%$ \\
\hline
$\mu$ & 83.0312 & 83.2369 & 83.2226 & 0.2478 & 0.2305 \\
\hline
$\sigma$ & 15.5418 & 14.4898 & 15.4278 & -6.7692 & -0.7340 \\
\hline
\end{tabular
\end{table}

Once we have the statistics of the PATC, a reasonable amount of TRM can be obtained from the CDF of PATC. Table \ref{tab:confid_levels_24} shows the TRM with different confidence levels and the resulting ATC values. For example, if 95\% confidence level is requested, i.e., $P(ATC_{actual} >= (\mu_{ATC} - TRM)) = 0.95$, the corresponding ATC is 61.8815 MW.
\begin{table}[]
\renewcommand{\arraystretch}{1.3}
\caption{The estimated TRM and resulting ATC for different confidence levels}
\label{tab:confid_levels_24}
\centering
\begin{tabular}{>{\centering\arraybackslash}m{1.6cm}|>{\centering\arraybackslash}m{1.5cm}|>{\centering\arraybackslash}m{1.5cm}|>{\centering\arraybackslash}m{1.5cm}}
\hline
Confid. Level & $\mathbb{E}($PATC$)$ & TRM (MW) & ATC (MW) \\
\hline
${99.0\%}$ & 83.0312 & 26.2326 & 56.7986 \\
\hline
${98.0\%}$ & 83.0312 & 24.4386 & 58.5926 \\
\hline
${95.0\%}$ & 83.0312 & 21.1497 & 61.8815 \\
\hline
${90.0\%}$ & 83.0312 & 17.9604 & 65.0707 \\
\hline
${80.0\%}$ & 83.0312 & 12.9772 & 70.0540 \\
\hline
\end{tabular
\end{table}

\subsection{The Modified 118-bus System}
The IEEE-118 bus system is a simplified representation of the Midwest U.S. transmission system in 1962, which contains 19 generators, 35 synchronous condensers, 177 transmission lines, 9 transformers and 91 loads.
Six wind farms, each with $50$ MW, are connected to bus $\{10, 25, 26, 49, 65, 66\}$ using the parameters $\{W1, W2, W3, W4, W5, W6\}$ in Table \ref{tab:wind_param} respectively. Six solar PV parks, each with an installed capacity of $30$ MW, are connected to bus $\{12, 59, 61, 80, 89, 100\}$ using parameters $\{S1, S2, S3, S4, S5, S6\}$ in Table \ref{tab:solar_param}. Besides, there are 99 stochastic loads, leading to 111 random inputs in total.
The total load in the base case is 4242 MW and 1438 Mvar. The new transaction under study is to transfer 60 MW from bus 89 to bus 91. The contingency list contains six $N-1$ outages $\{L88-89, L89-90, L90-91, L89-92, L91-92, L92-94\}$ surrounding the source and sink buses.

Likewise, we first apply the CPF-based ATC solver to compute the maximum transfer without considering the uncertainties of RES and loads. The deterministic ATC is 24.1835 MW. Next, we exploit the proposed LRA-PATC method to assess the probabilistic characteristics of the PATC and compare the results with those of the LHS-based MCS and the sparse PCE method.

Applying the proposed algorithm, 556 simulations are needed in Step 4-5 to build up the LRAs \eqref{eq:lra_rank_r_pce} of the PATC, which consist of 1 rank-one function with the highest degree $p_{i}=2$. The total number of coefficients plus weighting factors is 334.
Eventually, the mean and standard deviation of the responses are computed in Step 7 and are compared with those from the LHS-based MCS and from the sparse PCE method as shown in Table \ref{tab:comp_atc_118}.
All these results and comparisons clearly demonstrate that the LRA method can provide accurate estimation for the probabilistic characteristics of the PATC solutions.
However, to get comparable accuracy, the LHS-based MCS needs to run 10000 simulations (i.e., solving (\ref{eq:atc_equation})), taking 9110s, while the LRA takes only 480s, i.e., about $\frac{1}{19}$ of the computational time required by MCS. Clearly, the LRA is much more efficient than the MCS.
\begin{table}[]
\renewcommand{\arraystretch}{1.3}
\caption{Comparison of the estimated statistics of the overall ATC by the MCS, PCE and LRA methods}
\label{tab:comp_atc_118}
\centering
\begin{tabular}{c|c|c|c|c|c}
\hline
Indices & MCS & PCE & LRA & ${\frac{\Delta PCE}{MCS}} \%$ & ${\frac{\Delta LRA}{MCS}} \%$ \\
\hline
$\mu$ & 23.9278 & 23.7476 & 23.7550 & -0.7530 & -0.7220 \\
\hline
$\sigma$ & 4.8779 & 4.5220 & 4.8617 & -7.2973 & -0.3327 \\
\hline
\end{tabular
\end{table}

The TRM and the resulting overall ATC for different confidence levels are shown in Table \ref{tab:confid_levels_118}, which empowers the system operator to determine a reasonable amount of TRM to make full use of the transmission assets.
\begin{table}[]
\renewcommand{\arraystretch}{1.3}
\caption{The estimated TRM and resulting ATC for different confidence levels}
\label{tab:confid_levels_118}
\centering
\begin{tabular}{>{\centering\arraybackslash}m{1.6cm}|>{\centering\arraybackslash}m{1.5cm}|>{\centering\arraybackslash}m{1.5cm}|>{\centering\arraybackslash}m{1.5cm}}
\hline
Confid. Level & $\mathbb{E}($PATC$)$ & TRM (MW) & ATC (MW) \\
\hline
${99.0\%}$ & 23.9278 & 10.1839 & 13.7439 \\
\hline
${98.0\%}$ & 23.9278 & 8.9942 & 14.9336 \\
\hline
${95.0\%}$ & 23.9278 & 7.4975 & 16.4303 \\
\hline
${90.0\%}$ & 23.9278 & 6.1159 & 17.8119 \\
\hline
${80.0\%}$ & 23.9278 & 4.3121 & 19.6157 \\
\hline
\end{tabular
\end{table}

\section{Conclusion and Perspectives}

In this paper, we have proposed a mathematical formulation for the probabilistic ATC (PATC) problem in which the uncertainties from RES and loads are incorporated. Moreover, we have proposed a novel LRA method which can assess the PATC accurately and efficiently by building up a statistically-equivalent low-rank approximation.
Numerical studies show that the proposed method can accurately estimate the probabilistic characteristics of the PATC with much less computational effort compared to the LHS-based MCS.

The PATC provides important insights into how the uncertainties may affect the transfer capability of the power system. More importantly, the proposed method can help determine the TRM and thus the ATC in a more efficient manner to fully utilize the transmission assets while maintaining the security and reliability of the grid.

In the future, we plan to develop control measures to reduce the variance of PATC to decrease the TRM for a full utilization of transmission assets.


\section{Introduction}
\label{sec:intro}

High dynamic range (HDR) images, which provide wider range of luminance variation draw more attention in recent years.
But it is still not accessible to most people who use traditional standard dynamic range (SDR) devices.
Tone mapping operators and multi-exposure fusion methods~\cite{HDRbook} provide alternative approaches to access HDR contents with SDR devices.
The resulting images, called tone-mapped HDR images, can preserve the details and the contrast of HDR contents in a narrower dynamic range. But these methods also bring complex distortions into the HDR contents in the practical applications~\cite{Attributes}.
Therefore, proper image quality assessment (IQA) methods are necessary.
Ideally, the best solution is subjective assessment carried by human beings, but subjective assessment is complicated and time-consuming.
Thus proper objective assessment methods are needed to predict tone-mapped image quality automatically.
Furthermore, the absence of reference HDR images calls for no-reference image quality assessment (NR-IQA) methods, since tone-mapped images are designed for situations where HDR devices are not accessible.

\begin{figure}
\centering
\subfigure[MOS=18.3126018946] { \label{fig:a}
\includegraphics[width=0.45\columnwidth]{Fig2Sub1.PNG}
}
\subfigure[MOS=51.3429624349] { \label{fig:b}
\includegraphics[width=0.45\columnwidth]{Fig2Sub2.PNG}
}
\subfigure[MOS=54.4365059326] { \label{fig:c}
\includegraphics[width=0.45\columnwidth]{Fig1Sub1.PNG}
}
\subfigure[MOS=23.6067858845] { \label{fig:d}
\includegraphics[width=0.45\columnwidth]{Fig1Sub2.PNG}
}
\caption{Four images from ESPL-LIVE HDR Database. Higher MOS indicates better subjective image quality.}
\label{fig:0}
\end{figure}

However, the aim to preserve HDR contents on SDR devices also leads to the complex standard in the quality assessment of tone-mapped HDR images.
Normally, basic attributes like brightness and contrast are important to IQA because they can be used to measure the preservation of details and structures.
Images with moderate brightness and higher local contrast usually have more information and better image quality.
It is shown in the comparison between Fig.~\ref{fig:0} (a) and Fig.~\ref{fig:0} (b).
Focusing on the details and contrast preservation, existing objective assessment methods mainly take structural similarity and natural scene statistics as tools of tone-mapped image quality assessment (TMIQA).
But due to the complexity of TMIQA,
more factors should be considered.
Decided by the aim to preserve HDR contents, tone-mapped images should have more information and higher contrast than general images.
But then they can never preserve all the information from the original HDR images because of the reduction of dynamic range.
The proper level of preservation depends on the overall scene and the contents in the image.
For example, we can see that compared with Fig.~\ref{fig:0} (c), Fig.~\ref{fig:0} (d) preserves more details and color information in the sky and in the shadow of the woods.
But preserving higher local contrast also breaks the balance
between the sky and the woods and affects the recognition of the scene.
Without high level information and global information,
the image quality cannot be predicted using only low-level information and local information.

In this paper we propose a new no-reference image quality assessment (NR-IQA) method to evaluate tone-mapped image quality which makes use of multi-scale and multi-layer information, adapted to the complex standard of TMIQA.
A multi-scale representation is adopted to combine local and global information.
Multi-Layer features are extracted from different layers of a pre-trained deep convolutional neural network (DCNN) model, to give consideration to both low and high level information.
Then these informative features are mapped to quality prediction by regression.
We conduct experiments on the largest public available tone-mapped image database. The experimental results prove that our proposed method outperforms existing NR-IQA methods on quality assessment for tone-mapped images.

The paper is organized as follow. Section 2 reviews related research. Section 3 introduces the ESPL-LIVE HDR database. The proposed method is described in Section 4. Experimental results are given in Section 5, and lastly conclusions are made in Section 6.

\section{Related Work}

The preservation of information in extra dynamic range, the balance between local and global contrast, and the natural looking are important aspects of tone-mapped image quality.
Therefore, existing methods mainly access tone-mapped image quality from detail preservation, contrast preservation, and image naturalness.
Compared to the employment of HDR and tone mapping operators, research history of TMIQA is relatively short.
As far as we concerned, only a few objective assessment methods are proposed.
According to the availability and usage of original HDR images, they can be further classified as full-reference (FR) methods and no-reference (NR) methods.

\subsection{Full-Reference Methods}
 Due to the importance of detail preservation, many researchers naturally take structural similarity between the original HDR images and the tone-mapped ones as a measure of tone-mapped image quality.
 Aydin et al.~\cite{08SIGGRAPH} first break the limitation of different dynamic range with proper transformation and measure the similarity between HDR and reproduced SDR images.
 Yeganeh et al.~\cite{13TIP} propose the popular FR framework of Tone-Mapped image Quality Index (TMQI), which contains two components: multi-scale local structural similarity to evaluate the information preservation, and statistical naturalness to evaluate the natural looking of tone-mapped images.
 Many modified methods are proposed on the basis of TMQI.
 Nasrinpour et al.~\cite{15ICIP} use visual saliency, while Xie et al.~\cite{16ISCAS} choose sparse domain instead, to modify the measurement of local structural similarity;
 Kundu et al.~\cite{16ICIP} use the statistical features of mean subtracted contrast normalized (MSCN) pixels to measure the natural statistics, and choose a pooling strategy from standard deviation of MSCN pixels, visual saliency and local information entropy for the pooling of local structural similarity.

\subsection{No-Reference Methods}
 Researches on the contributions of various image attributes to tone-mapped image quality are carried by {\v{C}}ad{\'\i}k et al.~\cite{Attributes} and show interesting results that inherent attributes such as brightness, contrast and color contribute even more to tone-mapped image quality than the information preservation.
 Thus TMIQA can be done using these attributes without information from the reference images.

 Since structural similarity is not applicable, no-reference methods of TMIQA focus on the design and extraction of various features to deal with the complex distortions in tone-mapped images, and use machine learning methods to assess the quality.
 For example, Kundu et al.~\cite{17TIP}
 extract 40 statistical features such as mean, standard deviation, skewness and kurtosis from spatial domain, gradient domain, standard deviation field of MSCN pixels and 6 types of log-derivative domains, and implement $\epsilon$-SVR  on these features to predict the quality of tone-mapped images.
 Yue et al.~\cite{17TIE} adopt a human visual system inspired method, where tone-mapped images are decomposed into four pairs of opponent color channels, and the features of local structures in the simulated responses of Single-Opponent/Double-Opponent cells are used to train models for TMIQA.

 However, the structural information in the tone-mapped images is still accessible without reference HDR images.
 Inspired by the goal to preserve details in the extra dynamic range,
 Gu et al.~\cite{16TMM} replace structural similarity component in the TMQI framework with the mean of Sobel map of each image,

 while the information entropy of the brightened and darkened versions of the image is used to quantify the volume of details.

Existing FR-TMIQA methods based on TMQI make use of multi-scale structural features but neglect the high-level information in the tone-mapped images.
Existing NR-TMIQA methods focus on the design of different features, but the diversity of information is limited by the number of features.
Unlike existing methods, our method use the pre-trained DCNN model to extract different levels of features for quality assessment.
Furthermore both local and global features are considered by using multi-scale representations.
These steps adapt our method to the characteristic of TMIQA.

\begin{figure*}[htb]
\centering\includegraphics[width=5.8in]{nnn2}
\caption{The framework of proposed method. Features with necessary information are guaranteed mainly by multi-scale representation and multi-layer feature aggregation.}\label{fig:1}
\end{figure*}

\section{Database}

 To investigate the best proposal of our method, and to test its performance, we take the newest public available tone-mapped HDR image database, ESPL-LIVE HDR Database~\cite{ESPL} for experiments.
 So far it is the largest database for TMIQA,
 and it covers a greater diversity of methods to obtain tone-mapped images:

 images in the database are derived from HDR illuminance maps and SDR image stacks through 11 different methods.
 Therefore it is suitable for experiments of our method, which is aimed at predicting tone-mapped HDR image quality affected by complex factors.

 There are 1811 tone-mapped HDR images in total,
 and the images fall into three types by the ways they are obtained:
 747 images are created

 through 4 different tone mapping operators;
 710 images are created directly from SDR image stacks through 5 multi-exposure fusion methods;
 354 images are obtained by Photomatix with 2 different post-processing settings.
 The images have been evaluated by 5462 participants on a crowdsourcing platform, and the subjective sores are provided in the form of mean opinion scores (MOS) ranging from 0 to 100.

 We randomly split the data into training and testing sets at 4:1 ratio in the objective experiments which will be introduced in  Section~\ref{sec:exp}.
 20\% of the training sets are further taken as validation sets in the validation experiments which will be introduced in Section~\ref{sec:method}, to investigate proposals of our method.
 Care was taken to ensure the data independency between training, validation and testing sets.

\section{Proposed Method}
\label{sec:method}
As stated in Section~\ref{sec:intro}, both low-level and high-level information, as well as global and local information should be considered in TMIQA.
To describe information of this large span with designed features is difficult.
Therefore we decide to extract features
with the aid of the DCNN model.
Databases for IQA are relatively small owing to the expense of subjective assessment.
It is difficult to train a model deep enough to extract high-level features from scratch.
Inspired by the work of Li et al.~\cite{exploiting}, we choose the ResNet-50 model~\cite{resnet} trained on ImageNet as the pre-trained model for feature extraction.
On one hand, due to the variety in image quality and the large number of images in ImageNet, the model which is pre-trained on ImageNet is a good choice to extract features for IQA.
On the other hand, ResNet-50 is deep enough to extract high-level features we need, while preserving some primary information related to image quality in residual images.

To give proper quality assessment to tone-mapped HDR images, the proposed method consists of four parts: multi-scale image representation, feature extraction, multi-layer feature aggregation, and quality prediction.
\textcolor{black}
{Firstly an image is represented with its original and down-sampled version and fed into the modified ResNet-50,
to preserve both local and global information.
Then feature maps with information of different levels are extracted from selected layers of the pre-trained model.
After being extracted, feature maps are aggregated through mean and standard deviation pooling in each channel, and concatenated according to the information carried, to form a feature structure to clearly preserve the information for the final mapping.}
Finally, the multi-scale and multi-layer features are mapped to one single quality score by partial least squares regression (PLSR)~\cite{PLSR}.

The whole framework is shown in Fig.~\ref{fig:1}. Details of multi-scale representation and multi-layer feature aggregation will be discussed below.

\subsection{Multi-scale image representation and feature extraction}

One important role of the tone-mapped image is to preserve the contrast and the details in dark and highlight areas.
Therefore both local information that can be used to evaluate the detail preservation, and global information that can be used to evaluate the contrast preservation, are important to TMIQA.
Features with local information can be easily extracted from the pre-trained DCNN model, since the input size of ResNet-50 is smaller than the size of images in ESPL-LIVE HDR Database.
But most of the global information may be lost, if the image is cropped to fit the model.
Therefore we remove the fully-connected layer and consequently remove the limitation of fixed input size.
And feature maps will be extracted when the size of the input image is larger than the original input size of the pre-trained model.
However, these feature maps still contains little global information without fully-connected layer.
Because the receptive field of one neuron in ResNet-50 is at most  $32 \times 32$ after modified and that is much smaller than the image size.
As a result, one element in the output feature maps is mainly decided by a small area.
To preserve global information within this fixed receptive field, we take the down-sampled version of an image as one component of the image representation, to allow larger areas and more information being involved in the first place.
The down-sampling operation is done along both height and width with a factor of 2.
Furthermore, to avoid the impact of down-sampling operation on image quality and to make use of local information, the original image is also preserved and fed into the DCNN model.

Feature maps extracted from images of different scales will be in different size.
For example, an image with size of $960 \times 540 \times 3$ is fed into the modified ResNet-50, the feature maps extracted from layer $res4f$ will be $30 \times 16 \times 1024$, while its down-sampled version gets feature maps sized $15 \times 8 \times 1024$.
We implement mean and standard deviation pooling on each channel of the feature maps, and the results from all the channels are taken as the representative statistical features of each scale.
After that, features of the two scales share the same size of $2 \times 1 \times 1024$.
They are then simply concatenated to form multi-scale features to preserve the full and clear structure of the multi-scale representation.
After down-sampling, the multi-scale features of an image can be described as $$O \oplus D,$$  where $O$ is features from the original image, $D$ is the features from the image after down-sampling operation, $\oplus$ is the concatenation operator.

To validate the gain of multi-scale representation, we compare the performance of the multi-scale representation with that of the single-scale representation where only the original image is involved.
We take features from layer $res2c, res3d$ and $res4f$, as the representative features of different levels,
and then test the performance of different representations.
The average performance of the two representations is also reported.
The comparison is done on the validation sets from ESPL-LIVE HDR Database~\cite{ESPL}, and the results are reported in Table~\ref{tab:multi-scale}.
It can be seen that the multi-scale representation achieves better performance.
Therefore, in our method, we choose the multi-scale representation to preserve both local and global information of the tone-mapped HDR images.

\begin{table}[!htb]
\begin{center}
\caption{Comparison between the performance of the multi-scale representation and the single-scale representation. Medians of SROCC in 1000 runs are reported.} \label{tab:multi-scale}
\vspace{2mm}
\resizebox{0.9\columnwidth}{!}
{
\begin{tabular}{ccccc}
\toprule
Features & res2c & res3d & res4f & overall \\
\midrule  $O$ & 0.709 & 0.763 & 0.787 & 0.753 \\
$O \oplus D$ & \textbf{0.712} & \textbf{0.767} & \textbf{0.796} & \textbf{0.758} \\
\bottomrule
\end{tabular}
}
\end{center}
\end{table}

\subsection{Multi-layer feature aggregation}
\label{subsec:mulilayer}
 As stated in Section~\ref{sec:intro}, tone-mapped images should preserve more details and contrast while not affecting the contents in the original images.
 When the reference images are not available, it is difficult to estimate the proper level of detail and contrast preservation without high-level information, while basic attributes such as brightness and contrast are important as well.
 Therefore it can be expected that tone-mapped image quality can be better predicted with features extracted from different layers, i.e. multi-layer features.

 To construct these multi-layer features, we combine features from different layers of ResNet-50.
 Based on the core algorithm of residual learning, the main part of ResNet-50 can be divided into 16 building blocks.
 Each building block consists of multiple layers and only the ending layer outputs the result of residual learning.
 Thus we give consideration to these 16 ending layers of the blocks as the layers to extract features from.
 Layers after $res4f$ are ignored, which are aimed at image classification but may have less relation to TMIQA.
 The rest 13 layers are further divided into four sets according to the depth: \{$res2a, res2b, res2c$\}, \{$res3a, res3b, res3c, res3d$\}, \{$res4a, res4b, res4c$\}, \{$res4d, res4e, res4f$\}.
 Three layers are chosen to extract low-level, mid-level and high-level features, where at most one layer is chosen for one set in order to achieve a balance of low dimensions and diversity of features.
 Features from the three selected layers are then concatenated as multi-layer features in the form of $$f_l \oplus f_m \oplus f_h,$$ where $f_l$, $f_m$, $f_h$ are the low-level, the mid-level and the high-level features separately.
 And $\oplus$ is the concatenation operator, adopted to preserve the clear structure of the multi-layer representation.

\begin{table}[!tb]
\begin{center}
\caption{Comparison between the performance of the single layer features and the multi-layer features, where TM stands for Tone Mapping, MEF stands for Multi-Exposure Fusion, and PP stands for Post Processing. Medians of SROCC in 1000 runs are reported.}
\label{tab:multi-layer}
{\begin{tabular}{ccccc}
\toprule
Features & Overall & TM & MEF & PP   \\
\midrule
$f_l$ (res2a) & 0.712 & 0.786 & 0.654 & 0.576\\
$f_m$ (res4b) & 0.800 & 0.845 & 0.755 & 0.729\\
$f_h$ (res4f) & 0.796 & 0.838 & 0.754 & 0.722\\
$f_l \oplus f_m \oplus f_h$ & \textbf{0.814} & \textbf{0.854} & \textbf{0.771} & \textbf{0.744}\\
\bottomrule
\end{tabular}}
\end{center}
\end{table}

 Experiments are done to choose the combination of layers with the best performance.
 We choose three sets from all the four sets, and each set has three or four layers to choose.
 So there are 135 combinations of layers in total.
 We compare the performance of all the combinations by running experiments 1000 times on the validation sets from ESPL-LIVE HDR Database, and find that the combination of layer $res2a, res4b, res4f$ achieves the best validation performance among all the 135 combinations.
 Furthermore, the same experiments are also done with the single layer features to validate the gain of the multi-layer features.
 Results of the comparison between the features from the selected combination and from the single layers in the combination are reported in Table~\ref{tab:multi-layer}.
 As what we expect, multi-layer features performs better than the single layer features.

 Features are adapted for TMIQA after aggregation.
 However, the length of the multi-scale and multi-layer features after aggregation is $$l_A = (l_{f_l} + l_{f_m} + l_{f_h}) \times 2,$$ where $l_A, l_{f_l}, l_{f_m}, l_{f_h}$ are the length of the aggregated features and features from the single layers.
 The dimension of features we adopted is 4608, which is much larger than the number of training samples.
 Features of such dimension may cause overfitting in regression.
 Therefore, we adopt partial least squares regression (PLSR)~\cite{PLSR} in our work, because it is suitable for regression tasks which have more variables than observations.
 PLSR reduces the dimension of input features to a few unrelated latent components, then maps these components to quality predictions.
 The number of latent components in PLSR is selected from the range $[10,20]$ by running experiments 1000 times on the validation sets,
 and the results show that the PLSR model with 15 latent components achieves the best performance.

 After investigations in this section, we obtain the best proposal for our method.
 Tone-mapped images are down-sampled for once, and the original and the down-sampled version of the images are fed into the pre-trained DCNN model of modified ResNet-50.
 Then feature maps are extracted from layer $res2a, res4b$ and $res4f$, then pooled by the mean and standard deviation of each channel, and concatenated to form the multi-scale and multi-layer features.
 The aggregated features are mapped to quality scores by the PLSR model with 15 latent components.

 \begin{table}[!tb]
\centering
\caption{Medians of SROCC, PLCC and RMSE between quality predictions and the MOS scores on the ESPL-LIVE HDR database. The best two values in each column are marked in boldface.}
\vspace{2mm}
\label{tab:result2}
\begin{small}
\begin{threeparttable}
\begin{tabular}{cp{16mm}p{13mm}p{13mm}}
\toprule	IQA & SROCC & PLCC & RMSE\\
\midrule
Proposed(100 runs)&\textbf{0.823}&\textbf{0.827}&\textbf{5.697}\\
Proposed(1000 runs)&\textbf{0.820}&\textbf{0.824}&\textbf{5.768}\\
HIGRADE-2~\cite{17TIP} & 0.730 & 0.728 & 6.992\\
BIBQA~\cite{17TIE}     & 0.702 & 0.692 & 7.133\\
DESIQUE~\cite{DESIQUE} & 0.570 & 0.568 & 8.296\\
GM-LOG~\cite{GMLOG}    & 0.556 & 0.557 & 8.357\\
BRISQUE~\cite{MSCN}    & 0.418 & 0.444 & 9.049\\
\bottomrule
\end{tabular}
\begin{tablenotes}
\footnotesize
\item[] Results of HIGRADE-2, DESIQUE, GM-LOG and BRISQUE are quoted from~\cite{17TIP}.
Results of BIBQA are quoted from~\cite{17TIE}.
\end{tablenotes}
\end{threeparttable}
\end{small}
\end{table}

\section{Experiments}
\label{sec:exp}

In this section, we conduct experiments on ESPL-LIVE HDR Database~\cite{ESPL} to compare the performance of our proposed method with some existing TMIQA methods and some state-of-the-art NR-IQA methods which are aimed at general SDR images.
Because most of the existing TMIQA methods are full-referenced and the original HDR images are not available in the ESPL-LIVE HDR Database, we compared two NR-TMIQA methods among them.
The compared NR-TMIQA methods are HIGRADE method proposed by Kundu et al.~\cite{17TIP}, and
the biologically inspired blind quality assessment (BIBQA) method proposed by Yue et al.~\cite{17TIE}.
Other general NR-IQA methods compared are DESIQUE~\cite{DESIQUE}, GM-LOG~\cite{GMLOG}, and BRISQUE~\cite{MSCN}.
The proposed method is implemented on MATLAB platform.
The results of HIGRADE, DESIQUE, GM-LOG, BRISQUE on the same database are quoted from~\cite{17TIP}, and the result of BIBQA method are quoted from~\cite{17TIE}.
However, the experimental settings are different in these two cited papers. In~\cite{17TIP}, the data is randomly split into training and testing sets at 4:1 ratio for 100 times; in~\cite{17TIE}, the data is split for 1000 times at the same ratio.
To make the comparison fair, we conduct experiments on both settings for our method.
\textcolor{black}
{
Three commonly used criteria are adopted to evaluate the performance:
Spearman's rank correlation coefficient (SROCC), which can demonstrate the prediction monotonicity,
Pearson's linear correlation coefficient (PLCC) and
root mean-squared error (RMSE), which are used for measuring prediction accuracy.
These criteria are calculated for every run,
and the medians of the three criteria in 100 runs and 1000 runs are reported separately.
}
Table~\ref{tab:result2} reports the comparison results.
Performance of the best proposal of HIGRADE method, HIGRADE-2 is reported as the representative performance of HIGRADE method.
And the medians in 1000 runs are quoted from \cite{17TIE} as the result of BIBQA, because the data is randomly split and there shouldn't be obvious difference between the medians in 100 or 1000 runs.
As reported in the table, comparison results prove that our proposed method achieves better performance.

\section{Conclusion}

In this paper, we propose a new NR-IQA method for tone-mapped HDR images. Multi-scale and multi-layer features are extracted from the pre-trained ResNet-50, and mapped to quality predictions by partial least squares regression. The performance of our proposed method is tested on the largest database for tone-mapped HDR image assessment, ESPL-LIVE HDR Database, and compared to the existing NR-IQA methods aimed at tone-mapped images and general images. The experimental results show that our method outperforms all the existing NR-IQA methods on tone-mapped image quality assessment, and also prove that multi-scale and multi-layer features lead to performance improvements in tone-mapped image quality assessment.

\section{Acknowledgment}
    This work was partially supported by National Basic Research Program of China (973 Program) under contract 2015CB351803 and the Natural Science Foundation of China under contracts 61572042, 61390514, 61527804.
    We also acknowledge the High-Performance Computing Platform of Peking University for providing computational resources.

\bibliographystyle{IEEEbib}