Research Article  Open Access
A New Approach to Nonsinusoidal SteadyState Power System Analysis
Abstract
A new analysis method using wavelet domain for steadystate operating condition of power system is developed and introduced. Based on waveletGalerkin theory, the system components such as resistor, inductor, capacitor, transmission lines, and switching devices are modeled in discrete wavelet domain for the purpose of steadystate analysis. To solve system equations, they are transferred to wavelet domain by forming algebraic matrixvector relations using the wavelet transform coefficients and the equivalent circuit is thus built for system simulation. After describing the new algorithm, twocase 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 steadystate 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 steadystate 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 functionbased 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 timefrequency 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. WaveletGalerkin 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. Twocase 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]:
where , the function has a companion, the scaling function , and these functions satisfy the following relations:
The coefficients and 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:
such that
By defining as an orthonormal complement of in :
with an MRA, one can use and as the basis functions for Galerkin method.
2.3. WaveletGalerkin Solution of a Periodic Problem
In the MRA space the numerical solution of a differential equation based on WaveletGalerkin method in the th level can be written in this matrix form [13]
The decomposition allows the operator to be split into four pieces ( is called the wavelet space, i.e., the detail or finescale component of ) which can be written as follows:
where
And , are the orthonormal projections of and onto and spaces. The projection is the coarsescale component of the solution , and is the finescale component. To solve (2.7),
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
where
The general forms of the other pieces of are also similar to . For a circulant matrix such as , the eigenvalues are [14]
and the corresponding orthonormal eigenvectors are
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
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:
The volume of calculations is decreased significantly using the above technique. In other words, instead of calculating the inverse of matrices with dimensions in (2.7); (2.16) is used for 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
Then, this relation is expressed in a wavelet expansion as
where is an identity matrix with dimensions, is signal length, and is resolutionlevel.
3.1.2. Inductor
The relationship between voltage and current of an inductor is
The point discretization of (3.3) leads to
where is the discrete form of derivative operator.
In the th level of MRA space, (3.4) can be written as
where
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:
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
Projecting the above equation onto discrete time domain leads to
where is the discrete form of integral operator in periodic conditions. Thus, in the wavelet domain, is expressed as
To compute th value of this relation is used
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 VI characteristic of a differential element of transmission line in continuous time domain is represented by
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
where
If (3.14) is transformed to the new domain for th element of th level, the following relation can be obtained:
where and are the th diagonal elements of and , respectively. Also, , and is the disceretized form of secondorder 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
where is a linear operator. As load is in series with the switch, the relation between current and voltage is
where is switching signal. Switching signal is a periodic function defined as follows:
Discretizing equation (3.18) leads to
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 crosscouplings 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
According to the modified nodal method that is used in harmonic analysis:
where is derivative operator, (3.23) could also be written as follows:
where
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
Now according to Section 2 these equations are obtained
The formula (3.27) could be written for any network directly. Using these matrices, the th value of response vectors can be computed
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 steadystate analysis are as follows.
(1)Determine number of levels () and number of samples (), where is the index of finest scale in MRA space.(2)Calculate the , , and matrices according to the modified nodal method.(3)Set . (4)Set , 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 .(8)Calculate ,, , and using (3.27). Then by using (3.28), calculate and .(9)If then set 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 , then set and go to step (5).(12)End.Figure 2 shows the flowchart of nonsinusoidal steadystate analysis in the proposed domain.
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
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, ()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:
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 twocase studies and the results are compared with the time domain simulation results.
5.1. Case Study 1
The periodic steadystate 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 49RECT 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:Util69 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 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.


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.


Figures 6 and 7 show the simulation results from time domain and the new suggested domain for TCR current and load voltage, respectively.
6. Conclusions
This paper describes a new approach based on MRA space for the Nonsinusoidal SteadyState 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.
References
 A. Semlyen, E. Acha, and J. Arrillaga, “Newtontype algorithms for the harmonic phasor analysis of nonlinear power circuits in periodical steady state with special refrence to magnetic nonlinearities,” IEEE Transactions on Power Systems, vol. 103, pp. 310–317, 1991. View at: Google Scholar
 G. Murere, S. Lefebvre, and X. Dai Do, “A generalized harmonic balance method for EMTP initialization,” IEEE Transactions on Power Delivery, vol. 10, no. 3, pp. 1353–1359, 1995. View at: Publisher Site  Google Scholar
 B. K. Perkins, J. R. Marti, and H. W. Dommel, “Nonlinear elements in the EMTP: steadystate initialization,” IEEE Transactions on Power Systems, vol. 10, no. 2, pp. 593–601, 1995. View at: Publisher Site  Google Scholar
 Q. Wang and J. R. Marti, “A waveform relaxation technique for steady state initialization of circuits with nonlinear elements and ideal diodes,” IEEE Transactions on Power Delivery, vol. 11, no. 3, pp. 1437–1443, 1996. View at: Publisher Site  Google Scholar
 A. Semlyen and A. Medina, “Computation of the periodic steady state in systems with nonlinear components using a hybrid time and frequency domain methodology,” IEEE Transactions on Power Systems, vol. 10, no. 3, pp. 1498–1504, 1995. View at: Publisher Site  Google Scholar
 A. Semlyen and M. Shlash, “Principles of modular harmonic power flow methodology,” IEE Proceedings: Generation, Transmission and Distribution, vol. 147, no. 1, pp. 1–5, 2000. View at: Publisher Site  Google Scholar
 A. P. S. Meliopoulos and C.H. Lee, “An alternative method for transient analysis via wavelets,” IEEE Transactions on Power Delivery, vol. 15, no. 1, pp. 114–121, 2000. View at: Publisher Site  Google Scholar
 T. Zheng, E. B. Makram, and A. A. Girgis, “Power system transient and harmonic studies using wavelet transform,” IEEE Transactions on Power Delivery, vol. 14, no. 4, pp. 1461–1468, 1999. View at: Publisher Site  Google Scholar
 D. C. Robertson, O. I. Camps, and J. S. Meyer, “Wavelets and electromagnetics power system transients,” IEEE Transactions on Power Delivery, vol. 11, no. 2, pp. 1050–1058, 1996. View at: Publisher Site  Google Scholar
 K. Amaratunga, J. R. Williams, S. Qian, and J. Weiss, “WaveletGalerkin solutions for onedimensional partial differential equations,” International Journal for Numerical Methods in Engineering, vol. 37, no. 16, pp. 2703–2716, 1994. View at: Publisher Site  Google Scholar  Zentralblatt MATH  MathSciNet
 I. Daubechies, “The wavelet transform, timefrequency localization and signal analysis,” IEEE Transactions on Information Theory, vol. 36, no. 5, pp. 961–1005, 1990. View at: Google Scholar  Zentralblatt MATH  MathSciNet
 G. Beylkin, “On the representation of operators in bases of compactly supported wavelets,” SIAM Journal on Numerical Analysis, vol. 29, no. 6, pp. 1716–1740, 1992. View at: Publisher Site  Google Scholar  Zentralblatt MATH  MathSciNet
 N. A . Coult, A multiresolution analysis for homogenization of partial differential equations, Ph.D. thesis, Department of Applied Mathematics, University of Colorado, 1997.
 R. M. Gray, Toeplitz and Circulant Matrices: A Review, chapter 3, 2005.
 W. S. Meyer and H. W. Dommel, “Numerical modeling of frequency dependent transmission line parameters in an electromagnetic transient program,” IEEE Transactions on Power Apparatus and Systems, vol. 93, pp. 1401–1409, 1974. View at: Publisher Site  Google Scholar
 A. Ametani, “A highly efficient method for calculating transmission line transients,” IEEE Transactions on Power Apparatus and Systems, vol. 95, no. 5, pp. 1545–1551, 1976. View at: Publisher Site  Google Scholar
 J. R. Marti, “Accurate modeling of frequencydependent transmission lines in electromagnetic transient simulations,” IEEE Transactions on Power Apparatus and Systems, vol. 101, no. 1, pp. 147–157, 1982. View at: Publisher Site  Google Scholar
 R. Yacamini and J. W. Resende, “Thyristor controlled reactors as harmonic sources in HVDC converters station and AC systems,” IEE Proceeding, vol. 133, no. 4, part B, pp. 263–269, 1986. View at: Google Scholar
 W. Xu and H. W. Dommel, “Computation of steady state harmonics of Static Var Compensators,” in Proceedings of the International Conference on Harmonics in Power Systems, pp. 239–245, Nashville, Ind, USA, October 1988. View at: Google Scholar
 L. J. Bohmann and R. H. Lasseter, “Harmonic interactions in thyristor controlled reactor circuits,” IEEE Transactions on Power Delivery, vol. 4, no. 3, pp. 1919–1926, 1989. View at: Publisher Site  Google Scholar
 A. Medina, J. Arrillaga, and E. Acha, “Sparsityoriented hybrid formulation of linear multiports and its applications to harmonic analysis,” IEEE Transactions on Power Delivery, vol. 5, no. 3, pp. 1453–1458, 1989. View at: Publisher Site  Google Scholar
 J. Rico, E. Acha, and T. Miller, “Harmonic domain modeling of three phase ThyristorControlled reactors by means of switching vectors and discrete convolutions,” IEEE Transactions on Power Delivery, vol. 11, no. 3, 1996. View at: Publisher Site  Google Scholar
 T. Tarasiuk, “Hybrid waveletFourier spectrum analysis,” IEEE Transactions on Power Delivery, vol. 19, no. 3, pp. 957–964, 2004. View at: Publisher Site  Google Scholar
 L. Eren and M. J. Devaney, “Calculation of power system harmonics via wavelet packet decomposition in real time metering,” in Proceedings of the IEEE Instrumentation and Measurement Technology Conference, vol. 2, pp. 1643–1647, 2002. View at: Google Scholar
 R. Abuhashim, R. Burch, G. Chang et al., “Test systems for harmonics modeling and simulation,” IEEE Transactions on Power Delivery, vol. 14, no. 2, pp. 579–583, 1999. View at: Publisher Site  Google Scholar
Copyright
Copyright © 2009 A. Ketabi et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.