RMS and EMT simulations in PowerFactory with Python: batch runs and plots

The method

1. Activate a study case with dynamic models, and the events you need (here, the case's own fault at Bus 31 and its clearing).
2. Make a results object in the study case and add what to record: res.AddVariable(obj, "c:firel") for a rotor angle, "m:fehz" for a busbar frequency, "m:u:A" for an EMT phase voltage. Point ComInc's p_resvar to it.
3. Set the simulation: iopt_sim "rms" or "ins" on ComInc, tstop on ComSim, in seconds after app.SetAttributeModeInternal(1).
4. For each run, change what the batch varies (the clearing event's time), run ComInc (it must return 0) and then ComSim. sim.GetSimulationTime() below tstop means the run stopped early.
5. Read the curves: res.Load(), FindColumn(obj, variable), GetValue(row, column) (column -1 is time), then res.Release().
6. Put the case back (settings, event times, p_resvar), reset the calculation, delete the results object, and plot.

The RMS batch script

rms_batch.py

"""RMS simulations in PowerFactory with Python: run one fault with several clearing times,
record the rotor angle of every generator, and plot all the runs on one graph.

Runs on the "39 Bus New England System" example project that comes with PowerFactory 2024.
"""
import csv
import sys

sys.path.append(r"C:\Program Files\DIgSILENT\PowerFactory 2024\Python\3.12")
import powerfactory

PROJECT = "39 Bus New England System"
STUDY_CASE = "2.3 Simulation Fault Bus 31 Stable"   # a bolted fault at Bus 31 at 0 s
CLEARING_TIMES = [0.10, 0.15, 0.20, 0.25]           # s
T_STOP = 5.0                                         # s

app = powerfactory.GetApplicationExt()               # inside PowerFactory: powerfactory.GetApplication()
app.ActivateProject(PROJECT)
case = app.GetProjectFolder("study").GetContents(STUDY_CASE + ".IntCase")[0]
case.Activate()

inc = app.GetFromStudyCase("ComInc")                 # initial conditions
sim = app.GetFromStudyCase("ComSim")                 # run simulation
clearing = inc.GetAttribute("p_event").GetContents("Short-Circuit Clearance.EvtShc")[0]
gens = app.GetCalcRelevantObjects("*.ElmSym")

# A results object of our own that records the rotor angle of every generator.
res = case.CreateObject("ElmRes", "Rotor angles")
for g in gens:
    res.AddVariable(g, "c:firel")                    # rotor angle to the reference machine, deg

app.SetAttributeModeInternal(1)                      # times in seconds, whatever the project shows
saved = (inc.GetAttribute("p_resvar"), inc.GetAttribute("iopt_sim"), sim.GetAttribute("tstop"), clearing.GetAttribute("time"))
inc.SetAttribute("p_resvar", res)
inc.SetAttribute("iopt_sim", "rms")                  # "ins" runs an EMT simulation instead
sim.SetAttribute("tstop", T_STOP)

runs = []
try:
    for tc in CLEARING_TIMES:
        clearing.SetAttribute("time", tc)            # move the case's own clearing event
        if inc.Execute() != 0:
            raise SystemExit("The initial conditions failed")
        sim.Execute()
        reached = sim.GetSimulationTime()            # less than T_STOP when the run stopped early

        res.Load()
        cols = [res.FindColumn(g, "c:firel") for g in gens]
        t, spread = [], []
        for row in range(res.GetNumberOfRows()):
            t.append(res.GetValue(row, -1)[1])       # column -1 is the time
            angles = [res.GetValue(row, c)[1] for c in cols]
            spread.append(max(angles) - min(angles))
        res.Release()
        runs.append({"Clearing time (s)": tc, "Peak spread (deg)": max(spread), "Reached (s)": reached,
                     "Stable": max(spread) < 180 and reached >= T_STOP - 1e-6, "t": t, "spread": spread})
finally:
    inc.SetAttribute("p_resvar", saved[0])           # put the case back as it was
    inc.SetAttribute("iopt_sim", saved[1])
    sim.SetAttribute("tstop", saved[2])
    clearing.SetAttribute("time", saved[3])
    app.SetAttributeModeInternal(0)
    app.ResetCalculation()
    res.Delete()

for r in runs:
    print(f"Cleared after {r['Clearing time (s)']:.2f} s: peak spread {r['Peak spread (deg)']:7.1f} deg, "
          f"{'stable' if r['Stable'] else 'loses synchronism'}")

with open("rms_runs.csv", "w", newline="") as f:
    out = csv.writer(f)
    out.writerow(["Clearing time (s)", "Peak spread (deg)", "Stable"])
    out.writerows([r["Clearing time (s)"], r["Peak spread (deg)"], r["Stable"]] for r in runs)

# One graph for all runs, in the Adjie Brotools colours (blue, gold, teal, rose).
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt

fig, ax = plt.subplots(figsize=(9, 4.8))
for r, colour in zip(runs, ["#3A71C8", "#B45309", "#1F7F63", "#C46A7E"]):
    ax.plot(r["t"], r["spread"], color=colour, lw=2, label=f"Cleared after {r['Clearing time (s)']:.2f} s")
ax.axhline(180, color="#E66A6A", ls="--", lw=1.2)
ax.set(xlabel="Time (s)", ylabel="Rotor angle spread (deg)", xlim=(0, T_STOP), ylim=(0, 400))
ax.grid(color="#DCE3F0")
ax.spines[["top", "right"]].set_visible(False)
ax.legend(frameon=False, loc="upper right")
fig.savefig("rotor_angle_spread.png", dpi=150, bbox_inches="tight")

Result: four clearing times on the 39-bus example

The bolted fault at Bus 31 (the terminal of generator G 02) is cleared after 0.10, 0.15, 0.20 and 0.25 s. The rotor angle spread is the largest minus the smallest rotor angle of the ten generators. Cleared after up to 0.20 s, the swing peaks at 131.8 degrees and dies away; cleared after 0.25 s, the spread runs up to 359 degrees, where the angle wraps round: a pole slip. The example's own study cases put the critical clearing time between 0.21 s (stable) and 0.23 s (unstable).

An EMT run

The same API runs an EMT simulation: set iopt_sim to "ins" and record instantaneous phase quantities. This runs the 39-bus example's own EMT case (a two-phase fault at Bus 03 at 0.2 s, three breakers opening at 0.3 s) for 0.5 s and reads the three phase voltages of Bus 03: 50013 samples, 0.000010 s apart, in 23 s.

emt_run.py

"""EMT simulation in PowerFactory with Python: run the study case's fault, record the three phase
voltages of the faulted busbar, and read the waveforms back.

Runs on the "39 Bus New England System" example project that comes with PowerFactory 2024.
"""
import sys

sys.path.append(r"C:\Program Files\DIgSILENT\PowerFactory 2024\Python\3.12")
import powerfactory

PROJECT = "39 Bus New England System"
STUDY_CASE = "4. EMT Simulation Fault Bus 03"     # a two-phase fault at Bus 03 at 0.2 s
T_STOP = 0.5                                       # s

app = powerfactory.GetApplicationExt()             # inside PowerFactory: powerfactory.GetApplication()
app.ActivateProject(PROJECT)
case = app.GetProjectFolder("study").GetContents(STUDY_CASE + ".IntCase")[0]
case.Activate()

inc = app.GetFromStudyCase("ComInc")
sim = app.GetFromStudyCase("ComSim")
bus = app.GetCalcRelevantObjects("Bus 03.ElmTerm")[0]

res = case.CreateObject("ElmRes", "EMT voltages")
for phase in "ABC":
    res.AddVariable(bus, f"m:u:{phase}")           # instantaneous phase voltage, p.u.

app.SetAttributeModeInternal(1)                    # times in seconds
saved = (inc.GetAttribute("p_resvar"), inc.GetAttribute("iopt_sim"), sim.GetAttribute("tstop"))
inc.SetAttribute("p_resvar", res)
inc.SetAttribute("iopt_sim", "ins")                # "ins": EMT (instantaneous values), "rms": RMS
sim.SetAttribute("tstop", T_STOP)
try:
    if inc.Execute() != 0:
        raise SystemExit("The initial conditions failed")
    sim.Execute()
    res.Load()
    cols = [res.FindColumn(bus, f"m:u:{p}") for p in "ABC"]
    t = [res.GetValue(row, -1)[1] for row in range(res.GetNumberOfRows())]
    wave = {p: [res.GetValue(row, c)[1] for row in range(len(t))] for p, c in zip("ABC", cols)}
    res.Release()
finally:
    inc.SetAttribute("p_resvar", saved[0])
    inc.SetAttribute("iopt_sim", saved[1])
    sim.SetAttribute("tstop", saved[2])
    app.SetAttributeModeInternal(0)
    app.ResetCalculation()
    res.Delete()

print(f"{len(t)} samples, {t[1] - t[0]:.6f} s apart; phase A peak before the fault "
      f"{max(abs(u) for u, tt in zip(wave['A'], t) if tt < 0.2):.3f} p.u.")

The full Python script from the Builder

The Builder's Quick Start makes the same batch from blocks: a list of fault durations, a short circuit added at the busbar, an RMS run, the peak spread and a check, with every run drawn on one graph. On the 39-bus example it gives the same peaks within 2 degrees and the same unstable run. This is its file for the example (the fault at Bus 31):

rms_runs_for_several_fault_durations.py as the Builder writes it: the settings and the workflow⬇ Download the full file (.py, 3,265 lines)

#!/usr/bin/env python3
# RMS runs for several fault durations
"""Made with the PowerFactory Workflow Builder, https://tool.adjiebrotots.com/powerfactory-scripter/builder/

BEFORE YOU RUN
  1. Check the SETTINGS below (lines marked REQUIRED).
  2. Close PowerFactory, or run this file from inside PowerFactory (see RUN_MODE).
  3. Run it with the Python that matches your PowerFactory, e.g.  py -3.12 rms_runs_for_several_fault_durations.py
Results and run.log go to a new folder next to this file (or in OUTPUT_DIR).
"""
# ==== IMPORTS ====
import bisect
import csv
import difflib
import fnmatch
import glob
import importlib.util
import math
import os
import re
import subprocess
import sys
import time
import traceback
import warnings

# ==== SETTINGS (edit here) ====
# The name of this workflow: used for the run folder and in the log.
WORKFLOW_NAME = "RMS runs for several fault durations"
# REQUIRED. The PowerFactory project to open, exactly as the Data Manager shows it (letter case
# counts).
PROJECT = "39 Bus New England System"
# OPTIONAL. The study case to activate first ("" = keep the active one).
STUDY_CASE = "2.3 Simulation Fault Bus 31 Stable"
# OPTIONAL. Where each run writes its folder of results and run.log ("" = the folder of this
# script).
OUTPUT_DIR = ""
# OPTIONAL. True: put every value the workflow changed back when it ends.
UNDO_CHANGES = True
# OPTIONAL. How the results tables are saved in the run folder when the run ends: "csv", "xlsx"
# (Excel, needs pandas and openpyxl) or "both".
TABLE_FORMAT = "csv"
# OPTIONAL. Nominal frequency in Hz for the frequency checks: 50, 60 or "auto" (read from the
# grid).
F_NOMINAL = "auto"
# OPTIONAL. "auto", "engine" (start PowerFactory in the background) or "inside" (run from inside
# PowerFactory).
RUN_MODE = "auto"
# OPTIONAL. True: show the PowerFactory window while running (engine mode only).
SHOW_PF = False
# OPTIONAL. True: no screen updates and fewer output window messages while running.
FAST_MODE = True
# OPTIONAL. Seconds to wait for PowerFactory to start.
CONNECT_TIMEOUT = 120
# OPTIONAL. The folder with the powerfactory Python module ("" = find it).
PF_API_PATH = ""
# OPTIONAL. Log in to PowerFactory as this user ("" = your own).
PF_USER = ""

# How every graph of this run looks: set in the Builder (select nothing, then Graph style).
# Colours are hex codes. A value that cannot be used falls back to the Adjie Brotools style.
GRAPH_STYLE = {
    "template": "brotools",                         # where this style started (a name only)
    "colors": ["#5A91E8", "#B45309", "#4A9E87", "#C46A7E", "#9B59B6"],  # in order; more lines repeat them, paler
    "ramp": "blue",                                 # blue, night, viridis, cividis, magma, coolwarm, greys, classic
    "background": "#FFFFFF",
    "plot_background": "#FFFFFF",
    "text": "#14203A",
    "muted": "#4C5B76",
    "grid": "#C8D4E8",
    "axis": "#C8D4E8",
    "font": "DM Sans",
    "font_size": 9.5,                               # points, for the axis numbers (the rest follows)
    "title_size": 13,                               # points
    "title_bold": True,                             # True or False
    "title_align": "center",                        # center, left
    "grid_axes": "both",                            # both, y, x, none
    "grid_style": "solid",                          # solid, dashed, dotted
    "frame": "left-bottom",                         # left-bottom, box, bottom, none
    "ticks": "out",                                 # out, in, none
    "thousands": True,                              # True: 1,200 rather than 1200
    "line_width": 2,                                # points
    "markers": "auto",                              # auto, always, never
    "line_styles": "solid",                         # solid, vary
    "legend": "bottom",                             # bottom, top, right, inside, none
    "legend_frame": False,                          # True or False
    "width": 10,                                    # inches
    "height": 5.6,                                  # inches
    "dpi": 150,                                     # the PNG is width x dpi pixels wide
    "transparent": False,                           # True: no background (for slides)
    "also_save": "",                                # "", svg, pdf
}

# OPTIONAL. Fault duration (s): Values.
FAULT_DURATION_VALUES = [0.1, 0.15, 0.2, 0.25]

# REQUIRED. Fault: Fault at. The name of the element (ElmTerm), exactly as in the Data Manager.
FAULT_TARGET = "Bus 31"

# OPTIONAL. Rotor angle spread (deg): Highest allowed.
ROTOR_ANGLE_SPREAD_HI = 180

# True: connect, find every element named in this workflow, say what was found, then stop
# without running anything.
CHECK_ONLY = False


# ==== WORKFLOW ====
# Note: For each fault duration this applies a three-phase fault at the busbar, runs an RMS simulation for 5 s and records the rotor angle spread of all synchronous generators (the largest minus the smalle...
def main():
    preflight(["matplotlib"])
    app = pf_connect()
    try:
        if CHECK_ONLY:
            check_elements(app, [FAULT_TARGET + ".ElmTerm", ("*.ElmSym", "all")])
            return
        # -- Step 1: Start (Project: 39 Bus New England System) --
        # -- Step 2: Fault duration (s) (4 values: 0.1, 0.15, 0.2, 0.25) --
        fault_duration_s_items = FAULT_DURATION_VALUES
        for fault_duration_s_pass_no, fault_duration_s in enumerate(fault_duration_s_items, 1):
            ctx_push("Fault duration (s)", fault_duration_s, fault_duration_s_pass_no, len(fault_duration_s_items))
            mark = pf_mark()
            try:
                # -- Step 3: Fault (three-phase at Bus 31 at 1 s, cleared after 0.1 s) --
                add_short_circuit(app, "n4", pf_find(app, FAULT_TARGET + ".ElmTerm", "one"), 1, fault_duration_s, 0, 0, None, "fault", "3ph", None)
                # -- Step 4: Run RMS simulation (To 5 s, recording 1 curve) --
                rec = pf_find(app, "*.ElmSym", "all")
                rms = run_rms(app, 5, None, None, monitors=[], monitor_lists=[(rec, "c:firel")], net="case", case_events="auto", fast=False, extend_to=0)
                rotor_angle_deg = rms.series_combined(rec, "c:firel", "spread", "Spread of * \u00b7 Rotor angle (deg)")
                # Did not initialise: skips to the next pass of Fault duration (s)
                if not ((rms.status == "ok") or (rms.status == "early")):
                    case_skipped("RMS simulation did not initialise")
                    continue
                # -- Step 5: Rotor angle spread (deg) (At most 180, goes on either way) --
                rotor_angle_spread_deg_result = check_limits(rotor_angle_deg, lo=None, hi=ROTOR_ANGLE_SPREAD_HI, after=None, before=None)
                if rotor_angle_spread_deg_result.ok is None:
                    log("Rotor angle spread (deg): NO DATA (" + str(rotor_angle_spread_deg_result.get("reason") or "nothing to check") + "), so it did not pass.", "warn")
                # -- Step 6: Add row to results table (Into "RMS runs", saved at the end automatically) --
                peak_spread = peak(rotor_angle_deg, kind="max", after=None, before=None)  # Peak spread
                table("RMS runs").add({
                    "Peak rotor angle spread (deg)": peak_spread.value,
                    "Stable": rotor_angle_spread_deg_result.text,
                }, status=True)
                # -- Step 7: Rotor angle spread (Every pass on one graph) --
                collect_series("n9", [rotor_angle_deg])
            except StopWorkflow:
                raise
            except Exception as err:
                case_failed(err)
            finally:
                try:
                    pf_undo(mark)
                except Exception as fin:
                    log(f"Could not put things back after this pass: {fin}", "error")
                ctx_pop(["RMS runs"])
        draw_collected("n9", title="Rotor angle spread", ylabel="Rotor angle spread", unit=None, limits=[], events=True, stem=None)
    except StopWorkflow as stop:
        log(f"Stopped: {stop}")
    finally:
        finish(app)

The file goes on with 3,103 lines of helper functions copied from the Builder's runtime (connecting, logging, tables, graphs, putting values back), only those this workflow needs. The download has them all.

How the script works

• Step 2 loops over FAULT_DURATION_VALUES (0.1, 0.15, 0.2 and 0.25 s).
• Step 3, add_short_circuit, puts a three-phase fault at FAULT_TARGET at 1 s and its clearing (an EvtShc with i_shc 4) after the pass's duration; while the workflow's own events run, the study case's events are switched off.
• Step 4, run_rms, makes a temporary results object, records the rotor angle spread of every ElmSym, runs ComInc and ComSim to 5 s, reads the curve with GetColumnValues and puts every setting back.
• Step 5 checks the spread against ROTOR_ANGLE_SPREAD_HI (180 degrees); a run that stops early goes on to the check too. Step 6 adds the row, with the curve's peak from the Peak block.
• Step 7 collects one curve per pass and draws them together as a PNG when the loop ends.

Pitfalls

• Times are in the unit the project shows. SetAttribute("tstop", 5) means 5 ms in a project that shows milliseconds. Call app.SetAttributeModeInternal(1) to write seconds, and set it back to 0.
• Events must be in the event list ComInc reads, inc.GetAttribute("p_event"). An event made anywhere else does nothing.
• A fault is cleared by a second EvtShc with i_shc 4 at the clearing time. Its i_clearShc option does not clear a fault in an RMS or EMT run.
• Run ComInc before every ComSim and check it returned 0: after a failed initialisation there is nothing valid to simulate, so record the run as failed and go on.
• A run can stop early when a machine loses synchronism or the solution fails. Compare sim.GetSimulationTime() with tstop and keep the run as a result.
• Unbalanced faults need an unbalanced run. A balanced RMS simulation applies a single-phase fault as a three-phase one; set iopt_net on ComInc to "rst".
• Reading is slow value by value. GetValue is fine for a few curves; for many, GetColumnValues into an IntVec reads a whole column at once.
• Clean up: ResetCalculation(), then delete the results object you made, and put p_resvar and the event times back.

Questions

How do I run an RMS simulation in PowerFactory with Python?

Get ComInc and ComSim with app.GetFromStudyCase, set iopt_sim to "rms" on ComInc and tstop on ComSim, then call inc.Execute() and sim.Execute().

How do I run an EMT simulation in PowerFactory with Python?

The same two commands with iopt_sim set to "ins" (instantaneous values); record phase variables such as m:u:A.

How do I read simulation results in Python?

Add variables to an ElmRes with AddVariable before the run. After it, call Load(), find a column with FindColumn(obj, variable), read rows with GetValue(row, column) (column -1 is time) and call Release().

How do I find the critical clearing time with a script?

Run the fault with longer and longer clearing times and check the rotor angle spread, or search the edge between stable and unstable by halving the gap. The Builder's "Critical clearing time" Quick Start does the search.

Related

• Automate PowerFactory studies with Python: every guide
• PowerFactory Python API: getting started
• Load flow and results to CSV or Excel
• N-1 contingency analysis
• Short-circuit automation
• Hosting capacity study
• Parameter sweep: PV and QV curves

Attributes used on this page: ElmRes, EvtShc, EvtSwitch, ElmSym, ElmTerm, ElmDsl. All attributes.