Tensor-Based Predictor–Corrector Algorithm for Power Generation and Transmission Reliability Assessment with Sequential Monte Carlo Simulation

Abstract

The reliability assesment of large power systems, particularly when considering both generation and transmission facilities, is a computationally demanding and complex problem. The sequential Monte Carlo simulation is arguably the most versatile approach for tackling this problem. However, assessing sampled states in the sequential Monte Carlo simulation is time-intensive, rendering its use less appealing, particularly if nonlinear network representation must be deployed. In this context, this paper introduces a tensor-based predictor–corrector approach to reduce the burden of state evaluations in power generation and transmission reliability assessments. The approach allows for searching for sequences of operation points which can be assigned as success states within the sequential Monte Carlo simulation. If required, failure states are evaluated using a cross-entropy optimization algorithm designed to minimize load curtailments taking into account discrete variables. Numerical results emphasize the applicability of the developed algorithms using a small test system and the IEEE-RTS79 test system.

1. Introduction

Reliability evaluation of large power systems, particularly when considering both generation and transmission facilities, is a complex and computationally intensive challenge [1]. The sequential Monte Carlo simulation (SMCS) can be considered the most flexible method to address this problem, permitting the representation of time-connected aspects of operation and the availability of resources. Nevertheless, evaluating sampled states in the SMCS is a time-consuming task, making its application unattractive, especially if nonlinear network representation must be deployed.

Several approaches to reduce the computational burden of state evaluations can be found in the literature, with the aim of replacing optimal power flow (OPF) analyses with alternative analyses, within the Monte Carlo simulation. In [2,3], multi-label classification approaches trained with OPF data are applied to substitute state evaluation procedures in the SMCS. In [4], the authors use multi-parametric linear programming to replace linearized OPF evaluations, utilizing a transmission line state dictionary to handle transmission line outages, while analyzing generation states and load variations. Other approaches based on the application of neural networks [5,6], such as self-organizing maps [7], support vector machines [8], convolutional neural networks (CNNs) [9,10,11,12,13], multilayer perceptrons [14], and a deep neural networks [15], have also been utilized, aiming at a rapid assessment of system states. A state space classification technique in conjunction with a dynamically directed particle swarm optimization (PSO) can be found in [16], where the first aims to classify success or failure states, whereas the binary PSO is used to identify failure states within the unclassified subspace. Strategies considering an approximated power flow formulation are provided in [1], focusing on acquiring more accurate results through a detailed network modeling.

The reported approaches aim to reduce the time required for state evaluation by replacing the OPF calculation with more computationally efficient methodologies. In contrast, some studies have been employing techniques that seek to mitigate the computational burden by reducing the number of OPF evaluations, focusing on identifying states which are deemed to be assigned as failure states. These approaches use techniques that aim to group samples with similar characteristics, with each of those labeled as either success states, failure states, or mixed, and the evaluation via OPF is triggered according to the classification acquired with the technique. In [17], unsupervised machine learning techniques, including the Kohonen self-organizing map, K-means clustering, and K-medoids clustering, are embedded into non-sequential Monte Carlo simulation for the purpose of classifying state groups. Self-organizing maps have been applied in [18] as a state filter for failure analysis, taking into account adverse weather and spatial conditions. Moreover, this concept can be found in other works that employ machine learning techniques, such as CNNs [19], the group method of data handling [20], and binary logistic regressions [21], where classifiers are built aiming at verifying whether system states can be classified as success states (i.e., states without load curtailment) or failure states (i.e., states without load curtailment).

Unlike other state-of-the-art works, this work directly exploits the fact that load transitions occur more often than component state transitions in the SMCS, aiming to reduce computational efforts in state evaluations. For this, the paper introduces a tensor-based predictor–corrector method to aid in the identification of success system states (i.e., states without load curtailment) within the SMCS, considering a nonlinear network representation. The method permits establishing a sufficient, but not necessary, condition to check if a visited state and its subsequent ones can be assigned as success states. The tensor-based predictor–corrector method utilizes a problem formulation inspired by the modified tensor power flow method [22], which is also used for state evaluations in this work. The novelties of this work can be summarized as follows:

• The design of a tensor-based predictor-correction method to be applied within the SMCS to reduce the computation burden of generation and transmission reliability assessments;
• The application of a cross-entropy optimization [23] approach to search for generation re-scheduling solutions corresponding to feasible operation points with a minimum loss of load.

Numerical results are provided for a small test system and the IEEE-RTS79 test system, with the aim of highlighting the relevant aspects of the proposed approach. For the IEEE-RTS79 test system [24], in a set of about 20 million states, around 8 million evaluations have been suppressed with simple system analyses based on modified tensor power flow assessments, while around 7 million evaluations have been suppressed using the tensor-based predictor–corrector approach. Results also illustrate the applicability of the cross-entropy approach in finding remedial actions in cases where the rescheduling of generation units is mandatory.

2. Designed Generation and Transmission Reliability Evaluation Using Nonlinear Network Modeling

Generation and transmission reliability assessment, also referred in the state-of-the-art research in the literature as composite reliability assessment, can be performed through SMCS. The SMCS is comprised of three core stages: state duration sampling, state evaluation, and the update of test functions for index estimation, as illustrated in Figure 1. Examples of reliability indices evaluated in composite assessments are the Loss of Load Probability (LOLP), Loss of Load Frequency (LOLF), Loss of Load Duration (LOLD), and Expected Power Not Supplied (EPNS). The state evaluation stage is the most time-consuming stage and where the application of the nonlinear network modeling has its major impact.

State evaluation in composite assessment can be attained using the following rules for robust power systems modeled with nonlinear network representation:

• If all system components are in the operating state, the system state can be considered successful, meaning that load curtailments are not needed, and events of a loss of load are not expected;
• If at least one component is in the failure state, power flow evaluations can be performed assuming standard dispatch profiles and control settings. In case the power flow solution corresponds to an operation point without violating system constraints, the system state can also be considered successful, meaning that load curtailments are not needed, and events of a loss of load are not expected. Otherwise, optimization problems must be solved aiming to verify and quantify the occurrence of load curtailments and corresponding remedial actions.

In our approach, power flow evaluations are executed using the modified tensor power flow method [22], which is shown to have improved performance in relation to its original counterpart. Remedial actions are deployed considering the solution of an optimal power flow problem targeting the minimization of load curtailments subject to operational limits. The objective function of the problem is given by the sum of the squares of the productions of fictitious generation units located at load buses. These fictitious generation units can produce power up to the amount of the load connected to each bus, with a specified cost of production that can be used to prioritize buses according to their importance. The formulation utilized in our approach is given by

$$\begin{array}{cc}\mathrm{minimize}\hfill & {\displaystyle \frac{1}{2}}{\mathbf{P}}_{\mathbf{g}}^{\mathrm{t}}{\mathbf{D}}_C{\mathbf{P}}_{\mathbf{g}}\hfill \\ \mathrm{subject}\mathrm{to}\hfill & {\mathbf{A}}_{\mathbf{g}}{\mathbf{P}}_{\mathbf{g}}-{\mathbf{P}}_{\mathbf{d}}-\mathbf{P}(\mathbf{e},\mathbf{f})=\mathbf{0}\hfill \\ \hfill & {\mathbf{A}}_{\mathbf{g}}^{\mathbf{r}}{\mathbf{Q}}_{\mathbf{g}}-{\mathbf{Q}}_{\mathbf{d}}+{\mathbf{A}}_{\mathbf{g}}^{\mathbf{f}}{\mathbf{P}}_{\mathbf{g}}-\mathbf{Q}(\mathbf{e},\mathbf{f})=\mathbf{0}\hfill \\ \hfill & {\mathbf{P}}_{\mathbf{g}}^{\mathbf{min}}\le {\mathbf{P}}_{\mathbf{g}}\le {\mathbf{P}}_{\mathbf{g}}^{\mathbf{max}}\hfill \\ \hfill & {\mathbf{Q}}_{\mathbf{g}}^{\mathbf{min}}\le {\mathbf{Q}}_{\mathbf{g}}\le {\mathbf{Q}}_{\mathbf{g}}^{\mathbf{max}}\hfill \\ \hfill & {\left({\mathbf{V}}^{\mathbf{min}}\right)}^{\mathbf{2}}\le ({\mathbf{e}}^{\mathbf{2}}+{\mathbf{f}}^{\mathbf{2}})\le {\left({\mathbf{V}}^{\mathbf{max}}\right)}^{\mathbf{2}}\hfill \\ \hfill & {\left({\mathbf{I}}^{\mathbf{min}}\right)}^{\mathbf{2}}\le {\mathbf{I}}^{\mathbf{2}}(\mathbf{e},\mathbf{f})\le {\left({\mathbf{I}}^{\mathbf{max}}\right)}^{\mathbf{2}}\hfill \end{array}$$

(1)

where ${\mathbf{P}}_{\mathbf{g}}$ is the vector of active power productions in each generation unit; ${\mathbf{D}}_C$ is a diagonal matrix with inputs given by generation costs; ${\mathbf{A}}_{\mathbf{g}}$, ${\mathbf{A}}_{\mathbf{g}}^{\mathbf{r}}$ e ${\mathbf{A}}_{\mathbf{g}}^{\mathbf{f}}$ are generation connection matrices, where ${\mathbf{A}}_{\mathbf{g}}(i,j)$ equals 1 if generator i (real or fictitious) is connected to bus j or 0 otherwise, ${\mathbf{A}}_{\mathbf{g}}^{\mathbf{r}}(i,j)$ equals 1 if generator i (real) is connected to bus j or 0 otherwise, ${\mathbf{A}}_{\mathbf{g}}^{\mathbf{f}}(i,j)$ equals the ratio of bus reactive and active load if generator i (fictitious) is connected to bus j or 0 otherwise; ${\mathbf{P}}_{\mathbf{d}}$ and ${\mathbf{Q}}_{\mathbf{d}}$ are vectors of bus active and reactive loads; $\mathbf{e}$ and $\mathbf{f}$ are the real and imaginary bus voltage vectors; ${\mathbf{P}}_{\mathbf{g}}^{\mathbf{min}}$ and ${\mathbf{P}}_{\mathbf{g}}^{\mathbf{max}}$ are the minimum and maximum active generator capacity vectors; ${\mathbf{Q}}_{\mathbf{g}}^{\mathbf{min}}$ and ${\mathbf{Q}}_{\mathbf{g}}^{\mathbf{max}}$ are the minimum and maximum reactive generator capacity vectors; ${\mathbf{V}}^{\mathbf{min}}$ and ${\mathbf{V}}^{\mathbf{max}}$ are the minimum and maximum bus magnitude voltage vectors; ${\mathbf{I}}^{\mathbf{min}}$ and ${\mathbf{I}}^{\mathbf{max}}$ are the minimum and maximum limits for currents in transmission elements. The active and reactive powers injected in each bus are written as

$$\begin{array}{cc}\hfill \mathbf{P}(\mathbf{e},\mathbf{f})& ={\mathbf{D}}_{\mathbf{e}}\mathbf{Ge}-{\mathbf{D}}_{\mathbf{e}}\mathbf{Bf}+{\mathbf{D}}_{\mathbf{f}}\mathbf{Be}+{\mathbf{D}}_{\mathbf{f}}\mathbf{Gf}\hfill \end{array}$$

(2)

$$\begin{array}{cc}\hfill \mathbf{Q}(\mathbf{e},\mathbf{f})& ={\mathbf{D}}_{\mathbf{f}}\mathbf{Ge}-{\mathbf{D}}_{\mathbf{f}}\mathbf{Bf}+{\mathbf{D}}_{\mathbf{e}}\mathbf{Be}+{\mathbf{D}}_{\mathbf{e}}\mathbf{Gf}\hfill \end{array}$$

(3)

where $\mathbf{e}$ and $\mathbf{f}$ are vectors with the real and imaginary parts of the bus voltages, respectively; while $\mathbf{G}$ and $\mathbf{B}$ denote the real and imaginary parts of the bus admittance matrix. The notation ${\mathbf{D}}_{\mathbf{u}}$ represents a diagonal matrix with diagonal values given by the entries of vector $\mathbf{u}$, while ${\mathbf{y}}^{\mathbf{2}}$ means a vector composed of the squared values of the entries in $\mathbf{y}$.

Additional constraints related to, for instance, active power flow and apparent power flow can be included straightforwardly in the formulation. Nevertheless, in case discrete variables associated with control functions or rescheduling of generation units must be modeled, the formulation in (1) is insufficient to allow a proper state evaluation. As an example, for states with a high availability of generation units, minimum generation limits may prevent OPF formulations from finding feasible solutions. Consequently, a rescheduling of generation units is required to search for feasible operating points and minimize load curtailments.

In our approach, discrete variables are addressed through the application of the CE method, where it is assumed that each discrete variable in a decision vector ${\mathbf{z}}_i$ follows a Bernoulli distribution, as

$$Be{r}_{n_h}({\mathbf{z}}_i,\mathbf{p})=\prod _{j=1}^{n_h}{p}_j^{z_j}{(1-p_j)}^{(1-z_j)}$$

(4)

where $Be{r}_{n_h}(·,\mathbf{p})$ is the $n_h$-dimensional Bernoulli probability mass function, $\mathbf{p}=[p_1,\cdots ,p_j,\cdots ,p_{n_h}{]}^{\mathrm{t}}$ denotes the probability vector associated to the Bernoulli distribution, and $z_j$ is the j-th component of the random variable ${\mathbf{z}}_i$.

Algorithm 1 summarizes the application of the CE method to minimize load curtailments with discrete control variables and rescheduling of generation units. The algorithm reaches an end if one of the stopping criteria is satisfied: (a) if a system state with zero load curtailment is found, or (b) if the probabilities $p_j^{t+1}$ have reached convergence. In both cases, the candidate solution with the minimal load curtailment is selected as the solution to the problem. The general algorithm as well as an example of the application of the CE method in combinatorial optimization problems can be found in [25] (pp. 34–36).

Algorithm 1 CE method for optimization of discrete control variables and rescheduling of generation units

Let N be the number of samples, $\varrho$ the rarity parameter, $N^e=\lceil \varrho N\rceil$ the number of elite samples, $t=1$, $U_t$ an empty set and $r_t$ an empty vector. 1: Define an initial probability vector ${\mathbf{p}}^{\mathrm{t}}$. 2: Draw N samples ${\mathbf{z}}_1,\cdots ,{\mathbf{z}}_N$ with $z_{ij}\sim Ber\left(p_j\right)$ and define $U_t=\left\{{\mathbf{z}}_1,\cdots ,{\mathbf{z}}_N\right\}$. 3: Calculate $S\left({\mathbf{z}}_i\right)$ using $S\left({\mathbf{z}}_i\right)$ as the objective function value acquired by solving (1), $\forall i$. 4: Set $r_t=\left\{S\left({\mathbf{z}}_1\right),\cdots ,S\left({\mathbf{z}}_N\right)\right\}$. Set $r_t$ in decreasing order and define ${\widehat{\gamma }}_t=r_t\left(N^e\right)$. 5: Calculate ${\widehat{p}}_j^{t+1}={\displaystyle \frac{{\sum }_{\mathbf{u}\in U_t}\left\{{\mathbf{I}}_{S\left(u\right)\ge {\widehat{\gamma }}_t}{u}_j\right\}}{N^e}}$ using the same sample ${\mathbf{z}}_1,\cdots ,{\mathbf{z}}_N$ and $\mathbf{I}$ is the indicator function; 6: Update the probability vector using the expression $p_j^{(t+1)}\leftarrow \alpha {\widehat{p}}_j^{t+1}+(1-\alpha )p_j^t$, for a given smoothing factor $\alpha$. 7: If stopping conditions are met, stop the simulation. Otherwise, set $t\leftarrow t+1$ and return to step 2.

The execution of the algorithm procedures outlined in this section renders time-consuming state evaluations in the SMCS. The tensor-based predictor–corrector approach seeks to leverage the fact that load transitions occur more frequently than component state transitions in the SMCS, making the estimation of maximum/minimum loadability factors of interest to assess subsequent states characterized by load transitions. These loadability factors represent a proportion of the load increase in a specified direction, given a set of operational constraints. Consequently, their estimation can be exploited to avoid power flow and OPF analysis, within the state evaluation procedures.

3. Tensor-Based Predictor-Corrector Approach to State Evaluation

The application of the tensor-based predictor–corrector algorithm within the SMCS for composite reliability assessments is summarized in Algorithm 2. State evaluation is described in Step 4 of the Algorithm 2, and is further summarized in Figure 2. The application of the tensor-based predictor–corrector approach can be outlined in the following steps:

(a)

Sort states according to load levels until the next component state transition;

(b)

Starting from state k, calculate the bus voltages via the predictor–corrector model until the operational limits described in (1) are violated, taking into account the visiting of load levels in ascending order;

(c)

Starting from state k, calculate the bus voltages via the predictor–corrector model until operational limits provided in (1) are violated, taking into account the visiting of load levels in descending order.

• The states corresponding to operation points without violated limits are marked with a passing mark, indicating that their evaluation using PF, OPF, or CE is not required.

Figure 2. Proposed state evaluation.

Figure 2. Proposed state evaluation.

The general formulation for estimating loadability factors for a given state is introduced in Section 3.1. This formulation can be solved either by using a direct or predictor–corrector method, with the latter described in Section 3.2.

3.1. Problem Formulation for Loadability Factor Estimation

The maximum loadability problem, in rectangular coordinates, considering system loads increased by the proportion ${\rho }^2$ in the direction $\Delta {\mathbf{P}}_{\mathbf{d}}$, can be described as follows.

$$\begin{array}{cc}\mathrm{minimize}\hfill & c(\mathbf{e},\mathbf{f},{\mathbf{P}}_{\mathbf{g}},{\mathbf{Q}}_{\mathbf{g}},\rho )\hfill \\ \mathrm{subjec}\mathrm{to}\hfill & {\mathbf{P}}_{\mathbf{g}}-({\mathbf{P}}_{\mathbf{d}}+{\rho }^2\Delta {\mathbf{P}}_{\mathbf{d}})-\mathbf{P}(\mathbf{e},\mathbf{f})=\mathbf{0}\hfill \\ \hfill & {\mathbf{Q}}_{\mathbf{g}}-({\mathbf{Q}}_{\mathbf{d}}+{\rho }^2\Delta {\mathbf{Q}}_{\mathbf{d}})-\mathbf{Q}(\mathbf{e},\mathbf{f})=\mathbf{0}\hfill \\ \hfill & {\mathbf{P}}_{\mathbf{g}}^{\mathbf{min}}\le {\mathbf{P}}_{\mathbf{g}}\le {\mathbf{P}}_{\mathbf{g}}^{\mathbf{max}}\hfill \\ \hfill & {\mathbf{Q}}_{\mathbf{g}}^{\mathbf{min}}\le {\mathbf{Q}}_{\mathbf{g}}\le {\mathbf{Q}}_{\mathbf{g}}^{\mathbf{max}}\hfill \\ \hfill & {\left({\mathbf{V}}^{\mathbf{min}}\right)}^{\mathbf{2}}\le ({\mathbf{e}}^{\mathbf{2}}+{\mathbf{f}}^{\mathbf{2}})\le {\left({\mathbf{V}}^{\mathbf{max}}\right)}^{\mathbf{2}}\hfill \\ \hfill & {\left({\mathbf{I}}^{\mathbf{min}}\right)}^{\mathbf{2}}\le {\mathbf{I}}^{\mathbf{2}}(\mathbf{e},\mathbf{f})\le {\left({\mathbf{I}}^{\mathbf{max}}\right)}^{\mathbf{2}}\hfill \end{array}$$

(5)

where $c(\mathbf{e},\mathbf{f},{\mathbf{P}}_{\mathbf{g}},{\mathbf{Q}}_{\mathbf{g}},\rho )=-\rho$ and the other variables have been previously defined. The formulation presented in (5) has as constraints the active and reactive power balance equations, where the weighted load increment is included. Other constraints are designed as in (1). The formulation of the minimum loadability problem is analogous to the one presented in (5), differing only by a change in the signal of the objective function.

Algorithm 2 SMCS for composite evaluation with the application of the tensor-based predictor–corrector approach

Let n, $t_n$ and $N_y$ be a state counter, time counter and period counter, respectively; Initialize counters as $n\leftarrow 1;t_1\leftarrow 0$; $N_y\leftarrow 1$. 1: Initialize component and load states to compose the system state ${\mathbf{x}}_n$; 2: Sample state residence times using stochastic or deterministic functions; 3: Identify the next state transition time $t_{n+1}$. If $t_{n+1}$ does not overpass the end of the period T, perform the transition to constitute the next system state ${\mathbf{x}}_{n+1}$. Otherwise, $t_{n+1}\leftarrow T$ and ${\mathbf{x}}_{n+1}\leftarrow {\mathbf{x}}_n$. 4: If system state ${\mathbf{x}}_n$ has at least one component in the failure state, execute the steps i. or ii. below. Otherwise, go to step iii. i. If ${\mathbf{x}}_n$ has been constituted by a component state transition (a) Perform state evaluation via modified tensor power flow method; if the state is considered successful, execute step iii and go to step (d); otherwise, go to the next step; (b) Evaluate state via OPF. If the state is considered successful, execute step iii and go to step (d); otherwise, store the load curtailment. In case of non-convergence, go to the next step; (c) Evaluate state via CE method. If the state is considered successful, execute step iii and go to step (d); otherwise, store the load curtailment and finish the state evaluation. In case of non-convergence, define the state as successful, save the state for future analysis, and finish the state evaluation. (d) Compute passing marks for system states characterized by load transition via tensor-based predictor–corrector method; ii. If ${\mathbf{x}}_n$ has been constituted by a load state transition: (a) If the state has a passing mark, go to step iii and finish state evaluation; Otherwise, go to the next step; (b) Execute step i.(a) without resourcing to step i.(d). Continue the process eventually reaching step i.(c) without resourcing step i.(d); iii. Classify the state as successful and set load curtailment as zero. 5: Update test functions $H\left({\left\{{\mathbf{x}}_n\right\}}_{n=1}^{s_y}\right)$, where ${\left\{{\mathbf{x}}_n\right\}}_{n=1}^{s_y}$ is the sequence of states ${\mathbf{x}}_n$, with $s_y$ states in period y; 6: If a period is completed, update reliability indices using the expected value equation $\mathrm{E}\left[H\left(\mathbf{X}\right)\right]\leftarrow 1/N_y{\sum }_{i=1}^{N_y}H\left({\left\{{\mathbf{x}}_n\right\}}_{n=1}^{s_i}\right)$. Otherwise, return to step 2; 7: Update coefficients of variation $\beta \leftarrow \sqrt{\mathrm{Var}\left[H\left(\mathbf{X}\right)\right]/N_y}/\mathrm{E}\left[H\right(\mathbf{X}\left)\right]$. If the coefficients of variation are greater than a specified threshold, perform $n\leftarrow 1$, $N_y\leftarrow N_y+1$ and return to step 2. Otherwise, stop simulation.

3.2. Loadability Factor Estimation Thought Predictor–Corrector Tensor Method (NLMCS$\widehat{\rho }$)

The identification of loadability factors for a given state is a sufficient but not necessary condition, to conclude whether the state and its subsequent states are successful states, during the execution of the SMCS method. Indeed, in order to suppress state evaluations conducted with PF, OPF, or CE, one may identify if the load level belongs to the interval defined by the loadability factors with successful states. In this context, the application of a continuation method is proposed in the range of possible variations of load levels to estimate minimum/maximum loadability factors, above/below which the states can be considered as successful.

Let $\rho$ be the loadability factor. The effect of the loadability factor in the system can be formulated as

$$\underset{\mathbf{h}\left(\mathbf{x}\right)}{\underbrace{\left[\begin{array}{c}\Delta \mathbf{P}(\mathbf{e},\mathbf{f},\rho )\\ \Delta \mathbf{Q}(\mathbf{e},\mathbf{f},\rho )\\ \Delta \mathbf{V}(\mathbf{e},\mathbf{f},\rho )\\ \Delta \sigma (\mathbf{e},\mathbf{f},\rho )\end{array}\right]}}=\underset{{\mathbf{h}}_{\mathbf{s}}}{\underbrace{\left[\begin{array}{c}\mathbf{0}\\ \mathbf{0}\\ {\mathbf{V}}^{{\mathbf{esp}}^{\mathbf{2}}}\\ {\sigma }^{esp}\end{array}\right]}}-\underset{{\mathbf{h}}_{\mathbf{c}}\left(\mathbf{x}\right)}{\underbrace{\left[\begin{array}{c}-\rho {\mathbf{P}}^{\mathbf{esp}}-\mathbf{P}\left(\mathbf{e},\mathbf{f}\right)\\ -\rho {\mathbf{Q}}^{\mathbf{esp}}-\mathbf{Q}\left(\mathbf{e},\mathbf{f}\right)\\ {\mathbf{e}}^{\mathbf{2}}+{\mathbf{f}}^{\mathbf{2}}\\ \rho \end{array}\right]}}$$

(6)

where $\Delta \mathbf{P}(\mathbf{e},\mathbf{f},\rho ),\Delta \mathbf{Q}(\mathbf{e},\mathbf{f},\rho ),\Delta \mathbf{V}(\mathbf{e},\mathbf{f},\rho )$ belong to the set of power balance equations and $\mathbf{x}$ is a variable redefined as $[\mathbf{e},\mathbf{f},\rho ]$. With the inclusion of $\rho$, the problem specified only by the power balance equations has infinity solutions, since the number of variables is greater than the number of equations. To overcome this issue, equation $\Delta \sigma (\mathbf{e},\mathbf{f},\rho )={\sigma }^{esp}-\rho$ has been added to the problem formulation, where ${\sigma }^{esp}$ is an additional loadability factor.

The Taylor series at point $\mathbf{x}$ is given by the general expression

$$\mathbf{h}(\mathbf{x}+{\mathbf{d}}_n)={\mathbf{h}}_{\mathbf{s}}+{\mathbf{J}}_{\mathrm{m}\mathrm{o}\mathrm{d}}\left(\mathbf{x}\right){\mathbf{d}}_n$$

(7)

such that a prediction of $\mathbf{x}$ using the Newton direction ${\mathbf{d}}_n$ is given by

$$\underset{{\mathbf{J}}_{\mathrm{m}\mathrm{o}\mathrm{d}}\left(\mathbf{x}\right)}{\underbrace{\left[\begin{array}{ccc}{\displaystyle \frac{\partial \Delta \mathbf{P}}{\partial \mathbf{e}}}& {\displaystyle \frac{\partial \Delta \mathbf{P}}{\partial \mathbf{f}}}& {\mathbf{P}}^{\mathbf{esp}}\\ {\displaystyle \frac{\partial \Delta \mathbf{Q}}{\partial \mathbf{e}}}& {\displaystyle \frac{\partial \Delta \mathbf{Q}}{\partial \mathbf{f}}}& {\mathbf{Q}}^{\mathbf{esp}}\\ {\displaystyle \frac{\partial \Delta \mathbf{V}}{\partial \mathbf{e}}}& {\displaystyle \frac{\partial \Delta \mathbf{V}}{\partial \mathbf{f}}}& \mathbf{0}\\ 0& 0& -1\end{array}\right]}}\underset{{\mathbf{d}}_n}{\underbrace{\left[\begin{array}{c}\Delta \mathbf{e}\\ \Delta \mathbf{f}\\ \Delta \rho \end{array}\right]}}=\underset{-\mathbf{h}\left(\mathbf{x}\right)}{\underbrace{\left[\begin{array}{c}\mathbf{0}\\ \mathbf{0}\\ \mathbf{0}\\ -{\sigma }^{esp}\end{array}\right]}}.$$

(8)

Moreover, $\mathbf{h}\left(\mathbf{x}\right)$ can be expanded considering the tensor term as follows

$$\mathbf{h}(\mathbf{x}+\mathbf{d})=\mathbf{h}\left(\mathbf{x}\right)+{\mathbf{J}}_{\mathrm{m}\mathrm{o}\mathrm{d}}\left(\mathbf{x}\right)\mathbf{d}+{\displaystyle \frac{\mathbf{1}}{\mathbf{2}}}{\mathbf{d}}^{\mathrm{t}}{\mathbf{T}}_{\mathrm{m}\mathrm{o}\mathrm{d}}\mathbf{d}$$

(9)

and using (6), we have

$$\mathbf{h}(\mathbf{x}+\mathbf{d})={\mathbf{h}}_{\mathbf{s}}-{\mathbf{h}}_{\mathbf{c}}(\mathbf{x}+\mathbf{d}).$$

(10)

Applying (9) at $\mathbf{x}=\mathbf{0}$, we can obtain

$$\mathbf{h}\left(\mathbf{d}\right)=\mathbf{h}\left(\mathbf{0}\right)+{\mathbf{J}}_{\mathrm{m}\mathrm{o}\mathrm{d}}\left(\mathbf{0}\right)\mathbf{d}+{\displaystyle \frac{1}{2}}{\mathbf{d}}^{\mathrm{t}}{\mathbf{T}}_{\mathrm{m}\mathrm{o}\mathrm{d}}\mathbf{d}={\mathbf{h}}_{\mathbf{s}}+{\mathbf{J}}_{\mathrm{m}\mathrm{o}\mathrm{d}}\left(\mathbf{0}\right)\mathbf{d}+{\displaystyle \frac{1}{2}}{\mathbf{d}}^{\mathrm{t}}{\mathbf{T}}_{\mathrm{m}\mathrm{o}\mathrm{d}}\mathbf{d}$$

(11)

and using (10), we have

$$\mathbf{h}\left(\mathbf{d}\right)={\mathbf{h}}_{\mathbf{s}}-{\mathbf{h}}_{\mathbf{c}}\left(\mathbf{d}\right).$$

(12)

Therefore,

$${\mathbf{h}}_{\mathbf{c}}\left(\mathbf{d}\right)=-{\mathbf{J}}_{\mathrm{m}\mathrm{o}\mathrm{d}}\left(\mathbf{0}\right)\mathbf{d}-{\displaystyle \frac{1}{2}}{\mathbf{d}}^{\mathrm{t}}{\mathbf{T}}_{\mathrm{m}\mathrm{o}\mathrm{d}}\mathbf{d}.$$

(13)

Suppose that the direction $\mathbf{d}$ can be written as $\mathbf{d}={\mathbf{d}}_n+{\mathbf{d}}_t$, where ${\mathbf{d}}_n$ is the Newton direction and ${\mathbf{d}}_t$ is the tensor direction component. Then,

$$\begin{array}{cc}\hfill \mathbf{h}(\mathbf{x}+\mathbf{d})& =\mathbf{h}(\mathbf{x}+{\mathbf{d}}_n+{\mathbf{d}}_t)\hfill \\ \hfill & =\mathbf{h}\left(\mathbf{x}\right)+{\mathbf{J}}_{\mathrm{m}\mathrm{o}\mathrm{d}}\left(\mathbf{x}\right)({\mathbf{d}}_n+{\mathbf{d}}_t)+{\displaystyle \frac{1}{2}}{({\mathbf{d}}_n+{\mathbf{d}}_t)}^{\mathrm{t}}{\mathbf{T}}_{\mathrm{m}\mathrm{o}\mathrm{d}}({\mathbf{d}}_n+{\mathbf{d}}_t)\hfill \\ \hfill & =\mathbf{h}\left(\mathbf{x}\right)+{\mathbf{J}}_{\mathrm{m}\mathrm{o}\mathrm{d}}\left(\mathbf{x}\right){\mathbf{d}}_n+{\mathbf{J}}_{\mathrm{m}\mathrm{o}\mathrm{d}}\left(\mathbf{x}\right){\mathbf{d}}_t+{\displaystyle \frac{1}{2}}{({\mathbf{d}}_n+{\mathbf{d}}_t)}^{\mathrm{t}}{\mathbf{T}}_{\mathrm{m}\mathrm{o}\mathrm{d}}({\mathbf{d}}_n+{\mathbf{d}}_t).\hfill \end{array}$$

(14)

Since ${\mathbf{d}}_n=-{\left[{\mathbf{J}}_{\mathrm{m}\mathrm{o}\mathrm{d}}\left(\mathbf{x}\right)\right]}^{-1}\mathbf{h}\left(\mathbf{x}\right)$, then,

$$\mathbf{h}(\mathbf{x}+\mathbf{d})={\mathbf{J}}_{\mathrm{m}\mathrm{o}\mathrm{d}}\left(\mathbf{x}\right){\mathbf{d}}_t+{\displaystyle \frac{1}{2}}{({\mathbf{d}}_n+{\mathbf{d}}_t)}^{\mathrm{t}}{\mathbf{T}}_{\mathrm{m}\mathrm{o}\mathrm{d}}({\mathbf{d}}_n+{\mathbf{d}}_t).$$

(15)

By solving (6), we obtain

$$\left[\begin{array}{c}\Delta \mathbf{P}(\mathbf{e}+\Delta \mathbf{e},\mathbf{f}+\Delta \mathbf{f},\rho +\Delta \rho )\\ \Delta \mathbf{Q}(\mathbf{e}+\Delta \mathbf{e},\mathbf{f}+\Delta \mathbf{f},\rho +\Delta \rho )\\ \Delta \mathbf{V}(\mathbf{e}+\Delta \mathbf{e},\mathbf{f}+\Delta \mathbf{f},\rho +\Delta \rho )\\ \Delta \sigma (\mathbf{e}+\Delta \mathbf{e},\mathbf{f}+\Delta \mathbf{f},\rho +\Delta \rho )\end{array}\right]=\mathbf{0}$$

which is equivalent to searching for step components $\Delta \mathbf{e}$, $\Delta \mathbf{f}$ e $\Delta \rho$ such that $\mathbf{h}(\mathbf{x}+\mathbf{d})$ is a null vector. In that case, we have

$$\begin{array}{cc}\hfill {\mathbf{d}}_t& =-{\displaystyle \frac{1}{2}}{\left[{\mathbf{J}}_{\mathrm{m}\mathrm{o}\mathrm{d}}\left(\mathbf{x}\right)\right]}^{-1}\left[{({\mathbf{d}}_n+{\mathbf{d}}_t)}^{\mathrm{t}}{\mathbf{T}}_{\mathrm{m}\mathrm{o}\mathrm{d}}({\mathbf{d}}_n+{\mathbf{d}}_t)\right]\hfill \\ \hfill & ={\left[{\mathbf{J}}_{\mathrm{m}\mathrm{o}\mathrm{d}}\left(\mathbf{x}\right)\right]}^{-1}\left[{\mathbf{h}}_{\mathbf{c}}({\mathbf{d}}_n+{\mathbf{d}}_t)+{\mathbf{J}}_{\mathrm{m}\mathrm{o}\mathrm{d}}\left(\mathbf{0}\right)({\mathbf{d}}_n+{\mathbf{d}}_t)\right].\hfill \end{array}$$

(16)

In general, when the iteration step is small, the behavior of the functions is approximately linear, and therefore the Newton step is considered an adequate approximation. Furthermore, it is the component with the greatest weight (importance) in the direction $\mathbf{d}={\mathbf{d}}_n+{\mathbf{d}}_t$. Therefore, the following approximation can be made.

$${\mathbf{h}}_{\mathbf{c}}({\mathbf{d}}_n+{\mathbf{d}}_t)+{\mathbf{J}}_{\mathrm{m}\mathrm{o}\mathrm{d}}\left(\mathbf{0}\right)({\mathbf{d}}_n+{\mathbf{d}}_t)\approx {\mathbf{h}}_{\mathbf{c}}\left({\mathbf{d}}_n\right)+{\mathbf{J}}_{\mathrm{m}\mathrm{o}\mathrm{d}}\left(\mathbf{0}\right)\left({\mathbf{d}}_n\right).$$

(17)

As a consequence, ${\mathbf{d}}_t$ can be estimated by solving the linear equation system

$${\mathbf{J}}_{\mathrm{m}\mathrm{o}\mathrm{d}}\left(\mathbf{x}\right){\mathbf{d}}_t={\mathbf{h}}_{\mathbf{c}}\left({\mathbf{d}}_n\right)+{\mathbf{J}}_{\mathrm{m}\mathrm{o}\mathrm{d}}\left(\mathbf{0}\right){\mathbf{d}}_n$$

(18)

where

$${\mathbf{x}}^{pred}=\mathbf{x}+\mathbf{d}$$

(19)

corresponds to a prediction step for the solution of the problem. Each prediction step is followed by a correction step, performed by one iteration of the modified tensor power flow method described, using the $\rho$ parameter retrieved after applying the prediction step.

During the SMCS, the states under evaluation are subject to a load curve, whose variation for two consecutive states is generally small. Thus, instead of carrying out the traditional prediction and correction steps for load factors that are not part of the load curve, bus voltages are computed for each load factor belonging to the set of load variations in between components’ failures. States corresponding to operation points without violated limits receive a passing mark which is utilized in step ii.(a) of Algorithm 2.

4. Result Analysis and Discussions

This section provides a set of numerical results illustrating the validation of the implementation of the proposed approach and a set of results using the IEEE-RTS79 test system for two SMCS applications. The first application, named as reference method (NLMCS$\tau$), consists of a SMCS where the state evaluation is performed using the PF, OPF, and CE method. The NLMCS$\widehat{\rho }$ method employs the tensor-based prediction-correction approach to avoid time-consuming state evaluations. All simulations have been conducted for a coefficient of variation $\beta$ of $5\%$, using a Notebook Acer Intel(R) Core(TM) i5-5200U (2.20 GHz, 8 GB RAM, Made in China).

4.1. Results for Methodology Validation

In order to validate the algorithms implemented, a two-bus system has been used, where three generation units are allocated in one bus, while the other bus has a constant load connected, as shown in Figure 3.

The generation units have minimum and maximum active generation capacities of 0 MW and 100 MW, respectively, as well as minimum and maximum reactive generation capacities of $-25$ MVAr and 60 MVAr. The system load is modeled by a curve comprising 8760 constant hourly steps corresponding to a load of 200 MW. The buses are connected by three transmission lines in parallel, with a 0.30 p.u ($\Omega$) reactance and a 1.10 p.u capacity (A). The voltage at the generation bus is specified at 1.00 p.u. and the minimum and maximum limits of voltage magnitude at the buses are 0.95 p.u. and 1.05 p.u., respectively. The voltage and power base of the system is 230 kV and 100 MVA, respectively. It is assumed that only two generators, named $G_1$ and $G_2$, and one of the transmission lines, named LT1, are subject to failures. The component failure rate and mean time to repair are 0.002 occ./year and 10 h, respectively, totalizing a probability of failure of 0.02. All other components are considered 100% reliable.

The enumeration and analysis of the states of this system are illustrated in

Table 1. Results for the analytical method considering a two-bus test system.

States			Load Cut [MW]	Prob	${\mathit{V}}^{\mathit{min}}$ [p.u.]
${\mathit{G}}_{\mathbf{1}}$	${\mathit{G}}_{\mathbf{2}}$	${\mathit{LT}}_{\mathbf{1}}$			PF	OPF
1	1	1	0	$0.941192$	$0.9789$	-
1	1	0	0	$0.019208$	$0.9487$	$0.9880$
1	0	1	0	$0.019208$	$0.9789$	-
1	0	0	0	$0.000392$	$0.9487$	$0.9878$
0	1	1	0	$0.019208$	$0.9789$	-
0	1	0	0	$0.000392$	$0.9487$	$0.9878$
0	0	1	100	$0.000392$	$0.9789$	$0.9875$
0	0	0	100	$0.000008$	$0.9487$	$0.9886$

.

In Table 1, columns 1–3 indicate the component states (1—in operation, 0—in failure), columns 4–5 indicate the load curtailment and probability of occurrence of each system state, and columns 6–7 indicate the voltage magnitude in each bus after solving the PF and OPF, where the OPF is executed only in case of operational limit violations.

For the two-bus system under evaluation, only in two states is the existence of load curtailment verified, assuming always the value of 100 MW. In addition, the probability of occurrence of these two states is $3.92\times {10}^{-4}$ if $L{T}_1$ is available, and $8\times {10}^{-6}$ if $L{T}_1$ is unavailable, totalizing a LOLP of $4\times {10}^{-4}$. The EPNS index is equal to $(0.000392+0.000008)\left(100\right)=4\times {10}^{-2}$ MW. Only the states with $G_1$, $G_2$, and $L{T}_1$ available, or where only one of the generators $G_1$ or $G_2$ is unavailable, do not require the execution of an OPF. The states where $G_1$ and $G_2$ are unavailable have required a solution through an OPF due to the insufficiency of the generation capacity.

The comparative results for the composite adequacy assessment using the analytical method (AM), NLMCS$\tau$, NLMCS${\tau }^{\prime }$, and NLMCS$\widehat{\rho }$ are presented in 

Table 2. Composite evaluation results for the AM, NLMCS$\tau$, NLMCS${\tau }^{\prime }$, and NLMCS$\widehat{\rho }$ methods.

Method	AM	NLMCS$\mathit{\tau }$	NLMCS${\mathit{\tau }}^{\prime }$	NLMCS$\widehat{\mathit{\rho }}$
#Simulated years	-	1313	1313	1313
Runtime [s]	-	$900.31$	$3983.69$	$651.03$
LOLP $[\times$${10}^{-4}]$	$4.0000$	$3.9097$	$3.9097$	$3.9097$
EPNS [MW]	$0.0400$	$0.0391$	$0.0391$	$0.0391$

,

Table 3. A 95% confidence interval of the reliability indices using NLMCS$\tau$, NLMCS${\tau }^{\prime }$, and NLMCS$\widehat{\rho }$ methods.

Method	LOLP $[\times$${10}^{-4}]$	EPNS [MW]
NLMCS$\tau$	$3.5268$	$0.0353$
	$4.2926$	$0.0429$
NLMCS${\tau }^{\prime }$	$3.5268$	$0.0353$
	$4.2926$	$0.0429$
NLMCS$\widehat{\rho }$	$3.5268$	$0.0353$
	$4.2926$	$0.0429$

and

Table 4. State evaluation for AM, NLMCS$\tau$, NLMCS${\tau }^{\prime }$, and NLMCS$\widehat{\rho }$ methods.

Counter	AM	NLMCS$\mathit{\tau }$	NLMCS${\mathit{\tau }}^{\prime }$	NLMCS$\widehat{\mathit{\rho }}$
#Visited states	8	$\mathrm{11,585,012}$	$\mathrm{11,585,012}$	$\mathrm{11,585,012}$
#Evaluated states	8	$\mathrm{501,360}$	$\mathrm{501,360}$	$\mathrm{501,360}$
#Success states (PF)	3	$\mathrm{496,002}$	0	$\mathrm{89,643}$
#OPF executions	5	5358	$\mathrm{501,360}$	5358
#Generation rescheduling	-	0	0	0
#States with $\rho$ estimation	-	0	0	$\mathrm{45,946}$
#Suppressed PF/OPF/ generation rescheduling	-	0	0	$\mathrm{406,359}$

. The NLMCS${\tau }^{\prime }$ is a variant of the NLMCS$\tau$ method, where states are evaluated using OPF exclusively. Result analysis allows us to conclude that the reliability indices acquired with NLMCS$\tau$, NLMCS${\tau }^{\prime }$, and NLMCS$\widehat{\rho }$ are equal and numerically close to indices obtained with the AM. Furthermore, with 95% confidence, retrieved confidence intervals are equal and contain the indices computed using AM.

As depicted in Table 2, the number of simulated years until convergence is also the same for all methods, although runtimes are distinct. In fact, the NLMCS$\widehat{\rho }$ method has converged in a reduced runtime, achieving a speed up of $1.38$ times in comparison to the NLMCS$\tau$ method. Moreover, the runtimes of methods NLMCS$\tau$ and NLMCS${\tau }^{\prime }$ support highlighting the importance of the application of a pre-evaluation with PF in the SMCS, since the NLMCS${\tau }^{\prime }$ is $4.42$ times slower than the NLMCS$\tau$.

Table 4 highlights, for the AM, NLMCS$\tau$, NLMCS${\tau }^{\prime }$, and NLMCS$\widehat{\rho }$ methods, the number of visited states; the number state evaluations via PF, OPF, or CE; the number of success states assessed straightforwardly via PF; the number of states evaluated as success or failure via OPF; the number of executed generation reschedulings; the number of states with loadability factor estimation; and the number of states with suppressed PF/OPF/CE evaluations using loadability factors.

As shown in Table 4, the number of visited states and evaluated states are the same for all SMCS methods and are equal to 11,585,012 and 501,360, respectively. The number of states evaluated only with PF is 496,002 using the NLMCS$\tau$ method, and 89,643 in the NLMCS$\widehat{\rho }$ method, which correspond approximately to $18.07$% of the number of PFs required by the reference method. The number of $\rho$ estimations is equal to 45,946, which led to 406,359 suppressed evaluations as well as runtime reduction, as depicted in Figure 4.

It is important to highlight that, since the number of executed OPFs remains the same for the NLMCS$\tau$ and NLMCS$\widehat{\rho }$ methods, only PF evaluations have been suppressed in this test.

4.2. Results for the IEEE-RTS79

The IEEE-RTS79 system [24] is a benchmark test system, composed of 24 buses, 32 generating units, 33 transmission lines and 5 transformers, as illustrated in the Figure 5. Numerical results of the composite reliability assessment for the IEEE-RTS79 system using the NLMCS$\tau$ and NLMCS$\widehat{\rho }$ methods are depicted in

Table 5. Composite evaluation results for the NLMCS$\tau$ and NLMCS$\widehat{\rho }$ methods.

Method	NLMCS$\mathit{\tau }$	NLMCS$\widehat{\mathit{\rho }}$
#Simulated years	2077	2102
Runtime [s]	$\mathrm{104,845.47}$	$\mathrm{77,550.33}$
LOLP [$\times {10}^{-3}$]	$1.3660$	$1.3118$
EPNS [MW]	$0.1488$	$0.1454$

,

Table 6. The 95% confidence intervals of the reliability indices using NLMCS$\tau$ and NLMCS$\widehat{\rho }$ methods.

Method	LOLP [$\times {10}^{-3}$]	EPNS [MW]
NLMCS$\tau$	$1.2698$	$0.1343$
	$1.4621$	$0.1633$
NLMCS$\widehat{\rho }$	$1.2230$	$0.1311$
	$1.4005$	$0.1596$

and

Table 7. State evaluation for NLMCS$\tau$ and NLMCS$\widehat{\rho }$ methods.

Counter	NLMCS$\mathit{\tau }$	NLMCS$\widehat{\mathit{\rho }}$
#Visited states	$\mathrm{19,202,348}$	$\mathrm{19,433,516}$
#Evaluated states	$\mathrm{14,862,995}$	$\mathrm{15,040,817}$
#Success states (PF)	$\mathrm{14,778,582}$	$\mathrm{8,059,994}$
#OPF executions	$\mathrm{63,819}$	$\mathrm{51,303}$
#Generation rescheduling	$\mathrm{20,594}$	$\mathrm{15,857}$
#States with $\rho$ estimations	0	$\mathrm{904,273}$
#Suppressed PF/OPF/ generation rescheduling	0	$\mathrm{6,913,663}$

.

Table 5 shows that the composite reliability indices obtained with NLMCS$\tau$ and NLMCS$\widehat{\rho }$ methods are similar. The NLMCS$\tau$ method requires 2077 years and approximately $29.12$ h to reach convergence, while the NLMCS$\widehat{\rho }$ requires 2102 years and approximately $21.54$ h, which represents a difference of 25 years and $7.58$ h. Furthermore, as depicted in Table 6, the LOLP and EPNS indices estimated using the NLMCS$\tau$/NLMCS$\widehat{\rho }$ method are within the corresponding confidence intervals estimated using the NLMCS$\widehat{\rho }$/NLMCS$\tau$ method.

The number of states visited and evaluated in each method are considerably similar. In the reference method, it has sampled 19,202,348 states and evaluated 14,862,995 states, while in the NLMCS$\widehat{\rho }$ method, it has sampled 19,433,516 states and evaluated 15,040,817 states. Furthermore, using the NLMCS$\tau$ method, 14,778,582 have been evaluated via PF, 63,819 via OPF, and 20,594 via generation rescheduling, which correspond to 99.43%, 0.43%, and 0.14% of the evaluated states, respectively. Regarding the NLMCS$\widehat{\rho }$ method, 8,059,994 states have been evaluated via PF, 51,303 via OPF, and 15,857 via generation rescheduling, which represents 53.59%, 0.34%, and 0.11% of the evaluated states, respectively. The remainder of states (45.96%) had their evaluations via PF/OPF/generation rescheduling suppressed by using loadability factors. In fact, 904,273 loadability factors have been estimated, leading to 6,913,663 suppressed PF/OPF/generation rescheduling evaluations. It has been identified that the slight difference between the indices can be further reduced by aggregating intermediate steps of the increase/reduction of load in the application of the predictor-correction method, improving the estimates of the voltages according the load variation. Nonetheless, for electrical system planning purposes, the index estimates are considered adequate.

Results of the application of the NLMCS$\tau$ and NLMCS$\widehat{\rho }$ methods are shown in Figure 6.

Despite the corrective actions adopted, and for both methods of NLMCS$\tau$ and NLMCS$\widehat{\rho }$, 40 sampled states could not be characterized as successful or failure states due to the lack of convergence of the CE method. It is important to note that the number of inconclusive states correspond to $2.083\times {10}^{-4}$ of the visited states in the NLMCS$\tau$ method, and $2.058\times {10}^{-4}$ of visited states in the NLMCS$\widehat{\rho }$ method, which lead to a minor influence in the reliability index estimations. Finally, the NLMCS$\widehat{\rho }$ achieves convergence in 77,550.33 s (approximately 22 h), while the NLMCS$\tau$ requires more than 29 h to achieve convergence, which represents a speed-up equal to $1.35$.

5. Discussion and Final Remarks

A tensor-based predictor–corrector method has been proposed with the aim of reducing the computational burden of state evaluations in SMCS approaches for composite reliability assessments. Moreover, a cross-entropy optimization approach is provided to search for generation rescheduling solutions corresponding to feasible operation points with minimum losses of load. The tensor-based predictor–corrector method is executed, where a state is considered successful and at least one component failure is sampled. In this case, passing marks are stamped on subsequent states characterized by load variations within loadability factors.

As a final point of discussion, it is worth mentioning that the loadability factor estimation problem is addressed in the proposed approach using a prediction-correction method instead of a direct method. The application of a direct method, for instance by using the interior-point method to solve (5), would ultimately involve solving an OPF with the loadability factor as the objective function, eventually increasing the computational effort required in the state evaluation stage. Preliminary implementations of direct methods aiming applications to composite reliability assessments have encouraged the authors to design the tensor-based predictor–corrector method.

Future works are envisioned to address the adaptation of the proposed approach to the composite evaluation of power systems with high-level integrations of renewable sources.