Application of Cubic Spline Interpolation Technique in Power Systems: A Review

Abstract

In this chapter, a comprehensive review is made on the application of cubic spline interpolation techniques in the field of power systems. Domains like available transfer capability (ATC), electric arc furnace modeling, static var. compensation, voltage stability margin, and market power determination in deregulated electricity market are taken as samples to illustrate the significance of cubic spline interpolation.

1. Introduction

With electricity becoming an inevitable part of all spheres of human life, it is imminent that the increasing demand for electricity be met. The realization of this necessity has manifested in extensive research in the field of power systems, which has brought to light the complexity of power system. The power system involves the continuous variation in connected loads and increases in power requirement which demands a corresponding increase in the generation. This is met via several means most recently characterized by the penetration of renewable energy sources. The unpredictability of these events entails the deployment of probabilistic-based load flow techniques, estimation of unknown variables and load models which includes uncertainties for their analysis.

All the analysis and research aimed at providing solutions or mitigating the above-mentioned problems require vast amounts of data which is discrete in nature or the analysis techniques involve responses that are discrete in nature. Regardless of the origin of this non-continuity of data, obtaining a continuous response is imperative because of the desired accuracy and the continuity of real-time operation.

The challenges set the stage for interpolation techniques to play an active role in mitigating the problems prevailing in the current power system studies like load forecasting, power system reactive power planning, transmission network expansion, available transfer capability (ATC) determination and market power/clearing price forecasting. Among the various interpolation techniques available, the cubic spline method has been found to be a popular method. Cubic spline also has the desirable characteristics of continuous derivatives at data points which make the design of controllers around these regions possible, and its employment has been seconded by the high accuracy obtained.

The measurements we take for analysis in everyday life are wrought with noise that may be caused by the surrounding environment. Especially in electrical measurements, such random noises may be caused by the magnetic field produced by the current, the presence of stray charges, heating caused by the flow of eddy currents, and so on. One of the possible solutions to the problems caused by noise is to take a large number of measurements. This ensures that the random noise gets canceled out on an average, and hence the integrity of the data is maintained. Other complicated methods of tackling the errors caused by noise are available, but the detailed and in-depth analysis of these methods falls beyond the scope of this chapter.

2. Determination of available transfer capability

The first application of cubic splines is in finding the available transfer capability, which is an important parameter in power system operation. Ever since the advent of the deregulated power system, the computation of transfer capability has been a priority. Two quantities that require special attention in these computations are the total transfer capability (TTC) and the available transfer capability (ATC). The TTC of a power system is defined as the maximum amount of 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 [1]. ATC is defined as the measure of the additional amount of power that may flow over and above the base case flows without jeopardizing power system security [1]. The system operators of the deregulated power system normally can obtain previously calculated ATC values through an open-access same-time information system. However, the need for a quick calculation of ATC poses a challenge. Some of the methods employed include DC power flow, AC power flow, optimal power flow and sensitivity. Most of the methods mentioned earlier are not sufficiently agile as far as computation speed is concerned or trades accuracy for speed. For example, DC power flow yields quick results at the expense of accuracy, whereas AC power flow compromises speed for accuracy. A solution to the abovementioned problem is found by employing curve-fitting techniques especially the cubic spline interpolation technique.

In order to compute ATC, the cubic spline is employed to trace the curve of the variation of voltage magnitude and power flow with respect to the real power transfer. The ATC is then determined by the point where the limits of voltage magnitude or power flow intersect the curves. The computation of ATC takes place in two different forms—the point-to-point ATC and area-to-area ATC. The area-to-area ATC refers to the additional power that can be mobilized from the seller area to the buyer area, whereas the point-to-point ATC refers to bus-to-bus transfer (usually from a generating bus to a load bus). Another factor considered of paramount importance for the computation of ATC is the effect of contingencies like line outages. Individual consideration of line outages for a large-scale power system is an impractical approach, and hence contingency analysis is carried out using contingency ranking which helps select the critical lines.

The mathematical definition of ATC is shown in Eq. (1):

$$\mathrm{ATC}=\mathrm{TTC}-\mathrm{TRM}-\mathrm{ETC}-\mathrm{CBM}$$E1

The definition of the abovementioned terms is as follows:

Transmission reliability margin (TRM): the amount of transmission capability necessary to ensure the security of the interconnected system under a reasonable range of uncertainties in the system.

Capacity benefit margin (CBM): the amount of transmission capability reserved by load-serving entities to ensure access to generation from interconnected systems to meet generation reliability requirements.

Existing transmission commitments (ETC): normal transmission flows.

Traditionally, ATC computation involves the recursive application of AC power flow with increasing power transfers, thereby tracing voltage magnitudes and MVA power flow variations with respect to real power formally called the P-V and the P-S curves as shown in Figure 1. These curves are then employed in conjunction with the limits imposed by acceptable voltage magnitudes and power flows to calculate the ATC. In the cubic spline-based approach, the cubic spline interpolation technique is used to trace the P-V and P-S curves which are then employed to calculate the ATC. This is done by first determining four known points on the curve using AC power flow and then using cubic spline interpolation to trace the curves in between those points. The four known points are denoted as V_i (P_n) and S_ij (P_n), where V_i is for voltage of bus ‘I’, S_ij is the MVA power flow between bus ‘i’ and ‘j’, P_n is the real power transfer and n is the index for the four points and has values of n = 1, 2, 3 and 4. The incremental steps for tracing the P-V or the P-S curve will be of 1 MW each.

Figure 1.

Variation of voltage with real power transfer.

The first and foremost step in the process of curve tracing consists of the determination of the four points on the P-V and P-S curves. This is of paramount importance because the power system undergoes voltage collapse once the real power transfer crosses a certain limit. In the voltage collapse process, the particular bus faces a continuous drop in bus voltage once the critical point load is exceeded. Therefore, the points must be selected carefully making sure that the critical point is not exceeded.

2.1. Procedure to calculate ATC

1. Step 1. Perform the line contingency ranking using line-loading performance index (PI_MW) and bus voltage performance (PI_v). This is done to identify the critical lines. The critical lines will have a PI value greater than PI_{base case}.

2. Step 2. Find a base case by solving the AC load flow. This will act as the first point.

3. Step 3. Perform the simulation for line outage for one of the critical lines found in step 1.

4. Step 4. Specify the point or area of transfers. For the point-to-point transfer, generally a generator bus is taken as the selling bus and a load bus as the buying bus. On the other hand, for area-to-area transfer, all the generator busses in a particular area called the selling are and all the load busses in a specific buying area are considered.

5. Step 5. Now, the next three points are determined. The first step is to determine the fourth and final point P4. In order to do so, the sensitivity method is employed. This method is based on the limiting point of system constraints which is given in Eqs. (2)–(4)

   $${\mathrm{PT}}_{\mathrm{i},{\mathrm{V}}_{\mathrm{L}}}=\left|\frac{\mathrm{\partial \lambda }}{\partial {\mathrm{V}}_{\mathrm{i}}}\times \left({\mathrm{V}}_{\mathrm{L}}-{\mathrm{V}}_{\mathrm{i}}^0\right)\right|$$E2
   $${\mathrm{PT}}_{\mathrm{i},{\mathrm{V}}_{\mathrm{U}}}=\left|\frac{\mathrm{\partial \lambda }}{\partial {\mathrm{V}}_{\mathrm{i}}}\times \left({\mathrm{V}}_{\mathrm{U}}-{\mathrm{V}}_{\mathrm{i}}^0\right)\right|$$E3
   $${\mathrm{PT}}_{\mathrm{ij},\mathrm{S}}=\left|\frac{\mathrm{\partial \lambda }}{\partial {\mathrm{S}}_{\mathrm{ij}}}\times \left({\mathrm{S}}_{\mathrm{ij}}^{\text{limit}}-{\mathrm{S}}_{\mathrm{ij}}^0\right)\right|$$E4

   In Eqs. (2)–(4), PT_i,VL, PT_i,VU and PT_ij,S are the estimates of P₄ which are arrived at by using the linear estimates based on the lower voltage limit, the upper voltage limit and the thermal limit of each line, respectively. The lower limit of voltage V_L is taken as 0.9 p.u. and the upper limit V_u is taken as 1.1 p.u. S^limit_ij is the thermal limit of each line. V_i⁰ and S⁰_ij are the base case values computed in step (ii). The derivative terms are the reciprocal of the rate of change of voltage and MVA power flow with real power. This can be obtained by solving the AC power flow starting from the base case and making a small change in the transferred real power and noting the change in voltage and MVA transfer. The fourth known point P₄ is selected as the minimum of all three values computed above. Then, based on this, P₂ and P₃ are chosen using the formulae given in Eqs. (5) and (6)

   $${\mathrm{P}}_2=\frac{{\mathrm{P}}_4}{3}$$E5
   $${\mathrm{P}}_3={\mathrm{P}}_2\times 2$$E6

6. Step 6. Now, find the MVA and bus voltage Vi for each of the newly computed real power values—P2, P3 and P4.

7. Step 7. Now, use the cubic spline to trace the curve in between these points and obtain the P-V and the P-S curves.

Now, the area-to-area and point-to-point ATCs are calculated by obtaining the point where the voltage limit or the MVA limit line intersects the P-V and P-S curves, respectively.

Inference: The author of [1] has employed the proposed cubic spline-based method to compute the ATC of the Malaysian power system. The performance of the method is measured in terms of its accuracy and speed (which are the main reasons for the employment of the method). For testing the Malaysian power system, it has been simplified into a 241-bus system and further classified into three regions namely north, east, central, south and PUB. The lower and upper voltage limits are taken as 0.9 and 1.1 p.u., respectively. The results obtained have been compared with the traditional recursive AC power flow results.

It is observed from the results obtained in [1] that all the ATC values obtained are due to the MVA limitation. It is further observed that the performance of the cubic spline interpolation method is comparable to the AC power flow method in terms of accuracy. The observation of paramount importance is that while the proposed method obtains the high accuracy as found in the recursive AC power flow method, the time it requires to compute the ATC is much smaller than the AC recursive power flow method—to the extent that in some cases, the time required is up to 30 times lesser. Therefore, it may be concluded that the cubic spline method is superior to AC power flow in terms of speed and superior to DC power flow in terms of accuracy.

3. Electric arc furnace modeling

In this section, cubic spline is applied in modeling a very complex load—electric arc furnace (EAF). This example presents the true power of cubic splines in tracking extremely complex trajectories. The electric arc furnace’s (EAF) ability to efficiently smelt scrap iron raw materials has made it the backbone of the steel-making industries. The EAF employs high temperatures produced by low-voltage and high-current electric arcs to smelt scrap iron raw materials. The increase in productivity requirements has led to EAFs being designed for high-power applications. The operation of EAFs introduces a significant amount of harmonics, inter-harmonics and flickers in the supply system. Therefore, it is mandatory that the operators pay much attention to the power-quality considerations. In order to study these problems, it is of paramount importance to understand the nonlinear load characteristics that the EAF present to the power system. Also, in order to mitigate the power-quality issues and to further study the impact of EAFs on the power system, it is required that the EAF be modeled after obtaining its time response.

The entire operation of the EAF involves three stages namely striking, smelting and refining. In the striking process, the electric arc is built up by lowering the electrodes of all three phases, the melting process involves the melting of the material and the process ends with the stable refining. Due to the complexities involved in modeling the EAF operation in the striking and the melting processes, most of the research has been directed towards a steady-state modeling of the EAF in the refining stage of operation. The modeling of an EAF requires several parameters such as arc voltage, arc length and arc current, which are determined by the position of electrodes. Therefore, in order to accurately model the EAF, we need to know the field measurements of the electric response which involves the variation of voltage and current. The measured responses are then employed to develop an EAF conductance model using the cubic spline interpolation method.

The EAF is modeled as a function of nonlinear conductance using the cubic spline interpolation technique which is called the cubic spline interpolation model (CSIM). In this method, a set of cubic polynomials are obtained which helps understand the voltage–current characteristic of the EAF.

The steps involved are as follows:

First, a set of ‘n’ measured data points of conductance is obtained for one fundamental cycle of operation. The measurements lying in the interval [a, b] such that a = x₀ < x₁ < … < x_n = b. For the interval between two adjacent points, a cubic function is defined as shown in Eq. (7).

$${\mathrm{G}}_{\mathrm{i}}\left(\mathrm{x}\right)={\mathrm{a}}_{\mathrm{i}}+{\mathrm{b}}_{\mathrm{i}}\left(\mathrm{x}–{\mathrm{x}}_{\mathrm{i}}\right)+{\mathrm{c}}_{\mathrm{i}}{\left(\mathrm{x}–{\mathrm{x}}_{\mathrm{i}}\right)}^2+{\mathrm{d}}_{\mathrm{i}}{\left(\mathrm{x}–{\mathrm{x}}_{\mathrm{i}}\right)}^3$$E7

where i = 0, 1, 2,…., n-1. The coefficients a_i, b_i, c_i and d_i are unknown. These coefficients need to be determined based upon the following constraints:

1. Step 1. Each spline must pass through the given data points y_i.

   G_i(x_i) = y_i and G_n-1(x_n) = y_n

2. Step 2. Interior data points between each spline must be continuous.

   G_i + 1(x_i + 1) = G_i(x_i + 1)

3. Step 3. The first and the second derivatives of the splines must be continuous across the interior data points. Therefore, the spline forms a smooth function.

   G’_i + 1(x_i + 1) = G’_i (x_i + 1).

   G^”_i + 1(x_i + 1) = G^”_i(x_i + 1)

4. Step 4. In addition to the conditions mentioned in steps 1–3, another boundary condition must be satisfied which concerns the derivative of the functions at the boundaries (at x₀ = a and x_n = b). There are two types of boundary conditions that may be required to satisfy:

   Natural boundary condition: G^”₀(x₀) = G^”_n-1(x_n) = 0.

   Clamped boundary condition: G^’₀(x₀) = D₀ and G^’_n-1(x_n) = D_n.

where D₀ and D_n are the values of the first derivatives of the unknown functions.

It is generally found that the natural boundary conditions give less accurate results than the clamped boundary conditions. Alternately, one could possibly apply a boundary condition called not-a-knot condition, which in addition to the natural boundary condition also incorporates another condition that the third derivative of the function must be continuous at x₁ and x_n-1.

In order to find the coefficients—a_i, b_i, c_i and d_i—we follow the following steps and equations:

1. Step 1. Set $a$_i = y_i for i = 0, 1, 2, …, n.

2. Step 2. Solve

   

   E8

   and

   

   E9

   and

   

   E10

   where ${\mathrm{h}}_i=x_{i+1}-x_i,$ for i = 1, 2, …, n-1.

3. Step 3. Set

   

   E11
   

   E12

4. Step 4. Set

   

   E13

5. Step 5. Set

   

   E14

The coefficients $a$_i can be found in step 1. The coefficients c_i can be calculated by solving the linear equations developed in Eqs. (8)–(10). Then, finally coefficients b_i and d_i may be obtained from steps 4 and 5 (Eqs. (13) and (14), respectively).

Inference: The author in [2] has analyzed the performance of CSIM on an actual power system model developed in MATLAB/SIMULINK whose online diagram is also presented in [2]. The results obtained via employing the CSIM (shown in Figure 4) are then compared with the results obtained via employing two traditional methods namely (i) harmonic current injection model (HCIM), shown in Figure 2, and (ii) harmonic voltage source model (HVSM), shown in Figure 3. The results described in [2] have been obtained for three different parts—early, middle and later stages—of EAF operation during the refining period. This is done primarily because EAF does not possess a steady-state behaviour that lasts long. Since the HCIM and HVSM are unable to model dynamic behaviour, these methods are compared with the CSIM only during the latter part of the operation because of the dynamic behaviour of the EAF in the early and middle parts (Figure 4).

Figure 2.

HVSM V-I characteristic [2].

Figure 3.

HCIM V-I characteristic [2].

Figure 4.

CSIM V-I characteristic (later stage) [2].

Upon employing the CSIM for EAF modeling in the early and middle stages, it is found that the performance of the EAF during these stages is highly nonlinear. The V-I curve presents a multi-valued function which makes it difficult to model the EAF. The results shown in [2] suggest that upon the comparison, the HCIM method yields the largest errors in terms of both the EAF voltage and current determined. Then, HVSM, which albeit an improvement on HCIM, still contains errors. The proposed CSIM method performs better than HCIM and HVSM and provides a better fit for the voltage and current characteristics. It is also explained in the results that the error encountered during modeling using CSIM in the early and middle stages is larger when compared to that encountered when modeling the later stages of the refining period because the number of sampling points used for the EAF model is not sufficient for modeling in the early and middle stages of refining. As an extension, it is also proposed in [2] that the cubic spline interpolation technique can also be used for modeling other nonlinear loads.

4. Var compensator with thyristor-controlled reactor

The problem of maintaining a good power factor (greater than 0.85) is a challenge faced by most of the industries. Some of these bulk consumers are penalized for operation under a low power factor. The application of cubic splines is presented to topic address this challenge and provide a commendable solution. Rolling mills and electric arc furnaces constitute very large loads in the power system. When the problem of power factor is considered, these loads (particularly electric arc furnaces) may be termed as a necessary evil because of their extremely low power factor. The variation in the arc length during operation results in the introduction of severe and rapid fluctuations in the reactive power and the voltage, and when the short circuit occurs, the power factor drops to values as low as 0.1. The large impact current and reactive power generated result in significant waste of energy and may also cause the power system to lose its stability. This may cause decadence in the quality of load and endanger the users. The solution to the abovementioned problem is cited in a process called reactive power consumption. Currently, the reactive power compensation is achieved dynamically by the placement of a dynamic reactive power compensation device at access points of such interference loads. These devices are usually represented as a fixed capacitor and thyristor-controlled reactor (FCT). The FCT enables smooth control of the reactive power and also has the desirable feature of maintaining its voltage unchanged. In addition, the FCT can effectively suppress voltage fluctuation and solve the voltage distortion and flicker problem and improve the power quality. The continuous and smooth variation in the reactive power is obtained by the variation of the thyristor conduction angle. For the FCT to effectively carry out its task, it is imperative that it calculates the control angle quickly and accurately in a real-time environment. However, the existence of a nonlinear relation between the control angle and the reactor amplification factor makes the real-time calculation a cumbersome process.

The remainder of this section reviews the calculation and the corresponding control of the control angle via employing cubic spline interpolation technique. The calculation of the control angle, α, is based on the reactor amplification factor for each phase. Cubic spline interpolation is employed because it provides a low-order polynomial interpolation polynomial and also increases the smoothness of the interpolation function.

The basic thyristor-controlled reactor (TCR) consists of a pair of anti-parallel thyristors in series with an inductor as shown in Figure 5. The thyristor delay angle varies between 90° and 180°. As a result, the fundamental current is completely reactive. An increase in the delay angle leads to a decrease in the fundamental reactive current, which is equivalent to increasing the reactance or reducing the susceptance and hence results in a decrease in the reactive power. Hence, the TCR can be seen as being equivalent to a variable susceptance which can be controlled using the delay angle ‘α’. This is because the AC voltage remains constant but the value of the fundamental current changes, which in turn results in the variation in reactive power. The equivalent fundamental susceptance of the TCR is given in Eq. (15)

$${\mathrm{B}}_{\mathrm{L}}\left(\mathrm{\alpha }\right)=\frac{2\mathrm{\pi }-2\mathrm{\alpha }+sin2\mathrm{\alpha }}{{\mathrm{\pi X}}_{\mathrm{L}}}$$E15

Figure 5.

Thyristor controlled reactor.

The relationship between the delay angle, α, and the amplification factor of the equivalent impedance is shown in Eq. (16)

$$\mathrm{k}=\frac{\mathrm{\pi }}{2\left(\mathrm{\pi }-\mathrm{\alpha }\right)+sin2\mathrm{\alpha }}$$E16

As can be seen, there exists a nonlinear relation between the delay angle and the amplification factor. The application of cubic splines makes the control of delay angle with reactive power variation easier, quicker and smoother. It must be observed that the inductance current in the TCR does not depend on the inductance but is rather governed by the thyristor conduction angle. The accurate control of the thyristor angle makes the accurate control of inductor current and hence reactive power possible.

The problem is formulated as follows. Eq. (9) suggests that in order to control the delay angle using the amplification factor ‘k’, we need to find a solution to the nonlinear equation. Albeit possible, it is a cumbersome task to solve nonlinear equations, not to mention the considerable amount of time the solving takes which make real-time application and control a tedious task. Moreover, of paramount importance is the accuracy of the controlling action since it determines the reactive power compensation. The cubic spline interpolation is employed as a means to calculate the control angle ‘α’ having known the amplification factor ‘k’. Cubic spline is used as the method of interpolation because of the advantages it provides in terms of simplicity of calculation, numerical stability and smoothness of the interpolated curve.

As previously stated, the existence of a nonlinear relation between the control angle and the amplification factor makes the real-time computation a lengthy process. In order to achieve the aim using cubic spline interpolation, the control angle is regarded as the dependent variable and the amplification factor ‘k’ as the independent variable. According to [3], the author generated 158 data points using Eq. (2). The aim is now to generate 157 interpolating cubic polynomials that fit in between these points. The form of these cubic polynomials is shown in Eq. (17)

$$\mathrm{\alpha }={\mathrm{a}}_0{\left(\mathrm{k}-{\mathrm{k}}_{\mathrm{i}}\right)}^3+{\mathrm{a}}_1{\left(\mathrm{k}-{\mathrm{k}}_{\mathrm{i}}\right)}^2+{\mathrm{a}}_2\left(\mathrm{k}-{\mathrm{k}}_{\mathrm{i}}\right)+{\mathrm{a}}_3$$E17

In Eq. (10), i = 1, 2, …., 157 and k_i ≤ k ≤ k_i + =1.

The fitting coefficients a₀, a₁, a₂ and a₃ need to be computed for each interval which gives a unique cubic polynomial for each interval. When carrying out the control procedure first, the amplification is determined, and then the corresponding cubic polynomial is used to arrive at the required conduction angle value.

Inference: The results shown in [3] suggest that the conduction angle determined by using cubic splines shows significant match with the values obtained by simulation. Also, the use of cubic splines yields quicker results—a trait which would be beneficial for real-time applications. Hence, it can be said that the cubic spline interpolation method is able to solve the problem of quick computation of the conduction angle with accuracies relevant to engineering applications.

5. Estimation of voltage stability margin

In this section, the application of cubic splines is of paramount importance since it marks the limit of a power system—the maximum amount of power that can be supplied before the continuous uncontrolled drop of voltage. The impact of voltage instability on the power system has been so damaging that significant research has been conducted in this direction. The phenomenon of voltage instability is characterized by a sudden and uncontrollable drop in voltage as a response to a disturbance that has occurred on the power system. This disturbance could be anything from the variation of load to loss of a line or lightning strike, and so on. Most of the studies have been focused on the steady-state aspect of voltage stability. In order to determine the proximity of the system to voltage collapse, we need to estimate or find the voltage stability margin (VSM). Many methods have been proposed in the past to estimate the steady-state voltage stability margin. One of the popular methods is the continuous power flow method. A major disadvantage of the continuous power flow is that it requires a considerable amount of time and hence it cannot be employed in real-time applications.

An alternative approach is to use the time-synchronized voltage and current phasor measurements obtained from phasor measurement units (PMUs). In the presented method, the author proposed a combination of a coupled single-port Thevenin equivalent model and cubic spline extrapolation in order to find the point of voltage collapse or the voltage stability margin (VSM). The concept is based on the fact that the voltage collapse point of the load impedance equals the Thevenin equivalent impedance.

Thevenins’ equivalent voltage and equivalent impedance as seen from load bus ‘i’ as given in Figure 6 can be written as shown in Eq. (18)

$${\mathrm{E}}_{\mathrm{thi}}={\left[\mathrm{K}{\mathbf{V}}_{\mathrm{G}}\right]}_{\mathrm{i}}-\sum _{\mathrm{j}=1,\mathrm{j}\ne \mathrm{i}}^{\mathrm{n}}{\mathrm{Z}}_{\mathrm{ij}}{\mathrm{I}}_{\mathrm{Lj}}$$E18

$${\mathrm{Z}}_{\mathrm{thi}}={\mathrm{Z}}_{\mathrm{ii}}$$

Figure 6.

Coupled single port Thevenins’ equivalent.

The definitions for the abovementioned variables can be found in [4].

The load impedance Z_Li of bus ‘i’ can be arrived by using Eq. (19)

$${\mathrm{Z}}_{\mathrm{Li}}=\frac{{\mathrm{V}}_{\mathrm{Li}}}{{\mathrm{I}}_{\mathrm{Li}}}$$E19

The load bus voltage V_Li and the load bus current I_Li are obtained through the PMU measurements.

Figure 7 shows the variation of Z_li and Z_th as a function of the load parameter λ. As can be seen, the Z_li equals the Z_thi at the point of maximum critical loading, and this point gives the maximum value of the load parameter—λ_max. This maximum value of load parameter can be arrived at by equating an approximate function that extrapolates the Z_li versus λ curve to the point that it meets the Z_thi line.

Figure 7.

Variation of Thevenin and load impedance with λ [4].

The choice of cubic spline extrapolation is justified by the author as the superior fitting to the impedance trajectory as evidenced by extensive simulation results. The cubic spline extrapolation proceeds by developing different cubic polynomials for the interval between measurements based on certain constraints.

For ‘m’ PMU-measured data points obtained, there will be ‘m-1’ intervals present, each of which is represented by Eq. (20)

$${\mathrm{f}}_{\mathrm{j}}\left(\mathrm{x}\right)={\mathrm{a}}_{\mathrm{j}}+{\mathrm{b}}_{\mathrm{j}}\left(\mathrm{x}-{\mathrm{x}}_{\mathrm{j}}\right)+{\mathrm{c}}_{\mathrm{j}}{\left(\mathrm{x}–{\mathrm{x}}_{\mathrm{j}}\right)}^2+{\mathrm{d}}_{\mathrm{j}}{\left(\mathrm{x}–{\mathrm{x}}_{\mathrm{j}}\right)}^3$$E20

The inputs for extrapolation are taken as three sets of |Z_Li| and λ and the point of |Z_li| where λ is to be extrapolated and found. The λ thus found corresponds to the point where Z_li equals Z_thi and hence is equal to λ_maxi. Each bus in the considered system will have its own corresponding λ_maxs. The λ_max = λ_sys for the entire system is found by taking the minimum of all λ_maxi obtained for each of the individual load busses. This value of λ_sys corresponds to the value of λ_maxi for the weakest load bus in the system. The proximity of the load bus to the voltage collapse point is given in Eq. (21)

$$\mathrm{VSM}=\frac{{\mathrm{\lambda }}_{\mathrm{sys}}-{\mathrm{\lambda }}_0}{{\mathrm{\lambda }}_0}\times 100$$E21

After the computation of VSM, the author of [4] proceeds to find an index which helps in determining whether the load increase results in the violation of the reactive power limits at the generator busses. The author begins with the simple power equation and arrives at the equation of a surface given by Eq. (22)

$${\mathrm{Q}}_{\mathrm{G}}^2-\frac{{\left|{\mathrm{V}}_{\mathrm{G}}\right|}^2}{\mathrm{X}}{\mathrm{Q}}_{\mathrm{G}}+\frac{{\left|{\mathrm{V}}_{\mathrm{G}}\right|}^2}{\mathrm{X}}{\mathrm{Q}}_{\mathrm{L}}+{\mathrm{P}}_{\mathrm{L}}^2=0$$E22

where P_L, load real power; Q_L, load reactive power; Q_G, generator reactive power.

The surface defined lies in the (P_L, Q_L, Q_G) space. Upon cutting the surface with the constant power factor planes, we get the P_L–Q_G curves which are similar to PV curves or λV curves. The existence of an approximate quadratic relation between λ and V is extended to find a quadratic relation between λ and Q_G. This approximate quadratic relation is modeled and given in Eq. (23)

$$\mathrm{\lambda }={\mathrm{\alpha }}_{\mathrm{i}}{\mathrm{Q}}_{\mathrm{Gi}}^2+{\mathrm{\beta }}_{\mathrm{i}}{\mathrm{Q}}_{\mathrm{Gi}}+{\mathrm{\gamma }}_{\mathrm{i}}$$E23

In Eq. (16), α_i, β_i and γ_i are parameters that are different for each generator bus and need to be determined. The parameters can be determined using three sets of PMU readings. Then, using the computed values of the parameters, an estimate for the extreme of Q_Gi is arrived using the condition that at the extreme limit dλ/dQ_Gi = 0 holds. Enforcing this condition gives us the following value as the estimate shown in Eq. (24)

$${\mathrm{Q}}_{\mathrm{Gi}}^{\mathrm{ex}}=-\frac{{\mathrm{\beta }}_{\mathrm{i}}}{2{\mathrm{\alpha }}_{\mathrm{i}}}$$E24

The index i stands for all the generator busses of the system. At each load step, it is ensured that the estimated value of generator reactive power—Q _Gi^ex—remains within the bounds of the reactive power limits that are specified for a generator bus.

The algorithm for proceeding with the described method in order to compute the VSM is as follows:

1. Step 1. Obtain three sets of PMU measurements and use them to compute three sets of |Z_Li|, λ and Q_Gi. Also, we possess the knowledge of the admittance matrix.

2. Step 2. Compute the parameters α_i, β_i and γ_i for each of the generator bus of the system.

3. Step 3. Compute the estimate of the extreme of the reactive power for each generator ${\mathrm{Q}}_{\mathrm{Gi}}^{\mathrm{ex}}$

4. Step 4. If there is no violation of the reactive power limits of any generator, then proceed to step 6. Otherwise, proceed to step 5.

5. Step 5. Change the bus type from PV to PQ bus for the bus whose reactive power limit has been violated.

6. Step 6. Compute the impedance matrix as shown in [4].

7. Step 7. Compute the Thevenin equivalent impedance.

8. Step 8. Estimate the value of λ_maxi by using the cubic spline extrapolation technique for each load bus of the system.

9. Step 9. Find the value of λ_sys which is the minimum of all values computed in step 8.

10. Step 10. Determine the VSM.

Inference: The author of [4] has performed test runs of the proposed algorithm on the following test system—IEEE 30 bus system, IEEE 118 bus system and IEEE 300 bus system. Further, the obtained results have been compared with two other previously available methods. The obtained results indicate the superior performance of the proposed cubic spline technique when compared to the other methods. Upon comparison of the percentage error, it may be observed from [4] that the cubic spline method is almost 10 times more accurate than the other methods.

6. Estimation of market power in deregulated electricity market

In this final section, the application of cubic spline technique has been illustrated in the field of deregulated electricity market. In this market, market power issues predominantly spoil the basic idea of maintaining equilibrium within the market players.

Market power is the ability of showing one’s monopolistic nature on the price of the commodities in the market. This has become a challenging issue in the context of the present electricity market and will become more challenging and play a significant role when private generation companies start participating in buying/selling the power [5]. Due to the increase in demand and the regulatory policies, private parties have started investing in the power sector, especially in the renewable energy sources. Thus, it is inevitable for the independent system operator (ISO) to estimate market power for taking crucial decisions [6].

6.1. Market power indices

6.1.1. Herfindahl-Hirschman index

The Herfindahl-Hirschman index (HHI) is used to measure the market concentration that will reflect the number of players in the market and also the inequality in their market shares. The HHI is defined as the sum of the squares of market shares of all the players as given in Eq. (25)

$$\mathrm{HHI}=\sum _{i=1}^N{S}_i^2\mathrm{HHI}=\sum _{\mathrm{i}=1}^{\mathrm{N}}{{\mathrm{S}}_{\mathrm{i}}}^2$$E25

where N is the number of players and S_i is the i^th player market share in percentage [7].

6.2. Lerner index

It measures or indicates the proportional deviation of the price at the firm’s profit-maximizing output from the firm’s marginal cost at that output. It is defined as shown in Eq. (26)

$${\mathrm{LI}}_{\mathrm{i}}=\frac{\left({\mathrm{\rho }}_{\mathrm{i}}-{\mathrm{mc}}_{\mathrm{i}}\right)}{{\mathrm{\rho }}_{\mathrm{i}}}=\frac{1}{{\mathrm{\epsilon }}_{\mathrm{i}}^{\mathrm{d}}}{\mathrm{LI}}_{\mathrm{i}}=\frac{{\mathrm{\rho }}_{\mathrm{i}}-{\mathrm{mc}}_{\mathrm{i}}}{{\mathrm{\rho }}_{\mathrm{i}}}=\frac{1}{{\mathrm{\epsilon }}_{\mathrm{i}}^{\mathrm{d}}}$$E26

where LI_i is the Lerner index for a given firm i, ρ_i and mc_i are the price and marginal cost, respectively, and ${\mathrm{\epsilon }}_{\mathrm{i}}^{\mathrm{d}}$ is the elasticity of demand felt by the firm.

6.3. Nodal must run share

This index reflects the impact of load variation on market power and geographic difference of market due to network constraints. Must run share (MRS) represents the effect of load variation and nodal must run share (NMRS) represents the geographical difference of market powers. The equation for NMRS is shown in Eq. (27)

$${\text{NMRS}}_{\mathrm{k},\mathrm{i}}=\frac{{\mathrm{Pg}}_{\mathrm{k},\mathrm{i}}^{\text{must}}}{{\mathrm{Pd}}_{\mathrm{i}}}{\text{NMRS}}_{\mathrm{k}.\mathrm{i}}=\frac{{\text{Pgmust}}_{\mathrm{k},\mathrm{i}}}{{\mathrm{Pd}}_{\mathrm{i}}}\mathrm{i}=1,2,..\mathrm{N}$$E27

where $\mathrm{N}$ is the number of busses in a power system, ${\mathrm{Pd}}_{\mathrm{i}}$ is the load at bus i, and $\frac{{\mathrm{Pg}}_{\mathrm{k},\mathrm{i}}^{\text{must}}}{{\mathrm{Pd}}_{\mathrm{i}}}{\text{Pgmust}}_{\mathrm{k},\mathrm{i}}$ is the contribution of the must run generator k to ${\mathrm{Pd}}_{\mathrm{i}}$. The background calculation of NMRS is available in Ref. [7].

Steps involved in estimation of market power using NMRS [5, 8].

1. Step 1. Define the number of generators and their active power limits.

2. Step 2. Determine ${\mathit{Pg}}_k^{\mathit{\text{must}}}$ of generator k

3. Step 3. Calculate distribution matrix [M⁻¹]

4. Step 4. Calculate NMRS of generator 1 on load 1.

5. Step 5. Repeat step 4 for calculating NMRS for the remaining generators on each load.

6. Step 6. Repeat steps 4 and 5 for various cases.

7. Step 7. End.

8. Step 8. Calculate the NMRS for all the generators using the conventional method as discussed in steps–7.

9. Step 9. Plot NMRS of a generator against the maximum generation of other generators.

10. Step 10. Connect two points at a time using cubic spline interpolation technique (piecewise polynomial) using MATLAB built-in function, that is, for example, a = spline (b, c, de) and ff = spline (b, c) where ‘a’ gives the interpolated values which correspond to the query points in de [9].

The application of cubic’s spline interpolation for calculating the respective market power of generation companies on all busses for any given load or operating condition is implemented on a sample IEEE 14 bus test system. This system consists of four generator busses and an additional slack generator at bus number 1 in addition to the nine load busses. The interconnection of the system is accomplished with 20 transmission lines.

In Figure 8, the operating condition is such that the total load is increased by 10% from the base case. The fast and accurate calculation of the market power at any given operating point of the generator may serve or help the principal purposes of various bodies like the GENCOs, DISCOs and ISOs. The knowledge of market power for any given operating condition will enable these above bodies to take suitable actions to change their own operating points in order to fulfill their objectives. Figure 8 represents the interpolation plot for the case when load is increased by 10%.

Figure 8.

Normal case versus increase in load by 10%.

The interpolation is done between the five known operating points which are the maximum generation levels of Gen₃, starting from 60 MW in steps of 10 MW. For these operating points, the NMRS values are calculated. As we increase the maximum generation level of a particular generator, it affects the NMRS of its neighboring generators. In this case, we have monitored the NMRS of Gen₂ on bus 3 with the increase in the working range of Gen₃ to demonstrate this fact. As a result, the base case plot for Gen₂ tends to drop or go down as we keep increasing the generation level values for Gen₃. As the load is increased by 10%, the NMRS plot shifts up. Thus, it represents an important observation that as the load goes up at a certain time of day in a particular region, then the GENCOs which lie in such a subsystem will have an increased market power. Due to the piecewise polynomials which are attained due to the application of cubic’s spline interpolation between the operating points, GENCOs, DISCOs and ISOs can easily take suitable actions to not let anyone take undue advantage of the varying market power due to changes in operating conditions. Thus, it helps in attaining a zero market power. This work has been extended under various system conditions and an elaborate study is made in [10].

Inference: Market power reflects the amount of influence that a company has on the system in which it operates but in power systems Market power is the ability to maintain prices above the competitive levels for a significant period of time. Hence, it is of utmost importance to find the market power of the system under normal and abnormal condition that a system has to face.

7. Conclusions

In this chapter, a complete literature review is made on how cubic spline interpolation technique has been widely used in power systems application. ATC calculation, electric arc furnace modeling, static var compensation, voltage stability margin and market power estimation are some of the areas where cubic spline interpolation techniques are extensively used. The independent system operator (ISO), utilities and consumers can utilize this tool to predict the behaviour of the system/generation companies/utilities and bring back the system towards economic/system stability.