powerio-dist 0.3.0

Multiconductor distribution network model and lossless converters for OpenDSS, PMD JSON, and BMOPF JSON.
Documentation 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180

"""Re-solve emitted .dss cases against their originals.

Usage:
    cargo test -p powerio-dist --test matrix -- --ignored emit_for_physics_check
    <python-with-opendssdirect> powerio-dist/tools/physics_check.py

For every dss sourced fixture the harness writes three regenerated cases
under target/physics (canonical, via BMOPF, via PMD). This script solves
each against the original and reports the maximum per node voltage
deviation in per unit of the original node magnitude (nodes below 1 volt
are compared absolutely, in volts). The conversion contract bound is 1e-8.
"""

import glob
import os
import sys

ORIGINALS = {
    "opendss_ieee13_IEEE13Nodeckt": "tests/data/dist/opendss/ieee13/IEEE13Nodeckt.dss",
    "opendss_ieee34_ieee34Mod1": "tests/data/dist/opendss/ieee34/ieee34Mod1.dss",
    "opendss_ieee123_IEEE123Master": "tests/data/dist/opendss/ieee123/IEEE123Master.dss",
    "micro_xfmr_single_phase": "tests/data/dist/micro/xfmr_single_phase.dss",
    "micro_xfmr_center_tap": "tests/data/dist/micro/xfmr_center_tap.dss",
    "micro_xfmr_wye_delta": "tests/data/dist/micro/xfmr_wye_delta.dss",
    "micro_xfmr_delta_wye": "tests/data/dist/micro/xfmr_delta_wye.dss",
    "micro_xfmr_open_wye_open_delta": "tests/data/dist/micro/xfmr_open_wye_open_delta.dss",
    "micro_xfmr_1ph_delta_wye": "tests/data/dist/micro/xfmr_1ph_delta_wye.dss",
    "micro_switch": "tests/data/dist/micro/switch.dss",
    "micro_fourwire_linecode": "tests/data/dist/micro/fourwire_linecode.dss",
    "micro_defaults_degenerate": "tests/data/dist/micro/defaults_degenerate.dss",
    "micro_linecode_10x10": "tests/data/dist/micro/linecode_10x10.dss",
}


def solve(path):
    import opendssdirect as dss

    # The converter drops voltage regulator controls by documented policy
    # (RegControl becomes a fixed tap transformer). The cases run their own
    # Solve while loading, so control actions must be off before that:
    # inject the option right after the circuit line on both sides.
    with open(path, encoding="utf-8", errors="replace") as f:
        text = f.read()
    lines = text.splitlines()
    injected = False
    for i, line in enumerate(lines):
        head = line.lower().lstrip()
        # Both circuit spellings appear in the vendored masters: the writer
        # emits "New Circuit.x", ieee34/ieee123 use "New object=circuit.x".
        if head.startswith("new circuit") or head.startswith("new object=circuit"):
            # Tight solver tolerance: the default 1e-4 pu would swamp the
            # 1e-8 conversion bound with convergence noise.
            lines.insert(i + 1, "Set Controlmode=OFF")
            lines.insert(i + 2, "Set tolerance=0.0000000001")
            injected = True
            break
    if not injected:
        raise SystemExit(f"{path}: no circuit definition found to stage")
    staged = os.path.join(os.path.dirname(os.path.abspath(path)), "_staged_" + os.path.basename(path))
    with open(staged, "w") as f:
        f.write("\n".join(lines) + "\n")

    try:
        dss.Text.Command("Clear")
        dss.Text.Command(f'Redirect "{os.path.abspath(staged)}"')
        dss.Text.Command("Set Controlmode=OFF")
        dss.Text.Command("Solve")
    finally:
        os.unlink(staged)
    if not dss.Solution.Converged():
        return None
    volts = {}
    for bus in dss.Circuit.AllBusNames():
        dss.Circuit.SetActiveBus(bus)
        nodes = dss.Bus.Nodes()
        raw = dss.Bus.Voltages()
        for k, node in enumerate(nodes):
            volts[f"{bus}.{node}"] = complex(raw[2 * k], raw[2 * k + 1])
    return volts


def compare(base, emitted):
    # Deviation in per unit of the bus's own voltage scale (the largest
    # node magnitude at the bus), so near zero neutral nodes compare
    # against the working voltage, not their own tiny magnitude.
    bus_base = {}
    for node, v0 in base.items():
        bus = node.rsplit(".", 1)[0]
        bus_base[bus] = max(bus_base.get(bus, 0.0), abs(v0))
    worst = 0.0
    worst_node = ""
    for node, v0 in base.items():
        v1 = emitted.get(node)
        if v1 is None:
            return None, f"missing node {node}"
        base_v = max(bus_base[node.rsplit(".", 1)[0]], 1.0)
        dev = abs(v1 - v0) / base_v
        if dev > worst:
            worst, worst_node = dev, node
    return worst, worst_node


# Cells whose deviation has a documented cause. Bounds above 1e-8 carry
# the reason; "loss" cells are format losses every conversion reports in
# its warnings (constant power only loads in BMOPF, the center tap
# collapse, an unsupported transformer shape, no vminpu field in the
# ENGINEERING model). The engine seeding entries cover OpenDSS treating
# written properties differently from untouched defaults (an untouched
# load seeds VBase 7200 V; writing kv=12.47 computes 12470/sqrt(3)),
# amplified near vminpu boundaries.
DOCUMENTED = {
    ("micro_defaults_degenerate", "canonical"): (1e-6, "engine seeding asymmetry"),
    ("micro_defaults_degenerate", "via_pmd"): (1e-6, "engine seeding asymmetry"),
    ("micro_defaults_degenerate", "via_bmopf"): (1e-2, "BMOPF: constant power loads only"),
    ("opendss_ieee13_IEEE13Nodeckt", "via_bmopf"): (1e-1, "BMOPF: constant power loads only"),
    ("opendss_ieee34_ieee34Mod1", "via_bmopf"): (1e-1, "BMOPF: constant power loads only"),
    ("opendss_ieee34_ieee34Mod1", "via_pmd"): (1e-1, "no vminpu field in ENGINEERING"),
    ("opendss_ieee123_IEEE123Master", "via_bmopf"): (None, "transformer shape outside the four BMOPF subtypes"),
    ("opendss_ieee123_IEEE123Master", "via_pmd"): (1e-2, "regulator bank restatement"),
    ("micro_xfmr_center_tap", "via_bmopf"): (2e-1, "BMOPF: center tap collapses to two windings"),
    # Open wye / open delta and the single phase delta-wye leg. BMOPF
    # single_phase carries the terminals and impedance but no wye/delta label,
    # so the open delta secondary's floating ground reference is not
    # reproduced (its line to line voltages still agree); the grounded delta-wye
    # leg only shifts by the single phase delta winding's impedance base. PMD
    # reads a single phase delta winding back as a two phase transformer, which
    # re-solves differently. The canonical dss regeneration is exact for both.
    ("micro_xfmr_open_wye_open_delta", "via_bmopf"): (2e-1, "BMOPF: open delta floating reference, no wye/delta label"),
    ("micro_xfmr_open_wye_open_delta", "via_pmd"): (1.5e0, "PMD: single phase delta winding read as two phase"),
    ("micro_xfmr_1ph_delta_wye", "via_bmopf"): (1e-3, "BMOPF: single phase delta winding impedance base, no wye/delta label"),
    ("micro_xfmr_1ph_delta_wye", "via_pmd"): (1e0, "PMD: single phase delta winding read as two phase"),
    ("micro_xfmr_single_phase", "via_pmd"): (1e-6, "engine Z1/Z0 vs MVAsc input path"),
    # PMD models a dss switch as a 1e-7 ohm series element while the engine's
    # switch dummy works out near 1e-3 ohm over the forced length.
    ("micro_switch", "via_pmd"): (1e-5, "ENGINEERING switch impedance convention"),
    ("micro_xfmr_center_tap", "via_pmd"): (1e-6, "engine Z1/Z0 vs MVAsc input path"),
}


def main():
    failures = 0
    for stem, original in ORIGINALS.items():
        emitted_paths = sorted(glob.glob(f"target/physics/{stem}.*.dss"))
        if not emitted_paths:
            # An empty glob must fail, or the gate silently checks nothing
            # (forgotten emit step, renamed fixture).
            print(f"{stem}: NO EMITTED CASES under target/physics (run the emit test first)")
            failures += 1
            continue
        base = solve(original)
        if base is None:
            print(f"{stem}: ORIGINAL DID NOT CONVERGE")
            failures += 1
            continue
        for emitted_path in emitted_paths:
            kind = emitted_path.rsplit(".", 2)[-2]
            bound, reason = DOCUMENTED.get((stem, kind), (1e-8, None))
            emitted = solve(emitted_path)
            if emitted is None:
                print(f"{stem} [{kind}]: DID NOT CONVERGE")
                failures += 1
                continue
            worst, where = compare(base, emitted)
            if worst is None:
                if bound is None:
                    print(f"{stem} [{kind}]: {where} (documented: {reason})")
                else:
                    print(f"{stem} [{kind}]: {where}")
                    failures += 1
            elif bound is not None and worst <= bound:
                note = f" (documented: {reason})" if reason else ""
                print(f"{stem} [{kind}]: max deviation {worst:.3e} at {where} ok{note}")
            else:
                print(f"{stem} [{kind}]: max deviation {worst:.3e} at {where} FAIL")
                failures += 1
    return 1 if failures else 0


if __name__ == "__main__":
    sys.exit(main())