Automatized short-circuit analysis with Python in DIgSILENT PowerFactory

Sometimes you want to know the short-circuit currents and short-circuit powers at different locations in the power system. This could be of importance for different reasons. For example, design of protection systems or identification of weak grid areas in terms of short-circuit power. Therefore, in this post, I present an approach for automatized short-circuit analysis with Python in DIgSILENT PowerFactory. The code works on any PowerFactory model and conveniently calculates the short-circuit current/power at all buses of the network.

Grid

To showcase the functionality of the analysis I use the good old 39-bus New England system as shown below. The New England system is available in the example projects in PowerFactory.

Automatized short-circuit analysis in DIgSILENT PowerFactory

The approach for the automatized short-circuit analysis makes use of the ‘ComShc’ object in PowerFactory. It is the short-circuit calculation module that offers many options in terms of the calculation method, fault type, fault impedance, etc.

In the provided script, you can change the following short-circuit parameters:

• calculation method (I have only tested the method according to VDE 0102 Part 0 / DIN EN 60909-0, no guarantee that the other methods work as well)
• fault type: 3-phase-fault, 2-phase-fault, etc.
• fault impedance in terms of resistance and reactance

Import libraries and PowerFactory API

As always, at the beginning of the script definition of the needed libraries and functions. Here we need pandas and matplotlib. Additionally, we load the PowerFactory API.

Important: don’t forget to replace the path (shown in red) with the path to the Python directory of your PowerFactory installation as indicated.

"""
michael@thesmartinsights.com
www.thesmartinsights.com
"""
import pandas as pd
import matplotlib.pyplot as plt

import sys
sys.path.append(r'C:\Program Files\DIgSILENT\PowerFactory 2021 SP1\Python\3.8')
if __name__ == "__main__":
    import powerfactory as pf
app = pf.GetApplication()
if app is None:
    raise Exception('getting Powerfactory application failed')

Core of the automatized short-circuit analyis

Below you can find the core of the analysis. First, you need to define the project name and the study case. Then you have the possibility to set a voltage range for the evaluation buses. Only buses within this voltage range (and in-service) are considered.

The code below contains detailed comments for each step. So, if you want to change something, it should be easy to find.

#%% main

#define project name and study case    
projName =   '_TSI_39_Bus_New_England_System_SC'
study_case = '01_StudyCase.IntCase'

#evaluate buses within line-to-line voltage range
Vbusmin, Vbusmax = 0, 1000 # (kV) 

#activate project
project = app.ActivateProject(projName)
proj = app.GetActiveProject()

#get the study case folder and activate project
oFolder_studycase = app.GetProjectFolder('study')
oCase = oFolder_studycase.GetContents(study_case)[0]
oCase.Activate()

#get all the buses
oBuses = app.GetCalcRelevantObjects('.ElmTerm')

#select buses that are within specified voltage range and in-service
oBuseseval = []
for oBus in oBuses:
    Vbus = getattr(oBus, 'uknom') #bus voltage
    if (Vbusmin <= Vbus <= Vbusmax) & (getattr(oBus, 'outserv') == 0):
        oBuseseval.append(oBus)
numbuseval = len(oBuseseval)       
print('#%i buses to evaluate' %numbuseval)

#variable pre-allocation
Bus_col = []
Skss_col = []
Ikss_col = []
Vbus_col = []

#loop through all buses (that are within specified voltage range)
for oBus in oBuseseval:
    # Short-circuit settings
    oSC = app.GetFromStudyCase('ComShc') #get SC object
    setattr(oSC, 'iopt_mde', 0) # SC calculation method: 0 = VDE 0102 Part 0 / DIN EN 60909-0 (six other methods avail. -> check in SC calculation object in PF)
    setattr(oSC, 'iopt_shc', '3psc') # SC type: 3psc = 3-phase, 2psc = 2-phase, spgf = single-phase-to-ground, 2pgf = 2-phase-to-ground (more options -> check in SC calculation object in PF) 
    setattr(oSC, 'Rf', 0) # (Ohm) fault resistance
    setattr(oSC, 'Xf', 0) # (Ohm) fault reactance
    setattr(oSC, 'iopt_allbus', 0) # fault location option: 0 = User Selection, 1 = Busbars and Junction Nodes, 2 = All Busbars
    setattr(oSC, 'shcobj', oBus) # fault location -> bus/terminal object
    oSC.Execute() #execute load flow
    
    #get variables of interest
    Skss = getattr(oBus, 'm:Skss') # (MVA) short-circuit power
    Ikss = getattr(oBus, 'm:Ikss') # (kA) initial symmetrical short-circuit current
    Vbus = getattr(oBus, 'uknom') # (kV) voltage level of bus
    
    #put together important data
    Bus_col.append(getattr(oBus, 'loc_name')) # bus name
    Skss_col.append(Skss) # (MVA) short-circuit power
    Ikss_col.append(Ikss) # (kA) initial symmetrical short-circuit current
    Vbus_col.append(Vbus) # (kV) voltage level of bus
    
#Results: bus names, SC power, SC current, bus voltage
RES = pd.DataFrame()
RES['name'] = Bus_col
RES['Skss'] = Skss_col
RES['Ikss'] = Ikss_col
RES['Vbus'] = Vbus_col
#sorted results
RES_hightolow_Sk = RES.sort_values('Skss', ascending=False, ignore_index=True) #by SC power - high to low
RES_lowtohigh_Sk = RES.sort_values('Skss', ascending=True, ignore_index=True) #by SC power - low to high
RES_hightolow_Ik = RES.sort_values('Ikss', ascending=False, ignore_index=True) #by SC current - high to low
RES_lowtohigh_Ik = RES.sort_values('Ikss', ascending=True, ignore_index=True) #by SC current - low to high

Plot the results

To visualize the results, I have prepared some bar plots that show the results of the analysis in sorted order (find the results of the New England system below the script).

There are two plots, each containing two subplots. The first one shows the buses sorted from highest to lowest short-circuit power and vice versa. The second one shows the buses sorted from highest to lowest short-circuit current and vice versa.

#%% plot: buses with highest and lowest SC power
#plot specifications
left = 0.01  # the left side of the subplots of the figure
right = 0.99   # the right side of the subplots of the figure
bottom = 0.01  # the bottom of the subplots of the figure
top = 0.99     # the top of the subplots of the figure
wspace = 0.2  # the amount of width reserved for space between subplots,
              # expressed as a fraction of the average axis width
hspace = 0.7  # the amount of height reserved for space between subplots,
              # expressed as a fraction of the average axis height

numhigh_Sk = 10 
if len(RES_hightolow_Sk) < numhigh_Sk:
    numhigh_Sk = len(RES_hightolow_Sk)
numlow_Sk = 10 
if len(RES_lowtohigh_Sk) < numlow_Sk:
    numlow_Sk = len(RES_lowtohigh_Sk)   
    
#plot highest SC powers
plt.figure()
plt.subplots_adjust(left=left, bottom=bottom, right=right, top=top, wspace=wspace, hspace=hspace)

plt.subplot(2,1,1)
plt.bar(range(numhigh_Sk), RES_hightolow_Sk['Skss'][0:numhigh_Sk] , tick_label=RES_hightolow_Sk['name'][0:numhigh_Sk], color='darkred')
plt.xticks(range(numhigh_Sk), rotation='vertical')
plt.ylabel('Skss (MVA)') 
plt.grid(linestyle='--')
plt.title('Buses with highest short-circuit power')

#plot lowest SC powers   
plt.subplot(2,1,2)
plt.bar(range(numlow_Sk), RES_lowtohigh_Sk['Skss'][0:numlow_Sk] , tick_label=RES_lowtohigh_Sk['name'][0:numlow_Sk])
plt.xticks(range(numlow_Sk), rotation='vertical')
plt.ylabel('Skss (MVA)') 
plt.grid(linestyle='--')
plt.title('Buses with lowest short-circuit power')

#%% plot: buses with highest and lowest SC currents
numhigh_Ik = 10 
if len(RES_hightolow_Ik) < numhigh_Ik:
    numhigh_Ik = len(RES_hightolow_Ik)
numlow_Ik = 10 
if len(RES_lowtohigh_Ik) < numlow_Ik:
    numlow_Ik = len(RES_lowtohigh_Sk)   
    
#plot highest SC currents
plt.figure()
plt.subplots_adjust(left=left, bottom=bottom, right=right, top=top, wspace=wspace, hspace=hspace)

plt.subplot(2,1,1)
plt.bar(range(numhigh_Ik), RES_hightolow_Ik['Ikss'][0:numhigh_Ik] , tick_label=RES_hightolow_Ik['name'][0:numhigh_Ik], color='red')
plt.xticks(range(numhigh_Ik), rotation='vertical')
plt.ylabel('Ikss (kA)') 
plt.grid(linestyle='--')
plt.title('Highest SC currents')

#plot lowest SC currents   
plt.subplot(2,1,2)
plt.bar(range(numlow_Ik), RES_lowtohigh_Ik['Ikss'][0:numlow_Ik] , tick_label=RES_lowtohigh_Ik['name'][0:numlow_Ik], color='blue')
plt.xticks(range(numlow_Ik), rotation='vertical')
plt.ylabel('Ikss (kA)') 
plt.grid(linestyle='--')
plt.title('Lowest SC currents')

Results

Below you can find the results of the anaysis of the New England system.

The first plot shows the buses with the highest and lowest short-circuit power. If you want to show more/fewer buses, you can change the parameter ‘numhigh_Sk’ or ‘numlow_Sk’ in the code above.

Clearly, Bus 39 has the highest short-circuit power within the network. This makes sense because a large generator (10000 MVA) is connected to Bus 39. The generator represents the external network – the rest of the USA and Canada.

The second plot shows the buses with the highest and lowest short-circuit current. If you want to show more/fewer buses, you can change the parameter ‘numhigh_Ik’ or ‘numlow_Ik’ in the code above.

Here, we see a different distribution. Why? Because the voltage level of the bus plays a significant role. Bus 39 is connected at a 345-kV-bus and the short-circuit power and current are related by Skss = √3 * Vbus * Ikss, where Vbus…line-to-line bus voltage. On the other hand, the impedance of the grid seen from the bus highly influences the short-circuit current, i.e. the smaller the impedance, the higher the short-circuit current and vice versa. That is why it is important to look at both, the short-circuit power and current.

If you find this post insightful, you might also be interested in the following: Automated bottleneck analysis with Python in PowerFactory

More insights within power systems are to follow in the next article.

Michael