notebooks/analysis_tsg_nonlinear.ipynb

{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "fdae8563",
   "metadata": {},
   "source": [
    "## Analysis of Controller on Nonlinear Simulation\n",
    "\n",
    "This notebook generates the following tables and figures in the paper\n",
    "\n",
    "> Yeh, Christopher, et al. \"Online learning for robust voltage control under uncertain grid topology.\" _arXiv preprint arXiv:2306.16674_ (2023).\n",
    "\n",
    "- Tables: 1 (bottom)\n",
    "- Figures: 2b, 4, 5b, 6"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%load_ext autoreload\n",
    "%autoreload 2\n",
    "%matplotlib inline\n",
    "%cd .."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from __future__ import annotations\n",
    "\n",
    "from glob import glob\n",
    "import os\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "from tqdm.auto import tqdm\n",
    "\n",
    "from analysis import calculate_violations, plot_error_and_etahat, plot_voltages, plot_bus\n",
    "from utils import load_pkl, savefig\n",
    "\n",
    "# hide top and right splines on plots\n",
    "plt.rcParams['axes.spines.right'] = False\n",
    "plt.rcParams['axes.spines.top'] = False"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "02043ebf",
   "metadata": {},
   "outputs": [],
   "source": [
    "plots_dir = 'plots/tsg_final/'\n",
    "os.makedirs(plots_dir, exist_ok=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6d6fe05e",
   "metadata": {},
   "outputs": [],
   "source": [
    "from network_utils import (\n",
    "    create_56bus,\n",
    "    create_RX_from_net,\n",
    "    read_load_data)\n",
    "\n",
    "net = create_56bus()\n",
    "R, X = create_RX_from_net(net, noise=0)  # true R and X\n",
    "p, qe = read_load_data()  # in MW and MVar\n",
    "T, n = p.shape\n",
    "\n",
    "v_nom = 12**2  # nominal squared voltage magnitude, units kV^2\n",
    "v_sub = v_nom  # fixed squared voltage magnitude at substation, units kV^2\n",
    "\n",
    "# saved nonlinear vpars includes substation. Here, we ignore the substation.\n",
    "vpars = np.load('data/nonlinear_voltage_baseline.npy')[:, 1:]  # shape [T, n]\n",
    "vpars = (vpars * 12.)**2\n",
    "Vpar_min = np.min(vpars, axis=0) - 0.5\n",
    "Vpar_max = np.max(vpars, axis=0) + 0.5\n",
    "Vpar = (Vpar_min, Vpar_max)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6716963c",
   "metadata": {},
   "outputs": [],
   "source": [
    "print('max-voltage node:', np.argmax(vpars.max(axis=0)))\n",
    "print('min-voltage node:', np.argmin(vpars.min(axis=0)))"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "74b2fdce",
   "metadata": {},
   "source": [
    "## Figure 2b\n",
    "\n",
    "nonlinear simulation, no controller"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5aa06737",
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_voltages(vpars, ylim=(11.0, 13.4), yticks=None,\n",
    "              plots_dir=plots_dir, filename='nonlinear_nocontrol')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b467f8ab",
   "metadata": {},
   "source": [
    "## Table 1, bottom\n",
    "\n",
    "nonlinear simulation, mistakes and violations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e809c894",
   "metadata": {},
   "outputs": [],
   "source": [
    "outdir = 'out/nonlinear/'\n",
    "pkl_paths = {\n",
    "    ('known', None): outdir + 'CBCconst_δ20_η10_20230810_130842.pkl',  # fixed X̂, learned etahat\n",
    "}\n",
    "for seed in [8, 9, 10, 11]:\n",
    "    pkl_paths |= {\n",
    "        # default: δ=20\n",
    "        ('unknown', seed): glob(outdir + f'CBCproj_δ20_η10_noise1.0_perm_norm1.0_seed{seed}_2*.pkl')[0],\n",
    "        ('topo-14', seed): glob(outdir + f'CBCproj_δ20_η10_noise1.0_perm_norm1.0_seed{seed}_knowntopo14_2*.pkl')[0],\n",
    "        ('lines-14', seed): glob(outdir + f'CBCproj_δ20_η10_noise1.0_perm_norm1.0_seed{seed}_knownlines14_2*.pkl')[0],\n",
    "    }\n",
    "\n",
    "pkls = {}\n",
    "for (name, seed), pkl_path in pkl_paths.items():\n",
    "    pkl = load_pkl(pkl_path)\n",
    "    pkls[(name, seed)] = pkl\n",
    "    print(f'{name: <15} {str(seed): <4}', pkl.keys())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d1739d97",
   "metadata": {},
   "outputs": [],
   "source": [
    "figs_and_axs = {}\n",
    "for seed in [8, 9, 10, 11]:\n",
    "    figs_and_axs[seed] = plt.subplots(1, 3, figsize=(9, 3))\n",
    "\n",
    "rows = []\n",
    "for key, pkl in pkls.items():\n",
    "    info, seed = key\n",
    "    ax = None\n",
    "    if seed in figs_and_axs:\n",
    "        fig, axs = figs_and_axs[seed]\n",
    "        if key[0] == 'unknown': ax = axs[0]\n",
    "        elif key[0] == 'topo-14': ax = axs[1]\n",
    "        elif key[0] == 'lines-14': ax = axs[2]\n",
    "\n",
    "    num_mistakes, avg_viol, max_viol = calculate_violations(\n",
    "        key=key, pkl=pkl, ax=ax, T=T, n=n)\n",
    "    rows.append((info, seed, num_mistakes, avg_viol, max_viol))\n",
    "\n",
    "# df.groupby().agg('std') gives sample standard deviation (ddof=1), which is what we want\n",
    "df = pd.DataFrame.from_records(rows, columns=['info', 'seed', 'mistakes', 'avg_viol', 'max_viol'])\n",
    "stats = df.groupby(['info'])[['mistakes', 'avg_viol', 'max_viol']].agg(['mean', 'std'])\n",
    "\n",
    "with pd.option_context('display.precision', 1):\n",
    "    display(stats)\n",
    "with pd.option_context('display.precision', 2):\n",
    "    display(stats)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7b5312b1",
   "metadata": {},
   "source": [
    "## Figure 4\n",
    "\n",
    "nonlinear simulation, voltage curves and model error"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cb99fb65",
   "metadata": {},
   "outputs": [],
   "source": [
    "outdir = 'out/nonlinear/'\n",
    "pkl_paths = {\n",
    "    ('known', None): outdir + 'CBCconst_δ20_η10_20230810_130842.pkl',  # fixed X̂, learned etahat\n",
    "}\n",
    "seeds = [8]  # [8, 9, 10, 11]\n",
    "for seed in seeds:\n",
    "    pkl_paths |= {\n",
    "        # default: δ=20\n",
    "        ('unknown', seed): glob(outdir + f'CBCproj_δ20_η10_noise1.0_perm_norm1.0_seed{seed}_2*.pkl')[0],\n",
    "        ('topo-14', seed): glob(outdir + f'CBCproj_δ20_η10_noise1.0_perm_norm1.0_seed{seed}_knowntopo14_2*.pkl')[0],\n",
    "        ('lines-14', seed): glob(outdir + f'CBCproj_δ20_η10_noise1.0_perm_norm1.0_seed{seed}_knownlines14_2*.pkl')[0],\n",
    "    }\n",
    "\n",
    "pkls = {}\n",
    "for (name, seed), pkl_path in pkl_paths.items():\n",
    "    pkl = load_pkl(pkl_path)\n",
    "    pkls[(name, seed)] = pkl\n",
    "    print(f'{name: <15} {str(seed): <4}', pkl.keys())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a389d265",
   "metadata": {},
   "outputs": [],
   "source": [
    "def fig4abcd(seed: int) -> None:\n",
    "    for name in ['unknown', 'topo-14', 'lines-14']:\n",
    "        data = pkls[(name, seed)]\n",
    "        filename = f'nonlinear_{name}_s{seed}'\n",
    "        plot_voltages(data['vs'], plots_dir=plots_dir, filename=filename)\n",
    "    \n",
    "    data = pkls[('known', None)]\n",
    "    filename = f'nonlinear_known'\n",
    "    plot_voltages(data['vs'], plots_dir=plots_dir, filename=filename)\n",
    "\n",
    "\n",
    "fig4abcd(seed=8)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "480044b1",
   "metadata": {},
   "outputs": [],
   "source": [
    "def fig4e(seed: int) -> None:\n",
    "    fig4e_pkls = {\n",
    "        name: pkls[(name, seed)]\n",
    "        for name in ['unknown', 'topo-14', 'lines-14']\n",
    "    }\n",
    "    fig4e_pkls['known'] = pkls[('known', None)]\n",
    "    plot_error_and_etahat(\n",
    "        fig4e_pkls, plots_dir=plots_dir, filename=f'nonlinear_error_s{seed}',\n",
    "        legend_loc='top')\n",
    "\n",
    "\n",
    "fig4e(seed=8)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5f8c1258",
   "metadata": {},
   "source": [
    "## Figure 5b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "73c0e4ff",
   "metadata": {},
   "outputs": [],
   "source": [
    "outdir = 'out/nonlinear/'\n",
    "pkl_paths = {}\n",
    "seeds = [8]  # [8, 9, 10, 11]\n",
    "for seed in seeds:\n",
    "    pkl_paths = {\n",
    "        (r'η* known', seed): glob(outdir + f'CBCproj_noise1.0_perm_norm1.0_seed{seed}_knownlines14_2*.pkl')[0],\n",
    "        (r'δ=1', seed): glob(outdir + f'CBCproj_δ1_η10_noise1.0_perm_norm1.0_seed{seed}_knownlines14_2*.pkl')[0],\n",
    "        (r'δ=20', seed): glob(outdir + f'CBCproj_δ20_η10_noise1.0_perm_norm1.0_seed{seed}_knownlines14_2*.pkl')[0],\n",
    "        (r'δ=100', seed): glob(outdir + f'CBCproj_δ100_η10_noise1.0_perm_norm1.0_seed{seed}_knownlines14_2*.pkl')[0],\n",
    "        (r'δ=500', seed): glob(outdir + f'CBCproj_δ500_η10_noise1.0_perm_norm1.0_seed{seed}_knownlines14_2*.pkl')[0],\n",
    "    }\n",
    "\n",
    "pkls = {}\n",
    "for (name, seed), pkl_path in pkl_paths.items():\n",
    "    pkl = load_pkl(pkl_path)\n",
    "    pkls[(name, seed)] = pkl\n",
    "    print(f'{name: <15} {str(seed): <4}', pkl.keys())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ac36054f",
   "metadata": {},
   "outputs": [],
   "source": [
    "def fig5b(seed: int, legend_loc: str | None) -> None:\n",
    "    pkls_by_delta = {\n",
    "        name: pkls[(name, seed)]\n",
    "        for name in ['η* known', 'δ=1', 'δ=20', 'δ=100', 'δ=500']\n",
    "    }\n",
    "    filename = f'nonlinear_error_by_delta_s{seed}'\n",
    "    plot_error_and_etahat(\n",
    "        pkls_by_delta, plots_dir=plots_dir, filename=filename,\n",
    "        legend_loc=legend_loc, etamax=10)\n",
    "\n",
    "\n",
    "fig5b(seed=8, legend_loc=None)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "a2149f35",
   "metadata": {},
   "source": [
    "## Figure 6\n",
    "\n",
    "nonlinear simulation, partial control"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1be04b6f",
   "metadata": {},
   "outputs": [],
   "source": [
    "pkl_paths = {\n",
    "    ('known', None): outdir + 'CBCconst_δ20_η10_partialctrl_20230811_194616.pkl',  # fixed X̂, learned etahat\n",
    "}\n",
    "for seed in [8, 9, 10, 11]:\n",
    "    pkl_paths |= {\n",
    "        # default: δ=20\n",
    "        ('unknown', seed): glob(outdir + f'CBCproj_δ20_η10_noise1.0_perm_norm1.0_seed{seed}_partialctrl_2*.pkl')[0],\n",
    "        ('lines-14', seed): glob(outdir + f'CBCproj_δ20_η10_noise1.0_perm_norm1.0_seed{seed}_partialctrl_knownlines14_2*.pkl')[0],\n",
    "    }\n",
    "\n",
    "pkls = {}\n",
    "for (name, seed), pkl_path in pkl_paths.items():\n",
    "    pkl = load_pkl(pkl_path)\n",
    "    pkls[(name, seed)] = pkl\n",
    "    print(f'{name: <15} {str(seed): <4}', pkl.keys())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c770b31e",
   "metadata": {},
   "outputs": [],
   "source": [
    "def fig6(seeds: list[int]) -> None:\n",
    "    for i, bus in enumerate([18, 30]):\n",
    "        pkls_by_label = {\n",
    "            'unknown': [pkls[('unknown', seed)] for seed in seeds],\n",
    "            'lines-14': [pkls[('lines-14', seed)] for seed in seeds],\n",
    "            'known': [pkls[('known', None)]],\n",
    "        }\n",
    "        filename_base = 'nonlinear_partialctrl'\n",
    "        plot_bus(pkls_by_label, bus=bus, plots_dir=plots_dir,\n",
    "                 legend=(i==0), filename_base=filename_base)\n",
    "\n",
    "fig6(seeds=[8, 9, 10, 11])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1528ed03",
   "metadata": {},
   "source": [
    "## Unused"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f00816ec",
   "metadata": {},
   "outputs": [],
   "source": [
    "outdir = 'out/nonlinear/'\n",
    "pkl_paths = {\n",
    "    # ('known', None): outdir + 'CBCconst_20230810_130611.pkl',  # fixed X̂, fixed etahat\n",
    "    ('known', None): outdir + 'CBCconst_δ20_η10_20230810_130842.pkl',  # fixed X̂, learned etahat\n",
    "}\n",
    "for seed in [8, 9, 10, 11]:\n",
    "    pkl_paths |= {\n",
    "        # default: δ=20\n",
    "        ('unknown', seed): glob(outdir + f'CBCproj_δ20_η10_noise1.0_perm_norm1.0_seed{seed}_2*.pkl')[0],\n",
    "        ('topo-14', seed): glob(outdir + f'CBCproj_δ20_η10_noise1.0_perm_norm1.0_seed{seed}_knowntopo14_2*.pkl')[0],\n",
    "        ('lines-14', seed): glob(outdir + f'CBCproj_δ20_η10_noise1.0_perm_norm1.0_seed{seed}_knownlines14_2*.pkl')[0],\n",
    "\n",
    "        # (r'η* known', seed): glob(outdir + f'CBCproj_noise1.0_perm_norm1.0_seed{seed}_knownlines14_2*.pkl')[0],\n",
    "        # (r'δ=1', seed): glob(outdir + f'CBCproj_δ1_η10_noise1.0_perm_norm1.0_seed{seed}_knownlines14_2*.pkl')[0],\n",
    "        # (r'δ=20', seed): glob(outdir + f'CBCproj_δ20_η10_noise1.0_perm_norm1.0_seed{seed}_knownlines14_2*.pkl')[0],\n",
    "        # (r'δ=100', seed): glob(outdir + f'CBCproj_δ100_η10_noise1.0_perm_norm1.0_seed{seed}_knownlines14_2*.pkl')[0],\n",
    "        # (r'δ=500', seed): glob(outdir + f'CBCproj_δ500_η10_noise1.0_perm_norm1.0_seed{seed}_knownlines14_2*.pkl')[0],\n",
    "    }\n",
    "\n",
    "pkls = {}\n",
    "for (name, seed), pkl_path in pkl_paths.items():\n",
    "    pkl = load_pkl(pkl_path)\n",
    "    pkls[(name, seed)] = pkl\n",
    "    print(f'{name: <15} {str(seed): <4}', pkl.keys())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9ca8168b",
   "metadata": {},
   "outputs": [],
   "source": [
    "def check_consistency(data: dict):\n",
    "    v = data['vs']\n",
    "    qc = data['qcs']\n",
    "    u = qc[1:] - qc[:-1]  # u[t] = u(t) = q^c(t+1) - q^c(t)\n",
    "    Δv = v[1:] - v[:-1]  # Δv[t] = v(t+1) - v(t)\n",
    "    count_inconsistent = []\n",
    "    for t in sorted(data['params'].keys()):\n",
    "        if data['config']['δ'] == 0:\n",
    "            X̂ = data['params'][t]\n",
    "            etahat = 8.65\n",
    "        else:\n",
    "            X̂, etahat = data['params'][t]\n",
    "        vpar_hat = v[1:t+1] - qc[1:t+1] @ X̂\n",
    "        w_hat = Δv[:t] - u[:t] @ X̂\n",
    "        consistent = (\n",
    "            (Vpar_min - 0.05 <= vpar_hat).all(axis=1)\n",
    "            | (vpar_hat <= Vpar_max + 0.05).all(axis=1)\n",
    "            | (np.max(np.abs(w_hat), axis=1) <= etahat)\n",
    "        )\n",
    "        num_inconsistent = t - consistent.sum()\n",
    "        count_inconsistent.append(num_inconsistent)\n",
    "    return count_inconsistent"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "affed81e",
   "metadata": {},
   "outputs": [],
   "source": [
    "count_inconsistent = {}\n",
    "for (name, seed), data in tqdm(pkls.items()):\n",
    "    count_inconsistent[(name, seed)] = check_consistency(data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "53bfb2d8",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axs = plt.subplots(1, 4, sharey=True, figsize=(15, 4))\n",
    "for (name, seed), counts in count_inconsistent.items():\n",
    "    if seed is None:\n",
    "        continue\n",
    "    ax = axs[seed % 8]\n",
    "    ax.plot(np.arange(len(counts)) * 100, counts, label=f'{name}')\n",
    "    ax.legend(loc='upper left')\n",
    "for ax, seed in zip(axs, [8, 9, 10, 11]):\n",
    "    ax.set(xlabel='time $t$', title=f'seed {seed}')\n",
    "    if seed == 8:\n",
    "        ax.set(ylabel='# of inconsistent data points')\n",
    "\n",
    "# savefig(fig, plots_dir=plots_dir, filename='violations.png')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "cms21",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.14"
  },
  "vscode": {
   "interpreter": {
    "hash": "10b6275cc9ddd3a78a2191f6baae49e6f52bcb59b992509c77c5b1f4a0eb0ef2"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}