\section{Methodology}\label{sec:Methodology}

The proposed methodology is illustrated in Figure \ref{fig:Overall} and consists of A) an HVCB (component-level) reliability analysis, B) a contingency analysis that provides estimates of the system consequences of component failures, and C) a system-level reliability analysis. The following subsections describe each of the three parts of the methodology, with section numbers indicated in Figure \ref{fig:Overall}.
For the FMECA-based health index model (A1), the concept behind the health index model, feature extraction algorithms and the degradation scoring algorithms has previously been described in \cite{grant2023health} and \cite{Grant_RAMS}. For the condition-dependent HVCB reliability parameters model (A2), \cite{Grant_RAMS} has described the usage of the Cox Proportional Hazards Model. The integration of these models in the overall methodology as illustrated in Figure \ref{fig:Overall} is, however, described here for the first time.

\begin{figure}[]
\centering
\includegraphics[width=\linewidth]{IMG/Overall.pdf}
\caption{Methodology for incorporating HVCB technical condition into power system reliability analysis.}
\label{fig:Overall}
\end{figure}

\subsection{HVCB reliability analysis}

\subsubsection{FMECA-based health index model}

The health index model is based on the detection of the physical characteristics or failure mechanisms that lead to the HVCB failure modes. It includes the two failure modes of the HVCB tripping without a command and the HVCB failing to trip on command. These two failure modes were chosen because they make up the two main aspects of protection system reliability that are related to the technical condition of HVCB \cite{forootani2012circuit} \cite{bian2014protection}. The failure mechanisms that lead to the failure mode of the HVCB tripping without a command are shown in Table \ref{tab:FM2}, and the failure mechanisms that lead to the failure mode of the HVCB failing to trip on command are shown in Table \ref{tab:FM1} \cite{8751175}.

\begin{table}[]
\caption{Failure mechanisms that lead to the HVCB tripping without a command.}
\label{tab:FM2}
\begin{tabular}{l}
\toprule
FM01. Trip latch not secure \\
FM02. Stray current in trip circuit \\
FM03. Ground issues on trip circuit \\
FM04. Self-protective feature of HVCB \\
FM05. Loss of voltage on under voltage trip \\ \bottomrule
\end{tabular}
\end{table}

\begin{table}[]
\caption{Failure mechanisms that lead to the HVCB failing to trip on command.}
\label{tab:FM1}
\begin{tabular}{l}
\toprule
FM06. Open or shorted trip coil \\
FM07. Improper lubrication of the trip latch \\
FM08. Loss of spring energy \\
FM09. Control circuit failure \\
FM10. HVCB operation blocked \\
FM11. Mechanism failure between operating mechanism and interrupters \\
FM12. Trip latch surface wear \\
FM13. Operating mechanism below rated temperature \\
FM14. External circuit failure \\ \bottomrule
\end{tabular}
\end{table}

Condition monitoring measurements such as the trip coil current, travel curve, dynamic resistance, energy, and vibration measurements \cite{razi2020condition} can be used to detect the deterioration of failure mechanisms within HVCBs. The trip coil current (TCC) measurement can detect failure mechanisms in both the control and auxiliary circuits and the operating mechanism of an HVCB \cite{duan2016condition}. A typical TCC measurement is presented in Figure \ref{fig:typical}, which shows the current in the trip circuit of Figure \ref{fig:Control circuit} and resistance across the HVCB as a function of time. When there is a fault on a transmission line, the protective relay detects the fault via the current transformer and closes the trip circuit. The current in the trip circuit increases over time until Point A in yellow. An electromagnetic force causes an armature in the trip coil to start moving until Point B in green, when the armature hits the HVCB latch and stops. The current profile increases again until it reaches saturation at Point C in red. When the latch is pressed, the stored energy in the operating mechanism is released, opening the HVCB contacts and isolating the fault at Point O in blue on the resistance curve. Finally, the protective relay opens the trip circuit and the current decays to zero at Point D in black on the current curve \cite{hosseini2018current} \cite{razi2016circuit}.

\begin{figure}[]
\centering
\includegraphics[width=\linewidth]{IMG/Test 1.pdf}
\caption{Typical TCC measurements.}
\label{fig:typical}
\end{figure}

\begin{figure}[]
\centering
\includegraphics[width=\linewidth]{IMG/Control circuit.pdf}
\caption{Simplified circuit for tripping an HVCB.}
\label{fig:Control circuit}
\end{figure}

Features in the condition monitoring data that indicate the onset of failure mechanisms were identified in collaboration with HVCB experts at the Icelandic TSO. An unsecured trip latch failure mechanism was assessed based on the latch time, which was determined by measuring the duration between Points A and B. Any latch times that were too quick potentially indicated an unsecured trip latch. The noisiness of the TCC measurement was used to assess the possibility of a stray current or grounding issues within the trip circuit. The noise was evaluated after Point D, when the TCC should be stable. The Icelandic Transmission System Operator (TSO) observed that some HVCBs can trip automatically as a self-protective feature when the input power to the HVCB is too low and the tripping time is too long. Therefore, to assess this failure mechanism, the average current during the test was calculated, and the tripping time was calculated by determining the time until the resistivity curve transitioned to a high value at Point O. HVCBs with a low average current during the test and a long tripping time were considered to have a high risk of tripping due to the self-protective feature. The resistance curve was used to assess the possibility of an open or short circuit by measuring the open and closed circuit resistances. Anomalies between Points A and B were used to assess the possibility of improper lubrication or surface wear on the trip latch. A linear fit was computed between Points A and B, and the maximum deviation from this ideal line to the real data was extracted to quantify the latch anomalies.

A feature extraction algorithm was developed to extract the points illustrated in Figure \ref{fig:typical} from the TCC measurements. To select Point A, the process involved smoothing the current data, applying a derivative, and selecting the first point where the derivative equalled zero. Point C served as a helper point for determining the location of other points. It was calculated as the maximum current value after Point A. Point B was then selected by taking the minimum current value between Points A, and C. Point D was selected by locating the next zero current value after the test had begun. Finally, Point O was selected by taking the derivative of the resistivity curve and selecting the location of the value that had the highest rate of change.

Diagnostic data, such as the latch time, noisiness of the TCC data, average current during the test, tripping time, open/closed circuit resistances and latch anomalies, were extracted from the TCC measurements and used to quantify the amount of degradation of each failure mechanism. The amount of degradation was quantified as a score from 4 to 1, where 4 represents no signs of the failure mechanism occurring, and 1 represents a high degree of degradation. The range of the diagnostic data was observed, and thresholds were set to facilitate the mapping of the diagnostic data values to degradation scores. Threshold $a$ was selected to map the diagnostic data values to a score of 4, and threshold $b$ was selected to map the diagnostic data values to a score of 1. The generic algorithm for mapping diagnostic data values to scores is shown in Table \ref{tab:EJS_alg}.

\begin{table}[]
\caption{Algorithm for scoring the degradation of failure mechanisms.}
\label{tab:EJS_alg}
\centering
\begin{tabular}{l}
\toprule
\textbf{Algorithm:} Degradation of failure mechanisms \\
\midrule
\textbf{Input:} Diagnostic data, Thresholds  \\
\textbf{Output:} Score \\
\midrule
1: $a \gets$ \Comment{Threshold a} \\
2: $b \gets$ \Comment{Threshold b} \\
3: $x \gets$ \Comment{Diagnostic data} \\
4: \textbf{if} $x \geq a$ \textbf{then} \\
5: \quad $\text{Score[test]} \gets 4$  \\
6: \textbf{else} \\
7: \quad $\text{ min\_v} \gets b$ \\
8: \quad $\text{ max\_v} \gets a$ \\
9: \quad $\text{ min\_s} \gets 1$ \\
10: \quad $\text{max\_s} \gets 4$ \\
11: \quad $\text{Score[test]} \gets \left(\frac{x - \text{min\_v}}{\text{max\_v} - \text{min\_v}}\right) \times (\text{max\_s} - \text{min\_s}) + \text{min\_s}$ \\
12: \textbf{end if} \\
\bottomrule
\end{tabular}
\end{table}

Similar to the FMECA described in \cite{grant2023health}, a Risk Priority Number (RPN) is used as a weighting factor for each failure mechanism in the HVCB health index model. RPN values describe the risk of failure associated with a failure mechanism and are calculated as shown in Equation \ref{eq:RPN}. 
\begin{equation}
RPN = S \cdot O \cdot D
\label{eq:RPN}
\end{equation}
where $S$ represents the severity and $O$ represents the occurrence rate of a failure. The severity and occurrence rate were measured between 1 and 10, with higher values indicating increased severity and likelihood, respectively. $D$ is the detectability of a failure mechanism. It is also measured between 1 and 10, but higher values indicate a lower detectability, increasing the overall risk \cite{rausand2003system}. The calculated RPN values based on discussions with TSO experts for the failure mechanisms are shown in Figures \ref{tab:RPN_Ls} and \ref{tab:RPN_Pm}.

TCC measurements can detect degradation related to failure mechanisms FM01, FM02, FM03, and FM04, which lead to the HVCB tripping without a command, and failure mechanisms FM06, FM07, and FM12, which lead to the HVCB failing to trip on command. The weight of the health index that TCC measurements account for can be calculated by adding the RPN values of the failure mechanisms that TCC measurements can detect and dividing by the total of the RPN values for the respective failure mode. Based on this calculation, it is determined that TCC measurements can detect up to 92\% (272/296) of the weight of the failure mechanisms that contribute to the failure mode of the HVCB tripping without a command and up to 39\% (234/606) of the weight of the failure mechanisms that lead to the HVCB failing to trip on command.

\begin{figure}[]
\centering
\includegraphics[width=\linewidth]{IMG/RPN_Ls.pdf}
\caption{Calculated RPN values for the failure mechanisms that lead to the HVCB tripping without a command.}
\label{tab:RPN_Ls}
\end{figure}

\begin{figure}[]
\centering
\includegraphics[width=\linewidth]{IMG/RPN_Pm.pdf}
\caption{Calculated RPN values for the failure mechanisms that lead to the HVCB failing to trip on command.}
\label{tab:RPN_Pm}
\end{figure}

Based on the degradation scores and RPN weighting factors, health indices were quantified for the two HVCB failure modes using the health index assessment algorithm depicted in Table \ref{tab:HI_alg}.

\begin{table}[]
\caption{Algorithm for determining the health index considering failure mechanism score and risk}
\label{tab:HI_alg}
\centering
\begin{tabular}{l}
\toprule
\textbf{Algorithm:} Health index assessment \\
\midrule
\textbf{Input:} Score, RPN \\
\textbf{Output:} Health index\\
\midrule
1: $CalcValue \gets 0$ \\
2: $MaxValue \gets 0$ \\
3: \textbf{for} $n \gets 1$ \textbf{to} $NumFailMech$ \textbf{do} \\
4: \quad $CalcValue \gets CalcValue + (\text{Score}[n] \times \text{RPN}[n])$ \\
5: \quad $MaxValue \gets MaxValue + (4 \times \text{RPN}[n])$ \\
6: \textbf{end for} \\
7: $\text{$HI$} \gets \frac{CalcValue}{MaxValue}$ \\
\bottomrule
\end{tabular}
\end{table}

\subsubsection{Condition-dependent HVCB reliability parameters}

The Cox Proportional Hazards Model (PHM) was used to capture the dependence of the HVCB health index on the reliability parameters of individual HVCBs. These parameters represent the rate at which an HVCB trips without a command ($\lambda^s$) and the probability of an HVCB not tripping on command ($P^m$) when considering HVCB condition. The Cox PHM models the relationship between the failure rate and one or multiple covariates as shown in Equation \ref{eq:PHM1}.

\begin{equation}
\lambda(t;z) = \lambda_0(t)e^{(z \beta)} 
\label{eq:PHM1}
\end{equation}

where $\lambda_0(t)$ is the baseline failure rate, $z$ are the covariates and $\beta$ are the regression coefficients \cite{meeker2022statistical} \cite{leemis1995reliability}.

The Cox PHM calculates the regression coefficients without specifying the baseline failure rate by using the log-partial likelihood function \cite{meeker2022statistical}. The derivative of the log-partial likelihood function is then taken with respect to the regression coefficients, equated to zero, and solved for the regression coefficients using numerical methods like the Newton-Raphson method \cite{glover2012power}. Convergence is reached when the approximation to the function is within a predefined tolerance. 

The baseline failure rate is calculated using Breslows method \cite{kumar1994proportional}, and the covariates and regression coefficients are substituted for the component residual health ($1-HI\_\lambda^s$) as shown in Equation \ref{eq:Ls}.

\begin{equation}
\lambda^s(HI\_\lambda^s) = \lambda_0e^{(1-HI\_\lambda^s)} 
\label{eq:Ls}
\end{equation}

where $HI\_\lambda^s$ is the HVCB tripping without a command health index. The assumption of an exponential relationship between component residual health and failure rate was based on two observations. Firstly, the regression coefficients quantify the magnitude of the effects of the covariates \cite{kumar1994proportional} similar to how the product of scores and RPN values create the health index. Secondly, Brown \textit{et al.} \cite{brown2004failure} found that equipment condition data had an exponential effect on failure rates.

Model calibration is required if the predicted failure rates are to be expected to match observed historical rates \cite{brown2004failure,toftaker2023accounting}. Assuming that statistics on historical failure rates are available, this can be achieved by multiplying the predicted failure rates from the model by a calibration coefficient (CC) shown in Equation \ref{eq:CC}.

\begin{equation}
CC = \frac{\text{mean observed historical } \lambda}{\text{mean $\lambda$ predicted by model} }
\label{eq:CC}
\end{equation}

To model the condition dependence of the probability of missing operation parameter $P^m$, similarly, as for the parameter $\lambda^s$, the Cox PHM was used:   
\begin{equation}
P^m(HI\_P^m) = P_0e^{(1-HI\_P^m)} 
\label{eq:Pm}
\end{equation}
where $HI\_P^m$ is the HVCB failing to trip on command health index. 

\subsubsection{Equivalent line failure rate}

An equivalent line failure rate, including protection system failures, was calculated using existing methods outlined in \cite{vadlamudi2013dependability} \cite{vadlamudi2014impact} \cite{gjerde2016enhanced} \cite{kiel2020reliability}. With an approach adapted from \cite{kiel2020reliability}, adjacency matrices were used to indicate which transmission lines were adjacent as defined in Equation \ref{eq:A}, where $s$ represents the side of the line the HVCB is connected to, and $i$ and $j$ are transmission line indices. If the lines were adjacent the corresponding matrix element was given a 1; otherwise, 0.

\begin{equation}
A^s = [a_{i,j}^s] 
\label{eq:A}
\end{equation}

Probabilities also adapted from \cite{kiel2020reliability} were calculated. PT3 is the probability of an HVCB not tripping on command and modifies the adjacency matrix as shown in Equation \ref{eq:PT3}. PT4 is the probability of a fault being correctly cleared on line $i$ but an unwanted non-selective tripping event occurring on line $j$, as shown in Equation \ref{eq:PT4}. The probability of an unwanted non-selective tripping event occurring on a line that has two HVCBs ($A$ and $B$) was calculated by Equation \ref{eq:PNS}, where $P^u$ represents the probability of an unwanted non-selective tripping event occurring.

\begin{equation}
PT3 = \sum_s [A^s \cdot P^m]
\label{eq:PT3}
\end{equation}

\begin{equation}
PT4 = P_{NS} \cdot \sum_s [A^s \cdot (1 - P^m)]
\label{eq:PT4}
\end{equation}

\begin{equation}
P_{NS} = P_A^u + P_B^u - (P_A^u \cdot P_B^u)
\label{eq:PNS}
\end{equation}

The failure rates of the four fault types for each line $i$ were calculated using Equations \ref{eq:FT1} to \ref{eq:FT4} as adapted from \cite{vadlamudi2013dependability} \cite{vadlamudi2014impact} \cite{gjerde2016enhanced}.

\begin{equation}
\lambda_{FT1(i)} = \lambda_i
\label{eq:FT1}
\end{equation}

\begin{equation}
\lambda_{FT2(i)} = \sum_s \lambda_{i}^s 
\label{eq:FT2}
\end{equation}

\begin{equation}
\lambda_{FT3(i)} = PT3 \cdot \lambda_i
\label{eq:FT3}
\end{equation}

\begin{equation}
\lambda_{FT4(i)} = PT4 \cdot \lambda_i
\label{eq:FT4}
\end{equation}

The equivalent line failure rate, including protection system failures, was obtained by summing the four fault types as shown in Equation \ref{eq:lam_eq}.  

\begin{equation}
\lambda{(i)} =\lambda_{FT1(i)} + \lambda_{FT2(i)} + \lambda_{FT3(i)} + \lambda_{FT4(i)}
\label{eq:lam_eq}
\end{equation}

HVCB condition was incorporated into the equivalent line failure rate through the two condition parameters, $\lambda^s$ and $P^m$. $P^u$ was kept constant since this parameter was assumed to stem from relay selectivity issues rather than the HVCB's physical condition.

\subsection{Contingency analysis}
\label{sec:cont_anal}

\subsubsection{Interrupted power}

To effectively evaluate the system reliability indices, it is essential to consider both component failure rates and the consequences that arise from these failures or contingencies. The aim of the contingency analysis was to determine the interrupted power that arose from the contingencies. The interrupted power \(P_{int,j,k}\) as defined in \cite{Kjolle-Gjerde_TR_A7175_OPAL} was determined for each contingency state \(j\) at each delivery point \(k\) using Equation \ref{eq:LOL}, where \(P\) represents the load demand and \(SAC\) is the system available capacity. The $SAC$ is the available capacity to supply a load upon the occurrence of a contingency.

\begin{equation}
P_{int,j,k} = P_k - SAC_{j,k} 
\label{eq:LOL}
\end{equation}

When evaluating a system state, the \(SAC\) was calculated for each delivery point by considering certain factors:

\begin{itemize}
    \item What was the total generation capacity connected to the load (delivery point)?
    \item How did constraints such as transmission line limits reduce the generation capacity available to supply the loads? 
    \item Did the available generation capacity need to be shared with other loads?
    \item If so, what load curtailment philosophy should be used?
\end{itemize}

For this paper, a proportional load curtailment philosophy was used, as shown in Equation \ref{eq:prop}, where the available generation capacity \(G_a\) shared between a set of loads was shared proportionally according to their load.

\begin{equation}
SAC_k = \frac{P_k}{\sum_{k'} P_{k'}} \cdot G_a
\label{eq:prop}
\end{equation}

\subsection{System reliability analysis}

\subsubsection{Approximate methods}

Approximate methods for system reliability evaluation were utilised to reduce a system composed of components connected in series and parallel to an equivalent component \cite{billinton1992reliability}. Approximate methods are analytical reliability analysis methods that can be utilised when the failure frequency or mean time between failures (MTBF) is approximately equal to the failure rate or mean time to failure (MTTF), respectively. The equivalent basic reliability indices for $n$ components connected in series can be calculated using Equations \ref{eq:series_lam} to \ref{eq:series_r}.

\begin{equation}
\lambda_s = \sum_{i=1}^n \lambda_i \text{ f/yr}
\label{eq:series_lam}
\end{equation}

\begin{equation}
U_s = \sum_{i=1}^n (\lambda_i \cdot r_i) \text{ hours/yr}
\label{eq:series_U}
\end{equation}

\begin{equation}
r_s = \frac{U_s}{\lambda_s} \text{ hours}
\label{eq:series_r}
\end{equation}

Here $\lambda_s$, $U_s$, and $r_s$ denote the equivalent annual expected interruption frequency, the annual expected unavailability, and the expected interruption duration, respectively. The equivalent basic reliability indices for two components connected in parallel can be calculated using Equations \ref{eq:par_lam} to \ref{eq:par_U}. These equations can be used multiple times to calculate the equivalent basic reliability indices of multiple components connected in parallel.

\begin{equation}
\lambda_s = \frac{\lambda_1 \lambda_2 (r_1 + r_2)}{8760} \text{ f/yr}
\label{eq:par_lam}
\end{equation}

\begin{equation}
r_s = \frac{r_1 \cdot r_2}{(r_1 + r_2)} \text{ hours}
\label{eq:par_r}
\end{equation}

\begin{equation}
U_s = \lambda_s \cdot r_s \text{ hours/yr}
\label{eq:par_U}
\end{equation}

A cut set is a set of components that, when in an outage state, cause power supply interruptions. Furthermore, a minimal cut set is when a cut set contains no proper subsets that are also a cut set \cite{singh2018electric}. In the analysis of cut sets, up to second-order cut sets were considered since the reliability indices are generally dominated by lower-order cut sets \cite{billinton1992reliability}.

Equations \ref{eq:series_lam} to \ref{eq:par_U} assume independent failures. However, components can fail more frequently due to common-mode failures. Common-mode failures are when two or more components fail together due to a common cause \cite{singh2018electric} \cite{billinton1992reliability}. For example, when two lines $i$ and $j$ are adjacent within a cut set, a line failure on line $i$ could propagate to line $j$ due to protection system failures. A four-step process is utilised from \cite{vadlamudi2013dependability} \cite{vadlamudi2014impact} \cite{gjerde2016enhanced} to calculate the equivalent basic reliability indices when considering common-mode failures for neighbouring transmission lines.

% Step 1
\subsubsection*{Step 1: Remove dependency effects}

The dependency effects that come from the adjacent line are eliminated from both lines $x$ and $y$ within the cut set, as shown in Equation \ref{eq:step1}.

\begin{equation}\label{eq:step1}
\begin{split}
\lambda'_{FT3,x} = \lambda_{FT3,x} - (\lambda_{FT1,y} \cdot PT3)\\
\lambda'_{FT3,y} = \lambda_{FT3,y} - (\lambda_{FT1,x} \cdot PT3)\\
\lambda'_{FT4,x} = \lambda_{FT4,x} - (\lambda_{FT1,y} \cdot PT4 \cdot P_{NS})\\
\lambda'_{FT4,y} = \lambda_{FT4,y} - (\lambda_{FT1,x} \cdot PT4 \cdot P_{NS})
\end{split}
\end{equation}

% Step 2
\subsubsection*{Step 2: Calculate independent line equivalent failure rate}
The line equivalent failure rates are then calculated with the dependency effects removed as shown in Equation \ref{eq:step2}. The associated unavailability and outage times are calculated using a series system.

\begin{equation}\label{eq:step2}
\begin{split}
\lambda_x =\lambda_{FT1(x)} + \lambda_{FT2(x)} + \lambda'_{FT3,x} + \lambda'_{FT4,x}\\
\lambda_y =\lambda_{FT1(y)} + \lambda_{FT2(y)} + \lambda'_{FT3,y} + \lambda'_{FT4,y}
\end{split}
\end{equation}

% Step 3
\subsubsection*{Step 3: Calculate the dependency rate}
The dependency rate is then calculated by summing the dependency components from both lines as shown in Equation \ref{eq:step3}.

\begin{equation}
\begin{split}
\lambda_d = & (\lambda_{FT1,y} \cdot PT3) + (\lambda_{FT1,x} \cdot PT3) + \\
            & (\lambda_{FT1,y} \cdot PT4 \cdot P_{NS}) + (\lambda_{FT1,x} \cdot PT4 \cdot P_{NS})
\end{split}
\label{eq:step3}
\end{equation}

% Step 4
\subsubsection*{Step 4: Calculate the equivalent basic reliability indices}
Equations \ref{eq:dep_lam} to \ref{eq:dep_r} are used to calculate the basic reliability indices of two neighbouring lines.

\begin{equation}
\lambda_s = \frac{\lambda_x \lambda_y (r_x + r_y)}{8760} + \lambda_d \text{ f/yr}
\label{eq:dep_lam}
\end{equation}

\begin{equation}
U_s = \frac{\lambda_x \lambda_y r_x r_y}{8760} + \lambda_d \cdot r_d \text{ hours/yr}
\label{eq:dep_U}
\end{equation}

\begin{equation}
r_s = \frac{U_s}{\lambda_s} \text{ hours}
\label{eq:dep_r}
\end{equation}

\subsubsection{System indices}

The system interrupted power (\(P_{int,j,k}\)) was quantified using Equation \ref{eq:Pint} and the Expected Energy Not Supplied ($EENS$) was calculated using Equation \ref{eq:ENS} \cite{Kjolle-Gjerde_TR_A7175_OPAL}.

\begin{equation}
P_{int} = \sum_{j=1}^j [\lambda_j \cdot (P - SAC)] \text{ MW/year}
\label{eq:Pint}
\end{equation}

\begin{equation}
EENS = \sum_{j=1}^j [\lambda_j \cdot (r_j \cdot P_{int,j})] \text{ MWh/year}
\label{eq:ENS}
\end{equation}