## 2 Proposed Approach
### Continuation Technique  
Although the proposed approach for determination of \(\lambda^{\max}\) is independent of the model-based tool selected, the continuation method proposed in [22] was applied in this work. One feature of this approach is that it remains well-conditioned at all possible loading conditions in power flow calculation.  
Consider the following steady-state power system models  
\[F\left(\textbf{y},\,\lambda\right)=\textbf{0} \tag{11}\]  
Figure 5: Illustration of the cubic spline extrapolation technique for estimation of \(\lambda^{\max}\). (**a**) Case without \(Q\) limit violations; (**b**) Case with a \(Q\) limit violation.  
where \(\mathbf{y}\) denotes the vector of bus voltage phasors; \(\lambda\) represents the load increase parameter.  
The continuation approach is based on a locally parameterized continuation technique. Indeed, it is an iterative calculation process which starts from a known operating point and is solved as \(F\left(\mathbf{y},\lambda\right)=\mathbf{0}\) as \(\lambda\) changes. In addition, the predictor-corrector scheme is utilized to find each solution of the manifold.  
A widely accepted voltage stability measure can be expressed by the load power margin, i.e., \(\lambda^{\max}-\lambda^{0}\), where \(\lambda^{0}\) represents the base case loading point. Hence, VSM defined by the percentage of \(\lambda^{\max}-\lambda^{0}\) is  
\[VSM=\frac{\lambda^{\max}-\lambda^{0}}{\lambda^{0}}\times 100\% \tag{12}\]  
In summary, the proposed method aims to quickly determine the proximity to voltage collapse. It can cope with the cases that cause generators reaching their \(Q\) limits. The computation process to determine VSM is summarized in Algorithm 2. The input of this algorithm is \(Y_{\text{bus}}\) and three sets of \(\left|Z^{\text{L}}\right|\), \(\lambda\), and \(Q_{\text{G}}\), and the output is VSM.  
The main steps in Algorithm 2 for finding VSM are described as follows:  
1. Compute \(Q_{\text{G}}^{\text{ex}}\) for \(i\in\{PV\}\) to create a list of \(Q\) limit violations \(L_{PV}^{PQ}\);
2. Employ a multi-node Thevenin equivalent network to model a power system. Next, calculate \(\left|Z_{j}^{\text{th}}\right|\) and \(\left|Z_{j}^{\text{L}}\right|\) for \(j\in\{PQ\}\) based on \(L_{PV}^{PQ}\);
3. Utilize cubic spline extrapolation technique to estimate \(\lambda^{\text{es}}\) based on \(\left|Z_{j}^{\text{th}}\right|\) and \(\left|Z_{j}^{\text{L}}\right|\);
4. Execute a continuation program to determine the exact \(\lambda^{\max}\) through the information of \(\lambda^{\text{es}}\);
5. Compute VSM via Equation (12).  
```
1:Input:\(Y_{\text{bus}}\), three sets of \(Q_{\text{G}}\), \(\lambda\), and \(\left|Z_{j}^{\text{L}}\right|\);
2:\(D\leftarrow\left\{Q_{\text{G}i}(k),\,\lambda(k)\right\}_{k=1}^{3}\);
3:for\(i\in\{PV\}\)do
4:\(L_{PV}^{PQ}\leftarrow\text{Qlimit($D$)}\);
5:endfor
6: Compute \(Z\) by Equation (4) based on \(L_{PV}^{PQ}\);
7: Compute \(\left|Z_{j}^{\text{th}}\right|\) for \(j\in\{PQ\}\) by Equation (8);
8: Estimate \(\lambda^{\text{es}}\) via cubic spline extrapolation technique;
9: Determine \(\lambda^{\max}\) based on \(\lambda^{\text{es}}\) via a continuation program;
10: Compute VSM by Equation (12);
11:return VSM
```  
**Algorithm 2** Determine VSM via the proposed method.