## 1. Introduction

Power quality is an important branch of power system engineering. It plays an important role to ensure the quality of power being delivered to the industry customer. The emergence of smart grids further distinguishes the importance of power quality. A single power quality event such as voltage sag caused by a fault in transmission or distribution level may cost the affected industries up to millions of monetary losses [1]. Power quality disturbances are categorized into voltage sag, voltage swell, transient, harmonic, voltage notch, and flicker. Power quality research is the study of various phenomena that cause power quality disturbance to occur and the development of mitigation strategy. To develop the right mitigation strategy for the power quality problem, the power quality disturbance phenomenon and its root cause must be fully understood. With the advancement of computer technology and software development, simulation of power quality disturbance is made possible. It allows the researcher to model and simulate a given power system to trace, analyze, and understand the root cause of power quality disturbance [2]. It also allows the researcher to carry out what-if scenarios by changing the model or simulation parameters to test their hypothesis [3]. It provides an insight of how power quality disturbance propagates from the source and through the entire power system network. With that in mind, modeling and simulation approach became one of the widely used research method to model and simulate various power quality disturbances. There are many power system simulation tools [4] available, and each simulation tool has its own merits. However, the widely used tools in the academic research community are PSCAD/EMTDC [5], ATP/EMTP [6], MATLAB with Power System Toolbox [7, 8], Power System Analysis Toolbox [9, 10], Simulink with SimPowerSystems blockset [11–13], and Power Analysis Toolbox [14]. In this book chapter, MATLAB/Simulink with SimPowerSystems is chosen as the simulation platform. A comprehensive set of basic models developed to simulate voltage sag, swell, transient, harmonic, and flicker power quality disturbance are presented in this chapter.

## 2. Modeling approach

The simulation models were developed using MATLAB/Simulink with SimPowerSystems. It is then used to simulate various power quality disturbances and observe how these disturbances distort the power system sinusoidal waveform. The models were developed with minimum number of blocks in mind and use their default settings whenever possible to reserve their simplicity and reproducibility for the reader. The developed models present in the book chapter also serve as basic building blocks to a larger power system. The distribution voltage level used in the models is based on the Malaysia grid code. Simulation models including line fault, induction motor starting, transformer energizing, capacitor bank energizing, lightning impulse, nonlinear load, and electric arc furnace models used to simulate various power quality disturbances are described in Sections 2.1 to 2.8.

### 2.1. Line fault model

The line fault model developed in Simulink is show in Figure 1. The line fault model is used to simulate voltage sag caused by line fault. The line fault model consists of 11 kV, 30 MVA, 50 Hz three-phase source block feeding through 11 kV/0.4 kV, 1 MVA delta/wye transformers to a 10 kW resistive and 100VAR inductive load. There are instantaneous waveform and RMS measurement scopes located at 11 kV and 0.4 kV buses. There are two fault blocks located at the 11 kV bus to simulate line fault and multistage fault. A 0.4 second simulation time is set and ode23tb solver is selected to run the simulation.

This line fault model is capable of simulating various line faults including single line to ground, double line to ground, line-to-line, three phase fault, and multistage. Figure 2 shows voltage sag waveforms caused by a line-to-line fault between phase A and phase B at 11 kV feeder lines at 0.1 to 0.168 seconds. It can be observed that the 11 kV bus experiences two voltage sags at phase A and phase B with different sag magnitude. This is due to high fault resistance of 8 Ω between the two faulted lines. As the voltage sag propagates downstream through the 11 kV/0.4 kV, 1 MVA transformer to the load and the fault type [15] is altered by the delta/wye configuration of the transformer. The unfaulted phase C at 0.4 kV bus experiences a slight voltage swell due to the absent of ground point in the line-to-line fault and high fault resistance.

In power quality studies, voltage sag waveform magnitude is commonly presented in RMS waveform and normalized for better visualization. Figure 3 shows the RMS analysis of line-to-line fault voltage sag waveforms in Figure 2. The sag magnitudes for each phase can be clearly visualized. The slight oscillation occurs at the pre- and post-sag and swell is due to the phase shift during fault [16].

The line fault model is also capable of simulating multistage line faults. The multistage voltage sag is typically due to multiple faults protection relay clearing mechanisms that are not synchronous with each other, thus changing the power system impedance and network configuration leading to multiple stages of voltage sag before fully recovering to its nominal level [17], or the changes of fault or ground impedance during fault. It is sometimes also defined as multiple faults that occur consecutively within a short interval of time and captured within one single record as one event. Figure 4 shows multistage voltage sag instantaneous waveforms caused by double line to ground fault. The fault block is set to simulate fault from 0.1 to 0.168 second with fault impedance of 1 Ω and the multistage fault block is set to simulate fault from 0.168 to 0.3 second with fault impedance of 0.1 Ω. The changes of fault impedance during a fault create a multistage voltage sag power quality disturbance.

Figure 5 shows the RMS waveform of multistage voltage sag to better visualize the multistage voltage sag.

### 2.2. Induction motor starting model

The induction motor starting model is used to simulate voltage sag caused by starting a high-power industry induction motor. The induction motor starting model developed in Simulink is shown in Figure 6. The induction motor starting model can be used to simulate voltage sag caused by induction motor starting, The model consists of a 11 kV, 30 MVA, 50 Hz three-phase source feeder block feeding through 11 kV/0.4 kV, 1 MVA delta/wye transformers, a three-phase breaker as motor starting contactor, a three-phase induction motor, and 10 kW resistive load. There are instantaneous and RMS waveform scopes located at 11 kV and 0.4 kV buses for measurement. Similarly, a 0.4 second simulation time is set and ode23tb solver is selected to run the simulation

Figure 7 shows a three-phase voltage sag instantaneous waveform caused by a 75 kW (100 hp) induction motor starting upon closing of motor starting contactor at 0.1 second. The speed of the induction motor during starting is set at 1 rad/sec using the constant block. The voltage sag at 0.4 kV bus propagates upstream through the transformer to the 11 kV feeder bus. The voltage sag can only be noticed at 0.4 kV bus waveform. The voltage sag magnitude reduces as it propagates upstream toward the 11kV feeder where the voltage sag becomes insignificant and not noticeable.

Three-phase induction motor starting voltage sag is typically balanced and has a shallow drop up to 15% from its nominal magnitude. The sag magnitude of the induction motor voltage sag is dependent on the induction motor power rating. A higher induction motor power rating leads to a lower sag magnitude. The voltage sag pattern can be visualized clearly in RMS waveform as shown in Figure 8.

### 2.3. Transformer energizing model

The transformer energizing model developed in Simulink is shown in Figure 9. It is used to simulate voltage sag caused by transformer inrush current and core saturation during energizing. The model consists of an 11 kV, 30 MVA, 50 Hz three-phase source block feeding through a three-phase breaker block to a 11 kV/0. 4kV, 1 MVA saturable core transformer block to a 10 kW resistive and 100 VAR inductive load. The measurements of instantaneous waveform and RMS are located at the 11 kV feeder bus. The simout block is used to store the simulated data for harmonic FFT analysis using the power gui block. This model enables the simulation of voltage sag caused by transformer energizing. The switchgear is set to open at initiate stage and close at 0.06 second during simulation to simulate voltage sag caused by transformer energizing. A 1 second simulation time is set and ode23tb solver is selected to run the simulation. The voltage sag usually takes more than 1 second to rise back to its nominal voltage level.

Figure 10 shows voltage sag due to the energizing of the transformer upon closing of the switchgear that feeds the transformer. All three phases experience unbalanced voltage sag and gradually rise to its nominal voltage level. The sag magnitude of the transformer energizing voltage sag is dependent on the source feeder power rating and transformer power rating. The higher the transformer power rating, the lower the sag magnitude.

The unique characteristics of voltage sag caused by transformer energizing are as follows: it is unbalanced, has a shallow voltage sag up to 15% from its nominal magnitude, and consists high even harmonic of 2nd order. The power gui block is used to analyze the harmonic content with the setting of up to 2,000 Hz and 10 cycles window, which is the standard measurement window size based on the IEC 61000-4-7 standard. It can analyze up to 40 orders of harmonic content. The analysis starts at 0.2 second up to 0.4 second to allow coverage of voltage sag waveform because the voltage sag begins at 0.06 second. Figure 11 clearly shows phase A consisting of unusual high 2nd, 6th, 12th, 18th, and 24th order of even harmonic, which is the unique characteristic of transformer energizing voltage sag.

### 2.4. Capacitor bank energizing model

The capacitor bank energizing model developed in Simulink is shown in Figure 12. It is used to simulate voltage oscillatory transient caused by capacitor banking energizing for power factor correction in the power system [18, 19]. The model consists of 11 kV, 30 MVA, 50 Hz three-phase source block feeding through 11 kV/0.4 kV, 1 MAV delta/wye transformers to 100 kW resistive and 100 kvar inductive load. There are instantaneous waveform scopes located at 11 kV and 0.4 kV buses for measurement. Each feeder bus consists of a capacitor bank connecting through a three-phase breaker as a switching contactor block. The capacitor bank at 0.4 kV bus has a capacity of 40 kvar, which can compensate power factor up to 0.857 for a 100 kvar inductive load. The capacitor bank at 11 kV bus usually has a higher capacity, in this model 100 kvar is used. The capacitor bank energizing model is capable of simulating oscillatory voltage transient caused by energizing of capacitor bank at 11 kV bus or at 0.4 kV bus. A 0.1 second simulation time is set and ode23tb solver is selected to run the simulation.

Figure 13 shows the energizing of capacitor bank upon closing of three-phase breaker at 0.4 kV feeder line causing voltage transient at 0.4 kV and 11 kV bus. The voltage transient magnitude decreases significantly as it propagates upstream toward the 11 kV bus feeder due to a strong source. The oscillatory voltage transient frequency is determined by the size of the capacitor bank. The larger the capacitor bank size, the lower the voltage transient frequency. The speed of the transient settles down its oscillation depending on the size of the real power load. Larger resistive load size leads to higher damping factor, thus faster settling of transient oscillation.

Figure 14 shows oscillatory voltage transient at 0.4 kV and 11 kV bus caused by the energizing of the capacitor bank at 11 kV feeder line while the capacitor bank at 0.4 kV remains closed. The voltage transient magnitude does not decrease as it propagates downstream toward the 0.4 kV bus feeder due to strong sources at 11 kV bus. It can also be clearly seen that the transient frequency is higher due to the higher capacitive power energized at 11 kV bus.

### 2.5. Lightning impulse model

The lightning impulse model developed in Simulink is shown in Figure 15. It is used to simulate impulsive transient caused by lightning near the transmission line. The model consists of 0.4 kV, 1 MVA, 50 Hz three-phase source block fed to a 10 kW resistive and 10 kvar inductive load. There are instantaneous waveform scopes located at 0.4 kV buses for measurement. The lightning block is connected to the feeder line to induce impulsive transient. Since there is no lightning block available in the MATLAB/Simulink blockset library, it can be built by using existing Simulink blocks.

The lightning block subsystem is shown in Figure 16. A controlled voltage source with resistive and inductive network is used to couple the generated lightning impulse to a given phase of the power system line. For a three-phase system, three sets of controlled voltage source, resistive and inductive networks are required. The lightning impulse that is fed to the voltage controlled source input is computed using Equation 1, where *A* is the impulse magnitude, *α* is the damping factor, *t*_{1} is the time when the impulse starts, *t* is the time function, and *u* is the impulse rise step function.

The lightning impulse model can be implemented using MATLAB Function block, ramp block, constant block, and step function block as shown in Figure 16. The ramp block is used as time function *t*, constant block for impulse magnitude *A* and impulse start time *t*_{1}, and lastly step function block *u* for step rising. Figure 17 shows how Equation 1 is coded in the MATLAB function block.

The standard lightning impulse characteristic is 1.2/50 μs, where the impulse rises to a peak of 1.2 μs and decay to 50% at 50 μs time as defined by IEEE 1159.1 2009 and IEEE C62-41.2 2002 standards [20]. Figure 18 shows the 1.2/50 μs lightning impulse waveform model computed by Equation 1 with damping factor of alpha = 14,000 in the MATLAB function block. The waveform is captured using the scope sampled at 1 MHz, which yields 1 μs per sample. At sample 10,000, which is 0.01 second the impulse rise to the peak magnitude of 1, and at 10,050 sample, which is 50 μs later, the impulse magnitude decays to 0.5036. This validates that the model is closely approximate to the 1.2/50 μs lightning impulse characteristic defined by the standard.

A 0.04 second simulation time is set and ode23tb solver is selected to run the simulation. The impulse magnitude is set at 1 kV, begins at 0.012 second, and the coupling network is set at 10 Ω and 1μH. Figure 19 shows the simulated impulsive transient waveform for all three phase. The coupling network impedance determines how close the lightning discharge is to the transmission line, which in turn determines the impulsive transient magnitude induced in the waveform. The lower the coupling network impedance, the closer the lightning to the transmission line.

### 2.6. Single phase nonlinear load model

The single-phase, nonlinear model developed in Simulink is shown in Figure 20. It is used to simulate triplen harmonic voltage disturbance caused by single-phase bridge rectifier with filter capacitor that is commonly found in domestic and commercial buildings [21]. The model consists of 11 kV, 30 MVA, 50 Hz three-phase source block feeding through a 11 kV/0.4 kV, 1 MVA delta/wye transformer to a 1 MW resistive load, single-phase bridge with 2,000 μF capacitive filter and 10 Ω resistive load for each phase. There are instantaneous waveform scopes located at 11 kV and 0.4 kV buses for measurement. A 0.1 second simulation time is set and ode23tb solver is selected to run the simulation. Figure 21 shows the harmonic waveforms at 11 kV and 0.4 kV bus. Harmonic distortions at 0.4 kV are still noticeable along 60° and 240° of each cycle waveform. However, at 11 kV bus harmonic distortion has be significantly suppressed by the delta/wye transformer.

To visualize the harmonic distortion, the simulation time is set at 0.2 second so that 10 cycles will be simulated. The power gui block is used to analyze the harmonic content with the setting of up to 2,000 Hz and 10 cycles window, which is the standard measurement window size based on the IEC 61000-4-7 standard. Figure 22 clearly shows that at 0.4 kV, phase A consists of odd harmonic order with high 3rd order zero sequence harmonic, which is the unique characteristic of triplen harmonic generated by the single phase nonlinear load. However, as the waveform propagates upstream to 11 kV it can be clearly see that all the triplen harmonics 3rd, 9th, 15th, and 21st generated from the single phase nonlinear load has been suppressed by the delta/wye transformer.

### 2.7. Three-phase nonlinear load model

The three-phase nonlinear load model developed in Simulink is show in Figure 23. This nonlinear load model is used to simulate voltage notches and harmonic caused by a 6-pulse three-phase rectifier [22]. The model consists of 11 kV, 30 MVA, 50Hz three-phase source block feeding through a 11 kV/0.4 kV, 1MAV delta/wye transformer to a 6-pulse controlled three-phase rectifier connected to a 600 V, 10 kW resistive and 1 kVA inductive load. A phase lock loop (PLL) block is used to provide synchronization to the pulse generator block to generate the required pulses to control the three-phase rectifier. A constant block set the firing angle for the pulse-controlled, three-phase rectifier. There is a 400 V, 10 kW resistive load connected in front of the three-phase rectifier. There is a scope to monitor the instantaneous waveform at 11 kV and 0.4 kV bus. This three-phase nonlinear load model is capable to simulate voltage notches and negative sequence harmonic caused by the pulse-controlled, three-phase rectifier.

Figure 24 shows voltage notches caused by pulse controlled three-phase rectifier with pulse firing angle of 30°. The voltage notches for all three phase clearly visible across the sinusoidal waveform at 0.4 kV bus, again the voltage notches are significantly suppressed after it propagates upstream to 11 kV bus through the transformer.

Since voltage notches seriously distort the sinusoidal waveform, therefore it also introduces harmonic distortion. To visualize the harmonic distortion, the power gui block is used to analyze the harmonic content with the setting of up to 2,000 Hz and 10 cycles window start at 1 second. Figure 25 clearly shows that at 0.4 kV phase A consists of high 5th, 7th, 11th, 13th, 17th, 19th, 23rd, 25th, 29th, 31st, 35th, and 37th harmonic order, which is governed by a harmonic composite law in Equation 2, where *H*_{o} is the harmonic order, *n* is the natural number, and *q* is the converter pulse number. As the waveform propagates upstream to 11 kV it can be clearly seen that all these harmonic magnitudes have been suppressed by the strong 11 kV source harmonic order but still remain there.

The voltage notch location and depth is dependent on the firing angle control of the three-phase rectifier and the voltage notch width is dependent on the inductive load. The larger the inductive load, the wider the voltage notch. The voltage notch usually does not propagate upstream because it will be damped by the feeder line and transformer.

### 2.8. Electric arc furnace model

The electric arc furnace model developed in Simulink is shown in Figure 26. It is used to simulate flicker disturbance caused by the electric arc furnace. The model consists of 0.4 kV, 1 MVA, 50 Hz three-phase source block fed directly to an electric arc furnace block. There are instantaneous waveform scopes located at 0.4 kV buses for monitoring. Since there is no electric arc furnace block available in MATLAB/Simulink blockset library, it can be built by using existing Simulink blocks.

The electric arc furnace subsystem block developed in Simulink is shown in Figure 27. A controlled voltage source with resistive and inductive network is used to couple the generated flicker disturbance to a given phase of the power system line. For a three-phase system, three sets of controlled voltage source and resistive and inductive networks are required. The electric arc furnace model uses a hyperbolic model [23, 24] defined in Equation 3, where *V*_{at} is the arc length threshold voltage, *i* is the phase current, C is the arc power, and D is the arc current.

The effect of voltage flicker is determined by the threshold voltage shown in Equation 4, where *V*_{at0} is the base reference voltage when there is no arc activity, *m* is the modulation index, and *ω*_{f} is the flicker frequency.

The electric arc furnace model can be implemented using the MATLAB function block and a sinusoidal block as shown in Figure 27. Each electric arc furnace MATLAB Function block represent each electrode at each phase; therefore, three MATLAB function blocks are required for each phase. The sinusoidal is used to model the flicker frequency and magnitude variation. Lastly, an XY graph scope is located at phase C to monitor the electric arc furnace voltage and current curve. Figure 28 shows how Equations 2 and 3 are coded in the MATLAB function block.

To run the simulation, the sinusoidal block frequency is set at 55.3 rad/sec, which is approximately 8.8 Hz. This is the frequency of interest that can cause the light flickering effect that is sensitive and causes discomfort to the human eye [25]. The arc power C is set at 19 kW, the arc current D is set at 5 kA, the base threshold voltage is set at 200 V, the modulation index is set at 0.2 and lastly the coupling network is set at 0.01 Ω and 1 mH. The simulation time is set at 0.04 second and ode23tb solver is selected. Figure 29 shows the electric arc furnace voltage and current curve through the XY graph scope at phase C. The curve resembles the characteristics of an electric arc furnace operation. Figure 30 shows the instantaneous waveform at 0.4 kV bus and the flicker disturbances are clearly visible for all three phases.

Since flicker waveform also distorts the sinusoidal waveform in some way, therefore it also introduces harmonic distortion. To visualize the harmonic distortion the power gui block is used to analyze the harmonic content with settings of up to 2,000 Hz and 10 cycles window start at 0.1 second. Figure 31 clearly shows that at 0.4 kV phase A consists of odd harmonic of 3rd, 5th, 7th, 9th, 11th, and 13th order with high magnitude of 3rd harmonic order.

## 3. Conclusion

The simulation approach provides the researcher the flexibility to create power system models to simulate power quality disturbance by connecting various functional building blocks in the simulation environment. It gives an insight on how power quality disturbance propagates and behaves within the simulated power system model. The limitation of the simulation approach is its dependency on the capability of the chosen simulation software, basic knowledge of power quality and the simulation software, and the availability of power system building blocks required to build the power system model to simulate the intended power quality disturbance. This book chapter presents simulation models that are able to simulate power quality disturbance, including voltage sag due to fault, induction motor starting, transformer energizing, voltage swell, oscillatory transient, impulsive transient, harmonic, voltage notch, and flicker. These simulation models contribute as basic models to the power quality study as well as the development of power quality education and learning curriculum.

All Simulink models presented in this book chapter has been uploaded to the official Mathworks MATLAB Central File Exchange by the author for reader to download the model. The links for each model are listed below.

Line Fault Model; http://www.mathworks.com/matlabcentral/fileexchange/51928

Induction Motor Starting Model; http://www.mathworks.com/matlabcentral/fileexchange/51929

Transformer Energizing Modelhttp://www.mathworks.com/matlabcentral/fileexchange/51931

Capacitor Bank Energizing Model; http://www.mathworks.com/matlabcentral/fileexchange/51933

Lightning Impulse Model; http://www.mathworks.com/matlabcentral/fileexchange/51934

Single Phase Non Linear Load Model; http://www.mathworks.com/matlabcentral/fileexchange/51935

Three Phase Non Linear Load Model; http://www.mathworks.com/matlabcentral/fileexchange/51936

Electric Arc Furnace Model; http://www.mathworks.com/matlabcentral/fileexchange/51937