Open access peer-reviewed chapter

Electric Power System Simulator Tool in MATLAB

Written By

Mohamad Arnaout, Rabih Rammal and Samih Abdulnabi

Submitted: October 4th, 2016 Reviewed: April 3rd, 2017 Published: October 11th, 2017

DOI: 10.5772/intechopen.68955

From the Edited Volume

Science Education

Edited by Antonio Vanderlei Dos Santos and Joao Carlos Krause

Chapter metrics overview

1,941 Chapter Downloads

View Full Metrics


An electric power system is a network of electrical components used to supply, transmit, and use electric power. An example of an electric power system is the network that supplies a region’s homes and industry with power. Due to the complexity and nonlinearity of the power system, hand calculations may be very complicated in some cases, especially when the number of buses or inputs is very large. Here comes the role of software for convergence, time saving, and accuracy. The “Electric Power System Simulator” focuses on three main concepts in power system analysis, the “Power Flow Calculation,” “Faults Calculation,” and “Economic Dispatch Calculation.”


  • GUI
  • Newton‐Raphson
  • unsymmetrical faults
  • economic dispatch

1. Introduction

There are two classes of power system simulation tools analysis: commercial and educational programs. Several commercial programs are available in the market (Power World Simulator, Power System Simulator …). These tools present efficient computational programs for analysis. However, they are inadequate for research and education intentions. This drawback is present because they do not authorize the modification of algorithms or adding new models. For research and education goals, flexibility, and simplicity are more important than computation. Due to several educational and research features, MATLAB becomes one of the efficient and adequate programs in many scientific domains and especially in power systems.

This chapter describes a new MATLAB power system analysis toolbox that uses capabilities of MATLAB in numeric computations to investigate fault calculation, power flow, and economic dispatch of a given power systems. This tool is developed to help students in their education and research studies. In addition, it was included in the curriculum for the graduate electrical engineering students in the Lebanese International University. Since its adoption, students show better understanding of these concepts. In addition, they were able to enhance their basic knowledge and improve their way of thinking.


2. Power system analysis methods

An electric power system is sometimes very complex to analyze using hand calculations especially, if there are nonlinear equations, and a high number of buses. Human can deal with little number of buses, and if the number of buses is high, the hand calculations are very complex. In the Newton‐Raphson method, computer software may solve up to 100,000 or 150,000 buses in very short time, and more accurate when converging to the final solution obeying a specified level of tolerance. For the unbalanced faults, one must calculate the sequence and phase of voltages and currents depending on the type of fault, but computer software will calculate these values within few milliseconds and very accurately. Finally, for the economic dispatch, the value of the incremental cost and the generated powers will change as the value of the demand changes. Thus, the software performs several calculations as the load changes.

The MATLAB tool we are preparing performs several objectives through many power systems methods including: (i) the unsymmetrical faults analysis including line‐to‐ground fault, line‐to‐line fault, and double line‐to‐ground fault, (ii) the Newton‐Raphson method, and (iii) the economic dispatch.

2.1. Unsymmetrical faults analysis

Short circuits occur in three‐phase power systems as follows, in order of frequency of occurrence: single line‐to‐ground, double line‐to‐ground, and balanced three‐phase faults. The path of the fault current may have either zero impedance, which is called bolted short circuit, or nonzero impedance.

When an unbalanced fault occurs in an otherwise balanced system, the sequence networks are interconnected only at the fault location (Figure 1). As such, the computation of fault currents is greatly simplified by the use of sequence networks.

Figure 1.

Schematic representation of unsymmetrical fault.

As in the case of balanced three‐phase faults, unsymmetrical faults have two components of fault current: an AC or symmetrical components including sub‐transient, transient, and steady‐state currents, and a dc component [1, 2].

2.1.1. Unbalanced faults analysis

  • Single line‐ground (SLG) and line‐line (LL) are the principle types of faults in a power system. In addition, other types of faults, such as double line‐ground (DLG), open conductor, and balanced three phases, could be studied.

  • The fault of an unbalanced system is estimated using the concept of symmetrical components [3, 4].

2.1.2. Single line‐to‐ground (SLG) faults

Unbalanced faults will disturb the balancing of the network at the fault location. Therefore, the sequence network should be combined together with respect to the type of fault. A detailed derivation of these relationships will be discussed through this paragraph [3].

The terminal voltage at phase “a” can be transformed into its sequence components as:

V a = V a 0 + V a + + V a E1
I a 0 = V a 3 Z f = V a 0 + V a + + V a 3 Z f E2

The only way that these two constraints can be satisfied is by coupling the sequence networks in series as shown in Figure 2.

Figure 2.

Coupling sequence network for line‐to‐ground fault.

2.1.3. Line‐to‐line (LL) faults

The second most common fault is line‐to‐line, which occurs when two of the conductors come in contact with each other [3].

V a + = V a + I a + Z f E3

To satisfy: I a = I a + , V a + = V a + I a + Z f , I a 0 = 0 , the positive and negative sequence networks must be connected in parallel (Figure 3).

Figure 3.

Coupling sequence network for line‐to‐line fault.

2.1.4. Double line‐to‐ground (DLG) faults

With a double line‐to‐ground (DLG) fault, two line conductors come in contact both with each other and ground [3] as shown in Figure 4.

Figure 4.

Coupling sequence network for a double line‐to‐ground fault.

V a 0 V a + = 3 I a 0 Z f E4

To satisfy: I a = I a 0 + I a + + I a = 0 , and V a + = V a , the three symmetrical circuits during a double line‐to‐ground fault are connected as follows:

2.2. Power flow problem

The estimation of the power flow problem can be expressed using an adequate series of nonlinear equations. These equations represent both Kirchhoff?s Voltage Law and network operation limits. The assessment of the power flow problem is based on four variables for each “i” bus (network node) [4]:

  • V i : voltage magnitude

  • δ i : voltage angle

  • P i : net active power

  • Q i : net reactive power

Depending on which of the above four variables are known (given) and which ones are unknown (to be calculated), two basic types of buses can be defined:

  • PQ bus: P i and Q i are specified; V i and δ i are calculated.

  • PV bus: P i and V i are specified; Q i and δ i are calculated.

PQ buses are normally used to represent load buses without voltage control, and PV buses are used to represent generation buses with voltage control in power flow calculations. A third bus is also needed:

  • Vδ bus: V i and δ i are specified; P i and Q i are calculated.

The V δ bus, also called reference bus or slack bus, has double functions in the basic formulation of the power flow problem:

  • It serves as the voltage angle reference.

  • Since the active power losses are unknown in advance, the active power generation of V δ bus is used to balance generation, load, and losses [5, 6].

The polar form of the power flow equations is given by:

P i = n = 1 N | Y i n V i V n | cos ( θ i n + δ n δ i ) E5
Q i = n = 1 N | Y i n V i V n | sin ( θ i n + δ n δ i ) E6

For each line, numerical values for the series impedance Z and the total line‐charging admittance Y are necessary so that the computer can determine all the elements of the N × N bus admittance matrix of which the typical element Y i j is:

Y i j = | Y i j | θ i j = | Y i j | cos θ i j + j | Y i j | sin θ i j = G i j + j B i j E7

The voltage at any bus of the system is given by:

| V i | = | V i | δ i = | V i | ( cos δ i + j sin δ i ) E8

The net current injected to bus i is given by:

I i = Y i 1 V 1 + Y i 2 V 2 + + Y i N V N = n = 1 N Y i n V n E9

The net scheduled power being injected into the network at bus i is:

P i , s c h e d = P g i P d i E10

where Pgi is the scheduled power being generated at bus i, and Pdi is the scheduled power demand.

The mismatch value of the power is given by:

Δ P i = P i , s c h e d P i , c a l c . E11

Similarly, for the reactive power at bus i:

Δ Q i = Q i , s c h e d Q i , c a l c . E12

Table 1 lists the general number of equations and the state variables in function of the number of buses.

Bus type Number of buses Quantities specified Number of available equations Number of δi, |Vi| state variables
Slack (i = 1) 1 δi, |Vi| 0 0
PV (i = 2,…, Ng + 1) Ng Pi, |Vi| Ng Ng
PQ (i = Ng + 2,…,N) N‐Ng‐1 Pi, Qi 2(N‐Ng‐1) 2(N‐Ng‐1)
Totals N 2N 2N‐Ng‐2 2N‐Ng‐2

Table 1.

The number of equations and state variables of power flow problem.

2.2.1. Newton‐Raphson method applied to power flow study

In all realistic cases, the power flow problem cannot be solved analytically, and hence iterative solutions implemented in computers must be used. Here, we are going to discuss the Newton‐Raphson method.

To apply the Newton‐Raphson method to the solution of the power flow equations, we express bus voltages and line admittances in polar form as follows:

P i = | V i | 2 G i i + n = 1 n i N | V i V n Y i n | cos ( θ i n + δ n δ i ) E13
Q i = | V i | 2 B i i + n = 1 n i N | V i V n Y i n | sin ( θ i n + δ n δ i ) E14

Collecting all the mismatch equations into vector‐matrix form yields:

[ ( P 2 δ 2 P 2 δ n J 11 P n δ 2 P n δ n ) ( | V 2 | P 2 | V 2 | | V n | P 2 | V n | J 12 | V 2 | P n | V 2 | | V n | P n | V n | ) ( Q 2 δ 2 Q 2 δ n J 21 Q n δ 2 Q n δ n ) ( | V 2 | Q 2 | V 2 | | V n | Q 2 | V n | J 22 | V 2 | Q n | V 2 | | V n | Q n | V n | ) ] Jacobian [ Δ δ 2 Δ δ n Δ | V 2 | | V 2 | Δ | V n | | V n | ] Corrections = [ Δ P 2 Δ P n Δ Q 2 Δ Q n ] Mismatches E15

The solution of the above equation is found by an iterative method as follows [46]:

  • Estimate values δ i ( 0 ) and | V i | ( 0 ) for the state variables.

  • Use the estimates to calculate: P i , c a l c . ( 0 ) and Q i , c a l c . ( 0 ) , from (5) and (6).

  • Mismatches Δ P i ( 0 ) and Δ Q i ( 0 ) from (11) and (12).

  • Partial derivatives elements of the Jacobian matrix.

  • Solve the above equation of the Jacobian matrix, the corrections, and the mismatches to find the initial corrections δ i ( 0 ) and Δ | V i | ( 0 ) / | V i | ( 0 ) .

  • Add the solved corrections to the initial estimates to obtain:

    δ i ( 1 ) = δ i ( 0 ) + Δ δ i ( 0 ) E16
    | V i | ( 1 ) = | V i | ( 0 ) + Δ | V i | ( 0 ) = | V i | ( 0 ) ( 1 + Δ | V i | ( 0 ) | V i | ( 0 ) ) E17

Use the new values δ i ( 1 ) and | V i | ( 1 ) as starting values for iteration and then continue.

In more general terms, the updated formulas for starting values of the state variables are:

δ i ( k + 1 ) = δ i ( k ) + Δ δ i ( k ) E18
| V i | ( k + 1 ) = | V i | ( k ) + Δ | V i | ( k ) = | V i | ( k ) ( 1 + Δ | V i | ( k ) | V i | ( k ) ) E19

2.3. Economic dispatch

This section is dedicated to study the economic dispatch concept. For this reason, we consider the system configuration shown is Figure 5. This configuration based on N thermal units serving as a source of generation that would deliver the suitable electric power to the load. Each unit has the cost rate F as an input and its electrical power generated as an output. Therefore, the total system cost is represented by FT, which is the sum of each unit cost rate. The fundamental condition of this system considers that the total output powers should be equal to total power demand.

Figure 5.

Thermal units committed to serve electrical load.

The main objective from the economic dispatch concept is to minimize FTwith respect to the considered constraints. Note that any transmission losses are neglected and any operating limits are not explicitly stated when formulating this problem [7, 8]. That is,

F T = F 1 + F 2 + F 3 + + F N E20
F T = i = 1 N g e n F i ( P i ) E21
= 0 = P l o a d i = 1 N g e n P i E22

This type of optimization system is solved using the Lagrange concept. The extreme value condition of the objective function is determined using the multiplication of the constraint by a constant and adding this factor to the objective function as shown below:

L = F T + λ E23

The paramount conditions needed to determine the highest value of the objective function are based on the derivative of the Lagrange function with respect to the independent variables of each unit. These derivatives should be equal to 0. Consequently, there will be N + l variables (value of Pi for each N units and λ). In addition, the constraint equation is obtained by the derivative of the Lagrange function by Pi with respect to λ [911]:

L P i = d F i ( P i ) d P i λ = 0 E24


0 = d F i d P i λ E25

With respect to the above‐mentioned condition, the minimum operating cost is established when all incremental unit cost are equal to λ. The final step for this procedure is pointed out by the addition of the power demand constraint and the limitation value (minimum and maximum) of each power output unit (inequality constraint) [12].

These constraints are summarized below:

d F i d P i = λ           N g e n   e q u a t i o n s E26
P i , min P i P i , max               2 N g e n   e q u a t i o n s E27
i = 1 N P i = P l o a d                 1   c o n s t r a int E28

When we recognize the inequality constraints, then the necessary conditions may be expanded slightly as shown in the set of equations:

d F i d P i = λ for P i , min P i P i , max   E29
d F i d P i λ for P i = P i , max E30
d F i d P i λ for P i = P i , min E31

3. Flow chart

3.1. Unsymmetrical faults case

The implementation of the unsymmetrical faults analysis in MATLAB is based on the following flow chart (Figure 6):

  • Define the zero, positive, and negative impedance matrices.

  • Define the pre‐fault voltage and the faulted impedance.

  • Select the type of fault.

  • Calculate the admittance and impedance matrices of the power system.

  • Calculate the sequence current and voltage for the selected fault.

  • Calculate the phase voltages and currents for all buses or the faulted buses.

  • Save all results in a text file.

Figure 6.

Flow chart for the unsymmetrical faults.

3.2. Power flow case

Power flow solution is estimated using the Newton‐Raphson method. The fulfillment of this method is achieved using an adequate flowchart (Figure 7):

  • Define the number of buses.

  • Determine the resistance, admittance, and the power specification of each bus.

  • Assign the initial values of the variables.

  • Find the mismatches and Jacobian matrix.

  • Find the unknown variables.

  • Verify the accuracy of the calculation.

  • Calculate the real and imaginary power in each bus and line.

  • Calculate the real and imaginary power losses in each bus and line.

  • Save the power flow solution report in a text file.

Figure 7.

Flow chart for the power flow calculation.

3.3. Economic dispatch case

The execution of economic dispatch procedure depends on several parameters. Figure 8 shows the flow chart of this phenomenon.

  • Define the number of units ( P 1 , P 2 , P 3 P n ) .

  • Precise the lower and upper bounds of each unit.

  • Determine the total demand load and the fuel cost ($/MBtu).

  • Introduce the cost function ($/h) or heat rate function (MBtu/h).

  • Calculate the incremental cost rate λ ($/MWhr).

  • Estimate the cost function ($/h) or heat rate function (MBtu/h).

  • Compute the economic operating point using Lagrange and the efficiency.

Figure 8.

Flow chart for the economic dispatch.


4. MATLAB implementation

4.1. Main display page

The main display page of the MATLAB tool gives the choice for the user to choose between one of the three methods as shown in Figure 9 [13].

Figure 9.

Main display page of the electric power simulator tool.

4.2. Unsymmetrical fault analysis implementation

As the user chooses the first method, which is the “Fault Calculation,” the interface shown in Figure 10 appears [13].

Figure 10.

Unsymmetrical faults interface.

The principle of unsymmetrical faults method throughout the “Power System Simulator” software can be made through two modes:

  • Mode 1: The user enters all the data manually.

  • Mode 2: The user loads the data from a specified file.

Results will show the following: The Ybus and Zbus, the figure of the faults, and the output data (current and voltage) in sequence and phase domain as shown in Figure 11.

Figure 11.

Unsymmetrical faults output interface.

Finally, each result (Ybus, Zbus, Vbus, and Ibus) for both sequence and phase configuration is saved in a separated text file.

4.3. Power flow implementation

As the user chooses the first method, which is the “Power Flow Calculation,” Figure 12 appears [13].

Figure 12.

Power flow solution interface.

The principle of the Newton‐Raphson method throughout the “Power System Simulator” software can be made through two modes:

  • Mode 1: The user enters all the data manually.

  • Mode 2: The user loads the data from a specified file.

Two table results are now filled, the first one is the load flow analysis, and the second is the line flow and losses as shown in Figure 13.

Figure 13.

Power flow solution output interface.

4.4. Economic dispatch implementation

As the user chooses the first method, which is the “Power Flow Calculation,” Figure 14 will appear [13].

Figure 14.

Economic dispatch interface.

The principle of the economic and optimal dispatch method throughout the “Power System Simulator” software can be made through two modes:

  • Mode 1: The user enters all the data manually.

  • Mode 2: The user loads the data from a specified data.

In case of over or under limit estimation, the tool will provide a notification and the number or this unit to the user as shown in Figure 15.

Figure 15.

Economic dispatch output interface.


5. Conclusions

This chapter considers an appropriate guide for electrical engineers students who specialized in power systems analysis and design. The above‐mentioned paragraphs give to them an adequate MATLAB tool in order to facilitate the comprehension of some concept.

These concepts can solve a lot of problems such as power flow solution, unsymmetrical faults analysis, and economic dispatch with or without constraints.

They could also be adopting this tool for their practical studies without any complexity because it is related with their own theoretical knowledge.


  1. 1. Stoft S. Power system economics. IEEE Press. 2002;1(1):11
  2. 2. Rustebakke HM. Electric Utility Systems Practice. 4th ed. New York: Wiley; 1983
  3. 3. Golver D, editor. Power System Analysis and Design. 5th ed. New York: Global Engineering; 2012. p. 853
  4. 4. Saadat H, editor. Power System Analysis. 3rd ed. New York: PSA; 2010. p. 775
  5. 5. Vijay ARB, editor. Power Systems Analysis. 2nd ed. Prentice Hall; 2000. p. 619
  6. 6. Andersson G, editor. Power System Analysis. 1st ed. Zurich: ETH; 2012. p. 185
  7. 7. Wood AJ, Wollenberg BF, editors. Power Generation Operation and Control. New York: Wiley; 2014
  8. 8. Grainger J, Stevenson W, editors. Power System Analysis. McGraw‐Hill Inc.; 1994
  9. 9. Dommel HW, Tinney, WF. Optimal power flow solutions. IEEE Transactions on Power Apparatus and Systems. 1968;PAS‐87:1866–1876
  10. 10. Happ, H. Optimal power dispatch—A comprehensive survey. IEEE Transactions on Power Apparatus and Systems. 1977;PAS‐96:841–854
  11. 11. Horton JS, Grigsby LL. Voltage optimization using combined linear programming and gradient techniques. IEEE Transactions on Power Apparatus and Systems. 1984;PAS‐103(7):1637–1643
  12. 12. Lee KY, Park YM, Ortiz JL. A united approach to optimal real and reactive power dispatch. IEEE Power Engineering Review. 1985; PER‐5(5):1147–1153
  13. 13. Arnaout M. An educational electric power simulator. In: International Conference on Advances in Computational Tools for Engineering Applications; 14 July; Louayzi. Beirut: IEEE; 2016. p. 231–235

Written By

Mohamad Arnaout, Rabih Rammal and Samih Abdulnabi

Submitted: October 4th, 2016 Reviewed: April 3rd, 2017 Published: October 11th, 2017