How to do a hosting capacity study, and automate it in PowerFactory

How to do a hosting capacity study

The method does not depend on the software. Every step below is one decision you write down in the study report.

1. Choose the operating cases. The limit depends on the state of the network: for new generation the hardest case is usually minimum demand with the existing generation high (voltage rise, reverse flow); for new load, maximum demand. Run each case.
2. Set the limits from the planning standard: the rating each branch may reach and the voltage band each busbar must keep.
3. Model the connection as a generator at the connection point, at the power factor it will run at (unity here).
4. Decide where its power goes. In a distribution network the upstream grid takes it (the slack). In a transmission network, other generators turn down: pro rata, as here, or in merit order. This choice changes the answer, so state it.
5. Search for the largest output that passes. Start from a passing value (0 MW) and a failing one (the most that could ever connect), try the middle, keep the half where the limit lies, and stop when the gap is within the tolerance. From 0 to 2,000 MW with a 5 MW tolerance that is 11 load flows (the two ends and 9 halvings), against 400 for a sweep in 5 MW steps.
6. Record the hosting capacity and what limited it (the binding constraint), then repeat at every candidate busbar: the result is the hosting capacity map.
7. Test the sensitivities: another power factor, other demand cases, an outage, the projects already in the queue.

In PowerFactory

The new generation is a static generator (ElmGenstat) connected through a new cubicle (StaCubic) on the busbar, at constant reactive power 0. Each try sets its active power pgini, turns the existing synchronous generators down by the same total, runs ComLdf, and checks c:loading on every line and transformer and m:u on every busbar. Everything is put back and the test generator deleted after each busbar.

hosting_capacity.py

"""Hosting capacity study in PowerFactory with Python: at each busbar, the largest output a new
generator can have before a line or transformer overloads or a voltage leaves its band.

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_MAX = 100.0              # % of rating, every line and transformer
VOLTAGE_BAND = (0.90, 1.10)      # p.u., every busbar
SEARCH_MAX = 2000.0              # MW, the largest output tried
TOLERANCE = 5.0                  # MW, when to stop halving

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")
# The new generation displaces the generators in service pro rata, so the slack only covers losses.
dispatch = [(g, g.GetAttribute("pgini")) for g in app.GetCalcRelevantObjects("*.ElmSym") if g.outserv == 0]
total = sum(p for g, p in dispatch)


def limits_hold():
    """Run the load flow. (True, "") when every limit holds, else (False, the first one broken)."""
    if ldf.Execute() != 0:
        return False, "load flow does not converge"
    loading, branch = max((o.GetAttribute("c:loading"), o.loc_name) for o in branches if o.outserv == 0)
    if loading > LOADING_MAX:
        return False, f"{branch} loaded to {loading:.0f} %"
    volts = [(b.GetAttribute("m:u"), b.loc_name) for b in buses if b.GetAttribute("m:u") > 0]
    if max(volts)[0] > VOLTAGE_BAND[1]:
        return False, f"{max(volts)[1]} at {max(volts)[0]:.3f} p.u."
    if min(volts)[0] < VOLTAGE_BAND[0]:
        return False, f"{min(volts)[1]} at {min(volts)[0]:.3f} p.u."
    return True, ""


def inject(gen, p_mw):
    """The test generator at p_mw, every other generator turned down by the same total."""
    gen.SetAttribute("pgini", p_mw)
    for g, p in dispatch:
        g.SetAttribute("pgini", p * (1 - p_mw / total))


results = []
for bus in [b for b in buses if b.uknom >= 100]:          # the 345, 230 and 138 kV busbars
    cubicle = bus.CreateObject("StaCubic", "HC cubicle")
    gen = bus.GetAttribute("cpGrid").CreateObject("ElmGenstat", "HC test generator")
    gen.SetAttribute("bus1", cubicle)
    gen.SetAttribute("sgn", 1.5 * SEARCH_MAX)           # MVA, never the limit
    gen.SetAttribute("av_mode", "constq")               # fixed reactive power...
    gen.SetAttribute("qgini", 0.0)                      # ...of zero: unity power factor
    try:
        inject(gen, 0.0)
        ok, why = limits_hold()
        if not ok:
            results.append({"Busbar": bus.loc_name, "Hosting capacity (MW)": 0.0, "Limited by": why})
            continue
        low, high = 0.0, SEARCH_MAX
        inject(gen, high)
        ok, why = limits_hold()
        if ok:
            results.append({"Busbar": bus.loc_name, "Hosting capacity (MW)": high, "Limited by": "nothing up to the search limit"})
            continue
        while high - low > TOLERANCE:                    # halve the gap between pass and fail
            mid = (low + high) / 2
            inject(gen, mid)
            ok, reason = limits_hold()
            if ok:
                low = mid
            else:
                high, why = mid, reason
        results.append({"Busbar": bus.loc_name, "Hosting capacity (MW)": low, "Limited by": why})
    finally:
        for g, p in dispatch:                            # put the dispatch back
            g.SetAttribute("pgini", p)
        app.ResetCalculation()
        gen.Delete()
        cubicle.Delete()

for r in sorted(results, key=lambda r: r["Hosting capacity (MW)"]):
    print(f"{r['Busbar']:7} {r['Hosting capacity (MW)']:7.0f} MW  limited by {r['Limited by']}")

with open("hosting_capacity.csv", "w", newline="") as f:
    out = csv.DictWriter(f, fieldnames=["Busbar", "Hosting capacity (MW)", "Limited by"])
    out.writeheader()
    out.writerows(results)

Result: the hosting capacity map of the 39-bus example

The 30 busbars of 100 kV and above take between 8 MW (Bus 22) and 1,742 MW (Bus 03), 4.6 s of load flows in all. Every one is limited by a line rating, not by voltage: at transmission level, a generator at unity power factor moves the voltage little. The base case already loads Line 21 - 22 to 99.1 %, so busbars that push power through it can take almost nothing (Bus 22 and Bus 23), and Line 05 - 06 sets the limit of 5 busbars (Bus 06, Bus 11, Bus 12, Bus 10, Bus 07).

The full Python script from the Builder

The Builder's Quick Start is the same study as blocks: a For each element loop over the busbars, a Find limit search inside it, a test generator, and three limit checks. Run on the 39-bus example with the busbars of 100 kV and above, it gives every busbar within 7.5 MW of the script above. This is its file for that example:

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

#!/usr/bin/env python3
# Hosting capacity at every busbar
"""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 hosting_capacity_at_every_busbar.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 = "Hosting capacity at every busbar"
# 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 = ""

# 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. Connection point: Elements. A name (ElmTerm), exactly as in the Data Manager, or a
# pattern with * and ?: every element that matches is used.
CONNECTION_POINT_EL = "*"

# OPTIONAL. Generation (MW): Lowest value.
GENERATION_LOWER = 0

# OPTIONAL. Generation (MW): Highest value.
GENERATION_UPPER = 2000

# OPTIONAL. Generation (MW): Tolerance.
GENERATION_TOL = 5

# OPTIONAL. Line thermal limit: Highest allowed.
LINE_THERMAL_LIMIT_HI = 100

# OPTIONAL. Voltage limits: Lowest allowed.
VOLTAGE_LIMITS_LO = 0.9

# OPTIONAL. Voltage limits: Highest allowed.
VOLTAGE_LIMITS_HI = 1.1

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


# ==== WORKFLOW ====
# Note: At each busbar of Connection point this connects a test generator at unity power factor and searches its largest output for which the load flow converges, no line or 2-winding transformer is above...
def main():
    preflight(["matplotlib"])
    app = pf_connect()
    try:
        if CHECK_ONLY:
            check_elements(app, [("*.ElmSym", "all"), ("*.ElmLne", "all"), ("*.ElmTr2", "all"), ("*.ElmTerm", "all")])
            return
        # -- Step 1: Start (Project: 39 Bus New England System) --
        # -- Step 2: Generator dispatch (pgini of all * (ElmSym)) --
        generator_dispatch = pf_get_all(pf_find(app, "*.ElmSym", "all"), "pgini", None, False, True)
        # -- Step 3: Connection point (Each of all * (ElmTerm) where uknom >= 100) --
        connection_point_elements = find_elements(app, "ElmTerm", CONNECTION_POINT_EL, "uknom", ">=", "100", True)
        if not connection_point_elements:
            log("Connection point: no element to loop over.", "warn")
        for connection_point_pass_no, connection_point in enumerate(connection_point_elements, 1):
            ctx_push("Connection point", connection_point, connection_point_pass_no, len(connection_point_elements))
            mark = pf_mark()
            try:
                # -- Step 4: Generation (MW) (0 to 2000, within 5) --
                def try_n4(generation_mw):
                    trial = Trial({"Generation (MW)": generation_mw}, "min")
                    ctx_push("Generation (MW)", generation_mw)
                    mark_n4 = pf_mark()
                    try:
                        # -- Step 5: Test generator (At the wired busbar: wired MW, 0 Mvar) --
                        add_test_generator(app, connection_point, generation_mw, 0, "n5", None)
                        # -- Step 6: Other generators (all * (ElmSym): pgini x the wired factor) --
                        total_dispatch = stat("sum", generator_dispatch.values)  # Total dispatch
                        share_kept = apply_formula(lambda p, total: 1 - p / total, [generation_mw, total_dispatch], False, "Share kept")  # Share kept
                        scale_attr(pf_find(app, "*.ElmSym", "all"), "pgini", share_kept)
                        # -- Step 7: Run load flow (Settings of the study case) --
                        ldf = run_load_flow(app, "case", "case", "case")
                        # Did not converge: ends this try of Generation (MW)
                        if not ldf.ok:
                            case_skipped("Load flow did not converge")
                            return trial
                        # -- Step 8: 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 9: 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 10: 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 11: Line thermal limit (At most 100, goes on either way) --
                        line_thermal_limit_result = check_limits(line_loading.values, lo=None, hi=LINE_THERMAL_LIMIT_HI, after=None, before=None)
                        if line_thermal_limit_result.ok is None:
                            log("Line thermal limit: NO DATA (" + str(line_thermal_limit_result.get("reason") or "nothing to check") + "), so it did not pass.", "warn")
                        # -- Step 12: Transformer thermal limit (At most 100, goes on either way) --
                        transformer_thermal_limit_result = check_limits(transformer_loading.values, lo=None, hi=100, after=None, before=None)
                        if transformer_thermal_limit_result.ok is None:
                            log("Transformer thermal limit: NO DATA (" + str(transformer_thermal_limit_result.get("reason") or "nothing to check") + "), so it did not pass.", "warn")
                        # -- Step 13: Voltage limits (Between 0.9 and 1.1, goes on either way) --
                        voltage_limits_result = check_limits(busbar_voltage.values, lo=VOLTAGE_LIMITS_LO, hi=VOLTAGE_LIMITS_HI, after=None, before=None)
                        if voltage_limits_result.ok is None:
                            log("Voltage limits: NO DATA (" + str(voltage_limits_result.get("reason") or "nothing to check") + "), so it did not pass.", "warn")
                        thermal_limits = logic("and", line_thermal_limit_result.ok, transformer_thermal_limit_result.ok)  # Thermal limits
                        all_limits = logic("and", thermal_limits, voltage_limits_result.ok)  # All limits
                        trial.report(ok=all_limits)
                    except StopWorkflow:
                        raise
                    except Exception as err_n4:
                        case_failed(err_n4, trial)
                    finally:
                        try:
                            pf_undo(mark_n4)
                        except Exception as fin_n4:
                            log(f"Could not put things back after this pass: {fin_n4}", "error")
                        ctx_pop()
                    return trial
                generation_mw_result = find_limit(try_n4, GENERATION_LOWER, GENERATION_UPPER, GENERATION_TOL, max_steps=20, passes="low", bracket="given", on_unknown="fail", label="Generation (MW)", integer=False)
                # -- Step 14: Add row to results table (Into "Hosting capacity", saved at the end automatically) --
                table("Hosting capacity").add({
                    "Hosting capacity (MW)": generation_mw_result.limit,
                    "Search": generation_mw_result.status,
                }, 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(["Hosting capacity"])
        # -- Step 15: Hosting capacity at each busbar (Hosting capacity (MW) against Connection point) --
        plot_table("Hosting capacity", "Connection point", ["Hosting capacity (MW)"], kind="bar", group=None, title="Hosting capacity at each busbar", ylabel="Hosting capacity", yunit="MW", limits=[], stem=None)
    except StopWorkflow as stop:
        log(f"Stopped: {stop}")
    finally:
        finish(app)

The file goes on with 3,517 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 reads pgini of every synchronous generator, and "Total dispatch" adds them up: the existing generation that will make room.
• Step 3 loops over the busbars that match CONNECTION_POINT_EL; for this example a filter keeps those with uknom of 100 kV or more.
• Step 4, Find limit, calls try_n4 with one output after another and halves the gap between GENERATION_LOWER and GENERATION_UPPER down to GENERATION_TOL. A try whose load flow does not converge counts as a fail.
• Inside each try, add_test_generator connects (or moves) the static generator, "Share kept" works out 1 minus the output over the total dispatch, and Scale sets every generator to that share of its original pgini; pf_undo puts them back after the try.
• The three checks (LINE_THERMAL_LIMIT_HI, the transformer limit, VOLTAGE_LIMITS_LO and _HI) are combined with And blocks into the pass or fail the search reads.
• After the search, the busbar's row gets the hosting capacity and the search status; after the loop, the bar chart is drawn and the test generator deleted.

Pitfalls

• One operating case is one answer. Hosting capacity at minimum demand and at maximum demand can differ several times over. Run the cases that stress each limit, and report the lowest.
• The balancing rule changes the result. With a single slack taking every megawatt, the slack machine's own transformer soon becomes the "limit" at every busbar. Turn existing generation down (pro rata, as here) or let an upstream grid take it, and say which.
• Bisection assumes one crossing. It finds the edge between pass and fail when output only makes things worse. If more output first relieves a line and then loads it the other way, check a few points either side of the answer.
• Check the ratings. c:loading is measured against the rating in the model, the type's rated current times the element's derating. An out-of-date rating gives an out-of-date hosting capacity.
• Unity power factor is a choice. A connection that absorbs reactive power, or controls voltage, moves the voltage limit a long way: on LV and MV feeders it is often the difference that matters.
• Leave the network as you found it. Put the dispatch back and delete the test generator and its cubicle after each busbar (after ResetCalculation()), as the script does in finally.
• Load flow limits are not the whole story. Fault level, protection, voltage step and power quality limits can bind first; treat this result as an upper bound until they are checked.

Questions

What is hosting capacity?

The largest amount of new generation or load that can connect at a point of the network without breaking a limit, such as a thermal rating or a voltage limit, in the operating cases studied.

How is hosting capacity calculated?

Connect the new generation at the point, raise its output, and run a load flow at each step until the first limit breaks. A bisection search finds the edge in about ten load flows. The largest output that passes is the hosting capacity.

What limits hosting capacity?

Usually thermal ratings of lines and transformers and voltage rise; also voltage step on connection, fault level, protection and power quality.

What is a hosting capacity map?

The hosting capacity of every connection point in a network, shown on a map or a chart, often with the limit that binds at each point.

How do I do a hosting capacity study in PowerFactory?

Add a static generator (ElmGenstat) at the busbar, set its pgini, run ComLdf, check c:loading and m:u, and search the largest pgini that passes. The script and the Builder Quick Start on this page do it for every busbar.

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
• Parameter sweep: PV and QV curves
• RMS and EMT simulation batch runs

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