A New Approach to Nonsinusoidal Steady-State Power System Analysis

Abstract

A new analysis method using wavelet domain for steady-state operating condition of power system is developed and introduced. Based on wavelet-Galerkin theory, the system components such as resistor, inductor, capacitor, transmission lines, and switching devices are modeled in discrete wavelet domain for the purpose of steady-state analysis. To solve system equations, they are transferred to wavelet domain by forming algebraic matrix-vector relations using the wavelet transform coefficients and the equivalent circuit is thus built for system simulation. After describing the new algorithm, two-case studies are presented and compared with the simulations in the time domain to verify the accuracy and computational performance.

1. Introduction

New power systems include nonlinear, switching, and frequency dependent elements. An algorithm is required to calculate the periodic steady-state solutions of such systems. Different algorithms to this end have been developed by many researchers. These algorithms are classified, in terms of their formulation methodologies, into three categories: harmonic domain methods [1, 2], time domain methods [3, 4], and hybrid methods [5, 6]. In time domain method the nonlinearity and switching are modeled with ease, but frequency dependency of elements is a complicated concern. On the other hand, for the case of steady-state solution, the transient response of the system must be eliminated by adjusting the initial conditions [4]. Also when switching devices are modeled, especially in large systems with distributed elements, the step time of simulation must be decreased which leads to reduction of simulation speed. To solve this problem and to consider frequency dependency, the harmonic domain that relates to Fourier space is introduced and developed by many of researchers. In harmonic domain the system is studied discretely in different frequencies. On the other hand, the harmonic content of the system must be identified before any simulation is carried out. This means that not only a part of frequency content that is not known is eliminated, but also noninteger harmonics is neglected. Other problem that may occur is the weak homogenization of solution specially for switching devices.

To extract the complete spectrum of a physical signal, that contains harmonics of integer and noninteger orders, the main form of Fourier transform cannot be used. Hence, a modified form of this transform called Gabor’s transform is employed from which FFT algorithm is constructed. For signal analysis in steady-state condition, FFT method is adequate since it can extract the complete spectrum accurately. However, the kernel function of Gabor’s transform cannot construct a functional basis for power system simulation. Therefore, Fourier function-based methods cannot be used for power system analysis.

The wavelet analysis neither need to use a single window function in all frequency components, nor has linear resolution in the whole frequency domain, while these are essential and week points for Fourier analysis. Much of the interest regarding wavelet concentrates on time-frequency analysis. Power system analysis in multiresolution analysis (MRA) space is introduced by other authors in [7–9]. The speeds of solution described in these papers make them inefficient for numerical simulation and are hence considered impractical for a real power system with relevant large dimensions.

In this paper MRA space is used for power system simulation in nonsinusoidal and periodic conditions. Wavelet-Galerkin method guaranties the validation of this study from the mathematical point of view [10]. The FFT method can only be employed for signal analysis, but the proposed method can be used for spectrum analysis of a power system in less time as compared to other methods.

The paper is organized as follows. In Section 2 a brief description of mathematical theory of MRA is presented. Section 3 describes modeling of power system in the new suggested domain. In Section 4 the relationship between the new domain and spectral analysis is illustrated. Two-case studies are simulated in the new domain and the results of these simulations are compared with a time domain simulation in Section 5.

2. Mathematical Theory

2.1. Galerkin Method

The Galerkin method is one of the most reliable methods for finding numerical solution to differential equations [10]. Its simplicities make it perfect for many applications. The Galerkin approach is based on finding a functional basis for the solution space of the equation. It then projects the solution on the functional basis and minimizes the residual with respect to it. Standard polynomial basis or trigonometric basis is used for Galerkin method. However wavelets used to describe MRA space provide both time and frequency localization.

2.2. Multiresolution Analysis

In this section the orthonormal basis of compactly supported wavelets is reviewed briefly. The orthonormal basis of compactly supported wavelets of 𝐿²(ℜ) is formed by the dilation and translation of single function 𝜓(𝑥) [11, 12]:

𝜓_𝑗,𝑘(𝑥)=2^−𝑗/2𝜓(2^−𝑗𝑥−𝑘), (2.1) where 𝑗,𝑘∈𝑍, the function 𝜓(𝑥) has a companion, the scaling function 𝜑(𝑥), and these functions satisfy the following relations:

𝜑(𝑥)=√2^𝐽−1_𝑘=0ℎ_𝑘√𝜑(2𝑥−𝑘),𝜓(𝑥)=2^𝐽−1_𝑘=0𝑔_𝑘𝜑(2𝑥−𝑘). (2.2)

The coefficients 𝐻={ℎ_𝑘}^{𝑘=𝐽−1}_𝑘=0 and 𝐿={𝑔_𝑘}^{𝑘=𝐽−1}_𝑘=0 in (2.2) are quadrature mirror filters. The number 𝐽 of coefficients in (2.2) is related to the number of vanishing moments 𝑀. The wavelet basis induces an MRA on 𝐿²(ℜ), that is, the decomposition of Hilbert space into a chain of closed spaces:

⋯⊂𝑉₋₂⊂𝑉₋₁⊂𝑉₀⊂𝑉₁⊂𝑉₂⊂⋯, (2.3) such that

_𝑗∈𝑍𝑉_𝑗={0},^___𝑗∈𝑍𝑉_𝑗=𝐿²(ℜ). (2.4)

By defining 𝑊_𝑗 as an orthonormal complement of 𝑉_𝑗 in 𝑉_𝑗+1:

𝑉_𝑗+1=𝑉_𝑗⊕𝑊_𝑗, (2.5) with an MRA, one can use 𝜑_𝑛,𝑘(𝑥) and 𝜓_𝑛,𝑘(𝑥) as the basis functions for Galerkin method.

2.3. Wavelet-Galerkin Solution of a Periodic Problem

In the MRA space the numerical solution of a differential equation based on Wavelet-Galerkin method in the 𝑗th level can be written in this matrix form [13]

𝑆_𝑗𝑥=𝑓. (2.6)

The decomposition 𝑉_𝑗+1=𝑉_𝑗⊕𝑊_𝑗 allows the operator 𝑆_𝑗 to be split into four pieces (𝑊_𝑗 is called the wavelet space, i.e., the detail or fine-scale component of 𝑉_𝑗+1) which can be written as follows:

𝐴_𝑠𝑗𝐵_𝑠𝑗𝐶_𝑠𝑗𝑇_𝑠𝑗𝑑_𝑥𝑗𝑠_𝑥𝑗=𝑑_𝑓𝑗𝑠_𝑓𝑗, (2.7) where

𝐴_𝑆𝑗∶𝑊_𝑗→𝑊_𝑗,𝐵_𝑆𝑗∶𝑉_𝑗→𝑊_𝑗,𝐶_𝑆𝑗∶𝑊_𝑗→𝑉_𝑗,𝑇_𝑆𝑗∶𝑉_𝑗→𝑉_𝑗. (2.8)

And 𝑑_𝑥𝑗,𝑑_𝑓𝑗∈𝑊_𝑗, 𝑠_𝑥𝑗,𝑠_𝑓𝑗∈𝑉_𝑗 are the 𝐿²-orthonormal projections of 𝑥 and 𝑓 onto 𝑊_𝑗 and 𝑉_𝑗 spaces. The projection 𝑠_𝑥𝑗 is the coarse-scale component of the solution 𝑥, and 𝑑_𝑥𝑗 is the fine-scale component. To solve (2.7),

𝑅_𝑆𝑗=𝑇_𝑠𝑗−𝐶_𝑠𝑗𝐴⁻¹_𝑠𝑗𝐵_𝑠𝑗,𝑠_𝑥=𝑅⁻¹_𝑠𝑗(𝑠_𝑓−𝐶_𝑠𝑗𝐴⁻¹_𝑠𝑗𝑑_𝑓)𝑑_𝑥=𝐴⁻¹_𝑠𝑗(𝑑_𝑓−𝐵_𝑠𝑗𝑠_𝑥)., (2.9)

At this stage, 𝑇_𝑆𝑗 is selected and investigated. As the problem described above is periodic and supposing that the differential operator is equal to 𝑑^𝑚/𝑑𝑥^𝑚, the general form of 𝑇_𝑆𝑗 is

𝑇_𝑆𝑗=Ω₀^(𝑚)⋯Ω^(𝑚)_−𝑁+2⋯Ω^(𝑚)_𝑁−2⋯Ω₁^(𝑚)Ω₁^(𝑚)Ω₀^(𝑚)⋯Ω^(𝑚)_−𝑁+3⋯0⋯Ω₂^(𝑚)Ω⋮⋮⋯⋮⋯⋮⋮⋮_(𝑚)−1Ω^(𝑚)₋₂⋯0⋯Ω^(𝑚)_𝑁−3⋯Ω₀^(𝑚), (2.10) where

Ω^(𝑚)_𝑘−𝑙=^+∞_−∞𝜑_𝑗,𝑘𝑑(𝑥)^𝑚𝑑𝑡^𝑚(𝜑_𝑗,𝑘(𝑡))𝑑𝑥. (2.11)

The general forms of the other pieces of 𝑆_𝑗 are also similar to 𝑇_𝑆𝑗. For a circulant matrix such as 𝑇_𝑆𝑗, the eigenvalues 𝜆_𝛼 are [14]

𝜆_𝛼=^𝑁−2_{𝑘=−𝑁+2}Ω_𝑘^(𝑚)exp(−2𝜋𝑖𝛼𝑘/𝑛),𝛼=0,1,…,𝑛−1, (2.12) and the corresponding orthonormal eigenvectors 𝑣_𝛼 are

(𝑣_𝛼)_𝑘=(−1)^𝛼√𝑛exp(−2𝜋𝑖𝛼𝑘/𝑛),𝑘=0,1,…,𝑛−1. (2.13)

These relations lead to provision of quasidiagonal form of represented operators in MRA space without using conventional methods in a lesser time. Using diagonal form offers several advantages that are explained in the next parts. Using (2.12) and (2.13), (2.7) can be rewritten as

𝐴_𝑠𝑗𝐵_𝑠𝑗𝐶_𝑠𝑗𝑇_𝑠𝑗𝑑_𝑥𝑗𝑠̂_𝑥𝑗=𝑑_𝑓𝑗𝑠̂_𝑓𝑗, (2.14) where

𝑑_𝑥𝑗=Γ⁻¹𝑑_𝑥𝑗,𝑠̂_𝑥𝑗=Γ⁻¹𝑠_𝑥𝑗,𝑑_𝑓𝑗=Γ⁻¹𝑑_𝑓𝑗,𝑠̂_𝑓𝑗=Γ⁻¹𝑠_𝑓𝑗. (2.15)

In these equations Γ is the modal matrix. The columns of Γ are calculated using (2.13). 𝐴_𝑠𝑗, 𝐵_𝑠𝑗, 𝐶_𝑠𝑗, and 𝑇_𝑠𝑗 are diagonal matrices and their elements calculated by (2.12). So to calculate 𝑑_{𝑥𝑗,𝑖} and 𝑠̂_{𝑥𝑗,𝑖} (the 𝑖th values of 𝑑_𝑥𝑗 and 𝑠̂_𝑥𝑗), the following equation must be solved:

𝑎^𝑖_𝑠𝑗𝑏^𝑖_𝑠𝑗𝑐^𝑖_𝑠𝑗𝑡^𝑖_𝑠𝑗𝑑_{𝑥𝑗,𝑖}𝑠̂_{𝑥𝑗,𝑖}=𝑑_{𝑓𝑗,𝑖}𝑠̂_{𝑓𝑗,𝑖}. (2.16)

The volume of calculations is decreased significantly using the above technique. In other words, instead of calculating the inverse of matrices with 𝑁/2×𝑁/2 dimensions in (2.7); (2.16) is used for 𝑁/2 iterations. For problems with small dimensions this method seems not to be beneficial. However, as shown in the following sections, this approach could be useful for solving differential equations of large power systems. This is because in such systems the dimension of S_j in (2.6) is obtained by the multiplication of the system dimension and the number of considered samples (𝑁).

3. Power System Representation in the New Domain

3.1. Linear Elements Representation

The aim of this part is to obtain the expression for linear elements using mathematical operator representation in the new suggested domain. In this work, modeling in the MRA space has been carried out on the same basis as suggested by other researchers [7, 8]. For the purpose of wavelet domain study, the resistor, inductor, and capacitor models are set up in the following section. ‎

3.1.1. Resistor

The relationship between voltage and current of a resistor 𝑟 in the time domain is described as 𝑣(𝑡)=𝑟𝑖(𝑡). (3.1)

Then, this relation is expressed in a wavelet expansion as

𝑉_𝑑𝑗𝑉_𝑠𝑗𝐼=𝑟𝑈00𝑈𝐼^𝑑𝑗_𝑠𝑗, (3.2) where 𝑈‎ is‎ an identity ‎matrix ‎with 𝑁/2^{(𝐽max−𝑗+1)} dimensions, 𝑁 is signal length, and 𝐽 is resolution‎‎level‎.

3.1.2. Inductor

The relationship between voltage and current of an inductor ℓ is

𝑣(𝑡)=ℓ𝑑𝑖(𝑡)𝑑𝑡. (3.3)

The 𝑁 point discretization of (3.3) leads to

𝑉=ℓ[𝐷_𝑇]𝐼, (3.4) where 𝐷_𝑇 is the discrete form of derivative operator.

In the 𝑗th level of MRA space, (3.4) can be written as

𝑉_𝑑𝑗𝑉_𝑠𝑗=ℓ⋅𝑊𝐷_𝑇𝐼_𝑑𝑗𝐼_𝑠𝑗, (3.5) where

𝑊𝐷_𝑇=𝐻𝐷_𝑇𝐻𝐻𝐷_𝑇𝐿𝐿𝐷_𝑇𝐻𝐿𝐷_𝑇𝐿=𝐴_𝑠𝑗𝐵_𝑠𝑗𝐶_𝑠𝑗𝑇_𝑠𝑗. (3.6)

To transfer (3.5) from the highest level (the finest scale) to the next lower level (coarser scale) and, respectively, in a hierarchical form to other levels (scales) of MRA space, 𝐷_𝑇 is substituted with 𝐿𝐷_𝑇𝐿 of the higher resolution level. Of course in each subsequent level the dimensions of matrices will be different from previous ones and its magnitude is divided by 2. The submatrices of 𝑊𝐷_𝑇 have a circulant form and this feature is specific to all orders of derivative operator in MRA space. Also, the Γ matrix and eigenvectors are the same for all orders in each level of MRA space. Rewriting (3.5) using (2.12) and (2.13) leads to obtain a quasidiagonal form as follows:

𝑉_𝑑𝑗𝑉_𝑠𝑗=ℓ⋅𝑊𝐷_𝑇𝐼̂_𝑑𝑗𝐼̂_𝑠𝑗,𝑊𝐷_𝑇=𝐻𝐻𝐻𝐿𝐿𝐻𝐿𝐿, (3.7) where ℓ⋅𝑊𝐷_𝑇 represents‎  ‎the  wavelet  ‎domain‎  ‎impedance  ‎of  inductor. There are four submatrices for impedance definition of inductor, the first submatrix (ℓ𝐻𝐻 deals with 𝑊_𝑗→𝑊_𝑗) belongs to high frequency part of level 𝑗. Also the fourth submatrix (ℓ𝐿𝐿 relates to 𝑉_𝑗→𝑉_𝑗) represents the impedance in low frequency part.

3.1.3. Capacitor

There is a time domain relationship between the voltage and current of a capacitor c as represented by

1𝑣(𝑡)=𝑐𝑖(𝑡)𝑑𝑡. (3.8)

Projecting the above equation onto discrete time domain leads to

𝑉=1𝑐[𝐷_𝑇]⁻¹𝐼, (3.9) where 𝐷_𝑇⁻¹ is the discrete form of integral operator in periodic conditions. Thus, in the wavelet domain, 𝐷_𝑇⁻¹ is expressed as‎

𝑊𝐷^′_𝑇=𝐻𝐻^′𝐻𝐿^′𝐿𝐻^′𝐿𝐿^′. (3.10)

To compute 𝑖th value of 𝐻𝐻^′,𝐻𝐿^′,… this relation is used

ℎℎ^′_𝑖ℎ𝑙^′_𝑖𝑙ℎ^′_𝑖𝑙𝑙^′_𝑖=ℎℎ_𝑖ℎ𝑙_𝑖𝑙ℎ_𝑖𝑙𝑙_𝑖⁻¹, (3.11) where ℎℎ_𝑖, ℎ𝑙_𝑖, 𝑙ℎ_𝑖, and 𝑙𝑙_𝑖 are the 𝑖th values of 𝐻𝐻,𝐻𝐿,….

3.2. Transmission Line Modeling in the New Domain

There are many papers about transmission line modeling for transient studies [15–17]. In this section, the distributed modeling of single phase line for power system studies in the new suggested domain is discussed briefly.

The V-I characteristic of a differential element of transmission line in continuous time domain is represented by‎

−𝜕𝑣𝜕𝑥=𝑟𝑖+ℓ𝜕𝑖−𝜕𝑡,𝜕𝑖𝜕𝑥=𝑔𝑣+𝑐𝜕𝑣𝜕𝑡, (3.12)(3.13) where 𝑟, ℓ, 𝑔, and 𝑐 are resistance, inductance, conductance, and capacitance of the differential element of transmission line, respectively. If (3.12) is replaced in (3.13), we have

𝜕²𝑣𝜕𝑥²=𝑟₁⋅𝑣+𝑟₂⋅𝜕𝑣𝜕𝑡+𝑟₃⋅𝜕²𝑣𝜕𝑡², (3.14) where

𝑟₁=𝑟⋅𝑔,𝑟₂=𝑟⋅𝑐+ℓ⋅g,𝑟₃=ℓ⋅c. (3.15) If (3.14) is transformed to the new domain for 𝑖th element of 𝑗th level, the following relation can be obtained:

𝜕²𝑉_{𝑑𝑗,𝑖}𝜕𝑥²𝜕²𝑉_{𝑠𝑗,𝑖}𝜕𝑥²=𝑟₁+𝑟₂ℎℎ_𝑖+𝑟₃ℎℎ²_𝑖𝑟₂ℎ𝑙_𝑖+𝑟₃ℎ𝑙²_𝑖𝑟₂𝑙ℎ_𝑖+𝑟₃𝑙ℎ²_𝑖𝑟₁+𝑟₂𝑙𝑙_𝑖+𝑟₃𝑙𝑙²_𝑖𝑉_{𝑑𝑗,𝑖}𝑉_{𝑠𝑗,𝑖}, (3.16) where ℎℎ_𝑖,…,𝑙𝑙_𝑖 and ℎℎ²_𝑖,…,𝑙𝑙²_𝑖 are the 𝑖th diagonal elements of 𝐻𝐻,…,𝐿𝐿 and 𝐻𝐻²,…,𝐿𝐿², respectively. Also, 𝐻𝐻²=𝐻⋅𝐷²_𝑇⋅𝐻^𝑇,…, 𝐿𝐿²=𝐿⋅𝐷²_𝑇⋅𝐿^𝑇 and 𝐷²_𝑇 is the disceretized form of second-order derivative operator. In the above method, only the distribution of the line parameters is considered. In modeling of transmission line with frequency dependency, 𝑟₁, 𝑟₂, and 𝑟₃ in (3.16) are not scalars, as they are of matrix form. In the new domain they are diagonal matrices. Each element of these matrices belongs to a special frequency whose details are expressed in Section 4. Based on the relation between new domain and spectral analysis, the parameter adjustments for these frequency dependent matrices are performed.

3.3. Switching Devices Modeling

In this part, modeling method for switching devices is investigated. These devices are the main sources of harmonics in power network. Modeling of these devices is explained by many authors for harmonic studies [18–22]. Since wavelet makes a local analysis instead of a general analysis, modeling of switching devices in the new domain can be facilitated. Assume that a linear load is connected to network in series with a power electronic switch. The relation between voltage and current of load without considering the switch is

𝑖(𝑡)=𝑓(𝑣(𝑡)), (3.17) where 𝑓 is a linear operator. As load is in series with the switch, the relation between current and voltage is

𝑖(𝑡)=𝑝(𝑡)⋅𝑓(𝑝(𝑡)⋅𝑣(𝑡)), (3.18) where 𝑝(𝑡) is switching signal. Switching signal is a periodic function defined as follows:

𝑝(𝑡)=1∶switchison,0∶switchisoff. (3.19)

Discretizing equation (3.18) leads to

𝐼=[𝑆]⋅𝑉,=[𝑃]⋅[𝐹]⋅[𝑃]. (3.20)(3.21)

This relation is obtained based on this fact that mathematical operator 𝑓 is linear. It is not necessary to suppose that the switching device is synchronized with power system frequency. Transferring (3.21) to the new suggested domain does not result a diagonal matrix. This refers to existence of cross-couplings between harmonics. As the transferred matrix is not diagonal, using this matrix directly in the network equation increases the computational volume which leads to reduced efficiency in the numerical solution. To avoid this, the matrix is not considered in admittance matrix and the switching device is modeled as a voltage dependent current source. Therefore, simulation at each level is carried out without considering the admittance of switching device. Then, using the voltage of switching device node that is obtained from the simulation and admittance matrix obtained from (3.20), the current of switching device branch is calculated. For the first iteration this current is not exact. To have an exact solution this current is used for next iteration and the switching device will be modeled as a current source. This process is repeated until the solution homogenizes to a certain value.

3.4. Network Representation

To develop this method for power network simulation, a simple circuit is considered (see Figure 1). Applying the KCL relation yields

Figure 1 

Parallel RLC circuit.

−𝑖_𝑥+𝑣_in𝑟𝑑+𝑐𝑣𝑑𝑡ⁱⁿ+1ℓ𝑣_in𝑑𝑡=0. (3.22)

According to the modified nodal method that is used in harmonic analysis:

1𝑟1+𝑐𝑝+1𝑣ℓ𝑝10𝑖¹_𝑥=0𝑣_in, (3.23) where 𝑝 is derivative operator, (3.23) could also be written as follows:

(𝐶_𝑅+𝑝𝐶_𝐶+1𝑝𝐶_𝐿)𝑣₁𝑖_𝑥=0𝑣_in, (3.24) where

𝐶_𝑅=1𝑟110,𝐶_𝐿=1ℓ000,𝐶_𝑐=𝑐000. (3.25)

C_R, 𝐶_𝑐, and C_L are resistive coefficients matrix, capacitive coefficients matrix, and inductive coefficients matrix, respectively. These matrices can be defined according to the modified nodal method. In the new domain at the 𝑗th level, the general form of (3.24) is the same as (2.14), where

𝑑_𝑓𝑗=0Γ⁻¹⋅𝑉_in,𝑑𝑗,𝑠̂_𝑓𝑗=0Γ⁻¹⋅𝑉_in,𝑠𝑗,𝑑_𝑥𝑗=Γ⁻¹⋅𝑉_1,𝑑𝑗Γ⁻¹⋅𝐼_in,𝑑𝑗,𝑠̂_𝑥𝑗=Γ⁻¹⋅𝑉_1,𝑠𝑗Γ⁻¹⋅𝐼_in,𝑠𝑗. (3.26)

Now according to Section 2 these equations are obtained

𝑎_{𝑠𝑗,𝑖}=𝐶_𝑅+ℎℎ_𝑖⋅𝐶_𝐶+ℎℎ^′_𝑖⋅𝐶_𝐿,𝑏_{𝑠𝑗,𝑖}=ℎ𝑙_𝑖⋅𝐶_𝐶+ℎ𝑙^′_𝑖⋅𝐶_𝐿,𝑐_{𝑠𝑗,𝑖}=𝑙ℎ_𝑖⋅𝐶_𝐶+𝑙ℎ^′_𝑖⋅𝐶_𝐿,𝑡_{𝑠𝑗,𝑖}=𝐶_𝑅+𝑙𝑙_𝑖⋅𝐶_𝐶+𝑙𝑙^′_𝑖⋅𝐶_𝐿. (3.27)

The formula (3.27) could be written for any network directly. Using these matrices, the 𝑖th value of response vectors can be computed

𝑟_{𝑠𝑗,𝑖}=𝑡_{𝑠𝑗,𝑖}−𝑐_{𝑠𝑗,𝑖}𝑎⁻¹_{𝑠𝑗,𝑖}𝑏_{𝑠𝑗,𝑖},𝑠̂_{𝑥𝑗,𝑖}=𝑟⁻¹_{𝑠𝑗,𝑖}(𝑠̂_{𝑓𝑗,𝑖}−𝑐_{𝑠𝑗,𝑖}𝑎⁻¹_{𝑠𝑗,𝑖}𝑑_{𝑓𝑗,𝑖}𝑑),_{𝑥𝑗,𝑖}=𝑎⁻¹_{𝑠𝑗,𝑖}(𝑑_{𝑓𝑗,𝑖}−𝑏_{𝑠𝑗,𝑖}𝑠̂_{𝑥𝑗,𝑖}), (3.28) where 𝑑_{𝑥𝑗,𝑖} and 𝑠̂_{𝑥𝑗,𝑖} are the 𝑖th values of 𝑑_𝑥𝑗 and 𝑠̂_𝑥𝑗. To obtain 𝑗th level of response vector in MRA space, that is, 𝑑_𝑥𝑗 and 𝑑_𝑥𝑗, Γ is multiplied to 𝑑_𝑥𝑗 and 𝑠̂_𝑥𝑗, respectively.

The steps for Nonsinusoidal steady-state analysis are as follows.

(1)Determine number of levels (𝐽_max) and number of samples (𝑁), where 𝐽_max is the index of finest scale in MRA space.(2)Calculate the 𝐶_𝑅, 𝐶_𝐶, and 𝐶_𝐿 matrices according to the modified nodal method.(3)Set 𝐽=1. (4)Set 𝑗=𝐽_max−𝐽, where 𝑗 is the index of current resolution level. (5)Compute 𝐷_𝑇, 𝐻𝐷_𝑇𝐻,𝐻𝐷_𝑇𝐿,… for 𝑗th level. (6)Compute the Γ matrix using (2.13) and then transfer the input vector to the new domain. (7)Set 𝑖=1.(8)Calculate 𝑎_{𝑠𝑗,𝑖},𝑏_{𝑠𝑗,𝑖}, 𝑐_{𝑠𝑗,𝑖}, and 𝑡_{𝑠𝑗,𝑖} using (3.27). Then by using (3.28), calculate 𝑑_{𝑥𝑗,𝑖} and 𝑠̂_{𝑥𝑗,𝑖}.(9)If 𝑖=𝑁/2^𝐽 then set 𝑖=𝑖+1 and go to step (8).(10)Calculate the response vector in MRA space for 𝑗th resolution level, (i.e., 𝑑_𝑥𝑗 and 𝑠_𝑥𝑗). (11)If 𝐽 is not equal to 𝐽_max, then set 𝐽=𝐽+1 and go to step (5).(12)End.

Figure 2 shows the flowchart of nonsinusoidal steady-state analysis in the proposed domain.

Figure 2 

Flowchart of nonsinusoidal steady-state analysis.

4. The New Domain and Spectral Analysis

Spectral analysis is the most significant aim of power quality estimation in electrical power networks [23, 24]. Transferring the simulation results from the new proposed domain to time domain seems to be unnecessary step for obtaining harmonic information. Required information such as THD and harmonic amplitudes could be extracted directly from the results in the new suggested domain which has the advantage of less time being consumed for simulation.

If a column of Γ matrix is multiplied by a vector of wavelet coefficients of a signal such as 𝑑_𝑓𝑗, then

[(𝑣_𝛼)_𝑘]^𝑇⋅𝑑_𝑓𝑗=(−1)^𝛼√𝑛^𝑁−1_𝑘=0𝑑_{𝑓𝑗,𝑘}exp(−𝑗⋅2⋅𝜋⋅𝛼⋅𝑘𝑛), (4.1) where (𝑣_𝛼)_𝑘 is the 𝛼th column of Γ introduced in (2.13), and 𝑑_{𝑓𝑗,𝑘} is the 𝑘th element of 𝑑_𝑓𝑗. Equation (4.1) shows that multiplication by a column of Γ returns the DFT of 𝑑_𝑓𝑗. In extracting harmonic contents of a typical signal, in the coarsest scale, the first element of 𝑠̂_𝑓𝑗 represents the DC content of the signal which could easily be proven mathematically. For other harmonic contents, the number of considered periods has an important role in identifying the element representing the value of a special harmonic. If only one period of the original signal is considered, the second element of 𝑠̂_𝑓𝑗 from the coarsest scale (lowest level) represents the content of the main harmonic. With respect to the frequency band, (𝑛+1)th element of 𝑠̂_𝑓𝑗 and 𝑑_𝑓𝑗 in each level represents the contents of 𝑛th harmonic order.

To calculate THD value of a signal from its multiplication by Γ matrix, we define these relations as follows:

𝐴_𝑑,𝑗=1𝑁^𝑁/𝑐_𝑘=2|2⋅𝑑_𝑗,𝑘|²,𝑗=0,…,𝐽_max;𝑐=2^{(𝐽max−𝑗+2)},𝐴_𝑐,0=1𝑁^𝑁/𝑐_𝑘=2|2⋅𝑠̂_𝑗,𝑘|²𝑗=0,THD=(∑^𝐽max_𝑗=0𝐴_𝑑,𝑗²)+𝐴²_𝑐,0𝑆₁, (4.2) where 𝑁 is the number of original signal samples, and 𝑆₁ is the amplitude value of the main harmonic.

5. Case Studies

To verify the accuracy and computational performance, the proposed method is demonstrated for two-case studies and the results are compared with the time domain simulation results.

5.1. Case Study 1

The periodic steady-state solution of the test network shown in Figure 3 is calculated by the proposed algorithm. This network is a test system for Harmonics Modeling and Simulation [25], where a harmonic source is located on the bus 49-RECT as ASD load. Disregarding the frequency dependency, generally a simple model is used for lines and transformers. The test system is connected to a larger plant from 100:Util-69 bus. An equivalent model is hence obtained for the larger plant from the fault MVA level. The system consists of a local generator, modeled by a voltage source in series with a subtransient impedance. The whole system is therefore transformed into the new suggested domain where THDs are computed directly from results. Figure 4 shows the 𝐼_Util waveforms obtained from the time and new domain simulations for which there is a perfect overlap. The consumed time for simulation of test network was found to be 4.7 seconds. Simulation was coded with MATLAB version 7, using 512 samples, 3 levels for MRA space via a personal computer with Pentium 4 CPU (2.8 GHz) and 512 RAM. Tables 1 and 2 compare THDs for the new proposed method and those obtained from SIMULINK time domain simulations.

Table 1 

THD percentage values.

Table 2 

THD percentage values.

Figure 3 

The 13-bus test system.

Figure 4 

𝐼Util waveform comparison of the time and new domain simulations.

5.2. Case Study 2

To examine the validation of the proposed models for the distributed transmission line and switching devices mentioned earlier, a 132 KV test system is selected and simulated using the new suggested domain. This simple system consists of two transmission lines with 50 Km length and a switching load (see Figure 5). The associated test system parameters are listed in Table 3. The switching load is a SVC which has a thyristor controlled reactor (TCR) with the firing angle of 130^˚ and a fixed capacitor (FC). In Table 4 the THD percentages and the consumed times for this simulation are compared. The THDs are calculated directly from the new proposed domain. db4 wavelet function is used in this simulation. As far as the simulation time is concerned, the new proposed domain reduced this almost by half from 3.44 seconds to 1.85 seconds. This can further be reduced using fast numerical algorithms, especially in transferring input vectors and operators to the new domain which was not considered in this work.

Table 3 

Test system parameters.

Table 4 

Comparing of THD percentages and simulation times.

Figure 5 

Test system.

Figures 6 and 7 show the simulation results from time domain and the new suggested domain for TCR current and load voltage, respectively.

Figure 6 

TCR current: new proposed domain and time domain.

Figure 7 

Load Voltage (𝑉3): new proposed domain and time domain.

6. Conclusions

This paper describes a new approach based on MRA space for the Nonsinusoidal Steady-‎State Power System Analysis.‎ By applying operator representation theory, the system components ‎such as resistor, ‎inductor, capacitor, transmission lines, and switching devices are modeled in the wavelet domain. The model of switching device is based on switching signal while the interaction between network and switching device is also considered. Discrete nature of the model, easy adoptability ‎for nonlinear, and frequency dependent components are the main advantages of the proposed modeling technique.

Simulation results confirm the effectiveness and accuracy of the proposed system model and analysis scheme. The proposed method might well be applied to several fields including power quality analysis and power system protection.‎