Automated bottleneck (critical line) analysis with Python in DIgSILENT PowerFactory

On the 8^th of January 2021, the system separation incident in Continental Europe started with the tripping of a line in Croatia due to overcurrent protection. The subsequent shift of the power flow overloaded the neighboring lines and caused them to trip as well. These cascading events are undoubtedly very critical for power systems.

This shows how important it is to analyze the bottlenecks of the power system and under which conditions they appear.

In this post, I will provide Python code to conduct a basic bottleneck analysis with (probably) any DIgSILENT PowerFactory model. The basic analysis takes the actual generator and load setpoints and increases the total load (and so generation) until a pre-defined number of lines is overloaded. The distribution factors for the increase of the load and generator setpoints are determined by their initial setpoints.

If you want to do a more advanced bottleneck analysis, click on the banner below to download the advanced script. The advanced version includes a Monte Carlo simulation where the distribution of load and generator setpoints is varied to get a more complete picture of the possible critical lines aka bottlenecks. You can find more details about the advanced analysis at the end of the post.

If you are new to Python-PowerFactory simulation check this out: Run PowerFactory through Python

Before we continue…

For the demonstration, 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.

Basic bottleneck analysis in PowerFactory with Python

In the following section I will show you the different steps of the analysis. They include:

1. Import of libraries and definition of functions
2. The core of the bottleneck analysis
3. Print the summary of results
4. Restore the initial setpoints of generators and loads

Import libraries and define functions

As always, at the beginning of the script definition of the needed libraries and functions. Here we need NumPy and pandas. Additionally, I have also defined two functions to query the load and generation setpoints.

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 numpy as np
import pandas as pd

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')
#------------------------------------------------------------------------------

#get active/reactive power setpoints of all loads in service
def getLoadSpecs():
    Loads = app.GetCalcRelevantObjects('*.ElmLod')
    load_name, load_P, load_Q = [], [], []
    for load in Loads: #loop through list
        if getattr(load, 'outserv') == 0: #only loads in service
            load_name.append(getattr(load, 'loc_name')) # get name of the generator
            load_P.append(getattr(load, 'plini')) #active power
            load_Q.append(getattr(load, 'qlini'))  #reactive power
    SpecsLoad = pd.DataFrame(zip(load_name,load_P,load_Q), columns=['name','P','Q'])
    return SpecsLoad


#get active/reactive power setpoints, operation mode and min/max capacity of all generators in service
def getGenSpecs():
    Generators = app.GetCalcRelevantObjects('*.ElmSym')
    name, P, Q, avr_mode, Pmin, Pmax = [], [], [], [], [], []
    for gen in Generators: #loop through list    
        if getattr(gen, 'outserv') == 0: #only SGs in service
            name.append(getattr(gen, 'loc_name')) # get name of the generator
            P.append(getattr(gen, 'pgini')) #active power
            Q.append(getattr(gen, 'qgini')) #reactive power
            avr_mode.append(getattr(gen, 'av_mode')) #AVR mode
            Pmin.append(getattr(gen, 'Pmin_uc'))
            Pmax.append(getattr(gen, 'Pmax_uc'))  
    SpecsGen = pd.DataFrame(zip(name,P,Pmin,Pmax,Q,avr_mode), columns=['name','P','Pmin','Pmax','Q','avr_mode'])
    return SpecsGen

The core of the bottleneck analysis

After you have replaced the path to your PowerFactory installation (see above) you only have to specify the following points and you are ready to run the analysis.

• project name and study case
• number of permissible overloads
• threshold (in %) above which the loading counts as overload

The analysis runs until the line loadings exceed the defined number of permissible overloads, e.g. maxOverloads = 1 and OL_threshold = 100 % – the analysis stops when the loading of one line is above 100 %. More information on each step in the analysis can be directly found in the comments of the code (please see below).

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

maxOverloads = 1 # number of permissible overloaded lines
OL_threshold = 100 # (%) threshold for loading to count as overload

#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 initial setpoints of generators and loads by load flow with actual settings
oLoadflow = app.GetFromStudyCase('ComLdf') #get load flow object
setattr(oLoadflow, 'iopt_net', 0) # AC LF balanced
setattr(oLoadflow, 'iopt_apdist', 0) # P as dispatched
setattr(oLoadflow, 'iopt_plim', 1) # consider active power limits
setattr(oLoadflow, 'iopt_lim', 1) # consider reactive power limits
# setattr(oLoadflow, 'iPbalancing', 1) # slack: reference machine
setattr(oLoadflow, 'iopt_pq', 0) # consider voltage dependency of loads
oLoadflow.Execute() #execute load flow

#get the generators and their active power operational limits
Generators = app.GetCalcRelevantObjects('*.ElmSym')
name, P, Q, avr_mode, Pmin, Pmax = [], [], [], [], [], []
for gen in Generators: #loop through list    
    if getattr(gen, 'outserv') == 0: #only SGs in service
        name.append(getattr(gen, 'loc_name')) # get name of the generator
        # P.append(getattr(gen, 'pgini'))
        P.append(getattr(gen, 'm:P:bus1'))
        Q.append(getattr(gen, 'qgini')) #reactive power
        avr_mode.append(getattr(gen, 'av_mode')) #AVR mode
        Pmin.append(getattr(gen, 'Pmin_uc'))
        Pmax.append(getattr(gen, 'Pmax_uc'))  
SpecsGen_init = pd.DataFrame(zip(name,P,Pmin,Pmax,Q,avr_mode), columns=['name','P','Pmin','Pmax','Q','avr_mode'])
# SpecsGen_init = pd.DataFrame(zip(name,P,Pmin,Pmax), columns=['name','P','Pmin','Pmax'])

#get init load specs
SpecsLoad_init = getLoadSpecs()

#get initial line loading
Lines=app.GetCalcRelevantObjects('*.ElmLne')
LineName, LineLoading_init = [], []
for line in Lines: #loop through list
    if getattr(line, 'outserv') == 0: #only line is in service
        LineName.append(getattr(line, 'loc_name')) # get name of the line
        LineLoading_init.append(getattr(line, 'c:loading')) #get loading
LineLoading_init = pd.DataFrame(zip(LineName,LineLoading_init), columns=['name','loading'])
LineLoading_init = LineLoading_init.sort_values('loading',ascending=False)  

#total generation, generation capacity and load
TotGen_init = np.sum(SpecsGen_init['P'])
TotGenCap = np.sum(SpecsGen_init['Pmax'])
TotLoad_init = np.sum(SpecsLoad_init['P'])

#tangent phi = P / Q (assumed constant)
Loadtanphi = SpecsLoad_init['P']/SpecsLoad_init['Q']

#generation shift key / load shift key - based on total initial generaton and load, respectively
GSK = np.array(SpecsGen_init['P']) / TotGen_init
LSK= np.array(SpecsLoad_init['P']) / TotLoad_init

#load variation: starts at initial load and increases by 1 % after each step until permissible overloads are reached
TotLoad_MinToMax = np.arange(TotLoad_init, TotGenCap, step=TotLoad_init/100)

# loops until the number of overloaded lines (numOverloads) >= defined number of overloaded lines (maxOverloads)
i, numOverloads = -1, 0
while numOverloads < maxOverloads:
    i = i + 1
    #set new SG setpoints according to GSK
    Generators = app.GetCalcRelevantObjects('*.ElmSym')
    for ig, gen in enumerate(Generators): #loop through list    
        if getattr(gen, 'outserv') == 0: #only SGs in service
            setattr(gen, 'pgini', TotLoad_MinToMax[i] * GSK[ig]) # generator setpoints acc. to GSK     
    
    #set new load setpoints according to LSK
    Loads = app.GetCalcRelevantObjects('*.ElmLod')
    for il, load in enumerate(Loads): #loop through list
        if getattr(load, 'outserv') == 0: #only loads in service
            setattr(load, 'plini', TotLoad_MinToMax[i] * LSK[il]) # load setpoints acc. to LSK 
            setattr(load, 'qlini', TotLoad_MinToMax[i] * LSK[il] / Loadtanphi[il]) #const. tanphi (cosphi) for Q of loads
    #save scenario
    if app.GetActiveScenario() != None:
        actScen = app.GetActiveScenario()
        actScen.Save()
            
    #load flow settings / now distributed slack by SGs to avoid too high loading of reference machine
    oLoadflow = app.GetFromStudyCase('ComLdf') #get load flow object
    setattr(oLoadflow, 'iopt_net', 0) # AC LF balanced
    setattr(oLoadflow, 'iopt_apdist', 0) # P as dispatched
    setattr(oLoadflow, 'iopt_plim', 1) # consider active power limits
    setattr(oLoadflow, 'iopt_lim', 1) # consider reactive power limits
    setattr(oLoadflow, 'iPbalancing', 4) # distributed slack by SGs
    setattr(oLoadflow, 'iopt_pq', 1) # consider voltage dependency of loads
    oLoadflow.Execute() #execute load flow
    
    #get line loading
    Lines=app.GetCalcRelevantObjects('*.ElmLne')
    LineName, LineLoading = [], []
    for line in Lines: #loop through list
        if getattr(line, 'outserv') == 0: #only line is in service
            LineName.append(getattr(line, 'loc_name')) # get name of the line
            LineLoading.append(getattr(line, 'c:loading')) #get loading
    LineLoading = pd.DataFrame(zip(LineName,LineLoading), columns=['name','loading'])
    
    #get number of overloaded lines
    numOverloads = np.sum(LineLoading['loading'] > OL_threshold)
    
#get SG and load setpoints after calculation  
SpecsGen_after = getGenSpecs()
SpecsLoad_after = getLoadSpecs()  

Print summary of results

The code below prints a summary of the results in the Python console. The summary includes the active power of the generators, active/reactive power of the loads, and the line loadings before and after the analysis. The lines are sorted according to their loading from highest to lowest. Additionally, the total load initially and at the end of the analysis is given as well as the number of overloaded lines.

Results_SG = pd.DataFrame()
Results_SG['name'] = SpecsGen_init['name']
Results_SG['P_SG_init'] = SpecsGen_init['P']
Results_SG['P_SG_after'] = SpecsGen_after['P']
Results_SG['Delta Generation'] = SpecsGen_after['P'] - SpecsGen_init['P']

Results_Loads = pd.DataFrame()
Results_Loads['name'] = SpecsLoad_init['name']
Results_Loads['P_Load_init'] = SpecsLoad_init['P']
Results_Loads['P_Load_after'] = SpecsLoad_after['P']
Results_Loads['Delta Load'] = SpecsLoad_after['P'] - SpecsLoad_init['P']

Results_Lines = pd.DataFrame()
Results_Lines['name'] = LineLoading_init['name']
Results_Lines['Loading_init'] = LineLoading_init['loading']
Results_Lines['Loading_after'] = LineLoading['loading']
Results_Lines['Delta Loading'] = LineLoading['loading'] - LineLoading_init['loading']
Results_Lines = Results_Lines.sort_values('Loading_after',ascending=False)  

print('\n')
print('Summary------------------------------------------------')
print('\n')
print('Synchronous Generators----------------------')
print(Results_SG.to_string(index=False))
print('\n')
print('Loads---------------------------------------')
print(Results_Loads.to_string(index=False))
print('\n')
print('Lines---------------------------------------')
print(Results_Lines.to_string(index=False))
print('\n')

TotLoad_after = np.sum(Results_Loads['P_Load_after'])
LoadIncrease = (TotLoad_after-TotLoad_init) / TotLoad_init * 100
print('The total load increased from %.2f MW to %.2f MW by %i percent' %(TotLoad_init,TotLoad_after,LoadIncrease))
print('Number of overloaded lines: %i' %(numOverloads))

Restore initial setpoints of generators and loads

If you want to restore the initial setpoints of the generators and loads, you can use the following code. If you don’t use it, the final setpoints at the end of the analysis will remain in the PowerFactory model.

for ig, gen in enumerate(Generators): #loop through list    
    if getattr(gen, 'outserv') == 0: #only SGs in service
        setattr(gen, 'pgini', SpecsGen_init['P'][ig])         

for il, load in enumerate(Loads): #loop through list
    if getattr(load, 'outserv') == 0: #only loads in service
        setattr(load, 'plini', SpecsLoad_init['P'][il])
        setattr(load, 'qlini', SpecsLoad_init['Q'][il])
#save scenario
if app.GetActiveScenario() != None:
    actScen = app.GetActiveScenario()
    actScen.Save()

Advanced bottleneck analysis using Monte Carlo simulations

In the following, I will provide a brief description of the capabilities of the advanced bottleneck analysis. In order to get a more complete picture of the possible bottlenecks of a power system, different load flow situations need to be analyzed. Here, I am doing that by generating load samples within a defined minimum and maximum total load. For each load sample, I generate different generator setpoints that meet the total load to achieve various load flow scenarios. That basically means different load and generation combinations.

User-specific inputs for the advanced analysis

Basically, you can get started with the advanced analysis by specifying the PowerFactory project name, study case, and the total load range. The rest is already pre-defined as shown in the screenshot below.

The advanced bottleneck analysis has the following user-specific inputs:

• load samples: number of different load samples to analyze
• generation samples: number of generation samples analyzed for each load sample
• min. load: the minimum total load of the system
• max. load: the maximum total load of the system
• standard deviation of individual loads: defines how much the individual loads are varied around their mean (the initial load setpoints serve as the mean for each load individually)
• overload threshold: defines the loading level (in %) above which the loading counts as overload

Results of the advanced analysis

Below you will find an example of the results from the advanced analysis of the New England system. The most important output of the analysis is the ranking of the overloaded lines by the number of occurrences in the load flow simulations. By default, the 10 most overloaded lines are shown.

However, the advanced analysis also allows to dig deeper into the load flow situations. For example, if you want to know how were the setpoints of a specific generator distributed when a specific line was overloaded. Below an example of the generator G 01 setpoints when Line 02 – 03 was overloaded.

I hope you find this article insightful. Consider leaving a comment to let me know your thoughts on the topic.

More insights within energy are to follow in the next article.

Btw, have you also checked how to do peak shaving with energy storage? If not, find it here: Peak shaving with energy storage

Michael