A New Approach to Nonsinusoidal Steady-State Power System Analysis

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.


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

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.

Multiresolution Analysis
In this section the orthonormal basis of compactly supported wavelets is reviewed briefly.The orthonormal basis of compactly supported wavelets of L 2 R is formed by the dilation and translation of single function ψ x 11, 12 : where j, k ∈ Z, the function ψ x has a companion, the scaling function ϕ x , and these functions satisfy the following relations:

2.2
The coefficients 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 L 2 R , that is, the decomposition of Hilbert space into a chain of closed spaces: By defining W j as an orthonormal complement of V j in V j 1 : with an MRA, one can use ϕ n,k x and ψ n,k x as the basis functions for Galerkin method.

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 jth level can be written in this matrix form 13 S j x f.

2.6
The decomposition V j 1 V j ⊕ W j allows the operator S j to be split into four pieces W j is called the wavelet space, i.e., the detail or fine-scale component of V j 1 which can be written as follows: where And d xj , d fj ∈ W j , s xj , s fj ∈ V j are the L 2 -orthonormal projections of x and f onto W j and V j spaces.The projection s xj is the coarse-scale component of the solution x, and d xj is the fine-scale component.To solve 2.7 , , 2.9 At this stage, T Sj is selected and investigated.As the problem described above is periodic and supposing that the differential operator is equal to d m /dx m , the general form of T Sj is where The general forms of the other pieces of S j are also similar to T Sj .For a circulant matrix such as T Sj , the eigenvalues λ α are 14 12 and the corresponding orthonormal eigenvectors v α are

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 where

2.15
In these equations Γ is the modal matrix.The columns of Γ are calculated using 2.13 .A sj , B sj , C sj , and T sj 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: The volume of calculations is decreased significantly using the above technique.In other words, instead of calculating the inverse of matrices with N/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 S j in 2.6 is obtained by the multiplication of the system dimension and the number of considered samples N .

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.

Resistor
The relationship between voltage and current of a resistor r in the time domain is described as v t ri t .

3.1
Then, this relation is expressed in a wavelet expansion as where U is an identity matrix with N/2 J max −j 1 dimensions, N is signal length, and J is resolutionlevel.

Inductor
The relationship between voltage and current of an inductor is The N point discretization of 3.3 leads to where D T is the discrete form of derivative operator.
In the jth level of MRA space, 3.4 can be written as where

WD T HD T H HD T L LD T H LD T L
A sj B sj C sj T sj .

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, D T is substituted with LD T L 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 WD T 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: where • WD T represents the wavelet domain impedance of inductor.There are four submatrices for impedance definition of inductor, the first submatrix HH deals with W j → W j belongs to high frequency part of level j.Also the fourth submatrix LL relates to V j → V j represents the impedance in low frequency part.

Capacitor
There is a time domain relationship between the voltage and current of a capacitor c as represented by Projecting the above equation onto discrete time domain leads to where D −1 T is the discrete form of integral operator in periodic conditions.Thus, in the wavelet domain, D −1 T is expressed as To compute ith value of HH , HL , . . .this relation is used where hh i , hl i , lh i , and ll i are the ith values of HH, HL, . . . .

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 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 where If 3.14 is transformed to the new domain for ith element of jth level, the following relation can be obtained: where hh i , . . ., ll i and hh 2 i , . . ., ll 2 i are the ith diagonal elements of HH, . . ., LL and HH 2 , . . ., LL 2 , respectively.Also, T 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, r 1 , r 2 , and r 3 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.

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 where f is a linear operator.As load is in series with the switch, the relation between current and voltage is where p t is switching signal.Switching signal is a periodic function defined as follows: ⎩ 1 : switch is on, 0 : switch is off.

3.19
Discretizing equation 3.18 leads to This relation is obtained based on this fact that mathematical operator f 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.

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

Mathematical Problems in Engineering
According to the modified nodal method that is used in harmonic analysis: where p is derivative operator, 3.23 could also be written as follows: where C R , C c , 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 jth level, the general form of 3.24 is the same as 2.14 , where

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 ith value of response vectors can be computed The steps for Nonsinusoidal steady-state analysis are as follows.
1 Determine number of levels J max and number of samples N , where J max is the index of finest scale in MRA space.
2 Calculate the C R , C C , and C L matrices according to the modified nodal method.
4 Set j J max − J, where j is the index of current resolution level.
6 Compute the Γ matrix using 2.13 and then transfer the input vector to the new domain.11 If J is not equal to J max , then set J J 1 and go to step 5 .
Figure 2 shows the flowchart of nonsinusoidal steady-state analysis in the proposed domain.

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 d fj , then where v α k is the αth column of Γ introduced in 2. 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: where N is the number of original signal samples, and S 1 is the amplitude value of the main harmonic.

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.

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

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.Figures 6 and 7 show the simulation results from time domain and the new suggested domain for TCR current and load voltage, respectively.

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.

Figure 4 :
Figure 4: I Util waveform comparison of the time and new domain simulations.

Figure 6 :Figure 7 :
Figure 6: TCR current: new proposed domain and time domain.
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, d xj and d xj , Γ is multiplied to d xj and s xj , respectively.

7
Set i 1. 8 Calculate a sj,i , b sj,i , c sj,i , and t sj,i using 3.27 .Then by using 3.28 , calculate d xj,i and s xj,i .
9 If i N/2 J then set i i 1 and go to step 8 .10 Calculate the response vector in MRA space for jth resolution level, i.e., d xj and s xj .
13 , and d fj,k is the kth element of d fj .Equation 4.1 shows that multiplication by a column of Γ returns the DFT of d fj .In extracting harmonic contents of a typical signal, in the coarsest scale, the first element of s fj represents Start Determine level no.J max and sample no.N Calculation of C R , C C , and C L

Table 4 :
Comparing of THD percentages and simulation times.