pypower_test_case9.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import pypower.api as pp"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## PF benchmark"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Manual solution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pypower.api import ext2int, bustypes, makeBdc, makeSbus, dcpf, int2ext, case9\n",
    "\n",
    "from numpy import r_, c_, ix_, zeros, pi, ones, exp, argmax, union1d\n",
    "from numpy import flatnonzero as find\n",
    "\n",
    "from pypower.idx_bus import PD, QD, VM, VA, GS, BUS_TYPE, PV, PQ, REF\n",
    "from pypower.idx_brch import PF, PT, QF, QT\n",
    "from pypower.idx_gen import PG, QG, VG, QMAX, QMIN, GEN_BUS, GEN_STATUS\n",
    "\n",
    "# Read PyPower case\n",
    "ppc = case9()\n",
    "\n",
    "## add zero columns to branch for flows if needed \n",
    "if ppc[\"branch\"].shape[1] < QT: \n",
    "    ppc[\"branch\"] = c_[ppc[\"branch\"], \n",
    "                    zeros((ppc[\"branch\"].shape[0], \n",
    "                            QT - ppc[\"branch\"].shape[1] + 1))] \n",
    "\n",
    "## convert to internal indexing\n",
    "ppc = ext2int(ppc)\n",
    "baseMVA, bus, gen, branch = \\\n",
    "    ppc[\"baseMVA\"], ppc[\"bus\"], ppc[\"gen\"], ppc[\"branch\"]\n",
    "\n",
    "## get bus index lists of each type of bus\n",
    "ref, pv, pq = bustypes(bus, gen)\n",
    "\n",
    "## generator info\n",
    "on = find(gen[:, GEN_STATUS] > 0)      ## which generators are on?\n",
    "gbus = gen[on, GEN_BUS].astype(int)    ## what buses are they at?\n",
    "\n",
    "##-----  run the power flow  -----\n",
    "## initial state\n",
    "Va0 = bus[:, VA] * (pi / 180)\n",
    "\n",
    "## build B matrices and phase shift injections\n",
    "B, Bf, _, _ = makeBdc(baseMVA, bus, branch)\n",
    "\n",
    "## compute complex bus power injections [generation - load]\n",
    "Pbus = makeSbus(baseMVA, bus, gen).real\n",
    "\n",
    "## \"run\" the power flow\n",
    "pvpq = np.matrix(r_[pv, pq]) \n",
    "\n",
    "# initialize result vector \n",
    "# Va = np.copy(Va0)\n",
    "Va = np.zeros_like(Va0)\n",
    "\n",
    "B_dense = B.todense()\n",
    "B_bar = B_dense[pvpq.T, pvpq]\n",
    "P_bar = Pbus[pvpq]\n",
    "\n",
    "## update angles for non-reference buses\n",
    "Va[pvpq] = np.dot(np.linalg.inv(B_bar), P_bar.T).reshape(1,8)\n",
    " \n",
    "# update data matrices with solution\n",
    "branch[:, [QF, QT]] = zeros((branch.shape[0], 2))\n",
    "branch[:, PF] = (Bf * Va) * baseMVA\n",
    "branch[:, PT] = -branch[:, PF]\n",
    "bus[:, VM] = ones(bus.shape[0])\n",
    "bus[:, VA] = Va * (180 / pi)\n",
    "## update Pg for slack generator (1st gen at ref bus)\n",
    "## (note: other gens at ref bus are accounted for in Pbus)\n",
    "##      Pg = Pinj + Pload\n",
    "##      newPg = oldPg + newPinj - oldPinj\n",
    "refgen = find(gbus == ref)\n",
    "gen[on[refgen[0]], PG] = gen[on[refgen[0]], PG] + (B[ref, :] * Va - Pbus[ref])* baseMVA\n",
    "\n",
    "\n",
    "##-----  output results  -----\n",
    "## convert back to original bus numbering & print results\n",
    "ppc[\"bus\"], ppc[\"gen\"], ppc[\"branch\"] = bus, gen, branch\n",
    "results = int2ext(ppc)\n",
    "r_bus = results[\"bus\"]\n",
    "r_branch = results[\"branch\"]\n",
    "r_gen = results[\"gen\"]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### PyPower solution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "PYPOWER Version 5.1.4, 27-June-2018 -- DC Power Flow\n",
      "\n",
      "Converged in 0.00 seconds\n",
      "================================================================================\n",
      "|     System Summary                                                           |\n",
      "================================================================================\n",
      "\n",
      "How many?                How much?              P (MW)            Q (MVAr)\n",
      "---------------------    -------------------  -------------  -----------------\n",
      "Buses              9     Total Gen Capacity     820.0           0.0 to 0.0\n",
      "Generators         3     On-line Capacity       820.0           0.0 to 0.0\n",
      "Committed Gens     3     Generation (actual)    314.0               0.0\n",
      "Loads              3     Load                   315.0               0.0\n",
      "  Fixed            3       Fixed                315.0               0.0\n",
      "  Dispatchable     0       Dispatchable           0.0 of 0.0        0.0\n",
      "Shunts             0     Shunt (inj)              0.0               0.0\n",
      "Branches           9     Losses (I^2 * Z)         0.00              0.00\n",
      "Transformers       0     Branch Charging (inj)     -                0.0\n",
      "Inter-ties         0     Total Inter-tie Flow     0.0               0.0\n",
      "Areas              1\n",
      "\n",
      "                          Minimum                      Maximum\n",
      "                 -------------------------  --------------------------------\n",
      "Voltage Magnitude   1.000 p.u. @ bus 1          1.000 p.u. @ bus 1   \n",
      "Voltage Angle      -4.06 deg   @ bus 9          9.80 deg   @ bus 2   \n",
      "\n",
      "================================================================================\n",
      "|     Bus Data                                                                 |\n",
      "================================================================================\n",
      " Bus      Voltage          Generation             Load        \n",
      "  #   Mag(pu) Ang(deg)   P (MW)   Q (MVAr)   P (MW)   Q (MVAr)\n",
      "----- ------- --------  --------  --------  --------  --------\n",
      "    1  1.000    0.000*    66.00      0.00       -         -   \n",
      "    2  1.000    9.796    163.00      0.00       -         -   \n",
      "    3  1.000    5.061     85.00      0.00       -         -   \n",
      "    4  1.000   -2.211       -         -         -         -   \n",
      "    5  1.000   -3.738       -         -       90.00      0.00 \n",
      "    6  1.000    2.207       -         -         -         -   \n",
      "    7  1.000    0.822       -         -      100.00      0.00 \n",
      "    8  1.000    3.959       -         -         -         -   \n",
      "    9  1.000   -4.063       -         -      125.00      0.00 \n",
      "                        --------  --------  --------  --------\n",
      "               Total:    314.00      0.00    315.00      0.00\n",
      "\n",
      "================================================================================\n",
      "|     Branch Data                                                              |\n",
      "================================================================================\n",
      "Brnch   From   To    From Bus Injection   To Bus Injection     Loss (I^2 * Z)  \n",
      "  #     Bus    Bus    P (MW)   Q (MVAr)   P (MW)   Q (MVAr)   P (MW)   Q (MVAr)\n",
      "-----  -----  -----  --------  --------  --------  --------  --------  --------\n",
      "   0      1      4     67.00      0.00    -67.00      0.00     0.000      0.00\n",
      "   1      4      5     28.97      0.00    -28.97      0.00     0.000      0.00\n",
      "   2      5      6    -61.03      0.00     61.03      0.00     0.000      0.00\n",
      "   3      3      6     85.00      0.00    -85.00      0.00     0.000      0.00\n",
      "   4      6      7     23.97      0.00    -23.97      0.00     0.000      0.00\n",
      "   5      7      8    -76.03      0.00     76.03      0.00     0.000      0.00\n",
      "   6      8      2   -163.00      0.00    163.00      0.00     0.000      0.00\n",
      "   7      8      9     86.97      0.00    -86.97      0.00     0.000      0.00\n",
      "   8      9      4    -38.03      0.00     38.03      0.00     0.000      0.00\n",
      "                                                             --------  --------\n",
      "                                                    Total:     0.000      0.00\n"
     ]
    }
   ],
   "source": [
    "ppc = case9()\n",
    "pp_dcpf, success = pp.rundcpf(ppc)\n",
    "pp_bus = pp_dcpf[\"bus\"]\n",
    "pp_branch = pp_dcpf[\"branch\"]\n",
    "pp_gen = pp_dcpf[\"gen\"]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Comparison"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Bus voltage magitude:\n",
      "Manual:  [1. 1. 1. 1. 1. 1. 1. 1. 1.]\n",
      "Pypower:  [1. 1. 1. 1. 1. 1. 1. 1. 1.]\n",
      "Bus voltage angle:\n",
      "Manual:  [ 0.          9.79601886  5.06056005 -2.21115872 -3.73809125  2.20665727\n",
      "  0.82244106  3.95901132 -4.06340049]\n",
      "Pypower:  [ 0.          9.79601886  5.06056005 -2.21115872 -3.73809125  2.20665727\n",
      "  0.82244106  3.95901132 -4.06340049]\n",
      "Real power generation:\n",
      "Manual:  [ 66 163  85]\n",
      "PyPower:  [ 66 163  85]\n",
      "Real power flow on branch:\n",
      "Manual:  [  67.          28.9673913  -61.0326087   85.          23.9673913\n",
      "  -76.0326087 -163.          86.9673913  -38.0326087]\n",
      "PyPower:  [  67.          28.9673913  -61.0326087   85.          23.9673913\n",
      "  -76.0326087 -163.          86.9673913  -38.0326087]\n"
     ]
    }
   ],
   "source": [
    "print(\"Bus voltage magitude:\")\n",
    "print(\"Manual: \", r_bus[:,VM])\n",
    "print(\"Pypower: \", pp_bus[:,VM])\n",
    "\n",
    "print(\"Bus voltage angle:\")\n",
    "print(\"Manual: \", r_bus[:,VA])\n",
    "print(\"Pypower: \", pp_bus[:,VA])\n",
    "\n",
    "print(\"Real power generation:\")\n",
    "print(\"Manual: \", r_gen[:,PG])\n",
    "print(\"PyPower: \", pp_gen[:,PG])\n",
    "\n",
    "print(\"Real power flow on branch:\")\n",
    "print(\"Manual: \", r_branch[:,PF])\n",
    "print(\"PyPower: \", pp_branch[:,PF])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## OPF benchmark"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Manual solution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Converged in 0.00 seconds\n",
      "================================================================================\n",
      "|     System Summary                                                           |\n",
      "================================================================================\n",
      "\n",
      "How many?                How much?              P (MW)            Q (MVAr)\n",
      "---------------------    -------------------  -------------  -----------------\n",
      "Buses              9     Total Gen Capacity     820.0           0.0 to 0.0\n",
      "Generators         3     On-line Capacity       820.0           0.0 to 0.0\n",
      "Committed Gens     3     Generation (actual)    314.0               0.0\n",
      "Loads              3     Load                   315.0               0.0\n",
      "  Fixed            3       Fixed                315.0               0.0\n",
      "  Dispatchable     0       Dispatchable           0.0 of 0.0        0.0\n",
      "Shunts             0     Shunt (inj)              0.0               0.0\n",
      "Branches           9     Losses (I^2 * Z)         0.00              0.00\n",
      "Transformers       0     Branch Charging (inj)     -                0.0\n",
      "Inter-ties         0     Total Inter-tie Flow     0.0               0.0\n",
      "Areas              1\n",
      "\n",
      "                          Minimum                      Maximum\n",
      "                 -------------------------  --------------------------------\n",
      "Voltage Magnitude   1.000 p.u. @ bus 1          1.000 p.u. @ bus 1   \n",
      "Voltage Angle      -4.06 deg   @ bus 9          9.80 deg   @ bus 2   \n",
      "\n",
      "================================================================================\n",
      "|     Bus Data                                                                 |\n",
      "================================================================================\n",
      " Bus      Voltage          Generation             Load        \n",
      "  #   Mag(pu) Ang(deg)   P (MW)   Q (MVAr)   P (MW)   Q (MVAr)\n",
      "----- ------- --------  --------  --------  --------  --------\n",
      "    1  1.000    0.000*    66.00      0.00       -         -   \n",
      "    2  1.000    9.796    163.00      0.00       -         -   \n",
      "    3  1.000    5.061     85.00      0.00       -         -   \n",
      "    4  1.000   -2.211       -         -         -         -   \n",
      "    5  1.000   -3.738       -         -       90.00      0.00 \n",
      "    6  1.000    2.207       -         -         -         -   \n",
      "    7  1.000    0.822       -         -      100.00      0.00 \n",
      "    8  1.000    3.959       -         -         -         -   \n",
      "    9  1.000   -4.063       -         -      125.00      0.00 \n",
      "                        --------  --------  --------  --------\n",
      "               Total:    314.00      0.00    315.00      0.00\n",
      "\n",
      "================================================================================\n",
      "|     Branch Data                                                              |\n",
      "================================================================================\n",
      "Brnch   From   To    From Bus Injection   To Bus Injection     Loss (I^2 * Z)  \n",
      "  #     Bus    Bus    P (MW)   Q (MVAr)   P (MW)   Q (MVAr)   P (MW)   Q (MVAr)\n",
      "-----  -----  -----  --------  --------  --------  --------  --------  --------\n",
      "   0      1      4     67.00      0.00    -67.00      0.00     0.000      0.00\n",
      "   1      4      5     28.97      0.00    -28.97      0.00     0.000      0.00\n",
      "   2      5      6    -61.03      0.00     61.03      0.00     0.000      0.00\n",
      "   3      3      6     85.00      0.00    -85.00      0.00     0.000      0.00\n",
      "   4      6      7     23.97      0.00    -23.97      0.00     0.000      0.00\n",
      "   5      7      8    -76.03      0.00     76.03      0.00     0.000      0.00\n",
      "   6      8      2   -163.00      0.00    163.00      0.00     0.000      0.00\n",
      "   7      8      9     86.97      0.00    -86.97      0.00     0.000      0.00\n",
      "   8      9      4    -38.03      0.00     38.03      0.00     0.000      0.00\n",
      "                                                             --------  --------\n",
      "                                                    Total:     0.000      0.00\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "from numpy import flatnonzero as find\n",
    "from pypower.api import case9, ext2int, bustypes, makeBdc, rundcpf, ppoption, case39\n",
    "from pypower.idx_bus import BUS_TYPE, REF, VA, PD, LAM_P, LAM_Q, MU_VMAX, MU_VMIN\n",
    "from pypower.idx_gen import PG, MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN, PMAX, PMIN, GEN_BUS\n",
    "from pypower.idx_brch import PF, PT, QF, QT, RATE_A, MU_SF, MU_ST\n",
    "from pypower.idx_cost import COST\n",
    "\n",
    "# Run PyPower case to get initial condition\n",
    "ppc = case9()\n",
    "ppopt = ppoption(VERBOSE=0)\n",
    "ppopt = ppoption(ppopt, OPF_ALG_DC=200)\n",
    "pp_dcpf, success = rundcpf(ppc, ppopt)\n",
    "pp_bus = pp_dcpf[\"bus\"]\n",
    "pp_branch = pp_dcpf[\"branch\"]\n",
    "pp_gen = pp_dcpf[\"gen\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "## convert to internal indexing\n",
    "ppc = ext2int(ppc)\n",
    "baseMVA, bus, gen, branch, gencost = \\\n",
    "    ppc[\"baseMVA\"], ppc[\"bus\"], ppc[\"gen\"], ppc[\"branch\"], ppc[\"gencost\"]\n",
    "\n",
    "## get bus index lists of each type of bus\n",
    "ref, pv, pq = bustypes(bus, gen)\n",
    "pvpq = np.matrix(np.r_[pv, pq])\n",
    "\n",
    "## generator info\n",
    "gbus = gen[:, GEN_BUS].astype(int)    ## what buses are they at?\n",
    "refgen = find(gbus == ref)\n",
    "\n",
    "## build B matrices and phase shift injections\n",
    "B, Bf, _, _ = makeBdc(baseMVA, bus, branch)\n",
    "B = B.todense()\n",
    "Bf = Bf.todense()\n",
    "\n",
    "# Problem dimensions\n",
    "NG = gen.shape[0] # Number of generators\n",
    "NB = bus.shape[0] # Number of buses\n",
    "NBr = branch.shape[0] # Number of lines\n",
    "NL = 3 # Number of loads\n",
    "\n",
    "CG = np.zeros((NB,NG))\n",
    "CG[gbus,range(NG)] = 1\n",
    "\n",
    "# Generator capacity limit p.u.\n",
    "Pmax = gen[:, PMAX]/baseMVA\n",
    "Pmin = gen[:, PMIN]/baseMVA\n",
    "\n",
    "# Line flow limit p.u.\n",
    "Lmax = branch[:, RATE_A]/baseMVA\n",
    "Lmin = - branch[:, RATE_A]/baseMVA\n",
    "\n",
    "# Quadratic cost coefficients\n",
    "# Convert to p.u.\n",
    "Cg2 = gencost[:, COST]*baseMVA**2\n",
    "Cg1 = gencost[:, COST+1]*baseMVA\n",
    "Cg0 = gencost[:, COST+2]\n",
    "\n",
    "# Generator and demand set points\n",
    "# Power in p.u., Va in rad\n",
    "Pg0 = pp_gen[:, PG]/baseMVA\n",
    "Pd0 = pp_bus[:, PD]/baseMVA\n",
    "Pinj0 = (np.matmul(CG, Pg0) - Pd0)\n",
    "Va0 = pp_bus[:, VA]/180*np.pi"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pyomo.environ import *\n",
    "\n",
    "model = ConcreteModel(name=\"OPF\")\n",
    "\n",
    "# Define sets\n",
    "model.buses = Set(initialize=range(NB))\n",
    "model.lines = Set(initialize=range(NBr))\n",
    "model.gens = Set(initialize=range(NG))\n",
    "model.loads = Set(initialize=range(NL))\n",
    "\n",
    "# Define parameters\n",
    "def Pd_init(model, i):\n",
    "    return Pd0[i]\n",
    "model.Pd = Param(model.buses, initialize=Pd_init)\n",
    "\n",
    "def B_init(model, i, j):\n",
    "    return B[i,j]\n",
    "model.B = Param(model.buses, model.buses, initialize=B_init)\n",
    "\n",
    "def Bf_init(model, i, j):\n",
    "    return Bf[i,j]\n",
    "model.Bf = Param(model.lines, model.lines, initialize=Bf_init)\n",
    "\n",
    "def Pmax_init(model, i):\n",
    "    return Pmax[i]\n",
    "model.Pmax = Param(model.gens, initialize=Pmax_init)\n",
    "\n",
    "def Pmin_init(model, i):\n",
    "    return Pmin[i]\n",
    "model.Pmin = Param(model.gens, initialize=Pmin_init)\n",
    "\n",
    "def Lmax_init(model, i):\n",
    "    return Lmax[i]\n",
    "model.Lmax = Param(model.lines, initialize=Lmax_init)\n",
    "\n",
    "def Lmin_init(model, i):\n",
    "    return Lmin[i]\n",
    "model.Lmin = Param(model.lines, initialize=Lmin_init)\n",
    "\n",
    "def Cg0_init(model, i):\n",
    "    return Cg0[i]\n",
    "model.Cg0 = Param(model.gens, initialize=Cg0_init)\n",
    "\n",
    "def Cg1_init(model, i):\n",
    "    return Cg1[i]\n",
    "model.Cg1 = Param(model.gens, initialize=Cg1_init)\n",
    "\n",
    "def Cg2_init(model, i):\n",
    "    return Cg2[i]\n",
    "model.Cg2 = Param(model.gens, initialize=Cg2_init)\n",
    "\n",
    "def CG_init(model, i, j):\n",
    "    return CG[i,j]\n",
    "model.CG = Param(model.buses, model.gens, initialize=CG_init)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define variables\n",
    "def Pg_init(model, i):\n",
    "    return Pg0[i]\n",
    "model.Pg = Var(range(NG), within=PositiveReals, initialize=Pg_init, bounds=(0,3))\n",
    "\n",
    "def Pinj_init(model, i):\n",
    "    return Pinj0[i]\n",
    "model.Pinj = Var(range(NB), within=Reals, initialize=Pinj_init, bounds=(-2,2))\n",
    "\n",
    "def Va_init(model, i):\n",
    "    return Va0[i]\n",
    "model.Va = Var(range(NB), within=Reals, initialize=Va_init, bounds=(-1,1))\n",
    "\n",
    "model.dual = Suffix(direction=Suffix.IMPORT)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define objectives\n",
    "def Gen_cost(model):\n",
    "    cost0 = sum(model.Cg0[i] for i in model.gens)\n",
    "    cost1 = sum(model.Cg1[i] * model.Pg[i] for i in model.gens)\n",
    "    cost2 = sum(model.Cg2[i] * model.Pg[i]**2 for i in model.gens)\n",
    "    return (cost0 + cost1 + cost2)\n",
    "model.obj = Objective(rule=Gen_cost, sense=minimize)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define Constraints\n",
    "def Pmax_rule(model, i):\n",
    "    return model.Pg[i] <= model.Pmax[i]\n",
    "model.c_Pmax = Constraint(model.gens, rule=Pmax_rule)\n",
    "\n",
    "def Pmin_rule(model, i):\n",
    "    return model.Pg[i] >= model.Pmin[i]\n",
    "model.c_Pmin = Constraint(model.gens, rule=Pmin_rule)\n",
    "\n",
    "def Lmax_rule(model, i):\n",
    "    return sum(model.Bf[i,j]*model.Va[j] for j in model.buses) <= model.Lmax[i]\n",
    "model.c_Lmax = Constraint(model.lines, rule=Lmax_rule)\n",
    "\n",
    "def Lmin_rule(model, i):\n",
    "    return sum(model.Bf[i,j]*model.Va[j] for j in model.buses) >= model.Lmin[i]\n",
    "model.c_Lmin = Constraint(model.lines, rule=Lmin_rule)\n",
    "\n",
    "def Pinj_rule(model, i):\n",
    "    return model.Pinj[i] == (sum(model.CG[i,j]*model.Pg[j] for j in model.gens) - model.Pd[i])\n",
    "model.c_Pinj = Constraint(model.buses, rule=Pinj_rule)\n",
    "\n",
    "def PF_rule(model, i):\n",
    "    return model.Pinj[i] == sum(model.B[i,j]*model.Va[j] for j in model.buses)\n",
    "model.c_PF = Constraint(model.buses, rule=PF_rule)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solve the problem\n",
    "solver = SolverFactory(\"mindtpy\")\n",
    "solver.solve(model)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Model OPF\n",
      "\n",
      "  Variables:\n",
      "    Pg : Size=3, Index=Pg_index\n",
      "        Key : Lower : Value              : Upper : Fixed : Stale : Domain\n",
      "          0 :     0 : 0.8656449793509359 :     3 : False : False : PositiveReals\n",
      "          1 :     0 :  1.343775855535078 :     3 : False : False : PositiveReals\n",
      "          2 :     0 : 0.9405791651139859 :     3 : False : False : PositiveReals\n",
      "    Pinj : Size=9, Index=Pinj_index\n",
      "        Key : Lower : Value              : Upper : Fixed : Stale : Domain\n",
      "          0 :    -2 : 0.8656449793509359 :     2 : False : False :  Reals\n",
      "          1 :    -2 :  1.343775855535078 :     2 : False : False :  Reals\n",
      "          2 :    -2 : 0.9405791651139859 :     2 : False : False :  Reals\n",
      "          3 :    -2 :                0.0 :     2 : False : False :  Reals\n",
      "          4 :    -2 :               -0.9 :     2 : False : False :  Reals\n",
      "          5 :    -2 :                0.0 :     2 : False : False :  Reals\n",
      "          6 :    -2 :               -1.0 :     2 : False : False :  Reals\n",
      "          7 :    -2 :                0.0 :     2 : False : False :  Reals\n",
      "          8 :    -2 :              -1.25 :     2 : False : False :  Reals\n",
      "    Va : Size=9, Index=Va_index\n",
      "        Key : Lower : Value                 : Upper : Fixed : Stale : Domain\n",
      "          0 :    -1 :  0.017936769940642017 :     1 : False : False :  Reals\n",
      "          1 :    -1 :    0.1233578041060278 :     1 : False : False :  Reals\n",
      "          2 :    -1 :    0.0878006596754169 :     1 : False : False :  Reals\n",
      "          3 :    -1 :  -0.03192438086997189 :     1 : False : False :  Reals\n",
      "          4 :    -1 :  -0.06296310859816559 :     1 : False : False :  Reals\n",
      "          5 :    -1 :   0.03268272059973734 :     1 : False : False :  Reals\n",
      "          6 :    -1 : -0.005415308754642947 :     1 : False : False :  Reals\n",
      "          7 :    -1 :   0.03937181313508542 :     1 : False : False :  Reals\n",
      "          8 :    -1 :  -0.07682711871375289 :     1 : False : False :  Reals\n",
      "\n",
      "  Objectives:\n",
      "    obj : Size=1, Index=None, Active=True\n",
      "        Key  : Active : Value\n",
      "        None :   True : 5216.026607747273\n",
      "\n",
      "  Constraints:\n",
      "    c_Pmax : Size=3\n",
      "        Key : Lower : Body               : Upper\n",
      "          0 :  None : 0.8656449793509359 :   2.5\n",
      "          1 :  None :  1.343775855535078 :   3.0\n",
      "          2 :  None : 0.9405791651139859 :   2.7\n",
      "    c_Pmin : Size=3\n",
      "        Key : Lower : Body               : Upper\n",
      "          0 :   0.1 : 0.8656449793509359 :  None\n",
      "          1 :   0.1 :  1.343775855535078 :  None\n",
      "          2 :   0.1 : 0.9405791651139859 :  None\n",
      "    c_Lmax : Size=9\n",
      "        Key : Lower : Body                : Upper\n",
      "          0 :  None :  0.8656449793509359 :   2.5\n",
      "          1 :  None :  0.3373774753064533 :   2.5\n",
      "          2 :  None : -0.5626225246935467 :   1.5\n",
      "          3 :  None :  0.9405791651139856 :   3.0\n",
      "          4 :  None : 0.37795664042043936 :   1.5\n",
      "          5 :  None : -0.6220433595795607 :   2.5\n",
      "          6 :  None :  -1.343775855535078 :   2.5\n",
      "          7 :  None :  0.7217324959555175 :   2.5\n",
      "          8 :  None : -0.5282675040444824 :   2.5\n",
      "    c_Lmin : Size=9\n",
      "        Key : Lower : Body                : Upper\n",
      "          0 :  -2.5 :  0.8656449793509359 :  None\n",
      "          1 :  -2.5 :  0.3373774753064533 :  None\n",
      "          2 :  -1.5 : -0.5626225246935467 :  None\n",
      "          3 :  -3.0 :  0.9405791651139856 :  None\n",
      "          4 :  -1.5 : 0.37795664042043936 :  None\n",
      "          5 :  -2.5 : -0.6220433595795607 :  None\n",
      "          6 :  -2.5 :  -1.343775855535078 :  None\n",
      "          7 :  -2.5 :  0.7217324959555175 :  None\n",
      "          8 :  -2.5 : -0.5282675040444824 :  None\n",
      "    c_Pinj : Size=9\n",
      "        Key : Lower : Body  : Upper\n",
      "          0 :   0.0 :   0.0 :   0.0\n",
      "          1 :   0.0 :   0.0 :   0.0\n",
      "          2 :   0.0 :   0.0 :   0.0\n",
      "          3 :   0.0 :   0.0 :   0.0\n",
      "          4 :  -0.9 :  -0.9 :  -0.9\n",
      "          5 :   0.0 :   0.0 :   0.0\n",
      "          6 :  -1.0 :  -1.0 :  -1.0\n",
      "          7 :   0.0 :   0.0 :   0.0\n",
      "          8 : -1.25 : -1.25 : -1.25\n",
      "    c_PF : Size=9\n",
      "        Key : Lower : Body                    : Upper\n",
      "          0 :   0.0 :                     0.0 :   0.0\n",
      "          1 :   0.0 :                     0.0 :   0.0\n",
      "          2 :   0.0 :   2.220446049250313e-16 :   0.0\n",
      "          3 :   0.0 :  3.3306690738754696e-16 :   0.0\n",
      "          4 :   0.0 : -1.1102230246251565e-16 :   0.0\n",
      "          5 :   0.0 : -2.2898349882893854e-16 :   0.0\n",
      "          6 :   0.0 :                     0.0 :   0.0\n",
      "          7 :   0.0 :  1.6653345369377348e-16 :   0.0\n",
      "          8 :   0.0 :                     0.0 :   0.0\n"
     ]
    }
   ],
   "source": [
    "model.display()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "dual : Direction=Suffix.IMPORT, Datatype=Suffix.FLOAT\n",
      "    Key       : Value\n",
      "    c_Lmax[0] : -4.4325521001797886e-08\n",
      "    c_Lmax[1] :  -3.349818914657373e-08\n",
      "    c_Lmax[2] : -3.5122278666048176e-08\n",
      "    c_Lmax[3] : -3.5176834692668496e-08\n",
      "    c_Lmax[4] :  -6.456438359566623e-08\n",
      "    c_Lmax[5] :   -2.32037503681614e-08\n",
      "    c_Lmax[6] :  -1.884667120279402e-08\n",
      "    c_Lmax[7] :  -4.073906823229338e-08\n",
      "    c_Lmax[8] :  -2.392252900558152e-08\n",
      "    c_Lmin[0] :  2.1524446543098006e-08\n",
      "    c_Lmin[1] :   2.553193674238414e-08\n",
      "    c_Lmin[2] :    7.72851668465853e-08\n",
      "    c_Lmin[3] :   1.838383547929198e-08\n",
      "    c_Lmin[4] :   3.857639702718264e-08\n",
      "    c_Lmin[5] :   3.857639682176126e-08\n",
      "    c_Lmin[6] :     6.2657849865837e-08\n",
      "    c_Lmin[7] :   2.248573341263829e-08\n",
      "    c_Lmin[8] :     3.6741150727149e-08\n",
      "      c_PF[0] :       2404.418954510596\n",
      "      c_PF[1] :       2404.418954473698\n",
      "      c_PF[2] :      2404.4189544861547\n",
      "      c_PF[3] :       2404.418954533431\n",
      "      c_PF[4] :      2404.4189545435906\n",
      "      c_PF[5] :      2404.4189545035792\n",
      "      c_PF[6] :      2404.4189545322033\n",
      "      c_PF[7] :      2404.4189545185127\n",
      "      c_PF[8] :      2404.4189545437093\n",
      "    c_Pinj[0] :      -2404.418954472013\n",
      "    c_Pinj[1] :      -2404.418954384961\n",
      "    c_Pinj[2] :     -2404.4189544424094\n",
      "    c_Pinj[3] :      -2404.418954533431\n",
      "    c_Pinj[4] :     -2404.4189545844692\n",
      "    c_Pinj[5] :     -2404.4189545035792\n",
      "    c_Pinj[6] :        -2404.4189545805\n",
      "    c_Pinj[7] :     -2404.4189545185127\n",
      "    c_Pinj[8] :     -2404.4189546180123\n",
      "    c_Pmax[0] : -4.4325521001791765e-08\n",
      "    c_Pmax[1] :  -4.374145131051652e-08\n",
      "    c_Pmax[2] : -4.1174956913322165e-08\n",
      "    c_Pmin[0] :   9.462209130645108e-08\n",
      "    c_Pmin[1] :   5.824458089764794e-08\n",
      "    c_Pmin[2] :   8.618549340059997e-08\n"
     ]
    }
   ],
   "source": [
    "# Print dual variables (shadow prices)\n",
    "model.dual.pprint()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Average price = 24.044\n"
     ]
    }
   ],
   "source": [
    "# Calculate average price\n",
    "price = 0\n",
    "for i in range(NB):\n",
    "    price = price + model.dual[model.c_PF[i]]\n",
    "price = price/NB/baseMVA\n",
    "print(\"Average price = %.3f\" %price)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Pg = \n",
      " [[0.86564498]\n",
      " [1.34377586]\n",
      " [0.94057917]]\n",
      "Cost = \n",
      " [[1407.09784298]\n",
      " [2296.1266201 ]\n",
      " [1512.80214467]]\n",
      "MC = \n",
      " [[14.52209477]\n",
      " [12.62209477]\n",
      " [12.52209477]]\n"
     ]
    }
   ],
   "source": [
    "# Generation and cost\n",
    "r_pg = np.zeros((NG, 1))\n",
    "r_cost = np.zeros((NG, 1))\n",
    "r_mc = np.zeros((NG, 1))\n",
    "\n",
    "for i in range(NG):\n",
    "    r_pg[i] = model.Pg[i]()\n",
    "    r_cost[i] = Cg0[i] + Cg1[i]*r_pg[i] + Cg2[i]*r_pg[i]**2\n",
    "    r_mc[i] = Cg1[i]/baseMVA + Cg2[i]*r_pg[i]/baseMVA\n",
    "\n",
    "print(\"Pg = \\n\", r_pg)\n",
    "print(\"Cost = \\n\", r_cost)\n",
    "print(\"MC = \\n\", r_mc)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "5216.026607747273\n",
      "5216.026607747273\n"
     ]
    }
   ],
   "source": [
    "print(model.obj())\n",
    "print(np.sum(r_cost))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### PyPower solution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Converged in 0.06 seconds\n",
      "Objective Function Value = 5216.03 $/hr\n",
      "================================================================================\n",
      "|     System Summary                                                           |\n",
      "================================================================================\n",
      "\n",
      "How many?                How much?              P (MW)            Q (MVAr)\n",
      "---------------------    -------------------  -------------  -----------------\n",
      "Buses              9     Total Gen Capacity     820.0           0.0 to 0.0\n",
      "Generators         3     On-line Capacity       820.0           0.0 to 0.0\n",
      "Committed Gens     3     Generation (actual)    315.0               0.0\n",
      "Loads              3     Load                   315.0               0.0\n",
      "  Fixed            3       Fixed                315.0               0.0\n",
      "  Dispatchable     0       Dispatchable           0.0 of 0.0        0.0\n",
      "Shunts             0     Shunt (inj)              0.0               0.0\n",
      "Branches           9     Losses (I^2 * Z)         0.00              0.00\n",
      "Transformers       0     Branch Charging (inj)     -                0.0\n",
      "Inter-ties         0     Total Inter-tie Flow     0.0               0.0\n",
      "Areas              1\n",
      "\n",
      "                          Minimum                      Maximum\n",
      "                 -------------------------  --------------------------------\n",
      "Voltage Magnitude   1.000 p.u. @ bus 1          1.000 p.u. @ bus 1   \n",
      "Voltage Angle      -5.43 deg   @ bus 9          6.04 deg   @ bus 2   \n",
      "Lambda P           24.04 $/MWh @ bus 2         24.04 $/MWh @ bus 9   \n",
      "Lambda Q            0.00 $/MWh @ bus 1          0.00 $/MWh @ bus 1   \n",
      "\n",
      "================================================================================\n",
      "|     Bus Data                                                                 |\n",
      "================================================================================\n",
      " Bus      Voltage          Generation             Load          Lambda($/MVA-hr)\n",
      "  #   Mag(pu) Ang(deg)   P (MW)   Q (MVAr)   P (MW)   Q (MVAr)     P        Q   \n",
      "----- ------- --------  --------  --------  --------  --------  -------  -------\n",
      "    1  1.000    0.000*    86.56      0.00       -         -      24.044     -\n",
      "    2  1.000    6.040    134.38      0.00       -         -      24.044     -\n",
      "    3  1.000    4.003     94.06      0.00       -         -      24.044     -\n",
      "    4  1.000   -2.857       -         -         -         -      24.044     -\n",
      "    5  1.000   -4.635       -         -       90.00      0.00    24.044     -\n",
      "    6  1.000    0.845       -         -         -         -      24.044     -\n",
      "    7  1.000   -1.338       -         -      100.00      0.00    24.044     -\n",
      "    8  1.000    1.228       -         -         -         -      24.044     -\n",
      "    9  1.000   -5.430       -         -      125.00      0.00    24.044     -\n",
      "                        --------  --------  --------  --------\n",
      "               Total:    315.00      0.00    315.00      0.00\n",
      "\n",
      "================================================================================\n",
      "|     Branch Data                                                              |\n",
      "================================================================================\n",
      "Brnch   From   To    From Bus Injection   To Bus Injection     Loss (I^2 * Z)  \n",
      "  #     Bus    Bus    P (MW)   Q (MVAr)   P (MW)   Q (MVAr)   P (MW)   Q (MVAr)\n",
      "-----  -----  -----  --------  --------  --------  --------  --------  --------\n",
      "   0      1      4     86.56      0.00    -86.56      0.00     0.000      0.00\n",
      "   1      4      5     33.74      0.00    -33.74      0.00     0.000      0.00\n",
      "   2      5      6    -56.26      0.00     56.26      0.00     0.000      0.00\n",
      "   3      3      6     94.06      0.00    -94.06      0.00     0.000      0.00\n",
      "   4      6      7     37.80      0.00    -37.80      0.00     0.000      0.00\n",
      "   5      7      8    -62.20      0.00     62.20      0.00     0.000      0.00\n",
      "   6      8      2   -134.38      0.00    134.38      0.00     0.000      0.00\n",
      "   7      8      9     72.17      0.00    -72.17      0.00     0.000      0.00\n",
      "   8      9      4    -52.83      0.00     52.83      0.00     0.000      0.00\n",
      "                                                             --------  --------\n",
      "                                                    Total:     0.000      0.00\n"
     ]
    }
   ],
   "source": [
    "import pypower.api as pp\n",
    "ppc = case9()\n",
    "pp_dcopf = pp.rundcopf(ppc, ppopt)\n",
    "pp_bus = pp_dcopf[\"bus\"]\n",
    "pp_branch = pp_dcopf[\"branch\"]\n",
    "pp_gen = pp_dcopf[\"gen\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "pp_pg = pp_gen[:, PG]/baseMVA\n",
    "pp_va = pp_bus[:, VA]/180*np.pi\n",
    "pp_pd = pp_bus[:, PD]/baseMVA\n",
    "pp_pinj = (np.matmul(CG, pp_pg) - pp_pd)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 8.65644979e+01  1.34377586e+02  9.40579165e+01  3.49654948e-15\n",
      "  -9.00000000e+01  1.92557756e-14 -1.00000000e+02  2.77621649e-15\n",
      "  -1.25000000e+02]]\n",
      "[  86.56449793  134.37758556   94.05791651    0.          -90.\n",
      "    0.         -100.            0.         -125.        ]\n"
     ]
    }
   ],
   "source": [
    "# Power injection at buses\n",
    "pp_pbus = np.matmul(B, pp_va)*baseMVA\n",
    "print(pp_pbus)\n",
    "pp_pinj = (np.matmul(CG, pp_pg) - pp_pd)*baseMVA\n",
    "print(pp_pinj)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5216.0266077472725"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pp_cost = np.sum(Cg0) + np.sum(Cg1 * pp_pg) + np.sum(Cg2 * pp_pg**2)\n",
    "pp_cost"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.9.7 ('power')",
   "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.9.7"
  },
  "orig_nbformat": 4,
  "vscode": {
   "interpreter": {
    "hash": "f1e7b9e5b0ccd18e373e44bae5bf3ffa8873ffa9b0a195a69c8a44ccfc15a04d"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}