# Processing simulation data from power flow analysis

The intro is a bit long. If you're not interested, then please read the update, and have a look at the specific parts I'm highlighting.

I'm very interested in any improvements, but you can assume that all calls starting with psspy. are correct and can't be changed. Also, the calls to psspy that are inside a loop have to be inside a loop, I can't call the functions with several input values.

The code is written in Python 2.7. This is because the simulation tool doesn't work with newer versions.

Update:

In case the code is too long for anyone to bother reviewing it, I'll write up some specific lines that probably could need some reviewing:

• The way I create branch_string. First: nodes_from_string = [str(x) for x in nodes_from], followed by branch_string = [nodes_from_string[x] +'-' +nodes_to_string[x] for x in range(0, len(nodes_from))].
• Is this a good way of doing this. I'm not using nodes_from_string later in the code, but I do nf = nodes_from[branch] and nf_str = str(nf) inside the loop.
• What about the part where the plots are created? Starting with fig = plt.figure().
• Should I do this outside the loop instead? If so, how could this be done in a simple manner?
• Can I simplify the way I'm creating datasheets and writing to Excel?
• Do I use enumerate as it should be used in the last 3 rows? Izero_3 is a np.array.

Explanation:

The code below calls a simulation software (called Power System Simulator for Engineering), and runs power flow analysis using the module psspy (PSS/E to Python).

The code works fine, but unfortunately you can't run it as it stands now, since the software is not a free program, and you don't have the power system model. I've provided sample input variables that can be used instead of the function calls (there are only a few lines that must be changed). This way, the bottom half of the code can be tested properly. I'll also provide an overview of the different outputs the simulation tool gives, and how this works:

When importing psspy, a session is opened in the simulation tool. All calls to functions in that module (psspy.function_name) will run a command in that program, and that session. I have to reload the module in order to open a new session. The previous session (if any) is then closed.

Most of the functions give one or more outputs. The first one is always an error code. 0 means there are no errors. In general, I'm not interested in the error messages (since they are often just information, not actual errors). There are however some error messages that means the case I'm working with is bad. For those cases I've simply done sys.exit("Some message").

To test the function, simply substitute the line in question with a variable assignment. The substitutions that must be made are: Substitute the content of the inner loop with these lines. Note that the functions that are called in the inner loop must be done inside a loop, so I can't remove the loop in my actual code.

Izero_1[x] = (-0.079218775034+0.124398261309j)
Izero_2[x] = (1.59865610333e-16-2.18846208516e-16j)
Izero_dum[x] = (0.0790840312839-0.124262578785j)


The call to psspy.abuschar can be substituted with:

carray = [['Node from    ', 'Node to      ']]


ierr = psspy.scdone() can be removed.

Izero_1, I_zero_2 and Izero_dum contain 21 similar values. The output from print Izero_1 is included in the end. You can however just use the value above for all calls in the loop.

I believe everything after the following paragraph can be tested, if the substitutions above are made.

#==============================================================================
# Start doing the fault analysis.
#==============================================================================


The complete code:

Note that it might appear very long. That's mainly because I've tried to write good and descriptive comments.

# coding=utf-8
# The above comment is to ensure that æøå are allowed in the script

#==============================================================================
# This script is used to simulate faults in a power system. The resulting
# currents are calculated with a variable distance from  the stations on
# both ends.
#
# It start with a fault at a distance L from node 1, to a distance L from
# node 2, with N steps in between.
#==============================================================================

import os,sys
import xlwt
from datetime import datetime

#==============================================================================
# I'm using Spyder/Anaconda so both numpy and matplotlib.plot are
# pre-loaded. The two statements below are there to avoid the constant
# reminder that numpy and matplotlib aren't loaded.
#==============================================================================

import numpy as np
import matplotlib.pyplot as plt

#==============================================================================
# This is the file path to the simulation tool. The script must
# find this folder for this program to work
#==============================================================================

PSSEPATH = r'C:\Program Files (x86)\PTI\PSSE33\PSSBIN'
MODELFOLDER = r'C:\Program Files (x86)\PTI\PSSE33\EXAMPLE'

sys.path.append(PSSEPATH)
os.environ['PATH'] += ';' + PSSEPATH

#==============================================================================
# Import modules that are used in the script. psspy and redirect
# are part of the simulation tool.
# These import statements are not in the top of the script, since
# psspy and redirect can't be found until PSSPATH is part of the
# path
#==============================================================================

import psspy
import redirect
import excelpy
#==============================================================================
# The psspy module is reloaded to ensure it's empty. If it's not reloaded
# then the simulation tool will give back information that's already
# calculated / found last time the module was used.
#==============================================================================

#==============================================================================
# Redirect output from the simulation tool. All output will now be shown
# in the Python console window instead of the simulation tool.
# This way I don't need to have the other program open.
#==============================================================================

redirect.psse2py()
#==============================================================================
# CASENAME is a separate variable, so that I can easily change this
# without changing anything else
#==============================================================================

SCENARIO_FILE_PATH = r'\\server_name\model_folder\scenario_folder'
CASENAME = r'\case_name_1'
CASE = SCENARIO_FILE_PATH + CASENAME + r'.sav'

# Initiate the simulation tool with the maximum number of buses
# This is necessary and can't be skipped.
psspy.psseinit(12000)

# Open the case in the simulation tool:
ierr = psspy.case(CASE)

if ierr == 0:
print "Case loaded: " + CASE

if ierr != 0:
sys.exit("Case could not load. Check name and folder!")

#==============================================================================
# Run full load flow. (psspy.fnsl)

#==============================================================================

ierr = psspy.fnsl()
if ierr != 0:
print 'error code fnsl: ', ierr   # Don't do sys.exit here!

#==============================================================================
# Check if power flow solution was found. If it's not found then the
# model is corrupt and the simulations must stop. No calls to the program
# will work if it's not solved.
#==============================================================================

ierr = psspy.solved()
if ierr != 0:
sys.exit("Power flow did not converge!")
#==============================================================================
# Start doing the fault analysis.
#==============================================================================
#==============================================================================
#==============================================================================
# Decide how many steps the line shall be split into.
# I will not change this often, so I don't want it as a function argument
#==============================================================================

num_steps = 20          # 20 steps
min_dis = 1             # Minimum distance from node 1 (1/20 = 0.05
max_dis = num_steps     # Max distance: (num_steps-1)/num_steps.

dist = [0.001] + [float(x)/num_steps for x in range(min_dis, max_dis)] + [0.999]
#cmpval = [None] * len(dist)

#==============================================================================
# Datetime Will be used when presenting the output to document when the
# simulations were conducted
#==============================================================================

i = datetime.now()
time_str = i.strftime('%Y/%m/%d %H:%M:%S')

#==============================================================================
# Make a list of all lines that are of interest:
# A line goes from node_from[x] to node_to[x]. So 1-2, 2-3, 4-6, 6-8, 12-15

#==============================================================================

nodes_from = [1, 2, 4, 6, 12]
nodes_to =   [2, 3, 6, 8, 15]

nodes_from_string = [str(x) for x in nodes_from]
nodes_to_string = [str(x) for x in nodes_to]

branch_string = [nodes_from_string[x] +'-' +nodes_to_string[x] for x in range(0, len(nodes_from))]
#==============================================================================
# Create Excel-file and give names to sheets.
# Sheets will be called 1-2 for faults on lines 1-2, 2-3 for faults on line 2-3 etc.
#==============================================================================

EXCELPATH = SCENARIO_FILE_PATH

sheet_nr = 1
while os.path.exists(EXCELPATH + CASENAME + '_simulation_%s.xls' % sheet_nr) and sheet_nr < 100:
sheet_nr += 1

wbook_name = (EXCELPATH + CASENAME + '_simulation_%s.xls' % sheet_nr)
workbook = xlwt.Workbook(encoding = "UTF-8")
#==============================================================================
# Run the magical loop that does all the hard work
# Outer loop is used to iterate through all the lines
# Inner loop is used to vary the fault location (20 steps)
#==============================================================================
# This line can't be omitted. It has no output.
psspy.sequence_network_setup(1)

for branch in range(0,len(nodes_from)):

nf = nodes_from[branch]
nt = nodes_to[branch]

# Create a subsystem with the two nodes in question.
# This command has no output
psspy.bsys(1,0,[ 0.69, 400.],0,[],2,[nf,nt],0,[],0,[])

# Get names of each of the two nodes, and remove spaces.
# Output is ierr, and carray. I don't care about ierr.
ierr, carray = psspy.abuschar(1,2,'NAME')
names = [name.strip() for name in carray[0]]

nf_str = str(nf)   # Convert node number to string
nt_str = str(nt)   # Convert node number to string

# Make a name for the end points of the plots (x-axis)
# Format: Node name \n (Node number)
x0_tick = "\n".join([names[0], '(' + nf_str + ')'])
x1_tick = "\n".join([names[1], '(' + nt_str + ')'])

# Initiate list of current values
Izero_1 = len(dist) * [None]
Izero_2 = len(dist) * [None]
Izero_dum = len(dist) * [None]   # Dummy variable
# Dummy is actually a descriptive name

for x in range(0,len(dist)):

# The following commands calls the simulation tool. It has no
# outputs except ierr. I'm not interested in ierr unless I'm
# debugging the model
# I have to get one value at a time. It's impossible to fetch
# several results in one call to the program. The loop is therefore
# necessary.

ierr = psspy.scmu(1,[0,0,0,0,0,0,0],[0.0,0.0,0.0,0.0,0.0],"")
#print 'error code scmu1: ', ierr

ierr = psspy.scmu(2,[9,nf,nt,1,1,0,0],[0.0,0.0,0.0,0.0, dist[x]],r"""1""")
#print 'error code  scmu2: ', ierr

ierr = psspy.scmu(3,[9,nf,nt,1,1,0,0],[0.0,0.0,0.0,0.0, dist[x]],r"""1""")
#print 'error code  scmu3: ', ierr

ierr = psspy.scinit()
#print 'error code  scinit: ', ierr

if ierr != 0:
Izero_1[x] = 0
Izero_2[x] = 0
Izero_dum[x] = 0
else:
ierr, Izero_1[x] = psspy.scbrn2(nf, 999999, '1', 'IZERO')
#print 'error code  scbrn2: ', ierr

ierr, Izero_2[x] = psspy.scbrn2(nt, 999999, '1', 'IZERO')
#print 'error code  scbrn2: ', ierr

ierr, Izero_dum[x] = psspy.scbus2(999999,'FAULTZ')
#print 'error code  scbus2: ', ierr

ierr = psspy.scdone()
#print 'error code  scdone: ', ierr
#==========================================================================
# Convert currents to positive real numbers instead of complex values
#==========================================================================

Izero_1_3 = 3*np.absolute(Izero_1)      # 3I0 in per unit
Izero_2_3 = 3*np.absolute(Izero_2)      # 3I0 in per unit
Izero_dum_3 = 3*np.absolute(Izero_dum)  # 3I0 in per unit
I_base = 1000e6/130e3/np.sqrt(3)

Izero_1_3 = [I_base * x for x in Izero_1_3]     # 3I0 i Ampere
Izero_2_3 = [I_base * x for x in Izero_2_3]     # 3I0 i Ampere
Izero_dum_3 = [I_base * x for x in Izero_dum_3] # 3I0 i Ampere

Izero_3 = np.array([Izero_1_3, Izero_2_3, Izero_dum_3])
#==========================================================================
# Plot figures. Also add figure name, title etc.
#==========================================================================

distance = [100*x for x in dist]

fig = plt.figure()
fig.suptitle('Simulation of fault currents', fontsize=10)
line_1 = plt.plot(avstand, Izero_1_3, label = 'Current - '+names[0])
line_2 = plt.plot(avstand, Izero_2_3, label = 'Current - '+names[1])
line_3 = plt.plot(avstand, Izero_dum_3, label = 'Current - Fault location')
plt.ylabel('3I0 [A]')
plt.ylim(ymin = 0)
plt.legend(loc = 'center left', bbox_to_anchor=(1, 0.895), fontsize=10)
#plt.xticks([0, 20, 40, 60, 80, 100], [names[0], '20%', '40%' ,'60%', '80%', names[1]], fontsize=10)
plt.xticks([0, 20, 40, 60, 80, 100], [x0_tick, '20%', '40%' ,'60%', '80%', x1_tick], fontsize=10)
plot_text_1 = 'Current: Fault at different fault locations between {0} and {1}.'.format(*names)
plot_text_2 = '\nX-axis shows percentage distance from {0}.'.format(names[0])
fig.text(0.2, -0.05, plot_text_1 + plot_text_2, bbox=dict(facecolor='none'),fontsize=10)
fig.text(0.95, -0.0, time_str, fontsize=10)

#==========================================================================
# Write to Excel-file
#==========================================================================

worksheet.write(0, 0, label = 'Simulation of fault currents for line ' +names[0] + ' - ' + names[1])
worksheet.write(0, 5, label = 'Line: ' + branch_string[branch])
worksheet.write(1, 5, label = 'Simulation time:')
worksheet.write(2, 5, label = time_str)
worksheet.write(4, 5, label = 'Folder:')
worksheet.write(5, 5, label = SCENARIO_FILE_PATH)
worksheet.write(6, 5, label = 'File name:')
worksheet.write(7, 5, label = CASENAME[1:] + r'.sav')

headers = ['Distance', names[0], names[1], 'Fault location']

for r, row in enumerate(zip(dist)):
for c, col in enumerate(row):
worksheet.write(r+2, c, label = col)

for r, row in enumerate(zip(*Izero_3)):
for c, col in enumerate(row):
worksheet.write(r+2, c+1, label=col)

workbook.save(wbook_name)


The output from print Izero_1 is:

print Izero_1
[(-0.07938643544912338+0.12485483288764954j), (-0.0792187750339508+0.12439826130867004j), (-0.07901565730571747+0.12390241026878357j), (-0.07881283760070801+0.12340965121984482j), (-0.07861024886369705+0.12292001396417618j), (-0.0784088522195816+0.12243442237377167j), (-0.07820823043584824+0.12195241451263428j), (-0.07800843566656113+0.12147442251443863j), (-0.07780935615301132+0.12099983543157578j), (-0.07761109620332718+0.12052880972623825j), (-0.07741381973028183+0.12006136775016785j), (-0.07721710205078125+0.11959727853536606j), (-0.0770215317606926+0.11913706362247467j), (-0.07682641595602036+0.1186794713139534j), (-0.07663228362798691+0.11822568625211716j), (-0.07643887400627136+0.11777541041374207j), (-0.07624640315771103+0.11732850223779678j), (-0.0760544165968895+0.116884246468544j), (-0.07586311548948288+0.11644333600997925j), (-0.07567299902439117+0.11600597202777863j), (-0.07548230141401291+0.11555929481983185j)]


You're right,

nodes_from_string = [str(x) for x in nodes_from]
nodes_to_string = [str(x) for x in nodes_to]

branch_string = [nodes_from_string[x] +'-' +nodes_to_string[x] for x in range(0, len(nodes_from))]


could maybe better be written as:

branch_string = [ str(nodes_from[x]) + '-' + str(nodes_to[x]) for x in  range(len(nodes_from))


or even better, perhaps:

branch_string = [ "%d-%d" % (nodes_from[x], nodes_to[x]) for x in range(len(nodes_from)) ]


or a little less clear (but shorter):

branch_string = [ "%d-%d" % (x, nodes_to[i]) for i, x in enumerate(nodes_from) ]


Other tips:

• range(0,len(nodes_from)) can be range(len(nodes_from)) ; 0 is the default start if one isn't specified.
• range(len(...)) is a bit of a code smell; if possible, consider reorganizing your data, for instance, how about something like:

lines = [ (1, 2), (2, 3), (4,6), (6,8), (12,15) ]

Then your branch_string setter becomes:

branch_string = [ "%d-%d" % nodes for nodes in lines ]

and your for loop start becomes:

for nf, nt in lines:

• x0_tick = "\n".join([names[0], '(' + nf_str + ')']) can be more simply stated as:

x0_tick = "%s\n(%d)" % (names[0], nf)

which also lets you get rid of nf_str completely

• You did use enumerate correctly near the end, but you should use it more!

for x in range(0, len(headers)): worksheet.write(1, x, label = headers[x])

becomes:

for x, header in enumerate(headers): worksheet.write(1, x, label=header)

which is much cleaner looking.

There's only a single thing I can spot as a non-python programmer:

if ierr == 0:
print "Case loaded: " + CASE

if ierr != 0:
sys.exit("Case could not load. Check name and folder!")


You can use else here.

if ierr == 0:
print "Case loaded: " + CASE
else:
sys.exit("Case could not load. Check name and folder!")

• Going by the general rule that there are many possible ways to fail but usually just one way to succeed, I would just write if ierr != 0: sys.exit(…). Following that, print "Case loaded" would just be the normal flow of the program — no else needed. – 200_success Sep 7 '16 at 6:32