N-1 contingency analysis in PowerFactory with Python

The method

1. Run the base case load flow, and record it: an outage is judged against it.
2. List the outages: every line and network transformer in service (outserv == 0). A generator's own transformer is a loss of generation; study those separately (see Pitfalls).
3. For each outage, set outserv to 1 and run the load flow. A load flow that does not converge is a result too: record it and go on.
4. Record the worst values: the highest c:loading of the branches still in service and the lowest m:u of the energised busbars, with where they are.
5. Put the element back in a finally block, so an error never leaves the project with an outage in it.
6. Rank the outages by loading and save the table.

The script

contingency.py

"""N-1 contingency analysis in PowerFactory with Python: take each line and network transformer
out in turn, run the load flow, and keep the highest loading and the lowest voltage of each outage.

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 = "1. Power Flow"
LOADING_LIMIT = 100.0                               # %

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

ldf = app.GetFromStudyCase("ComLdf")
branches = app.GetCalcRelevantObjects("*.ElmLne") + app.GetCalcRelevantObjects("*.ElmTr2")
buses = app.GetCalcRelevantObjects("*.ElmTerm")
# A generator's own transformer is left out: losing it loses the generator, a generation
# contingency that needs governor response, not a single slack taking up all of it.
gen_terminals = {g.GetAttribute("bus1").GetParent().GetFullName() for g in app.GetCalcRelevantObjects("*.ElmSym")}
outages = [o for o in branches if o.outserv == 0 and not (
    o.GetClassName() == "ElmTr2" and o.GetAttribute("buslv").GetParent().GetFullName() in gen_terminals)]


def worst():
    """The highest branch loading and the lowest busbar voltage of the last load flow."""
    loading, branch = max((o.GetAttribute("c:loading"), o.loc_name) for o in branches if o.outserv == 0)
    # A busbar cut off by the outage reads exactly 0 p.u.: it is islanded, not low.
    u, bus = min((b.GetAttribute("m:u"), b.loc_name) for b in buses if b.GetAttribute("m:u") > 0)
    return {"Highest loading (%)": loading, "Most loaded": branch, "Lowest voltage (p.u.)": u, "Lowest at": bus}


rows = []
if ldf.Execute() != 0:
    raise SystemExit("The base case load flow does not converge")
rows.append({"Outage": "Base case", "Converged": "yes", **worst()})

for out in outages:
    out.SetAttribute("outserv", 1)                  # take it out of service
    try:
        if ldf.Execute() == 0:
            rows.append({"Outage": out.loc_name, "Converged": "yes", **worst()})
        else:
            rows.append({"Outage": out.loc_name, "Converged": "no"})
    finally:
        out.SetAttribute("outserv", 0)              # always put it back

overloads = [r for r in rows if r.get("Highest loading (%)", 0) > LOADING_LIMIT]
print(f"{len(rows) - 1} outages, {len(overloads)} with a branch above {LOADING_LIMIT:.0f} %")
for r in sorted(overloads, key=lambda r: -r["Highest loading (%)"])[:5]:
    print(f"  {r['Outage']:16} {r['Most loaded']:16} {r['Highest loading (%)']:6.1f} %")

with open("n-1_results.csv", "w", newline="") as f:
    out_csv = csv.DictWriter(f, fieldnames=["Outage", "Converged", "Highest loading (%)", "Most loaded",
                                            "Lowest voltage (p.u.)", "Lowest at"])
    out_csv.writeheader()
    out_csv.writerows(rows)

Result on the 39-bus example

The base case already loads Line 21 - 22 to 99.1 %. Of the 37 outages (34 lines and 3 network transformers), 18 push a branch above 100 %. The two worst are each other's parallel path: with Line 21 - 22 out, Line 23 - 24 carries 158.3 %, and with Line 23 - 24 out, Line 21 - 22 carries 156.9 %. Line 16 - 19 is the only link of generators G 04 and G 05 to the rest of the grid: without it their export is lost, the slack machine makes it up, and its transformer Trf 06 - 31 reaches 156.5 %.

The base case, then every outage that overloads a branch, worst first. Every outage converged.

The full Python script from the Builder

The Builder's N-1 Quick Start runs the same loop with its Contingency block, which also takes out generator transformers, and keeps line and transformer loading in separate columns. On this example it gives the same numbers: 158.3 % with Line 21 - 22 out, 108.0 % with Line 10 - 11 out. This is its file for the 39-bus example:

n-1_screening.py as the Builder writes it: the settings and the workflow⬇ Download the full file (.py, 1,695 lines)

#!/usr/bin/env python3
# N-1 screening
"""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 n_1_screening.py
Results and run.log go to a new folder next to this file (or in OUTPUT_DIR).
"""
# ==== IMPORTS ====
import csv
import difflib
import fnmatch
import glob
import importlib.util
import itertools
import math
import os
import re
import string
import subprocess
import sys
import time
import traceback

# ==== SETTINGS (edit here) ====
# The name of this workflow: used for the run folder and in the log.
WORKFLOW_NAME = "N-1 screening"
# 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 = "1. Power Flow"
# 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 = ""

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


# ==== WORKFLOW ====
# Note: After a base load flow this takes out each line and transformer in service in turn (the base case first), runs the load flow and reads the loading of every line and 2-winding transformer and the vo...
def main():
    preflight([])
    app = pf_connect()
    try:
        if CHECK_ONLY:
            check_elements(app, [("*.ElmLne", "all"), ("*.ElmTerm", "all"), ("*.ElmTr2", "all")])
            return
        # -- Step 1: Start (Project: 39 Bus New England System) --
        # -- Step 2: Base load flow (Settings of the study case) --
        base_load_flow = run_load_flow(app, "case", "case", "case")
        # Did not converge: ends the workflow
        if not base_load_flow.ok:
            log("Load flow did not converge. The workflow ends here.", "warn")
            return
        # -- Step 3: Contingency (Base case, then each N-1 outage of lines, transformers) --
        contingency_contingencies = define_contingencies(app, ["*.ElmLne", "*.ElmTr2", "*.ElmTr3"], 1, True)
        for contingency_pass_no, contingency in enumerate(each_contingency(app, contingency_contingencies, True), 1):
            ctx_push("Contingency", contingency, contingency_pass_no, len(contingency_contingencies) + (1 if True else 0))
            try:
                # -- Step 4: Run load flow (Settings of the study case) --
                ldf = run_load_flow(app, "case", "case", "case")
                # Did not converge: skips to the next pass of Contingency
                if not ldf.ok:
                    case_skipped("Load flow did not converge")
                    continue
                # -- Step 5: Line loading (c:loading of all * (ElmLne) (load flow results)) --
                line_loading = pf_get_all(pf_find(app, "*.ElmLne", "all"), "c:loading", None, False, True)
                # -- Step 6: Busbar voltage (m:u of all * (ElmTerm) (load flow results)) --
                busbar_voltage = pf_get_all(pf_find(app, "*.ElmTerm", "all"), "m:u", None, False, True)
                # -- Step 7: Transformer loading (c:loading of all * (ElmTr2) (load flow results)) --
                transformer_loading = pf_get_all(pf_find(app, "*.ElmTr2", "all"), "c:loading", None, False, True)
                # -- Step 8: Add row to results table (Into "N-1", saved at the end automatically) --
                table("N-1").add({
                    "Highest loading (%)": line_loading.max,
                    "Most loaded line": line_loading.at_max,
                    "Lowest voltage (p.u.)": busbar_voltage.min,
                    "Lowest voltage at": busbar_voltage.at_min,
                    "Highest transformer loading (%)": transformer_loading.max,
                    "Most loaded transformer": transformer_loading.at_max,
                }, status=True)
            except StopWorkflow:
                raise
            except Exception as err:
                case_failed(err)
            finally:
                ctx_pop(["N-1"])
        # -- Step 9: Print the worst outage (Worst outage: {worst}, with a line loaded to {loading:.1f} %) --
        worst_outage = table_pick("N-1", "Contingency", "max", "Highest loading (%)", "text")  # Worst outage
        worst_outage_highest = table_pick("N-1", "Highest loading (%)" or "Contingency", "max", "Highest loading (%)", "num")  # Worst outage
        log(text_format("Worst outage: {worst}, with a line loaded to {loading:.1f} %", {"worst": worst_outage, "loading": worst_outage_highest}))
    except StopWorkflow as stop:
        log(f"Stopped: {stop}")
    finally:
        finish(app)

The file goes on with 1,573 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 runs the base load flow; step 3, the Contingency block, gives the base case first and then one pass per line and transformer in service, putting each element back after its pass.
• Steps 4 to 7 run the load flow and read the highest line loading, the lowest busbar voltage and the highest transformer loading, each with where it is (busbars out of service or not energised are left out).
• Step 8 adds one row per outage; an outage whose load flow does not converge keeps its row with the reason.
• After all passes, "Table column" picks the outage with the highest line loading and "Print a message" writes it to the log.

Pitfalls

• Always put the element back. Set outserv to 0 in a finally block. A script stopped between the two calls leaves the project with the line out of service, and every later study wrong.
• A single slack takes up every lost megawatt. Taking out a generator's own transformer loses the generator, and on this example the slack machine G 02 then carries it all: its transformer Trf 06 - 31 reaches 211 % when Trf 29 - 38 (G 09) is out. That is the balancing rule speaking, not the network. Study generator outages with a distributed slack or as their own contingency set.
• Islanded busbars read 0 p.u. An outage that cuts part of the network off leaves its busbars at exactly 0, and IsEnergized() can still return 1 for them after a load flow. Leave out busbars at 0 p.u. when you look for the lowest voltage.
• A load flow that does not converge is a finding, not an error. Record it and go on: it often means voltage collapse or an island without a reference.
• Which rating is 100 %? c:loading uses the rating in the model (the type's rated current with the element's derating). Post-contingency short-term ratings are higher: use a different limit, or a separate check, for them.
• outserv is not switching. It removes the element from the calculation. To model the breakers that clear it, open them (on_off of ElmCoup, or the switch in the cubicle) instead.

Questions

How do I run a contingency analysis in PowerFactory with Python?

Loop over the lines and transformers: SetAttribute("outserv", 1), run app.GetFromStudyCase("ComLdf").Execute(), read c:loading and m:u, then set outserv back to 0 in a finally block.

How do I take an element out of service in a PowerFactory script?

Set its input parameter outserv to 1 with obj.SetAttribute("outserv", 1), and back to 0 to return it to service.

What is the N-1 criterion?

A planning and operating rule: the network must stay within its thermal and voltage limits after the loss of any single element, such as a line, transformer or generator.

Does PowerFactory have its own contingency analysis?

Yes, the Contingency Analysis command (ComSimoutage). A Python loop is a lighter way to the same screening that is easy to adapt; the Builder's N-1 Quick Start writes one for you.

Related

• Automate PowerFactory studies with Python: every guide
• PowerFactory Python API: getting started
• Load flow and results to CSV or Excel
• Short-circuit automation
• Hosting capacity study
• Parameter sweep: PV and QV curves
• RMS and EMT simulation batch runs

Attributes used on this page: ElmLne, ElmTr2, ElmTerm, ElmSym, EvtOutage. All attributes.