{
  "nbformat": 4,
  "nbformat_minor": 5,
  "metadata": {
    "kernelspec": {
      "display_name": "Python 3",
      "language": "python",
      "name": "python3"
    },
    "language_info": {
      "name": "python",
      "version": "3.12"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "# PandaPower-based Implementation\n",
        "\n",
        "400 V microgrid/distribution feeder · pandapower 3.2.1\n",
        "\n",
        "**学习目标 / Learning goals:** construct equipment tables, run balanced and three-phase AC power flow, check units and supply, and generate line-outage labels.\n",
        "\n",
        "The three-phase model uses sequence impedance and earth return; it does not reproduce the preceding chapter's explicit finite-neutral four-wire model. The limits below are teaching settings.\n"
      ],
      "id": "pp-lesson-00"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## 1. Install the pinned environment / 安装固定版本依赖\n",
        "\n",
        "Use a fresh Python 3.12 environment and run this cell once. Restart the kernel after installation if replacing packages that were already imported. Installation needs internet access.\n"
      ],
      "id": "pp-lesson-01"
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "%pip install pandapower==3.2.1 numpy==2.2.5 scipy==1.14.1 pandas==2.3.1 networkx==3.4.2 packaging==25.0 tqdm==4.67.1 deepdiff==8.6.1 geojson==3.2.0 typing_extensions==4.15.0 lxml==6.0.0\n"
      ],
      "execution_count": null,
      "outputs": [],
      "id": "pp-lesson-02"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## 2. Model source / 模型源码\n",
        "\n",
        "All helper functions are embedded here. `build_network` builds equipment tables, `run_network` solves, and `collect_results` checks the line-only feeder. Extend connectivity and balance checks when adding other equipment.\n"
      ],
      "id": "pp-lesson-03"
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "\"\"\"Real pandapower examples for the interactive teaching chapter (3.2.1).\n",
        "\n",
        "Run locally: python -m pip install -r requirements-pandapower.txt\n",
        "            python pandapower_implementation.py\n",
        "The three-phase model uses sequence impedances and earth return, not the\n",
        "explicit finite-neutral four-wire model of the preceding chapter.\n",
        "\"\"\"\n",
        "import json\n",
        "import math\n",
        "import pandapower as pp\n",
        "\n",
        "DEFAULT_CASE = dict(mode=\"balanced\", algorithm=\"nr\", load_scale=1.0,\n",
        "                    p2_kw=120.0, p3_a_kw=90.0, p3_b_kw=55.0, p3_c_kw=35.0,\n",
        "                    pf=0.95, dg_kw=50.0, dg_phase=\"balanced\", vm_pu=1.0,\n",
        "                    z_scale=1.0, zero_ratio=3.0, limit_a=600.0,\n",
        "                    tie=False, trip12=False, trip23=False, trip13=False)\n",
        "\n",
        "\n",
        "def settings(case=None):\n",
        "    c = dict(DEFAULT_CASE, **(case or {}))\n",
        "    if c[\"mode\"] not in (\"balanced\", \"three_phase\"):\n",
        "        raise ValueError(\"mode must be balanced or three_phase\")\n",
        "    if c[\"algorithm\"] not in (\"nr\", \"bfsw\"):\n",
        "        raise ValueError(\"algorithm must be nr or bfsw\")\n",
        "    if c[\"dg_phase\"] not in (\"balanced\", \"a\", \"b\", \"c\"):\n",
        "        raise ValueError(\"dg_phase must be balanced, a, b, or c\")\n",
        "    for key in (\"load_scale\", \"p2_kw\", \"p3_a_kw\", \"p3_b_kw\", \"p3_c_kw\", \"dg_kw\"):\n",
        "        if not math.isfinite(float(c[key])) or c[key] < 0:\n",
        "            raise ValueError(key + \" must be finite and nonnegative\")\n",
        "    for key in (\"z_scale\", \"zero_ratio\", \"limit_a\", \"vm_pu\"):\n",
        "        if not math.isfinite(float(c[key])) or c[key] <= 0:\n",
        "            raise ValueError(key + \" must be finite and positive\")\n",
        "    if not 0 < c[\"pf\"] <= 1:\n",
        "        raise ValueError(\"pf must be in (0, 1]\")\n",
        "    return c\n",
        "\n",
        "\n",
        "def build_network(case=None):\n",
        "    \"\"\"Build in physical units: kV, MW, Mvar, ohm/km, km, kA.\"\"\"\n",
        "    c = settings(case)\n",
        "    net = pp.create_empty_network(sn_mva=1.0, f_hz=60.0)\n",
        "    buses = [pp.create_bus(net, vn_kv=0.4, name=f\"Bus {i+1}\") for i in range(3)]\n",
        "    # Positive-sequence slack; 0/2-sequence source impedances follow these\n",
        "    # illustrative short-circuit parameters, not a physical cable spec.\n",
        "    pp.create_ext_grid(net, buses[0], vm_pu=c[\"vm_pu\"], va_degree=0,\n",
        "                       s_sc_max_mva=1000, s_sc_min_mva=1000,\n",
        "                       rx_max=0.1, rx_min=0.1, r0x0_max=0.1, x0x_max=1.0)\n",
        "    branches = [(0, 1, .012, .008, not c[\"trip12\"]),\n",
        "                (1, 2, .008, .006, not c[\"trip23\"]),\n",
        "                (0, 2, .022, .014, c[\"tie\"] and not c[\"trip13\"])]\n",
        "    for a, b, r, x, active in branches:\n",
        "        # Length = 1 km is an equivalent parameterization: length * r/km\n",
        "        # reproduces the previous balanced chapter's total branch resistance.\n",
        "        r, x = r * c[\"z_scale\"], x * c[\"z_scale\"]\n",
        "        pp.create_line_from_parameters(\n",
        "            net, buses[a], buses[b], length_km=1.0,\n",
        "            r_ohm_per_km=r, x_ohm_per_km=x, c_nf_per_km=0,\n",
        "            max_i_ka=c[\"limit_a\"] / 1000, name=f\"Line {a+1}-{b+1}\",\n",
        "            in_service=bool(active), r0_ohm_per_km=r*c[\"zero_ratio\"],\n",
        "            x0_ohm_per_km=x*c[\"zero_ratio\"], c0_nf_per_km=0)\n",
        "    q_ratio = math.tan(math.acos(c[\"pf\"]))\n",
        "    phases = [c[f\"p3_{phase}_kw\"] for phase in \"abc\"]\n",
        "    if c[\"mode\"] == \"balanced\":\n",
        "        for bus, kw in ((buses[1], c[\"p2_kw\"]), (buses[2], sum(phases))):\n",
        "            p = kw * c[\"load_scale\"] / 1000\n",
        "            pp.create_load(net, bus, p_mw=p, q_mvar=p*q_ratio)\n",
        "        # sgen is fixed PQ; it is not a voltage-controlled PV bus.\n",
        "        pp.create_sgen(net, buses[2], p_mw=c[\"dg_kw\"]/1000, q_mvar=0)\n",
        "    else:\n",
        "        for bus, powers in ((buses[1], [c[\"p2_kw\"]/3]*3), (buses[2], phases)):\n",
        "            kwargs = {}\n",
        "            for phase, kw in zip(\"abc\", powers):\n",
        "                p = kw * c[\"load_scale\"] / 1000\n",
        "                kwargs[f\"p_{phase}_mw\"] = p\n",
        "                kwargs[f\"q_{phase}_mvar\"] = p*q_ratio\n",
        "            pp.create_asymmetric_load(net, bus, type=\"wye\", **kwargs)\n",
        "        generation = [c[\"dg_kw\"]/3]*3 if c[\"dg_phase\"] == \"balanced\" else [\n",
        "            c[\"dg_kw\"] if phase == c[\"dg_phase\"] else 0 for phase in \"abc\"]\n",
        "        pp.create_asymmetric_sgen(net, buses[2], type=\"wye\",\n",
        "                                 **{f\"p_{p}_mw\": kw/1000 for p, kw in zip(\"abc\", generation)})\n",
        "    return net\n",
        "\n",
        "\n",
        "def run_network(net, case=None):\n",
        "    c = settings(case)\n",
        "    kwargs = dict(numba=False, init=\"flat\", max_iteration=100,\n",
        "                  tolerance_mva=1e-9, check_connectivity=True)\n",
        "    if c[\"mode\"] == \"balanced\":\n",
        "        pp.runpp(net, algorithm=c[\"algorithm\"], calculate_voltage_angles=True,\n",
        "                 voltage_depend_loads=False, **kwargs)\n",
        "    else:\n",
        "        pp.runpp_3ph(net, calculate_voltage_angles=True, **kwargs)\n",
        "    return net\n",
        "\n",
        "\n",
        "def supplied_buses(net):\n",
        "    \"\"\"Connectivity for this line-only feeder with ext_grid sources.\n",
        "\n",
        "    This deliberately small graph routine does not model transformer/switch\n",
        "    connectivity. Use pandapower.topology for richer network experiments.\n",
        "    \"\"\"\n",
        "    reached = {int(row.bus) for _, row in net.ext_grid.iterrows() if row.in_service}\n",
        "    while True:\n",
        "        previous = set(reached)\n",
        "        for _, line in net.line.iterrows():\n",
        "            if line.in_service and (int(line.from_bus) in reached or int(line.to_bus) in reached):\n",
        "                reached.update((int(line.from_bus), int(line.to_bus)))\n",
        "        if reached == previous:\n",
        "            return reached\n",
        "\n",
        "\n",
        "def finite(value):\n",
        "    return float(value) if math.isfinite(float(value)) else None\n",
        "\n",
        "\n",
        "def collect_results(net, case=None):\n",
        "    \"\"\"Read actual result DataFrames; convert NaN to JSON null, never zero.\"\"\"\n",
        "    c = settings(case)\n",
        "    three = c[\"mode\"] == \"three_phase\"\n",
        "    bus_table = net.res_bus_3ph if three else net.res_bus\n",
        "    line_table = net.res_line_3ph if three else net.res_line\n",
        "    source_table = net.res_ext_grid_3ph if three else net.res_ext_grid\n",
        "    reached = supplied_buses(net)\n",
        "    buses, lines = [], []\n",
        "    for index, row in bus_table.iterrows():\n",
        "        voltages = [finite(row[f\"vm_{p}_pu\"]) for p in \"abc\"] if three else [finite(row.vm_pu)]\n",
        "        angles = [finite(row[f\"va_{p}_degree\"]) for p in \"abc\"] if three else [finite(row.va_degree)]\n",
        "        # An unsupplied bus has no operating voltage; some 3ph result paths\n",
        "        # return zero instead of NaN. Preserve a missing value in the lesson.\n",
        "        if int(index) not in reached:\n",
        "            voltages, angles = [None]*len(voltages), [None]*len(angles)\n",
        "        buses.append(dict(id=int(index), name=str(net.bus.at[index, \"name\"]),\n",
        "                          supplied=int(index) in reached, vm_pu=voltages,\n",
        "                          va_degree=angles, vuf_pct=(finite(row.unbalance_percent) if three else 0.0) if int(index) in reached else None))\n",
        "    for index, row in line_table.iterrows():\n",
        "        currents = [finite(1000*max(row[f\"i_{p}_from_ka\"], row[f\"i_{p}_to_ka\"])) for p in \"abc\"] if three else [finite(1000*row.i_ka)]\n",
        "        loss = sum(row[f\"pl_{p}_mw\"] for p in \"abc\") if three else row.pl_mw\n",
        "        power = sum(row[f\"p_{p}_from_mw\"] for p in \"abc\") if three else row.p_from_mw\n",
        "        lines.append(dict(id=int(index), name=str(net.line.at[index, \"name\"]),\n",
        "                          from_bus=int(net.line.at[index, \"from_bus\"]),\n",
        "                          to_bus=int(net.line.at[index, \"to_bus\"]),\n",
        "                          in_service=bool(net.line.at[index, \"in_service\"]),\n",
        "                          loading_pct=finite(row.loading_percent), current_a=currents,\n",
        "                          loss_kw=finite(loss*1000), p_from_kw=finite(power*1000)))\n",
        "    source_columns = [f\"p_{p}_mw\" for p in \"abc\"] if three else [\"p_mw\"]\n",
        "    slack_kw = float(source_table[source_columns].sum().sum()*1000)\n",
        "    loss_kw = sum(line[\"loss_kw\"] or 0 for line in lines)\n",
        "    unsupplied = [b[\"id\"] for b in buses if not b[\"supplied\"]]\n",
        "    voltage_values = [v for b in buses for v in b[\"vm_pu\"] if v is not None]\n",
        "    max_loading = max((line[\"loading_pct\"] or 0 for line in lines), default=0)\n",
        "    max_vuf = max((b[\"vuf_pct\"] or 0 for b in buses), default=0)\n",
        "    violations = []\n",
        "    if any(v < .95 or v > 1.05 for v in voltage_values):\n",
        "        violations.append(\"voltage\")\n",
        "    if max_loading > 100:\n",
        "        violations.append(\"current\")\n",
        "    if three and max_vuf > 2:\n",
        "        violations.append(\"unbalance\")\n",
        "    if unsupplied:\n",
        "        violations.append(\"unsupplied\")\n",
        "    # Balance check uses only served load/generation result tables.\n",
        "    demand_kw, generation_kw = 0.0, 0.0\n",
        "    for name, is_generation in ((\"load\", False), (\"sgen\", True),\n",
        "                                (\"asymmetric_load\", False), (\"asymmetric_sgen\", True)):\n",
        "        result_name = \"res_\"+name+(\"_3ph\" if three else \"\")\n",
        "        table = net[result_name]\n",
        "        cols = [f\"p_{p}_mw\" for p in \"abc\"] if three else [\"p_mw\"]\n",
        "        if not table.empty:\n",
        "            power_kw = float(table[cols].sum().sum()*1000)\n",
        "            if is_generation:\n",
        "                generation_kw += power_kw\n",
        "            else:\n",
        "                demand_kw += power_kw\n",
        "    return dict(ok=bool(net.converged), mode=c[\"mode\"], version=pp.__version__,\n",
        "                buses=buses, lines=lines, unsupplied=unsupplied,\n",
        "                slack_kw=slack_kw, loss_kw=loss_kw, min_voltage=min(voltage_values, default=None),\n",
        "                max_loading=max_loading, max_vuf=max_vuf, violations=violations,\n",
        "                secure=bool(net.converged) and not violations,\n",
        "                balance_error_kw=slack_kw+generation_kw-demand_kw-loss_kw)\n",
        "\n",
        "\n",
        "def solve(case=None):\n",
        "    c = settings(case)\n",
        "    net = build_network(c)\n",
        "    if len(supplied_buses(net)) == 1:\n",
        "        # runpp_3ph cannot solve an empty set of supplied PQ buses in 3.2.1.\n",
        "        return dict(ok=False, mode=c[\"mode\"], reason=\"island\",\n",
        "                    message=\"No load bus is connected to the source.\", unsupplied=[1, 2])\n",
        "    try:\n",
        "        run_network(net, c)\n",
        "    except pp.LoadflowNotConverged as error:\n",
        "        return dict(ok=False, mode=c[\"mode\"], reason=\"nonconvergence\", message=str(error),\n",
        "                    unsupplied=sorted(set(map(int, net.bus.index))-supplied_buses(net)))\n",
        "    return collect_results(net, c)\n",
        "\n",
        "\n"
      ],
      "execution_count": null,
      "outputs": [],
      "id": "pp-lesson-04"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## 3. Balanced baseline / 平衡基准\n",
        "\n",
        "Bus 2: 120 kW; bus 3: 180 kW; PV: 50 kW; PF: 0.95. Voltage base: 0.4 kV line-to-line. Expect bus 3 ≈ 0.966565 pu and losses ≈ 6.838895 kW.\n"
      ],
      "id": "pp-lesson-05"
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "case = dict(DEFAULT_CASE)\n",
        "net = build_network(case)\n",
        "run_network(net, case)\n",
        "display(net.bus, net.line, net.load, net.sgen)\n",
        "display(net.res_bus, net.res_line)\n",
        "balanced = collect_results(net, case)\n",
        "print(\"Loss (kW):\", balanced[\"loss_kw\"])\n",
        "print(\"Power-balance residual (kW):\", balanced[\"balance_error_kw\"])\n"
      ],
      "execution_count": null,
      "outputs": [],
      "id": "pp-lesson-06"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## 4. Same total demand, unequal phases / 同总量下的不平衡\n",
        "\n",
        "Change bus 3 to 90/55/35 kW with equal PV. Compare 60/60/60 kW without changing total demand. The source sequence parameters and line zero-sequence impedances are illustrative.\n"
      ],
      "id": "pp-lesson-07"
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "import pandas as pd\n",
        "comparison = []\n",
        "for name, phases in ((\"unequal\", (90,55,35)), (\"equal\", (60,60,60))):\n",
        "    trial = dict(case, mode=\"three_phase\", **dict(zip((\"p3_a_kw\",\"p3_b_kw\",\"p3_c_kw\"), phases)))\n",
        "    solved = solve(trial)\n",
        "    comparison.append({\"case\": name, \"bus3_voltage\": solved[\"buses\"][2][\"vm_pu\"],\n",
        "                       \"loss_kw\": solved[\"loss_kw\"], \"vuf_pct\": solved[\"max_vuf\"], \"secure\": solved[\"secure\"]})\n",
        "display(pd.DataFrame(comparison))\n"
      ],
      "execution_count": null,
      "outputs": [],
      "id": "pp-lesson-08"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## 5. Edit and rerun / 修改后重新求解\n",
        "\n",
        "Solar PV here is fixed PQ, not a voltage-controlled PV bus. Changing a current rating changes loading, not an unconstrained power-flow operating point.\n"
      ],
      "id": "pp-lesson-09"
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "net.line[\"max_i_ka\"] = 0.3\n",
        "run_network(net, case)\n",
        "display(net.res_line)\n",
        "net.load.loc[net.load.bus == 2, \"p_mw\"] = .240\n",
        "net.load.loc[net.load.bus == 2, \"q_mvar\"] = .240 * math.tan(math.acos(.95))\n",
        "run_network(net, case)\n",
        "display(net.res_bus)\n"
      ],
      "execution_count": null,
      "outputs": [],
      "id": "pp-lesson-10"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## 6. N-1 screen / 线路 N-1 筛查\n",
        "\n",
        "Each trial starts from the same base. Supply and limits must be checked in addition to convergence. A fixed-PQ inverter cannot supply an isolated network without a modeled grid-forming source.\n"
      ],
      "id": "pp-lesson-11"
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "base = dict(case, tie=True, trip12=False, trip23=False, trip13=False)\n",
        "rows = []\n",
        "for outage in (None, \"trip12\", \"trip23\", \"trip13\"):\n",
        "    trial = dict(base)\n",
        "    if outage:\n",
        "        trial[outage] = True\n",
        "    solved = solve(trial)\n",
        "    rows.append({\"outage\": outage or \"base\", \"converged\": solved[\"ok\"],\n",
        "                 \"unsupplied\": solved.get(\"unsupplied\", []), \"min_v_pu\": solved.get(\"min_voltage\"),\n",
        "                 \"max_loading_pct\": solved.get(\"max_loading\"), \"secure\": solved.get(\"secure\", False)})\n",
        "display(pd.DataFrame(rows))\n"
      ],
      "execution_count": null,
      "outputs": [],
      "id": "pp-lesson-12"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## 7. Exercises / 练习\n",
        "\n",
        "- Compare `nr` and `bfsw` on the balanced feeder.\n",
        "- Try phase-A-only PV in three-phase mode.\n",
        "- Repeat outages with the tie open and record unsupplied buses separately from nonconvergence.\n",
        "- Export scenario features and labels while preserving bus IDs and phase order.\n",
        "\n",
        "Sources: [pandapower units](https://pandapower.readthedocs.io/en/v3.2.1/about/units.html), [three-phase assumptions](https://pandapower.readthedocs.io/en/v3.2.1/powerflow/ac_3ph.html).\n"
      ],
      "id": "pp-lesson-13"
    }
  ]
}