Parameter sweeps in PowerFactory with Python: PV and QV curves

The method

1. Keep the original values of everything the sweep changes, before the first step.
2. Set the next value from the original, not from the previous step (multiplying the previous value compounds).
3. Run and check: Execute() returns 0 when the load flow converged. For a PV curve the first failure is the end of the curve; for most sweeps it is a point to record and pass.
4. Read what you sweep for: busbar voltages m:u, the delivered load m:P:bus1, a condenser's reactive power m:Q:bus1.
5. Put everything back in a finally block, and delete any object the sweep made.

PV curve: scale every load to the nose

Each step multiplies the active and reactive power of the three loads by the same factor, so their power factor stays as it is. The extra power comes from the slack generator; generator reactive limits follow the study case (ignored in this example).

pv_curve.py

"""PV curve in PowerFactory with Python: scale every load up step by step until the load flow no
longer converges, recording the load busbar voltages on the way. The last point is the nose.

Runs on the "Nine-bus 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 = "Nine-bus System"
STUDY_CASE = "01- Load Flow"
WATCH = ["Bus 5", "Bus 6", "Bus 8"]            # the busbars of the three loads
STEP = 0.05                                    # load scaling step

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

ldf = app.GetFromStudyCase("ComLdf")
loads = app.GetCalcRelevantObjects("*.ElmLod")
base = [(ld, ld.GetAttribute("plini"), ld.GetAttribute("qlini")) for ld in loads]
watch = [app.GetCalcRelevantObjects(name + ".ElmTerm")[0] for name in WATCH]

curve = []
scale = 1.0
try:
    while scale <= 5.0:
        for ld, p, q in base:                  # same power factor: P and Q scale together
            ld.SetAttribute("plini", p * scale)
            ld.SetAttribute("qlini", q * scale)
        if ldf.Execute() != 0:                 # no solution any more: past the nose
            break
        point = {"Scale": scale, "Load (MW)": sum(ld.GetAttribute("m:P:bus1") for ld in loads)}
        point.update({b.loc_name: b.GetAttribute("m:u") for b in watch})
        curve.append(point)
        scale = round(scale + STEP, 6)
finally:
    for ld, p, q in base:                      # put the loads back
        ld.SetAttribute("plini", p)
        ld.SetAttribute("qlini", q)

nose = curve[-1]
weakest = min(WATCH, key=lambda name: nose[name])
print(f"Nose: {nose['Load (MW)']:.0f} MW ({nose['Scale']:.2f} x the base load), "
      f"weakest busbar {weakest} at {nose[weakest]:.3f} p.u.")

with open("pv_curve.csv", "w", newline="") as f:
    out = csv.DictWriter(f, fieldnames=list(curve[0]))
    out.writeheader()
    out.writerows(curve)

QV curve: a fictitious condenser at Bus 5

A static generator with no active power, in voltage control, holds Bus 5 at each voltage from 1.10 p.u. down; the reactive power it injects is the curve. It is rated far above what it will need, so it never limits, and deleted at the end.

qv_curve.py

"""QV curve in PowerFactory with Python: hold one busbar at a set voltage with a fictitious
voltage-controlling static generator, step the voltage down and record the reactive power it
needs. The bottom of the curve is the reactive power margin of the busbar.

Runs on the "Nine-bus 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 = "Nine-bus System"
STUDY_CASE = "01- Load Flow"
BUS = "Bus 5"
VOLTAGES = [round(1.10 - 0.02 * k, 2) for k in range(36)]      # 1.10 p.u. down to 0.40 p.u.

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

ldf = app.GetFromStudyCase("ComLdf")
bus = app.GetCalcRelevantObjects(BUS + ".ElmTerm")[0]
grid = bus.GetAttribute("cpGrid")                # the grid (ElmNet) the busbar belongs to

# The fictitious condenser: no active power, voltage control, rated well above what it will need.
cubicle = bus.CreateObject("StaCubic", "QV cubicle")
gen = grid.CreateObject("ElmGenstat", "QV condenser")
gen.SetAttribute("bus1", cubicle)
gen.SetAttribute("sgn", 2000.0)                  # MVA
gen.SetAttribute("pgini", 0.0)
gen.SetAttribute("av_mode", "constv")            # hold the voltage at usetp

curve = []
try:
    for v in VOLTAGES:
        gen.SetAttribute("usetp", v)
        if ldf.Execute() != 0:
            break
        curve.append({"Voltage (p.u.)": bus.GetAttribute("m:u"), "Q (Mvar)": gen.GetAttribute("m:Q:bus1")})
finally:
    app.ResetCalculation()                       # results first, then the objects can go
    gen.Delete()
    cubicle.Delete()

bottom = min(curve, key=lambda p: p["Q (Mvar)"])
if bottom is curve[-1]:
    print(f"{BUS}: Q still falls at {bottom['Voltage (p.u.)']:.2f} p.u.; extend VOLTAGES to find the bottom")
else:
    print(f"{BUS}: the curve bottoms out at {bottom['Q (Mvar)']:.0f} Mvar, at {bottom['Voltage (p.u.)']:.2f} p.u.")

with open("qv_curve.csv", "w", newline="") as f:
    out = csv.DictWriter(f, fieldnames=list(curve[0]))
    out.writeheader()
    out.writerows(curve)

The full Python script from the Builder

The Builder's PV curve Quick Start is the same sweep as a loop of blocks, ending at the first load flow that fails. On the Nine-bus example it gives the same 28 points and the same nose at 2.35 times the base load. This is its file for that example:

pv_curve_by_load_scaling.py as the Builder writes it: the settings and the workflow⬇ Download the full file (.py, 2,784 lines)

#!/usr/bin/env python3
# PV curve by load scaling
"""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 pv_curve_by_load_scaling.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 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 = "PV curve by load scaling"
# REQUIRED. The PowerFactory project to open, exactly as the Data Manager shows it (letter case
# counts).
PROJECT = "Nine-bus System"
# OPTIONAL. The study case to activate first ("" = keep the active one).
STUDY_CASE = "01- Load 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 = ""

# 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. Load scaling: First.
LOAD_SCALING_FROM = 1

# OPTIONAL. Load scaling: Last.
LOAD_SCALING_TO = 3

# OPTIONAL. Load scaling: Step.
LOAD_SCALING_STEP = 0.05

# 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 scaling factor this multiplies the active and reactive power of every load (or of the loads named on Loads) by the factor, from their original values, runs the load flow and reads the lowe...
def main():
    preflight(["matplotlib"])
    app = pf_connect()
    try:
        if CHECK_ONLY:
            check_elements(app, [("*.ElmLod", "all"), ("*.ElmTerm", "all")])
            return
        # -- Step 1: Start (Project: Nine-bus System) --
        # -- Step 2: Load scaling (1 to 3, step 0.05) --
        load_scaling_values = num_range(LOAD_SCALING_FROM, LOAD_SCALING_TO, LOAD_SCALING_STEP)
        for load_scaling_pass_no, load_scaling in enumerate(load_scaling_values, 1):
            ctx_push("Load scaling", load_scaling, load_scaling_pass_no, len(load_scaling_values))
            mark = pf_mark()
            try:
                # -- Step 3: Loads (all * (ElmLod): plini, qlini x the wired factor) --
                scaled = pf_find(app, "*.ElmLod", "all")
                factor = load_scaling
                scale_attr(scaled, "plini", factor)
                scale_attr(scaled, "qlini", factor)
                # -- Step 4: Run load flow (Settings of the study case) --
                ldf = run_load_flow(app, "case", "case", "case")
                # Run load flow: Did not converge
                if not ldf.ok:
                    # -- Step 5: Nose of the curve (The load flow no longer converges: the nose of the PV curve) --
                    log("The load flow no longer converges: the nose of the PV curve")
                    break
                # -- 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: Add row to results table (Into "PV curve", saved at the end automatically) --
                table("PV curve").add({
                    "Lowest voltage (p.u.)": busbar_voltage.min,
                    "Lowest at": busbar_voltage.at_min,
                }, status=True)
            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(["PV curve"])
        # -- Step 8: PV curve (Lowest voltage (p.u.) against Load scaling) --
        plot_table("PV curve", "Load scaling", ["Lowest voltage (p.u.)"], kind="line", group=None, title="PV curve", ylabel="Lowest busbar voltage", yunit="p.u.", limits=[], stem=None)
    except StopWorkflow as stop:
        log(f"Stopped: {stop}")
    finally:
        finish(app)

The file goes on with 2,625 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 is a "For a range of numbers" loop over the scaling factor, from LOAD_SCALING_FROM to LOAD_SCALING_TO in steps of LOAD_SCALING_STEP; it puts the loads back after every pass.
• Step 3, Scale, sets plini and qlini of every load to the factor times its original value.
• Step 4 runs the load flow. When it does not converge, "Stop the loop" ends the sweep with a message in the log: the nose.
• Step 6 reads the lowest busbar voltage and where it is; step 7 adds the row.
• After the loop, "Draw graph from results table" saves the PV curve as a PNG.

Pitfalls

• The last converged point is an estimate of the nose. A load flow can also fail for numerical reasons. Use a smaller step near the end, or the continuation method, when the exact margin matters.
• Scale P and Q together to keep each load's power factor; scaling only P sweeps the power factor too.
• Where the extra power comes from shapes the curve. Here the slack takes it all. A transfer from one area to another, or generator reactive limits (iopt_lim on ComLdf), move the nose.
• Voltage-dependent loads draw less as the voltage falls. Record the delivered power m:P:bus1, not the setpoint plini, on the x axis.
• The QV condenser must never limit. Rate it well above the reactive power it will need, run without its limits, and delete it (after ResetCalculation()) when the sweep ends.
• Put the inputs back from the stored originals in finally, so a stopped sweep never leaves the loads scaled.

Questions

What is a PV curve?

The voltage of a busbar plotted against the power delivered as load or transfer grows. Its tip, the nose, is the maximum loadability; beyond it the load flow has no solution.

What is a QV curve?

The reactive power a fictitious condenser must inject to hold a busbar at each voltage. The bottom of the curve marks the critical voltage, and its depth below zero is the reactive power margin.

How do I run a parameter sweep in PowerFactory with Python?

Store the original value, then for each step set it with SetAttribute, run the calculation, check that Execute() returned 0 and read the results. Put the original back in a finally block.

How do I scale all loads in PowerFactory with Python?

Loop over app.GetCalcRelevantObjects("*.ElmLod") and set plini and qlini to the factor times their original values. The load's scaling factor scale0 is another way.

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
• RMS and EMT simulation batch runs

Attributes used on this page: ElmLod, ElmTerm, ElmGenstat, TypLod. All attributes.