Transient stability analysis by evaluation of the critical clearing time with Python in DIgSILENT PowerFactory

Most power system engineers have faced some kind of stability-related questions in their careers. In my case, I have worked quite a bit on transient stability related matters during my PhD (the publications are available on Google Scholar). That is why I decided to write this article about transient stability analysis by evaluation of the critical clearing time (CCT). I have implemented the analysis in Python and it works with any DIgSILENT PowerFactory model that includes synchronous generators.

In the post, I will provide the Python code to calculate the CCT for a user-specific bus. The script delivers the CCT and the associated critical synchronous generator as a result.

If you want to do a more advanced transient stability analysis, you can download the advanced script. The advanced version evaluates the CCT for all buses and their associated synchronous generators. Additionally, it also provides a ranking of the most critical generators by occurrence throughout the analysis.

Definition of transient stability

“Transient stability is the ability of the power system to maintain synchronism when subjected to a severe transient disturbance.”

According to Kundur’s “Power System Stability and Control” book (in my opinion the bible for each power system engineer)

Definition of critical clearing time

The critical clearing time defines the maximum time a disturbance (here a 3-phase-fault) can be applied without the power system losing its stability. In other words, if you remove a 3-phase-fault before reaching the critical clearing time, the power system will remain stable. If you remove the fault after the critical clearing time, one or several synchronous generators will lose synchronism, thus the power system its stability.

Test grid for demonstration

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.

Transient stability analysis

In this section I will go through the different steps of the transient stability analysis. They include:

1. Preliminaries
2. Import of libraries and definition of functions
3. Critical clearing time assessment
4. Summary of results

Preliminaries

Before starting the analysis you need to define the out-of-step signal for all generators under “Results for Simulation RMS/EMT” as shown below. This signal is necessary for the assessment of the CCT.

Import libraries and definition of 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 defined a function to get the out-of-step signal (s:outofstep) of all generators from the RMS simulation.

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

"""
Created March 2021
@author: Michael Pertl, 
michael@thesmartinsights.com
www.thesmartinsights.com
------------------------------------------------------------------------------- 
"""
import numpy as np
import pandas as pd

import sys
sys.path.append(r'C:Program FilesDIgSILENTPowerFactory 2021 SP1Python3.8')
if __name__ == "__main__":
    import powerfactory as pf
app = pf.GetApplication()
if app is None:
    raise Exception('getting Powerfactory application failed')
#------------------------------------------------------------------------------

def getResults():    
    #get result file
    elmRes = app.GetFromStudyCase('*.ElmRes')    
    app.ResLoadData(elmRes)
  
    #Get number of rows and columns
    NrRow = app.ResGetValueCount(elmRes,0)    
    oSGall = app.GetCalcRelevantObjects('*.ElmSym')    
    SGname = []
    SGfirel = []
    SGoutofstep = []
    #get results for each SG in service
    for oSG in oSGall:        
        if getattr(oSG, 'outserv') == 0:
            ColIndex_firel = app.ResGetIndex(elmRes,oSG,'c:firel')        
            ColIndex_outofstep = app.ResGetIndex(elmRes,oSG,'s:outofstep')        
            result_firel = np.zeros((NrRow))  
            result_outofstep = np.zeros((NrRow))  
            #get results for each time step
            for i in range(NrRow):    
                result_firel[i] = app.ResGetData(elmRes,i,ColIndex_firel)[1]
                result_outofstep[i] = app.ResGetData(elmRes,i,ColIndex_outofstep)[1]
            
            SGname.append(getattr(oSG, 'loc_name'))
            SGfirel.append(result_firel)    
            SGoutofstep.append(result_outofstep)    
    
    #Get index of variable of interest    
    ColIndex_time = app.ResGetIndex(elmRes,elmRes,'b:tnow')    
    #pre-allocate result variables
    result_time = np.zeros((NrRow,))
    #get results for each time step
    for i in range(NrRow):    
        result_time[i] = app.ResGetData(elmRes,i,ColIndex_time)[1]
    
    #put together the needed result variables
    results = pd.DataFrame()
    results['time'] = result_time
    for i, SG in enumerate(SGname):
        results[SG + '_firel'] = SGfirel[i]
        results[SG + '_outofstep'] = SGoutofstep[i] 
    return results

Critical clearing time assessment

At the beginning of the script you can define the following parameters:

• Project name and study case
• Bus at which to evaluate the CCT
• Min/max CCT for evaluation. If CCT is below/above these limits, the assessment stops and provides a message.
• Decimals for CCT evaluation (CCT accuracy) – my recommendation is 3 decimals.

Then you can start the script and run the assessment.

How does the script find the CCT?

First, the script checks whether the CCT is within the min/max CCT limits by applying a 3-phase-fault with both of the limits one after another. If at least one generator is out of step with the lower limit, the CCT is below the minimum. If no generator is out of step with the upper limit, the CCT is above the maximum. In both cases, the assessment stops and provides a message saying whether the CCT is below or above the limit.

If the CCT is within the limits, in the next iteration the fault clearing time is set in the middle between the minimum and maximum CCT.

If at least one of the generators is out of step, the CCT must be smaller than the current fault duration. In this case, in the next iteration, the fault clearing time is set in the middle between the current fault clearing time and the minimum. Further, the current fault clearing time is set as the new maximum.

If no generator is out of step, the CCT must be smaller than the current fault duration. In this case, in the next iteration, the fault clearing time is set in the middle between the current fault clearing time and the maximum. Further, the current fault clearing time is set as the new minimum.

In this way, each iteration gets closer to the CCT. The script prints the results of each iteration in the console (more info in the next section).

#%%main

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

#define bus at which to analyze the CCT
bus = 'Bus 15'

CCT_min, CCT_max = 0.005, 2 #min/max CCT for evaluation (s)
decimals = 3 #decimals for CCT evaluation
#--------------------------------------------

#threshold to stop evaluation based on decimals for CCT evaluation
CCT_thresh = 10**-decimals 

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

 #delete all existing events
ofold_Events=app.GetFromStudyCase('IntEVt')
all_ParEvents=ofold_Events.GetContents('*.EvtParam') #parameter events
for tmpEvent in all_ParEvents: 
    tmpEvent.Delete()  
all_ShcEvents=ofold_Events.GetContents('*.EvtShc') #short-circuit events
for tmpEvent in all_ShcEvents:
    tmpEvent.Delete() 

#get the bus object
oBus = app.GetCalcRelevantObjects(bus + '.ElmTerm')[0]
    
#add fault events
#start fault
ofold_Events.CreateObject('EvtShc','SC_start')
oSCstart=ofold_Events.GetContents('SC_start.EvtShc')[0]
setattr(oSCstart, 'time', 0) #start time of fault
setattr(oSCstart, 'p_target', oBus) #bus of fault
setattr(oSCstart, 'i_shc', 0) #fault type (0 = 3-phase-fault)
setattr(oSCstart, 'ZfaultInp', 0) #fault impedance input type (0 = resistance, reactance)
setattr(oSCstart, 'ZfaultInp', 0) #fault impedance input type (0 = resistance, reactance)
setattr(oSCstart, 'R_f', 0) #fault resistance (Ohm)
setattr(oSCstart, 'X_f', 0) #fault reactance (Ohm)
#clear fault
ofold_Events.CreateObject('EvtShc','SC_clear')
oSCclear=ofold_Events.GetContents('SC_clear.EvtShc')[0]
setattr(oSCclear, 'time', 0) #clearing time of fault
setattr(oSCclear, 'p_target', oBus) #bus of fault
setattr(oSCclear, 'i_shc', 4) #fault type (4 = clear SC)

    
#calculate CCT
i = -1 #iteration variable
cont = 1 #stop when set to 0
log_CCT = [99] #log CCT through iterations
log_SGcrit = [['NA']] #log critical SG through iterations
while cont > 0:
    i = i + 1    
    #evaluate the min/max CCTs (to check if act. CCT is between the limits)
    if i == 0:
        #set fault clearing time
        setattr(oSCclear, 'time', CCT_min) #clearing time of fault    
        # calculate initial conditions
        oInit = app.GetFromStudyCase('ComInc') #get initial condition calculation object
        oInit.Execute()    
        #run RMS-simulation
        oRms = app.GetFromStudyCase('ComSim') #get RMS-simulation object
        setattr(oRms, 'tstop', CCT_max+1)
        oRms.Execute()   
        #retrieve results
        RES = getResults()       
        #check if any SG is out of step
        outofstep_cols = [col for col in RES.columns if 'outofstep' in col]
        outofstep_res = np.array(RES[outofstep_cols])
        outofstep = np.max(outofstep_res)
        
        if outofstep: #if any SG out of step CCT <= CCT_min
            cont = 0
            print('CCT <= %.3f sec' %CCT_min)    
            CCTfin = CCT_min
        else:
            #set fault clearing time
            setattr(oSCclear, 'time', CCT_max) #clearing time of fault    
            # calculate initial conditions
            oInit = app.GetFromStudyCase('ComInc') #get initial condition calculation object
            oInit.Execute()    
            #run RMS-simulation
            oRms = app.GetFromStudyCase('ComSim') #get RMS-simulation object
            oRms.Execute() 
            #retrieve results
            RES = getResults()        
            #check if any SG is out of step
            outofstep_cols = [col for col in RES.columns if 'outofstep' in col]
            outofstep_res = np.array(RES[outofstep_cols])
            outofstep = np.max(outofstep_res)
            
            if not outofstep: #if NO SG out of step CCT >= CCT_max
                cont = 0
                print('CCT >= %.3f sec' %CCT_max)
                CCTfin = CCT_max
    else:
        if i == 1: #second iteration tries CCT in the middle between CCT_min and CCT_max
            CCT_new = (CCT_max+CCT_min) / 2 #middle between min/max CCT
            #set fault clearing time
            setattr(oSCclear, 'time', CCT_new) #clearing time of fault    
            # calculate initial conditions
            oInit = app.GetFromStudyCase('ComInc') #get initial condition calculation object
            oInit.Execute()    
            #run RMS-simulation
            oRms = app.GetFromStudyCase('ComSim') #get RMS-simulation object
            oRms.Execute()  
            #retrieve results
            RES = getResults()        
            #check if any SG is out of step
            outofstep_cols = [col for col in RES.columns if 'outofstep' in col]
            outofstep_res = np.array(RES[outofstep_cols])
            outofstep = np.max(outofstep_res)
            #set temporary min/max CCTs
            CCT_min_tmp = CCT_min
            CCT_max_tmp = CCT_max
            
        #from third iteration CCT_new set in the middle between the previous value and 
        #temporary min/max value depending if any SG was out of step or not.
        #in this way step-wise approaching the real CCT
        else: 
            if outofstep: 
                CCT_max_tmp = CCT_new #if out of step -> tmp. CCT max is previous CCT
                CCT_new = (CCT_new+CCT_min_tmp) / 2 #middle between previous CCT and tmp min CCT        
                CCT_new = np.round(CCT_new, decimals)
            else:
                CCT_min_tmp = CCT_new #if NOT out of step -> tmp. CCT min is previous CCT
                CCT_new = (CCT_new+CCT_max_tmp) / 2 #middle between previous CCT and tmp max CCT 
                CCT_new = np.round(CCT_new, decimals)
            #set fault clearing time
            setattr(oSCclear, 'time', CCT_new) #clearing time of fault    
            # calculate initial conditions
            oInit = app.GetFromStudyCase('ComInc') #get initial condition calculation object
            oInit.Execute()    
            #run RMS-simulation
            oRms = app.GetFromStudyCase('ComSim') #get RMS-simulation object
            oRms.Execute()    
            #retrieve results
            RES = getResults()
            #check if any SG is out of step
            outofstep_cols = [col for col in RES.columns if 'outofstep' in col]
            outofstep_res = np.array(RES[outofstep_cols])
            outofstep = np.max(outofstep_res)
            
        log_CCT.append(CCT_new) #log the calculated CCTs 
        #get critical (out of step) SGs
        tmp_SGcrit = list(RES[outofstep_cols].columns[np.max(RES[outofstep_cols]) == 1])
        tmp_SGcrit = [string.replace('_outofstep', '') for string in tmp_SGcrit]
        log_SGcrit.append(tmp_SGcrit) #log critical (out of step) SGs
        
        #stop criterion: if change of iteration is smaller than threshold
        if np.abs(log_CCT[i] - log_CCT[i-1]) <= CCT_thresh:
            cont = 0
            #if any SG out of step in last iteration reduce CCT by threshold for final CCT
            if outofstep:
                CCTfin = CCT_new - CCT_thresh 
            else:
                CCTfin = CCT_new
        
        #print iteration results: current CCT | out of step (0/1=no/yes) | critical SGs
        print('CCT = %.3f s | out of step = %i | %s' %(CCT_new, outofstep, tmp_SGcrit))

Summary of results

Below the code to put together the results of the iterations and to print the final results. You can check the variable ‘Summary’ for info about all iterations.

#%% summary prints analyzed bus, final CCT and corresponding critical SG(s)
#summary of results during iterations       
Summary = pd.DataFrame()        
Summary['CCT of iteration (s)'] = log_CCT[1:]
Summary['critical SGs of iteration'] = log_SGcrit[1:]

#critical SG(s) which first lose synchronism above final CCT
lastcritSG = [SGs for SGs in log_SGcrit if SGs != []][-1]

print('n')
print('Results------------------------------------------------')
print('n')
print('CCT of %s: %.3f s | critical SG = %s' %(bus, CCTfin, lastcritSG))

The screenshot below shows the iteration and final results of the CCT assessment. Throughout the iterations, the CCT (fault clearing time), out-of-step indicator, and associated generators are printed.

In this example, I analyzed ‘Bus 15’ of the New England system and it shows that the CCT is 0.245 s and that ‘G 09’ is the critical generator. That means ‘G09’ is the first generator to lose synchronism.

Advanced transient stability analysis

In this section, I am explaining the additional capabilities of the advanced transient stability analysis. The advanced assessment provides a more complete picture of the level of transient stability of the investigated grid.

Advantages

• Automatically assess the CCT for all buses of the grid or a selection of buses within a specified voltage range that can be defined as shown in the screenshot below.
• Provides the lowest CCTs of the assessed buses in ascending order
• Provides the most critical synchronous generators

Results of the advanced analysis

As already mentioned, I did the analysis on the New England system. However, you can use the script to analyze any PowerFactory model. The script can be downloaded by clicking on the banner at the end of the post.

Below a screenshot of the results printed to the console using the advanced analysis for the New England system. The results include the bus name, CCT in seconds, and the associated critical SGs. As you can see, sometimes there is only one critical SG and sometimes several.

The plot below shows the 10 lowest CCTs of the New England system. Clearly, the results show that the CCTs for some buses are quite low. I consider CCTs below 120 – 150 ms critical because this time frame comes in the area of the reaction time of the protection systems. Generally, the CCT should always be higher than the reaction time of the protection systems to be on the safe side.

In order to improve the transient stability of a power system, you need to know which generators are critical. Therefore, the analysis provides the most critical generators by their occurrence throughout the analysis. That means, the more often a generator loses synchronism, the more critical it is. Here we see that ‘G 09’ is the most critical one.

In this post, I just focused on how to analyze transient stability. If you are interested in how to improve transient stability, check my follow-up post about Transient Stability Improvement by Re-Dispatching of Synchronous Generators.

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

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

You might also be interested in how to do automatized bottleneck analysis with Python in PowerFactory.

Michael