# A New Approach to Nonsinusoidal Steady-State Power System Analysis

Authors: A. Ketabi; M. Khoshkholgh; R. Feuillet
Journal: Mathematical Problems in Engineering
 (2009)

---

## 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.


---

## Body

## 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 Section2 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 ofL2(ℜ) is formed by the dilation and translation of single function ψ(x) [11, 12]:(2.1)ψj,k(x)=2-j/2ψ(2-jx-k),
where j,k∈Z, the function ψ(x) has a companion, the scaling function φ(x), and these functions satisfy the following relations:(2.2)φ(x)=2∑k=0J-1hkφ(2x-k),ψ(x)=2∑k=0J-1gkφ(2x-k).The coefficientsH={hk}k=0k=J-1 and L={gk}k=0k=J-1 in (2.2)   are quadrature mirror filters. The number J of coefficients in (2.2)    is related to the number of vanishing moments M. The wavelet basis induces an MRA on  L2(ℜ), that is, the decomposition of Hilbert space into a chain of closed spaces:(2.3)⋯⊂V-2⊂V-1⊂V0⊂V1⊂V2⊂⋯,
such that(2.4)⋂j∈ZVj={0},⋃j∈Z__Vj=L2(ℜ).By definingWj as an orthonormal complement of  Vj in Vj+1:(2.5)Vj+1=Vj⊕Wj,
with an MRA, one can use φn,k(x) and ψn,k(x) 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 thejth level can be written in this matrix form [13](2.6)Sjx=f.The decompositionVj+1=Vj⊕Wj allows the operator Sj to be split into four pieces (Wj is called the wavelet space, i.e., the detail or fine-scale component of Vj+1) which can be written as follows:(2.7)(AsjBsjCsjTsj)(dxjsxj)=(dfjsfj),
where(2.8)ASj:Wj→Wj,BSj:Vj→Wj,CSj:Wj→Vj,TSj:Vj→Vj.Anddxj,dfj∈Wj,  sxj,sfj∈Vj are the L2-orthonormal projections of  x and f onto Wj and Vj spaces. The projection sxj is the coarse-scale component of the solution x, and  dxj is the fine-scale component. To solve (2.7),(2.9)RSj=Tsj-CsjAsj-1Bsj,sx=Rsj-1(sf-CsjAsj-1df)dx=Asj-1(df-Bsjsx).,At this stage,TSj is selected and investigated. As the problem described above is periodic and supposing that the differential operator is equal to dm/dxm, the general form of  TSj is(2.10)TSj=(Ω0(m)⋯Ω-N+2(m)⋯ΩN-2(m)⋯Ω1(m)Ω1(m)Ω0(m)⋯Ω-N+3(m)⋯0⋯Ω2(m)⋮⋮⋯⋮⋯⋮⋮⋮Ω-1(m)Ω-2(m)⋯0⋯ΩN-3(m)⋯Ω0(m)),
where(2.11)Ωk-l(m)=∫-∞+∞φj,k(x)dmdtm(φj,k(t))dx.The general forms of the other pieces ofSj are also similar to TSj. For a circulant matrix such as TSj, the eigenvalues λα are [14](2.12)λα=∑k=-N+2N-2Ωk(m)exp(-2πiαk/n),α=0,1,…,n-1,
and the corresponding orthonormal eigenvectors vα are(2.13)(vα)k=(-1)αnexp(-2πiαk/n),k=0,1,…,n-1.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)(Asj¯Bsj¯Csj¯Tsj¯)(d̂xjŝxj)=(d̂fjŝfj),
where(2.15)d̂xj=Γ-1dxj,ŝxj=Γ-1sxj,d̂fj=Γ-1dfj,ŝfj=Γ-1sfj.In these equationsΓ is the modal matrix. The columns of  Γ are calculated using (2.13). Asj¯, Bsj¯, Csj¯, and Tsj¯ are diagonal matrices and their elements calculated by (2.12). So to calculate d̂xj,i and ŝxj,i (the ith values of d̂xj and ŝxj), the following equation must be solved:(2.16)(asjibsjicsjitsji)(d̂xj,iŝxj,i)=(d̂fj,iŝfj,i).The volume of calculations is decreased significantly using the above technique. In other words, instead of calculating the inverse of matrices withN/2×N/2 dimensions in (2.7);   (2.16) is used for N/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 Sj in (2.6) is obtained by the multiplication of the system dimension and the number of considered samples (N).


## 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 ofL2(ℜ) is formed by the dilation and translation of single function ψ(x) [11, 12]:(2.1)ψj,k(x)=2-j/2ψ(2-jx-k),
where j,k∈Z, the function ψ(x) has a companion, the scaling function φ(x), and these functions satisfy the following relations:(2.2)φ(x)=2∑k=0J-1hkφ(2x-k),ψ(x)=2∑k=0J-1gkφ(2x-k).The coefficientsH={hk}k=0k=J-1 and L={gk}k=0k=J-1 in (2.2)   are quadrature mirror filters. The number J of coefficients in (2.2)    is related to the number of vanishing moments M. The wavelet basis induces an MRA on  L2(ℜ), that is, the decomposition of Hilbert space into a chain of closed spaces:(2.3)⋯⊂V-2⊂V-1⊂V0⊂V1⊂V2⊂⋯,
such that(2.4)⋂j∈ZVj={0},⋃j∈Z__Vj=L2(ℜ).By definingWj as an orthonormal complement of  Vj in Vj+1:(2.5)Vj+1=Vj⊕Wj,
with an MRA, one can use φn,k(x) and ψn,k(x) 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 thejth level can be written in this matrix form [13](2.6)Sjx=f.The decompositionVj+1=Vj⊕Wj allows the operator Sj to be split into four pieces (Wj is called the wavelet space, i.e., the detail or fine-scale component of Vj+1) which can be written as follows:(2.7)(AsjBsjCsjTsj)(dxjsxj)=(dfjsfj),
where(2.8)ASj:Wj→Wj,BSj:Vj→Wj,CSj:Wj→Vj,TSj:Vj→Vj.Anddxj,dfj∈Wj,  sxj,sfj∈Vj are the L2-orthonormal projections of  x and f onto Wj and Vj spaces. The projection sxj is the coarse-scale component of the solution x, and  dxj is the fine-scale component. To solve (2.7),(2.9)RSj=Tsj-CsjAsj-1Bsj,sx=Rsj-1(sf-CsjAsj-1df)dx=Asj-1(df-Bsjsx).,At this stage,TSj is selected and investigated. As the problem described above is periodic and supposing that the differential operator is equal to dm/dxm, the general form of  TSj is(2.10)TSj=(Ω0(m)⋯Ω-N+2(m)⋯ΩN-2(m)⋯Ω1(m)Ω1(m)Ω0(m)⋯Ω-N+3(m)⋯0⋯Ω2(m)⋮⋮⋯⋮⋯⋮⋮⋮Ω-1(m)Ω-2(m)⋯0⋯ΩN-3(m)⋯Ω0(m)),
where(2.11)Ωk-l(m)=∫-∞+∞φj,k(x)dmdtm(φj,k(t))dx.The general forms of the other pieces ofSj are also similar to TSj. For a circulant matrix such as TSj, the eigenvalues λα are [14](2.12)λα=∑k=-N+2N-2Ωk(m)exp(-2πiαk/n),α=0,1,…,n-1,
and the corresponding orthonormal eigenvectors vα are(2.13)(vα)k=(-1)αnexp(-2πiαk/n),k=0,1,…,n-1.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)(Asj¯Bsj¯Csj¯Tsj¯)(d̂xjŝxj)=(d̂fjŝfj),
where(2.15)d̂xj=Γ-1dxj,ŝxj=Γ-1sxj,d̂fj=Γ-1dfj,ŝfj=Γ-1sfj.In these equationsΓ is the modal matrix. The columns of  Γ are calculated using (2.13). Asj¯, Bsj¯, Csj¯, and Tsj¯ are diagonal matrices and their elements calculated by (2.12). So to calculate d̂xj,i and ŝxj,i (the ith values of d̂xj and ŝxj), the following equation must be solved:(2.16)(asjibsjicsjitsji)(d̂xj,iŝxj,i)=(d̂fj,iŝfj,i).The volume of calculations is decreased significantly using the above technique. In other words, instead of calculating the inverse of matrices withN/2×N/2 dimensions in (2.7);   (2.16) is used for N/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 Sj in (2.6) is obtained by the multiplication of the system dimension and the number of considered samples (N).


## 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 resistorr in the time domain is described as
(3.1)v(t)=ri(t).Then, this relation is expressed in a wavelet expansion as(3.2)(VdjVsj)=r(U00U)(IdjIsj),
where U‎ is‎   an identity ‎matrix ‎with N/2(Jmax-j+1) dimensions,  N is signal length, and J is resolution‎‎level‎.


#### 3.1.2. Inductor

The relationship between voltage and current of an inductorℓ  is(3.3)v(t)=ℓdi(t)dt.TheN point discretization of (3.3) leads to(3.4)V̲=ℓ[DT]I̲,
where DT is the discrete form of derivative operator.In thejth level of MRA space, (3.4) can be written as(3.5)(VdjVsj)=ℓ·WDT(IdjIsj),
where(3.6)WDT=(HDTHHDTLLDTHLDTL)=(AsjBsjCsjTsj).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, DT is substituted with LDTL 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 WDT 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)(V̂djV̂sj)=ℓ·WDT¯(ÎdjÎsj),WDT¯=(HH¯HL¯LH¯LL¯),
where ℓ·WDT¯ represents‎  ‎the  wavelet  ‎domain‎  ‎impedance  ‎of  inductor. There are four submatrices for impedance definition of inductor, the first submatrix (ℓHH¯ deals with Wj→Wj) belongs to high frequency part of level j. Also the fourth submatrix (ℓLL¯ relates to Vj→Vj) represents the impedance in low frequency part.


#### 3.1.3. Capacitor

There is a time domain relationship between the voltage and current of a capacitorc as represented by(3.8)v(t)=1c∫i(t)dt.Projecting the above equation onto discrete time domain  leads to(3.9)V̲=1c[DT]-1I̲,
where DT-1  is the discrete form of integral operator in periodic conditions. Thus, in the wavelet domain, DT-1 is expressed as‎(3.10)WDT′̅=(HH′¯HL′¯LH′¯LL′¯).To computeith value of  HH′¯,HL′¯,…  this relation is used(3.11)(hhi′hli′lhi′lli′)=(hhihlilhilli)-1,
where hhi, hli, lhi, and lli are the ith values of HH¯,HL¯,….


### 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)-∂v∂x=ri+ℓ∂i∂t,(3.13)-∂i∂x=gv+c∂v∂t,
where r, ℓ, g, and c 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)∂2v∂x2=r1·v+r2·∂v∂t+r3·∂2v∂t2,
where(3.15)r1=r·g,r2=r·c+ℓ·g,r3=ℓ·c.
If (3.14) is transformed to the new domain for ith element of jth level, the following relation can be obtained:(3.16)(∂2V̂dj,i∂x2∂2V̂sj,i∂x2)=(r1+r2hhi+r3hhi2r2hli+r3hli2r2lhi+r3lhi2r1+r2lli+r3lli2)(Vdj,iVsj,i),
where hhi,…,lli and hhi2,…,lli2 are the ith diagonal elements of HH¯,…,LL¯ and HH2¯,…,LL2¯, respectively. Also, HH2=H·DT2·HT,…, LL2=L·DT2·LT and  DT2 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, r1, r2, and r3 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)i(t)=f(v(t)),
where f is a linear operator. As load is in series with the switch, the relation between current and voltage is(3.18)i(t)=p(t)·f(p(t)·v(t)),
where p(t) is switching signal. Switching signal is a periodic function defined as follows:(3.19)p(t)={1:switchison,0:switchisoff.Discretizing equation (3.18) leads to(3.20)I̲=[S]·V̲,(3.21)[S]=[P]·[F]·[P].This relation is obtained based on this fact that mathematical operatorf 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 Figure1). Applying the  KCL relation  yieldsFigure 1
Parallel RLC circuit.(3.22)-ix+vinr+cddtvin+1ℓ∫vindt=0.According to the modified nodal method that is used in harmonic analysis:(3.23)(1r+cp+1ℓp110)(v1ix)=(0vin),
where p is derivative operator,    (3.23) could also be written as follows:(3.24)(CR+pCC+1pCL)(v1ix)=(0vin),
where(3.25)CR=(1r110),CL=(1ℓ000),Cc=(c000).CR, Cc, and CL 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 jth level, the general form of (3.24) is the same as (2.14), where(3.26)d̂fj=(0Γ-1·Vin,dj),ŝfj=(0Γ-1·Vin,sj),d̂xj=(Γ-1·V1,djΓ-1·Iin,dj),ŝxj=(Γ-1·V1,sjΓ-1·Iin,sj).Now according to Section2 these equations are obtained(3.27)asj,i=CR+hhi·CC+hhi′·CL,bsj,i=hli·CC+hli′·CL,csj,i=lhi·CC+lhi′·CL,tsj,i=CR+lli·CC+lli′·CL.The formula (3.27) could be written for any network directly. Using these matrices, the ith value of response vectors can be computed(3.28)rsj,i=tsj,i-csj,iasj,i-1bsj,i,ŝxj,i=rsj,i-1(ŝfj,i-csj,iasj,i-1d̂fj,i),d̂xj,i=asj,i-1(d̂fj,i-bsj,iŝxj,i),
where d̂xj,i and ŝxj,i are the ith values of d̂xj and ŝxj. To obtain jth  level of response vector in MRA space, that is, dxj and dxj,  Γ is multiplied to d̂xj and ŝxj, respectively.The steps for Nonsinusoidal steady-state analysis are as follows.(1)
Determine number of levels (Jmax) and number of samples (N), where Jmax is the index of finest scale in MRA space.(2)
Calculate theCR, CC, and CL matrices according to the modified nodal method.(3)
SetJ=1.(4)
Setj=Jmax-J, where j is the index of current resolution level.(5)
ComputeDT, HDTH,HDTL,… for jth level.(6)
Compute theΓ matrix using (2.13) and then transfer the input vector to the new domain.(7)
Seti=1.(8)
Calculateasj,i,bsj,i, csj,i, and tsj,i using (3.27). Then by using (3.28), calculate d̂xj,i and  ŝxj,i.(9)
Ifi=N/2J then set i=i+1 and go to step (8).(10)
Calculate the response vector in MRA space forjth resolution level, (i.e., dxj  and sxj).(11)
IfJ is not equal to Jmax, then set J=J+1 and go to step (5).(12)
End.Figure2 shows the flowchart of nonsinusoidal steady-state analysis in the proposed domain.Figure 2
Flowchart of nonsinusoidal steady-state analysis.


## 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 resistorr in the time domain is described as
(3.1)v(t)=ri(t).Then, this relation is expressed in a wavelet expansion as(3.2)(VdjVsj)=r(U00U)(IdjIsj),
where U‎ is‎   an identity ‎matrix ‎with N/2(Jmax-j+1) dimensions,  N is signal length, and J is resolution‎‎level‎.


### 3.1.2. Inductor

The relationship between voltage and current of an inductorℓ  is(3.3)v(t)=ℓdi(t)dt.TheN point discretization of (3.3) leads to(3.4)V̲=ℓ[DT]I̲,
where DT is the discrete form of derivative operator.In thejth level of MRA space, (3.4) can be written as(3.5)(VdjVsj)=ℓ·WDT(IdjIsj),
where(3.6)WDT=(HDTHHDTLLDTHLDTL)=(AsjBsjCsjTsj).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, DT is substituted with LDTL 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 WDT 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)(V̂djV̂sj)=ℓ·WDT¯(ÎdjÎsj),WDT¯=(HH¯HL¯LH¯LL¯),
where ℓ·WDT¯ represents‎  ‎the  wavelet  ‎domain‎  ‎impedance  ‎of  inductor. There are four submatrices for impedance definition of inductor, the first submatrix (ℓHH¯ deals with Wj→Wj) belongs to high frequency part of level j. Also the fourth submatrix (ℓLL¯ relates to Vj→Vj) represents the impedance in low frequency part.


### 3.1.3. Capacitor

There is a time domain relationship between the voltage and current of a capacitorc as represented by(3.8)v(t)=1c∫i(t)dt.Projecting the above equation onto discrete time domain  leads to(3.9)V̲=1c[DT]-1I̲,
where DT-1  is the discrete form of integral operator in periodic conditions. Thus, in the wavelet domain, DT-1 is expressed as‎(3.10)WDT′̅=(HH′¯HL′¯LH′¯LL′¯).To computeith value of  HH′¯,HL′¯,…  this relation is used(3.11)(hhi′hli′lhi′lli′)=(hhihlilhilli)-1,
where hhi, hli, lhi, and lli are the ith values of HH¯,HL¯,….


## 3.1.1. Resistor

The relationship between voltage and current of a resistorr in the time domain is described as
(3.1)v(t)=ri(t).Then, this relation is expressed in a wavelet expansion as(3.2)(VdjVsj)=r(U00U)(IdjIsj),
where U‎ is‎   an identity ‎matrix ‎with N/2(Jmax-j+1) dimensions,  N is signal length, and J is resolution‎‎level‎.


## 3.1.2. Inductor

The relationship between voltage and current of an inductorℓ  is(3.3)v(t)=ℓdi(t)dt.TheN point discretization of (3.3) leads to(3.4)V̲=ℓ[DT]I̲,
where DT is the discrete form of derivative operator.In thejth level of MRA space, (3.4) can be written as(3.5)(VdjVsj)=ℓ·WDT(IdjIsj),
where(3.6)WDT=(HDTHHDTLLDTHLDTL)=(AsjBsjCsjTsj).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, DT is substituted with LDTL 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 WDT 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)(V̂djV̂sj)=ℓ·WDT¯(ÎdjÎsj),WDT¯=(HH¯HL¯LH¯LL¯),
where ℓ·WDT¯ represents‎  ‎the  wavelet  ‎domain‎  ‎impedance  ‎of  inductor. There are four submatrices for impedance definition of inductor, the first submatrix (ℓHH¯ deals with Wj→Wj) belongs to high frequency part of level j. Also the fourth submatrix (ℓLL¯ relates to Vj→Vj) represents the impedance in low frequency part.


## 3.1.3. Capacitor

There is a time domain relationship between the voltage and current of a capacitorc as represented by(3.8)v(t)=1c∫i(t)dt.Projecting the above equation onto discrete time domain  leads to(3.9)V̲=1c[DT]-1I̲,
where DT-1  is the discrete form of integral operator in periodic conditions. Thus, in the wavelet domain, DT-1 is expressed as‎(3.10)WDT′̅=(HH′¯HL′¯LH′¯LL′¯).To computeith value of  HH′¯,HL′¯,…  this relation is used(3.11)(hhi′hli′lhi′lli′)=(hhihlilhilli)-1,
where hhi, hli, lhi, and lli are the ith values of HH¯,HL¯,….


## 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)-∂v∂x=ri+ℓ∂i∂t,(3.13)-∂i∂x=gv+c∂v∂t,
where r, ℓ, g, and c 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)∂2v∂x2=r1·v+r2·∂v∂t+r3·∂2v∂t2,
where(3.15)r1=r·g,r2=r·c+ℓ·g,r3=ℓ·c.
If (3.14) is transformed to the new domain for ith element of jth level, the following relation can be obtained:(3.16)(∂2V̂dj,i∂x2∂2V̂sj,i∂x2)=(r1+r2hhi+r3hhi2r2hli+r3hli2r2lhi+r3lhi2r1+r2lli+r3lli2)(Vdj,iVsj,i),
where hhi,…,lli and hhi2,…,lli2 are the ith diagonal elements of HH¯,…,LL¯ and HH2¯,…,LL2¯, respectively. Also, HH2=H·DT2·HT,…, LL2=L·DT2·LT and  DT2 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, r1, r2, and r3 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)i(t)=f(v(t)),
where f is a linear operator. As load is in series with the switch, the relation between current and voltage is(3.18)i(t)=p(t)·f(p(t)·v(t)),
where p(t) is switching signal. Switching signal is a periodic function defined as follows:(3.19)p(t)={1:switchison,0:switchisoff.Discretizing equation (3.18) leads to(3.20)I̲=[S]·V̲,(3.21)[S]=[P]·[F]·[P].This relation is obtained based on this fact that mathematical operatorf 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 Figure1). Applying the  KCL relation  yieldsFigure 1
Parallel RLC circuit.(3.22)-ix+vinr+cddtvin+1ℓ∫vindt=0.According to the modified nodal method that is used in harmonic analysis:(3.23)(1r+cp+1ℓp110)(v1ix)=(0vin),
where p is derivative operator,    (3.23) could also be written as follows:(3.24)(CR+pCC+1pCL)(v1ix)=(0vin),
where(3.25)CR=(1r110),CL=(1ℓ000),Cc=(c000).CR, Cc, and CL 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 jth level, the general form of (3.24) is the same as (2.14), where(3.26)d̂fj=(0Γ-1·Vin,dj),ŝfj=(0Γ-1·Vin,sj),d̂xj=(Γ-1·V1,djΓ-1·Iin,dj),ŝxj=(Γ-1·V1,sjΓ-1·Iin,sj).Now according to Section2 these equations are obtained(3.27)asj,i=CR+hhi·CC+hhi′·CL,bsj,i=hli·CC+hli′·CL,csj,i=lhi·CC+lhi′·CL,tsj,i=CR+lli·CC+lli′·CL.The formula (3.27) could be written for any network directly. Using these matrices, the ith value of response vectors can be computed(3.28)rsj,i=tsj,i-csj,iasj,i-1bsj,i,ŝxj,i=rsj,i-1(ŝfj,i-csj,iasj,i-1d̂fj,i),d̂xj,i=asj,i-1(d̂fj,i-bsj,iŝxj,i),
where d̂xj,i and ŝxj,i are the ith values of d̂xj and ŝxj. To obtain jth  level of response vector in MRA space, that is, dxj and dxj,  Γ is multiplied to d̂xj and ŝxj, respectively.The steps for Nonsinusoidal steady-state analysis are as follows.(1)
Determine number of levels (Jmax) and number of samples (N), where Jmax is the index of finest scale in MRA space.(2)
Calculate theCR, CC, and CL matrices according to the modified nodal method.(3)
SetJ=1.(4)
Setj=Jmax-J, where j is the index of current resolution level.(5)
ComputeDT, HDTH,HDTL,… for jth level.(6)
Compute theΓ matrix using (2.13) and then transfer the input vector to the new domain.(7)
Seti=1.(8)
Calculateasj,i,bsj,i, csj,i, and tsj,i using (3.27). Then by using (3.28), calculate d̂xj,i and  ŝxj,i.(9)
Ifi=N/2J then set i=i+1 and go to step (8).(10)
Calculate the response vector in MRA space forjth resolution level, (i.e., dxj  and sxj).(11)
IfJ is not equal to Jmax, then set J=J+1 and go to step (5).(12)
End.Figure2 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 dfj,  then(4.1)[(vα)k]T·dfj=(-1)αn∑k=0N-1dfj,kexp(-j·2·π·α·kn),
where (vα)k is the αth column of Γ introduced in (2.13), and dfj,k is the kth element of dfj. Equation (4.1) shows that  multiplication by a column of Γ returns the DFT of dfj.  In extracting harmonic contents of a typical signal, in the coarsest scale, the first element of ŝfj 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 ŝfj from the coarsest scale (lowest level) represents the content of the main harmonic. With respect to the frequency band, (n+1)th element of ŝfj and d̂fj in each level represents the contents of nth harmonic order.To calculate THD value of a signal from its multiplication byΓ matrix, we define these relations as follows:(4.2)Ad,j=1N∑k=2N/c|2·d̂j,k|2,j=0,…,Jmax;c=2(Jmax-j+2),Ac,0=1N∑k=2N/c|2·ŝj,k|2j=0,THD=(∑j=0JmaxAd,j2)+Ac,02S1,
where N is the number of original signal samples, and S1 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 Figure3 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 IUtil 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.
Bus voltagesNew proposed methodTime domain methodError (%)49-RECT2.482.573.56805:FDR F2.412.524.40009:FDR H2.402.524.78903:MILL2.412.524.39650:GEN-12.322.434.497100:UTIL-690.170.183.111Table 2
THD percentage values.
CurrentNew proposed methodTime domain methodError (%)IUtil2.492.542.106IG20.5421.162.945Figure 3
The 13-bus test system.Figure 4
IUtil 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 Figure5). 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.
TL1 and TL2Rl=0.0955 [Ω/Km]Ll=2.134×10-3 [H/Km]Cl=12.371×10-9 [F/Km]SVCXc=-377j [Ω]XL=377j [Ω]RXL=10 [Ω]LoadsR1=R2=300 ΩXL1=XL2=264j [Ω]Table 4
Comparing of THD percentages and simulation times.
V2V3IinSimulation time (sec)Time domain (FFT)3.38066.51684.94023.44New domain (N=512,  J=3)3.36796.49004.89841.85Figure 5
Test system.Figures6 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 (V3): new proposed domain and time domain.


## 5.1. Case Study 1

The periodic steady-state solution of the test network shown in Figure3 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 IUtil 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.
Bus voltagesNew proposed methodTime domain methodError (%)49-RECT2.482.573.56805:FDR F2.412.524.40009:FDR H2.402.524.78903:MILL2.412.524.39650:GEN-12.322.434.497100:UTIL-690.170.183.111Table 2
THD percentage values.
CurrentNew proposed methodTime domain methodError (%)IUtil2.492.542.106IG20.5421.162.945Figure 3
The 13-bus test system.Figure 4
IUtil 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 Figure5). 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.
TL1 and TL2Rl=0.0955 [Ω/Km]Ll=2.134×10-3 [H/Km]Cl=12.371×10-9 [F/Km]SVCXc=-377j [Ω]XL=377j [Ω]RXL=10 [Ω]LoadsR1=R2=300 ΩXL1=XL2=264j [Ω]Table 4
Comparing of THD percentages and simulation times.
V2V3IinSimulation time (sec)Time domain (FFT)3.38066.51684.94023.44New domain (N=512,  J=3)3.36796.49004.89841.85Figure 5
Test system.Figures6 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 (V3): 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.‎



---

*Source: 584637-2009-09-01.xml*