Harmonics and THD from a PSCAD run

The method

1. Choose the plot step for the highest harmonic you need (below).
2. Run long enough to settle, and take the window at the end of the run.
3. Take whole cycles: one or two cycles of the fundamental is the usual window.
4. Transform and sum: each harmonic as a percentage of the fundamental, THD from harmonic 2 upwards.

How far the plot step lets the spectrum reach

The highest frequency is half the rate the channel is written at. Keep a margin below it: the harmonics close to that limit are the least accurate.

The calculation

v = run.series("Vbus")
h = thd(v, frequency=50, at=None, cycles=2, harmonics=50)   # at=None: the window ends at the end of the run
print(h.value)            # THD, % of the fundamental
print(h.worst, h.worst_order)   # the largest harmonic, %, and its order
print(h.fundamental)      # the magnitude of harmonic 1, in the channel's unit
spectrum = list(zip(h.orders, h.pct))   # every order, for a table or a bar chart

thd takes the cycles whole cycles that end at at, transforms them, and counts THD the way a grid code states it: the square root of the sum of the squares of harmonics 2 upwards, over harmonic 1. The DC component is left out unless include_dc=True.

The full Python script from the Builder

The Harmonics and THD Quick Start runs the case for 1 s with a 20 us step and a 50 us plot step, takes two cycles of a voltage at the end of the run, tables the THD, the worst harmonic and its order, and draws the spectrum as a bar chart. This is the file with the case set to C:/Studies/inverter/inverter.pscx:

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

#!/usr/bin/env python3
# Harmonics
"""Made with the PSCAD Workflow Builder, https://tool.adjiebrotots.com/pscad-scripter/builder/

BEFORE YOU RUN
  1. Check the SETTINGS below (lines marked REQUIRED).
  2. Have PSCAD installed and licensed; the script opens one, or uses the one already open (see RUN_MODE).
  3. Run it with Python 3.7 or newer and mhi-pscad installed, e.g.  py -3.12 harmonics.py
Results and run.log go to a new folder next to this file (or in OUTPUT_DIR).
"""
# ==== IMPORTS ====
import bisect
import csv
import fnmatch
import importlib
import math
import os
import re
import sys
import time
import warnings

# ==== SETTINGS (edit here) ====
# The name of this workflow: used for the run folder and in the log.
WORKFLOW_NAME = "Harmonics"
# REQUIRED. The PSCAD case to work on: a .pscx path, or the name of a case already in the
# workspace.
PROJECT = "C:/Studies/inverter/inverter.pscx"
# OPTIONAL. A .pswx workspace to load before the case ("" = keep the one that is open).
WORKSPACE = ""
# OPTIONAL. Extra .pslx / .pscx libraries to load before the case.
LIBRARIES = []
# OPTIONAL. Where each run writes its folder of results and run.log ("" = the folder of this
# script).
OUTPUT_DIR = ""
# OPTIONAL. True: put every parameter, layer state and substitution 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 the case's
# own global substitution).
F_NOMINAL = "auto"
# OPTIONAL. "auto", "launch" (start a new PSCAD), "connect" (use the one that is open) or
# "inside" (run from inside PSCAD).
RUN_MODE = "auto"
# OPTIONAL. True: show the PSCAD window while running.
SHOW_PSCAD = False
# OPTIONAL. The Fortran compiler that builds the case, as PSCAD names it ("" = the one PSCAD has
# selected).
FORTRAN_VERSION = ""
# OPTIONAL. Which installed PSCAD to start ("" = the newest).
PSCAD_VERSION = ""
# OPTIONAL. The port PSCAD listens on for automation (0 = let PSCAD choose).
PSCAD_PORT = 0
# OPTIONAL. Seconds to wait for PSCAD to start.
CONNECT_TIMEOUT = 120

# 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
}

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


# ==== WORKFLOW ====
# Note: PSCAD has no frequency-domain harmonic scan: harmonics come out of the time-domain result. Give the run long enough to settle and take the window at the end. The plot step sets how far up the spect...
def main():
    preflight(["matplotlib"])
    app = px_connect()
    try:
        if CHECK_ONLY:
            check_elements(app, [])
            return
        # -- Step 1: Start (Case: inverter.pscx) --
        # -- Step 2: Run simulation (1 s, 20 us step) --
        run = px_run(app, 1, 20, 50, None, events=(), build="auto", read="out", keep_files=False)
        # Did not build: ends the workflow
        if not run.ok and (run.status == "failed"):
            log("The case did not build. The workflow ends here.", "warn")
            return
        # Stopped early: ends the workflow
        if not (run.ok or (run.status == "failed")):
            log("The simulation stopped before the end. The workflow ends here.", "warn")
            return
        # -- Step 3: Add row to results table (Into "Harmonics", saved at the end automatically) --
        vbus = run.series("Vbus", "")  # Read a channel
        thd_result = thd(vbus, frequency=f_nominal(app), at=None, cycles=2, harmonics=50, include_dc=False)  # Harmonic distortion (THD)
        table("Harmonics").add({
            "THD (%)": thd_result.value,
            "Worst (%)": thd_result.worst,
            "Worst order": thd_result.worst_order,
        })
        # -- Step 4: Draw 2D graph (Bars) --
        plot_xy(thd_result.orders, [thd_result.pct], labels=["Harmonic"], kind="bar", title="", xlabel="Harmonic order", xunit="", ylabel="Percent of the fundamental", yunit="", limits=[], stem="spectrum")
    except StopWorkflow as stop:
        log(f"Stopped: {stop}")
    finally:
        finish(app)

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

How the script works

• The 50 us plot step holds frequencies up to 10 kHz, so the 50 harmonics asked for are well inside what the file has.
• thd(vbus, frequency=f_nominal(app), ...) reads the nominal frequency from the case's own global substitution when F_NOMINAL is "auto".
• The window ends at the end of the run. Make the run long enough that the case has settled well before that.
• plot_xy(thd_result.orders, [thd_result.pct], kind="bar") is the spectrum as a bar chart, each harmonic as a percentage of the fundamental.

Pitfalls

• Harmonics above the plot step's limit. They are not in the file, and the transform folds anything near the limit back down as a false lower harmonic. Lower the plot step, or ask for fewer harmonics.
• A window during a transient. The start-up of a case, a fault and a step all spread energy across the spectrum. Take the window once the waveform repeats cycle after cycle.
• A fundamental that is not what you said. The window is cycles of the frequency you give. A case running at 60 Hz analysed at 50 Hz gives a meaningless spectrum. F_NOMINAL = "auto" reads it from the case.
• A three phase channel is three channels. A channel titled Vbus on a three phase signal is recorded as Vbus:1, Vbus:2 and Vbus:3. Analyse each phase, or the one the requirement names.

Questions

How do I calculate THD from a PSCAD simulation?

Record the voltage or current with an output channel, run until the case has settled, then take a whole number of cycles at the end and work out the Fourier transform. THD is the root of the sum of the squares of harmonics 2 upwards, over the fundamental. The Builder's Harmonic distortion block does this from any recorded channel.

Does PSCAD have a harmonic load flow?

No. PSCAD is a time-domain simulator, so harmonics are read from the simulated waveform rather than calculated in the frequency domain. Its on-line frequency scanner does the same during the run.

What plot step do I need for harmonics in PSCAD?

At least twice the highest frequency you need, with a margin: a 250 us plot step reaches about 2 kHz, a 50 us plot step about 10 kHz.

Related

• Automate PSCAD studies with Python: every guide
• PSCAD automation with Python: getting started
• Build and run a PSCAD case from Python
• Fault studies in PSCAD with Python: ride through and critical clearing time
• Sweep a PSCAD parameter with Python
• Snapshots, simulation sets and batch runs
• Control tuning: overshoot, rise time and damping

Every parameter name a script writes is in the component reference, read out of PSCAD itself. The Workflow Builder writes the Python from blocks.