pypower_test_case39.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, case39\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 = case39()\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,38)\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             39     Total Gen Capacity    7367.0           0.0 to 0.0\n",
      "Generators        10     On-line Capacity      7367.0           0.0 to 0.0\n",
      "Committed Gens    10     Generation (actual)   6254.2               0.0\n",
      "Loads             21     Load                  6254.2               0.0\n",
      "  Fixed           21       Fixed               6254.2               0.0\n",
      "  Dispatchable     0       Dispatchable           0.0 of 0.0        0.0\n",
      "Shunts             0     Shunt (inj)              0.0               0.0\n",
      "Branches          46     Losses (I^2 * Z)         0.00              0.00\n",
      "Transformers      12     Branch Charging (inj)     -                0.0\n",
      "Inter-ties         6     Total Inter-tie Flow   736.4               0.0\n",
      "Areas              3\n",
      "\n",
      "                          Minimum                      Maximum\n",
      "                 -------------------------  --------------------------------\n",
      "Voltage Magnitude   1.000 p.u. @ bus 1          1.000 p.u. @ bus 1   \n",
      "Voltage Angle     -13.46 deg   @ bus 39         7.40 deg   @ bus 36  \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  -12.304       -         -       97.60      0.00 \n",
      "    2  1.000   -8.104       -         -         -         -   \n",
      "    3  1.000  -10.989       -         -      322.00      0.00 \n",
      "    4  1.000  -11.650       -         -      500.00      0.00 \n",
      "    5  1.000  -10.346       -         -         -         -   \n",
      "    6  1.000   -9.580       -         -         -         -   \n",
      "    7  1.000  -11.944       -         -      233.80      0.00 \n",
      "    8  1.000  -12.509       -         -      522.00      0.00 \n",
      "    9  1.000  -13.128       -         -        6.50      0.00 \n",
      "   10  1.000   -7.151       -         -         -         -   \n",
      "   11  1.000   -7.991       -         -         -         -   \n",
      "   12  1.000   -8.058       -         -        8.53      0.00 \n",
      "   13  1.000   -7.912       -         -         -         -   \n",
      "   14  1.000   -9.667       -         -         -         -   \n",
      "   15  1.000  -10.103       -         -      320.00      0.00 \n",
      "   16  1.000   -8.569       -         -      329.00      0.00 \n",
      "   17  1.000   -9.721       -         -         -         -   \n",
      "   18  1.000  -10.664       -         -      158.00      0.00 \n",
      "   19  1.000   -3.429       -         -         -         -   \n",
      "   20  1.000   -4.871       -         -      680.00      0.00 \n",
      "   21  1.000   -5.979       -         -      274.00      0.00 \n",
      "   22  1.000   -1.096       -         -         -         -   \n",
      "   23  1.000   -1.323       -         -      247.50      0.00 \n",
      "   24  1.000   -8.416       -         -      308.60      0.00 \n",
      "   25  1.000   -6.814       -         -      224.00      0.00 \n",
      "   26  1.000   -7.818       -         -      139.00      0.00 \n",
      "   27  1.000   -9.972       -         -      281.00      0.00 \n",
      "   28  1.000   -3.870       -         -      206.00      0.00 \n",
      "   29  1.000   -0.830       -         -      283.50      0.00 \n",
      "   30  1.000   -5.447    250.00      0.00       -         -   \n",
      "   31  1.000    0.000*   634.23      0.00      9.20      0.00 \n",
      "   32  1.000    0.819    650.00      0.00       -         -   \n",
      "   33  1.000    2.073    632.00      0.00       -         -   \n",
      "   34  1.000    0.415    508.00      0.00       -         -   \n",
      "   35  1.000    4.363    650.00      0.00       -         -   \n",
      "   36  1.000    7.405    560.00      0.00       -         -   \n",
      "   37  1.000    0.543    540.00      0.00       -         -   \n",
      "   38  1.000    6.774    830.00      0.00       -         -   \n",
      "   39  1.000  -13.461   1000.00      0.00   1104.00      0.00 \n",
      "                        --------  --------  --------  --------\n",
      "               Total:   6254.23      0.00   6254.23      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      2   -178.35      0.00    178.35      0.00     0.000      0.00\n",
      "   1      1     39     80.75      0.00    -80.75      0.00     0.000      0.00\n",
      "   2      2      3    333.43      0.00   -333.43      0.00     0.000      0.00\n",
      "   3      2     25   -261.78      0.00    261.78      0.00     0.000      0.00\n",
      "   4      2     30   -250.00      0.00    250.00      0.00     0.000      0.00\n",
      "   5      3      4     54.12      0.00    -54.12      0.00     0.000      0.00\n",
      "   6      3     18    -42.69      0.00     42.69      0.00     0.000      0.00\n",
      "   7      4      5   -177.69      0.00    177.69      0.00     0.000      0.00\n",
      "   8      4     14   -268.20      0.00    268.20      0.00     0.000      0.00\n",
      "   9      5      6   -514.75      0.00    514.75      0.00     0.000      0.00\n",
      "  10      5      8    337.07      0.00   -337.07      0.00     0.000      0.00\n",
      "  11      6      7    448.48      0.00   -448.48      0.00     0.000      0.00\n",
      "  12      6     11   -338.20      0.00    338.20      0.00     0.000      0.00\n",
      "  13      6     31   -625.03      0.00    625.03      0.00     0.000      0.00\n",
      "  14      7      8    214.68      0.00   -214.68      0.00     0.000      0.00\n",
      "  15      8      9     29.75      0.00    -29.75      0.00     0.000      0.00\n",
      "  16      9     39     23.25      0.00    -23.25      0.00     0.000      0.00\n",
      "  17     10     11    340.90      0.00   -340.90      0.00     0.000      0.00\n",
      "  18     10     13    309.10      0.00   -309.10      0.00     0.000      0.00\n",
      "  19     10     32   -650.00      0.00    650.00      0.00     0.000      0.00\n",
      "  20     12     11     -2.70      0.00      2.70      0.00     0.000      0.00\n",
      "  21     12     13     -5.83      0.00      5.83      0.00     0.000      0.00\n",
      "  22     13     14    303.27      0.00   -303.27      0.00     0.000      0.00\n",
      "  23     14     15     35.07      0.00    -35.07      0.00     0.000      0.00\n",
      "  24     15     16   -284.93      0.00    284.93      0.00     0.000      0.00\n",
      "  25     16     17    225.97      0.00   -225.97      0.00     0.000      0.00\n",
      "  26     16     19   -460.00      0.00    460.00      0.00     0.000      0.00\n",
      "  27     16     21   -334.78      0.00    334.78      0.00     0.000      0.00\n",
      "  28     16     24    -45.12      0.00     45.12      0.00     0.000      0.00\n",
      "  29     17     18    200.69      0.00   -200.69      0.00     0.000      0.00\n",
      "  30     17     27     25.28      0.00    -25.28      0.00     0.000      0.00\n",
      "  31     19     20    172.00      0.00   -172.00      0.00     0.000      0.00\n",
      "  32     19     33   -632.00      0.00    632.00      0.00     0.000      0.00\n",
      "  33     20     34   -508.00      0.00    508.00      0.00     0.000      0.00\n",
      "  34     21     22   -608.78      0.00    608.78      0.00     0.000      0.00\n",
      "  35     22     23     41.22      0.00    -41.22      0.00     0.000      0.00\n",
      "  36     22     35   -650.00      0.00    650.00      0.00     0.000      0.00\n",
      "  37     23     24    353.72      0.00   -353.72      0.00     0.000      0.00\n",
      "  38     23     36   -560.00      0.00    560.00      0.00     0.000      0.00\n",
      "  39     25     26     54.22      0.00    -54.22      0.00     0.000      0.00\n",
      "  40     25     37   -540.00      0.00    540.00      0.00     0.000      0.00\n",
      "  41     26     27    255.72      0.00   -255.72      0.00     0.000      0.00\n",
      "  42     26     28   -145.37      0.00    145.37      0.00     0.000      0.00\n",
      "  43     26     29   -195.13      0.00    195.13      0.00     0.000      0.00\n",
      "  44     28     29   -351.37      0.00    351.37      0.00     0.000      0.00\n",
      "  45     29     38   -830.00      0.00    830.00      0.00     0.000      0.00\n",
      "                                                             --------  --------\n",
      "                                                    Total:     0.000      0.00\n"
     ]
    }
   ],
   "source": [
    "ppc = case39()\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. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.\n",
      " 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]\n",
      "Pypower:  [1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.\n",
      " 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]\n",
      "Bus voltage angle:\n",
      "Manual:  [-12.3043699   -8.10439551 -10.98912011 -11.64954416 -10.34642157\n",
      "  -9.57959794 -11.94362211 -12.50942954 -13.12810348  -7.15074661\n",
      "  -7.99063881  -8.05839236  -7.91227145  -9.66724476 -10.10326532\n",
      "  -8.56868354  -9.72097327 -10.66384373  -3.42925212  -4.8708231\n",
      "  -5.97921632  -1.09597681  -1.32272635  -8.41614329  -6.81447228\n",
      "  -7.81782625  -9.97159067  -3.86996915  -0.83007595  -5.44694564\n",
      "   0.           0.81909632   2.07263699   0.41545511   4.36280697\n",
      "   7.40456678   0.54299335   6.77404802 -13.46108183]\n",
      "Pypower:  [-12.3043699   -8.10439551 -10.98912011 -11.64954416 -10.34642157\n",
      "  -9.57959794 -11.94362211 -12.50942954 -13.12810348  -7.15074661\n",
      "  -7.99063881  -8.05839236  -7.91227145  -9.66724476 -10.10326532\n",
      "  -8.56868354  -9.72097327 -10.66384373  -3.42925212  -4.8708231\n",
      "  -5.97921632  -1.09597681  -1.32272635  -8.41614329  -6.81447228\n",
      "  -7.81782625  -9.97159067  -3.86996915  -0.83007595  -5.44694564\n",
      "   0.           0.81909632   2.07263699   0.41545511   4.36280697\n",
      "   7.40456678   0.54299335   6.77404802 -13.46108183]\n",
      "Real power generation:\n",
      "Manual:  [ 250.    634.23  650.    632.    508.    650.    560.    540.    830.\n",
      " 1000.  ]\n",
      "PyPower:  [ 250.    634.23  650.    632.    508.    650.    560.    540.    830.\n",
      " 1000.  ]\n",
      "Real power flow on branch:\n",
      "Manual:  [-178.35372641   80.75372641  333.43008077 -261.78380719 -250.\n",
      "   54.11537211  -42.68529133 -177.68578059 -268.19884731 -514.75373515\n",
      "  337.06795456  448.47831902 -338.20205417 -625.03        214.67831902\n",
      "   29.74627359   23.24627359  340.90428326  309.09571674 -650.\n",
      "   -2.70222909   -5.82777091  303.26794583   35.06909852 -284.93090148\n",
      "  225.96909852 -460.         -334.77576923  -45.12423077  200.68529133\n",
      "   25.28380719  172.         -632.         -508.         -608.77576923\n",
      "   41.22423077 -650.          353.72423077 -560.           54.21619281\n",
      " -540.          255.71619281 -145.3652     -195.1348     -351.3652\n",
      " -830.        ]\n",
      "PyPower:  [-178.35372641   80.75372641  333.43008077 -261.78380719 -250.\n",
      "   54.11537211  -42.68529133 -177.68578059 -268.19884731 -514.75373515\n",
      "  337.06795456  448.47831902 -338.20205417 -625.03        214.67831902\n",
      "   29.74627359   23.24627359  340.90428326  309.09571674 -650.\n",
      "   -2.70222909   -5.82777091  303.26794583   35.06909852 -284.93090148\n",
      "  225.96909852 -460.         -334.77576923  -45.12423077  200.68529133\n",
      "   25.28380719  172.         -632.         -508.         -608.77576923\n",
      "   41.22423077 -650.          353.72423077 -560.           54.21619281\n",
      " -540.          255.71619281 -145.3652     -195.1348     -351.3652\n",
      " -830.        ]\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             39     Total Gen Capacity    7367.0           0.0 to 0.0\n",
      "Generators        10     On-line Capacity      7367.0           0.0 to 0.0\n",
      "Committed Gens    10     Generation (actual)   6254.2               0.0\n",
      "Loads             21     Load                  6254.2               0.0\n",
      "  Fixed           21       Fixed               6254.2               0.0\n",
      "  Dispatchable     0       Dispatchable           0.0 of 0.0        0.0\n",
      "Shunts             0     Shunt (inj)              0.0               0.0\n",
      "Branches          46     Losses (I^2 * Z)         0.00              0.00\n",
      "Transformers      12     Branch Charging (inj)     -                0.0\n",
      "Inter-ties         6     Total Inter-tie Flow   736.4               0.0\n",
      "Areas              3\n",
      "\n",
      "                          Minimum                      Maximum\n",
      "                 -------------------------  --------------------------------\n",
      "Voltage Magnitude   1.000 p.u. @ bus 1          1.000 p.u. @ bus 1   \n",
      "Voltage Angle     -13.46 deg   @ bus 39         7.40 deg   @ bus 36  \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  -12.304       -         -       97.60      0.00 \n",
      "    2  1.000   -8.104       -         -         -         -   \n",
      "    3  1.000  -10.989       -         -      322.00      0.00 \n",
      "    4  1.000  -11.650       -         -      500.00      0.00 \n",
      "    5  1.000  -10.346       -         -         -         -   \n",
      "    6  1.000   -9.580       -         -         -         -   \n",
      "    7  1.000  -11.944       -         -      233.80      0.00 \n",
      "    8  1.000  -12.509       -         -      522.00      0.00 \n",
      "    9  1.000  -13.128       -         -        6.50      0.00 \n",
      "   10  1.000   -7.151       -         -         -         -   \n",
      "   11  1.000   -7.991       -         -         -         -   \n",
      "   12  1.000   -8.058       -         -        8.53      0.00 \n",
      "   13  1.000   -7.912       -         -         -         -   \n",
      "   14  1.000   -9.667       -         -         -         -   \n",
      "   15  1.000  -10.103       -         -      320.00      0.00 \n",
      "   16  1.000   -8.569       -         -      329.00      0.00 \n",
      "   17  1.000   -9.721       -         -         -         -   \n",
      "   18  1.000  -10.664       -         -      158.00      0.00 \n",
      "   19  1.000   -3.429       -         -         -         -   \n",
      "   20  1.000   -4.871       -         -      680.00      0.00 \n",
      "   21  1.000   -5.979       -         -      274.00      0.00 \n",
      "   22  1.000   -1.096       -         -         -         -   \n",
      "   23  1.000   -1.323       -         -      247.50      0.00 \n",
      "   24  1.000   -8.416       -         -      308.60      0.00 \n",
      "   25  1.000   -6.814       -         -      224.00      0.00 \n",
      "   26  1.000   -7.818       -         -      139.00      0.00 \n",
      "   27  1.000   -9.972       -         -      281.00      0.00 \n",
      "   28  1.000   -3.870       -         -      206.00      0.00 \n",
      "   29  1.000   -0.830       -         -      283.50      0.00 \n",
      "   30  1.000   -5.447    250.00      0.00       -         -   \n",
      "   31  1.000    0.000*   634.23      0.00      9.20      0.00 \n",
      "   32  1.000    0.819    650.00      0.00       -         -   \n",
      "   33  1.000    2.073    632.00      0.00       -         -   \n",
      "   34  1.000    0.415    508.00      0.00       -         -   \n",
      "   35  1.000    4.363    650.00      0.00       -         -   \n",
      "   36  1.000    7.405    560.00      0.00       -         -   \n",
      "   37  1.000    0.543    540.00      0.00       -         -   \n",
      "   38  1.000    6.774    830.00      0.00       -         -   \n",
      "   39  1.000  -13.461   1000.00      0.00   1104.00      0.00 \n",
      "                        --------  --------  --------  --------\n",
      "               Total:   6254.23      0.00   6254.23      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      2   -178.35      0.00    178.35      0.00     0.000      0.00\n",
      "   1      1     39     80.75      0.00    -80.75      0.00     0.000      0.00\n",
      "   2      2      3    333.43      0.00   -333.43      0.00     0.000      0.00\n",
      "   3      2     25   -261.78      0.00    261.78      0.00     0.000      0.00\n",
      "   4      2     30   -250.00      0.00    250.00      0.00     0.000      0.00\n",
      "   5      3      4     54.12      0.00    -54.12      0.00     0.000      0.00\n",
      "   6      3     18    -42.69      0.00     42.69      0.00     0.000      0.00\n",
      "   7      4      5   -177.69      0.00    177.69      0.00     0.000      0.00\n",
      "   8      4     14   -268.20      0.00    268.20      0.00     0.000      0.00\n",
      "   9      5      6   -514.75      0.00    514.75      0.00     0.000      0.00\n",
      "  10      5      8    337.07      0.00   -337.07      0.00     0.000      0.00\n",
      "  11      6      7    448.48      0.00   -448.48      0.00     0.000      0.00\n",
      "  12      6     11   -338.20      0.00    338.20      0.00     0.000      0.00\n",
      "  13      6     31   -625.03      0.00    625.03      0.00     0.000      0.00\n",
      "  14      7      8    214.68      0.00   -214.68      0.00     0.000      0.00\n",
      "  15      8      9     29.75      0.00    -29.75      0.00     0.000      0.00\n",
      "  16      9     39     23.25      0.00    -23.25      0.00     0.000      0.00\n",
      "  17     10     11    340.90      0.00   -340.90      0.00     0.000      0.00\n",
      "  18     10     13    309.10      0.00   -309.10      0.00     0.000      0.00\n",
      "  19     10     32   -650.00      0.00    650.00      0.00     0.000      0.00\n",
      "  20     12     11     -2.70      0.00      2.70      0.00     0.000      0.00\n",
      "  21     12     13     -5.83      0.00      5.83      0.00     0.000      0.00\n",
      "  22     13     14    303.27      0.00   -303.27      0.00     0.000      0.00\n",
      "  23     14     15     35.07      0.00    -35.07      0.00     0.000      0.00\n",
      "  24     15     16   -284.93      0.00    284.93      0.00     0.000      0.00\n",
      "  25     16     17    225.97      0.00   -225.97      0.00     0.000      0.00\n",
      "  26     16     19   -460.00      0.00    460.00      0.00     0.000      0.00\n",
      "  27     16     21   -334.78      0.00    334.78      0.00     0.000      0.00\n",
      "  28     16     24    -45.12      0.00     45.12      0.00     0.000      0.00\n",
      "  29     17     18    200.69      0.00   -200.69      0.00     0.000      0.00\n",
      "  30     17     27     25.28      0.00    -25.28      0.00     0.000      0.00\n",
      "  31     19     20    172.00      0.00   -172.00      0.00     0.000      0.00\n",
      "  32     19     33   -632.00      0.00    632.00      0.00     0.000      0.00\n",
      "  33     20     34   -508.00      0.00    508.00      0.00     0.000      0.00\n",
      "  34     21     22   -608.78      0.00    608.78      0.00     0.000      0.00\n",
      "  35     22     23     41.22      0.00    -41.22      0.00     0.000      0.00\n",
      "  36     22     35   -650.00      0.00    650.00      0.00     0.000      0.00\n",
      "  37     23     24    353.72      0.00   -353.72      0.00     0.000      0.00\n",
      "  38     23     36   -560.00      0.00    560.00      0.00     0.000      0.00\n",
      "  39     25     26     54.22      0.00    -54.22      0.00     0.000      0.00\n",
      "  40     25     37   -540.00      0.00    540.00      0.00     0.000      0.00\n",
      "  41     26     27    255.72      0.00   -255.72      0.00     0.000      0.00\n",
      "  42     26     28   -145.37      0.00    145.37      0.00     0.000      0.00\n",
      "  43     26     29   -195.13      0.00    195.13      0.00     0.000      0.00\n",
      "  44     28     29   -351.37      0.00    351.37      0.00     0.000      0.00\n",
      "  45     29     38   -830.00      0.00    830.00      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 = case39()\n",
    "ppopt = ppoption(VERBOSE=0)\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.buses, 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)\n",
    "\n",
    "def Pinj_init(model, i):\n",
    "    return Pinj0[i]\n",
    "model.Pinj = Var(range(NB), within=Reals, initialize=Pinj_init)\n",
    "\n",
    "def Va_init(model, i):\n",
    "    return Va0[i]\n",
    "model.Va = Var(range(NB), within=Reals, initialize=Va_init)\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=10, Index=Pg_index\n",
      "        Key : Lower : Value              : Upper : Fixed : Stale : Domain\n",
      "          0 :     0 :  6.608459942695003 :  None : False : False : PositiveReals\n",
      "          1 :     0 :  6.460000062804398 :  None : False : False : PositiveReals\n",
      "          2 :     0 :  6.608459942471459 :  None : False : False : PositiveReals\n",
      "          3 :     0 :   6.52000005983922 :  None : False : False : PositiveReals\n",
      "          4 :     0 : 5.0800000506333065 :  None : False : False : PositiveReals\n",
      "          5 :     0 : 6.6084599418883965 :  None : False : False : PositiveReals\n",
      "          6 :     0 :  5.800000057684317 :  None : False : False : PositiveReals\n",
      "          7 :     0 :  5.640000056136645 :  None : False : False : PositiveReals\n",
      "          8 :     0 :  6.608459942793593 :  None : False : False : PositiveReals\n",
      "          9 :     0 :  6.608459943071543 :  None : False : False : PositiveReals\n",
      "    Pinj : Size=39, Index=Pinj_index\n",
      "        Key : Lower : Value                : Upper : Fixed : Stale : Domain\n",
      "          0 :  None :               -0.976 :  None : False : False :  Reals\n",
      "          1 :  None :                  0.0 :  None : False : False :  Reals\n",
      "          2 :  None :                -3.22 :  None : False : False :  Reals\n",
      "          3 :  None :                 -5.0 :  None : False : False :  Reals\n",
      "          4 :  None :                  0.0 :  None : False : False :  Reals\n",
      "          5 :  None :                  0.0 :  None : False : False :  Reals\n",
      "          6 :  None :               -2.338 :  None : False : False :  Reals\n",
      "          7 :  None :                -5.22 :  None : False : False :  Reals\n",
      "          8 :  None :               -0.065 :  None : False : False :  Reals\n",
      "          9 :  None :                  0.0 :  None : False : False :  Reals\n",
      "         10 :  None :                  0.0 :  None : False : False :  Reals\n",
      "         11 :  None : -0.08529999999999999 :  None : False : False :  Reals\n",
      "         12 :  None :                  0.0 :  None : False : False :  Reals\n",
      "         13 :  None :                  0.0 :  None : False : False :  Reals\n",
      "         14 :  None :                 -3.2 :  None : False : False :  Reals\n",
      "         15 :  None :                -3.29 :  None : False : False :  Reals\n",
      "         16 :  None :                  0.0 :  None : False : False :  Reals\n",
      "         17 :  None :                -1.58 :  None : False : False :  Reals\n",
      "         18 :  None :                  0.0 :  None : False : False :  Reals\n",
      "         19 :  None :                 -6.8 :  None : False : False :  Reals\n",
      "         20 :  None :                -2.74 :  None : False : False :  Reals\n",
      "         21 :  None :                  0.0 :  None : False : False :  Reals\n",
      "         22 :  None :               -2.475 :  None : False : False :  Reals\n",
      "         23 :  None :  -3.0860000000000003 :  None : False : False :  Reals\n",
      "         24 :  None :                -2.24 :  None : False : False :  Reals\n",
      "         25 :  None :                -1.39 :  None : False : False :  Reals\n",
      "         26 :  None :                -2.81 :  None : False : False :  Reals\n",
      "         27 :  None :                -2.06 :  None : False : False :  Reals\n",
      "         28 :  None :               -2.835 :  None : False : False :  Reals\n",
      "         29 :  None :    6.608459942695018 :  None : False : False :  Reals\n",
      "         30 :  None :   6.3680000628043745 :  None : False : False :  Reals\n",
      "         31 :  None :    6.608459942471462 :  None : False : False :  Reals\n",
      "         32 :  None :    6.520000059839218 :  None : False : False :  Reals\n",
      "         33 :  None :    5.080000050633311 :  None : False : False :  Reals\n",
      "         34 :  None :    6.608459941888401 :  None : False : False :  Reals\n",
      "         35 :  None :     5.80000005768433 :  None : False : False :  Reals\n",
      "         36 :  None :    5.640000056136649 :  None : False : False :  Reals\n",
      "         37 :  None :    6.608459942793603 :  None : False : False :  Reals\n",
      "         38 :  None :    -4.43154005692846 :  None : False : False :  Reals\n",
      "    Va : Size=39, Index=Va_index\n",
      "        Key : Lower : Value              : Upper : Fixed : Stale : Domain\n",
      "          0 :  None : 455.93166690709273 :  None : False : False :  Reals\n",
      "          1 :  None :  456.0891139676162 :  None : False : False :  Reals\n",
      "          2 :  None :  456.0210412895039 :  None : False : False :  Reals\n",
      "          3 :  None : 455.99270081549946 :  None : False : False :  Reals\n",
      "          4 :  None : 456.00541427170094 :  None : False : False :  Reals\n",
      "          5 :  None :  456.0186999123687 :  None : False : False :  Reals\n",
      "          6 :  None : 455.97193549996973 :  None : False : False :  Reals\n",
      "          7 :  None : 455.95930809377035 :  None : False : False :  Reals\n",
      "          8 :  None : 455.89971397235604 :  None : False : False :  Reals\n",
      "          9 :  None :  456.0664784434236 :  None : False : False :  Reals\n",
      "         10 :  None : 456.05006450266995 :  None : False : False :  Reals\n",
      "         11 :  None :  456.0504038478972 :  None : False : False :  Reals\n",
      "         12 :  None :  456.0544760064246 :  None : False : False :  Reals\n",
      "         13 :  None :  456.0272240863766 :  None : False : False :  Reals\n",
      "         14 :  None :  456.0267469607269 :  None : False : False :  Reals\n",
      "         15 :  None : 456.05662027956987 :  None : False : False :  Reals\n",
      "         16 :  None : 456.03490919711186 :  None : False : False :  Reals\n",
      "         17 :  None :  456.0216053997916 :  None : False : False :  Reals\n",
      "         18 :  None : 456.15022028172405 :  None : False : False :  Reals\n",
      "         19 :  None :  456.1250601224648 :  None : False : False :  Reals\n",
      "         20 :  None :  456.1041787596626 :  None : False : False :  Reals\n",
      "         21 :  None :   456.191858664944 :  None : False : False :  Reals\n",
      "         22 :  None : 456.18854081312344 :  None : False : False :  Reals\n",
      "         23 :  None :  456.0600694763417 :  None : False : False :  Reals\n",
      "         24 :  None : 456.10399620878167 :  None : False : False :  Reals\n",
      "         25 :  None : 456.05007113599686 :  None : False : False :  Reals\n",
      "         26 :  None :  456.0207745234465 :  None : False : False :  Reals\n",
      "         27 :  None :   456.078884741441 :  None : False : False :  Reals\n",
      "         28 :  None :  456.1191697592091 :  None : False : False :  Reals\n",
      "         29 :  None :   456.211717420703 :  None : False : False :  Reals\n",
      "         30 :  None :  456.1890439140486 :  None : False : False :  Reals\n",
      "         31 :  None :  456.2078994861925 :  None : False : False :  Reals\n",
      "         32 :  None :  456.2492851626332 :  None : False : False :  Reals\n",
      "         33 :  None : 456.21732308338443 :  None : False : False :  Reals\n",
      "         34 :  None :  456.2887221665423 :  None : False : False :  Reals\n",
      "         35 :  None :  456.3463008146924 :  None : False : False :  Reals\n",
      "         36 :  None : 456.23811541011656 :  None : False : False :  Reals\n",
      "         37 :  None :  456.2248390336944 :  None : False : False :  Reals\n",
      "         38 :  None :  455.8602961890128 :  None : False : False :  Reals\n",
      "\n",
      "  Objectives:\n",
      "    obj : Size=1, Index=None, Active=True\n",
      "        Key  : Active : Value\n",
      "        None :   True : 41263.9407472222\n",
      "\n",
      "  Constraints:\n",
      "    c_Pmax : Size=10\n",
      "        Key : Lower : Body               : Upper\n",
      "          0 :  None :  6.608459942695003 :  10.4\n",
      "          1 :  None :  6.460000062804398 :  6.46\n",
      "          2 :  None :  6.608459942471459 :  7.25\n",
      "          3 :  None :   6.52000005983922 :  6.52\n",
      "          4 :  None : 5.0800000506333065 :  5.08\n",
      "          5 :  None : 6.6084599418883965 :  6.87\n",
      "          6 :  None :  5.800000057684317 :   5.8\n",
      "          7 :  None :  5.640000056136645 :  5.64\n",
      "          8 :  None :  6.608459942793593 :  8.65\n",
      "          9 :  None :  6.608459943071543 :  11.0\n",
      "    c_Pmin : Size=10\n",
      "        Key : Lower : Body               : Upper\n",
      "          0 :   0.0 :  6.608459942695003 :  None\n",
      "          1 :   0.0 :  6.460000062804398 :  None\n",
      "          2 :   0.0 :  6.608459942471459 :  None\n",
      "          3 :   0.0 :   6.52000005983922 :  None\n",
      "          4 :   0.0 : 5.0800000506333065 :  None\n",
      "          5 :   0.0 : 6.6084599418883965 :  None\n",
      "          6 :   0.0 :  5.800000057684317 :  None\n",
      "          7 :   0.0 :  5.640000056136645 :  None\n",
      "          8 :   0.0 :  6.608459942793593 :  None\n",
      "          9 :   0.0 :  6.608459943071543 :  None\n",
      "    c_Lmax : Size=46\n",
      "        Key : Lower : Body                 : Upper\n",
      "          0 :  None :  -3.8308287231993745 :   6.0\n",
      "          1 :  None :    2.854828723200626 :  10.0\n",
      "          2 :  None :    4.508124378298817 :   5.0\n",
      "          3 :  None :  -1.7304931587714236 :   5.0\n",
      "          4 :  None :   -6.608459942694026 :   9.0\n",
      "          5 :  None :   1.3305386856554833 :   5.0\n",
      "          6 :  None : -0.04241430734691676 :   5.0\n",
      "          7 :  None :  -0.9932387657390791 :   6.0\n",
      "          8 :  None :  -2.6762225486090756 :   5.0\n",
      "          9 :  None :   -5.109861795295728 :  12.0\n",
      "         10 :  None :   4.1166230295202695 :   9.0\n",
      "         11 :  None :   5.0830883042290225 :   9.0\n",
      "         12 :  None :  -3.8249500367455767 :   4.8\n",
      "         13 :  None :   -6.368000062801002 :  18.0\n",
      "         14 :  None :   2.7450883042183705 :   9.0\n",
      "         15 :  None :   1.6417113337265619 :   9.0\n",
      "         16 :  None :   1.5767113337315095 :   9.0\n",
      "         17 :  None :   3.8171955240977695 :   6.0\n",
      "         18 :  None :   2.7912644183816155 :   6.0\n",
      "         19 :  None :   -6.608459942475747 :   9.0\n",
      "         20 :  None : 0.007754512631436228 :   5.0\n",
      "         21 :  None : -0.09305451263389841 :   5.0\n",
      "         22 :  None :    2.698209905749536 :   6.0\n",
      "         23 :  None : 0.021987357125908602 :   6.0\n",
      "         24 :  None :   -3.178012642871181 :   6.0\n",
      "         25 :  None :   2.4394474671935313 :   6.0\n",
      "         26 :  None :   -4.800000110470137 :   6.0\n",
      "         27 :  None :  -3.5228503772377735 :   6.0\n",
      "         28 :  None :  -0.5846096223540371 :   6.0\n",
      "         29 :  None :    1.622414307341387 :   6.0\n",
      "         30 :  None :   0.8170331598485063 :   6.0\n",
      "         31 :  None :   1.7199999493604992 :   9.0\n",
      "         32 :  None :   -6.520000059837912 :   9.0\n",
      "         33 :  None :   -5.080000050635135 :   9.0\n",
      "         34 :  None :   -6.262850377243012 :   9.0\n",
      "         35 :  None :   0.3456095646470203 :   6.0\n",
      "         36 :  None :  -6.6084599418936705 :   9.0\n",
      "         37 :  None :    3.670609622335178 :   6.0\n",
      "         38 :  None :   -5.800000057683064 :   9.0\n",
      "         39 :  None :   1.6695068973622256 :   6.0\n",
      "         40 :  None :   -5.640000056133431 :   9.0\n",
      "         41 :  None :   1.9929668401600793 :   6.0\n",
      "         42 :  None :  -0.6078819713966368 :   6.0\n",
      "         43 :  None :  -1.1055779713960874 :   6.0\n",
      "         44 :  None :  -2.6678819713997655 :   6.0\n",
      "         45 :  None :   -6.608459942795889 :  12.0\n",
      "    c_Lmin : Size=46\n",
      "        Key : Lower : Body                 : Upper\n",
      "          0 :  -6.0 :  -3.8308287231993745 :  None\n",
      "          1 : -10.0 :    2.854828723200626 :  None\n",
      "          2 :  -5.0 :    4.508124378298817 :  None\n",
      "          3 :  -5.0 :  -1.7304931587714236 :  None\n",
      "          4 :  -9.0 :   -6.608459942694026 :  None\n",
      "          5 :  -5.0 :   1.3305386856554833 :  None\n",
      "          6 :  -5.0 : -0.04241430734691676 :  None\n",
      "          7 :  -6.0 :  -0.9932387657390791 :  None\n",
      "          8 :  -5.0 :  -2.6762225486090756 :  None\n",
      "          9 : -12.0 :   -5.109861795295728 :  None\n",
      "         10 :  -9.0 :   4.1166230295202695 :  None\n",
      "         11 :  -9.0 :   5.0830883042290225 :  None\n",
      "         12 :  -4.8 :  -3.8249500367455767 :  None\n",
      "         13 : -18.0 :   -6.368000062801002 :  None\n",
      "         14 :  -9.0 :   2.7450883042183705 :  None\n",
      "         15 :  -9.0 :   1.6417113337265619 :  None\n",
      "         16 :  -9.0 :   1.5767113337315095 :  None\n",
      "         17 :  -6.0 :   3.8171955240977695 :  None\n",
      "         18 :  -6.0 :   2.7912644183816155 :  None\n",
      "         19 :  -9.0 :   -6.608459942475747 :  None\n",
      "         20 :  -5.0 : 0.007754512631436228 :  None\n",
      "         21 :  -5.0 : -0.09305451263389841 :  None\n",
      "         22 :  -6.0 :    2.698209905749536 :  None\n",
      "         23 :  -6.0 : 0.021987357125908602 :  None\n",
      "         24 :  -6.0 :   -3.178012642871181 :  None\n",
      "         25 :  -6.0 :   2.4394474671935313 :  None\n",
      "         26 :  -6.0 :   -4.800000110470137 :  None\n",
      "         27 :  -6.0 :  -3.5228503772377735 :  None\n",
      "         28 :  -6.0 :  -0.5846096223540371 :  None\n",
      "         29 :  -6.0 :    1.622414307341387 :  None\n",
      "         30 :  -6.0 :   0.8170331598485063 :  None\n",
      "         31 :  -9.0 :   1.7199999493604992 :  None\n",
      "         32 :  -9.0 :   -6.520000059837912 :  None\n",
      "         33 :  -9.0 :   -5.080000050635135 :  None\n",
      "         34 :  -9.0 :   -6.262850377243012 :  None\n",
      "         35 :  -6.0 :   0.3456095646470203 :  None\n",
      "         36 :  -9.0 :  -6.6084599418936705 :  None\n",
      "         37 :  -6.0 :    3.670609622335178 :  None\n",
      "         38 :  -9.0 :   -5.800000057683064 :  None\n",
      "         39 :  -6.0 :   1.6695068973622256 :  None\n",
      "         40 :  -9.0 :   -5.640000056133431 :  None\n",
      "         41 :  -6.0 :   1.9929668401600793 :  None\n",
      "         42 :  -6.0 :  -0.6078819713966368 :  None\n",
      "         43 :  -6.0 :  -1.1055779713960874 :  None\n",
      "         44 :  -6.0 :  -2.6678819713997655 :  None\n",
      "         45 : -12.0 :   -6.608459942795889 :  None\n",
      "    c_Pinj : Size=39\n",
      "        Key : Lower                : Body                    : Upper\n",
      "          0 :               -0.976 :                  -0.976 :               -0.976\n",
      "          1 :                  0.0 :                     0.0 :                  0.0\n",
      "          2 :                -3.22 :                   -3.22 :                -3.22\n",
      "          3 :                 -5.0 :                    -5.0 :                 -5.0\n",
      "          4 :                  0.0 :                     0.0 :                  0.0\n",
      "          5 :                  0.0 :                     0.0 :                  0.0\n",
      "          6 :               -2.338 :                  -2.338 :               -2.338\n",
      "          7 :                -5.22 :                   -5.22 :                -5.22\n",
      "          8 :               -0.065 :                  -0.065 :               -0.065\n",
      "          9 :                  0.0 :                     0.0 :                  0.0\n",
      "         10 :                  0.0 :                     0.0 :                  0.0\n",
      "         11 : -0.08529999999999999 :    -0.08529999999999999 : -0.08529999999999999\n",
      "         12 :                  0.0 :                     0.0 :                  0.0\n",
      "         13 :                  0.0 :                     0.0 :                  0.0\n",
      "         14 :                 -3.2 :                    -3.2 :                 -3.2\n",
      "         15 :                -3.29 :                   -3.29 :                -3.29\n",
      "         16 :                  0.0 :                     0.0 :                  0.0\n",
      "         17 :                -1.58 :                   -1.58 :                -1.58\n",
      "         18 :                  0.0 :                     0.0 :                  0.0\n",
      "         19 :                 -6.8 :                    -6.8 :                 -6.8\n",
      "         20 :                -2.74 :                   -2.74 :                -2.74\n",
      "         21 :                  0.0 :                     0.0 :                  0.0\n",
      "         22 :               -2.475 :                  -2.475 :               -2.475\n",
      "         23 :  -3.0860000000000003 :     -3.0860000000000003 :  -3.0860000000000003\n",
      "         24 :                -2.24 :                   -2.24 :                -2.24\n",
      "         25 :                -1.39 :                   -1.39 :                -1.39\n",
      "         26 :                -2.81 :                   -2.81 :                -2.81\n",
      "         27 :                -2.06 :                   -2.06 :                -2.06\n",
      "         28 :               -2.835 :                  -2.835 :               -2.835\n",
      "         29 :                  0.0 :   1.509903313490213e-14 :                  0.0\n",
      "         30 :                  0.0 :  -2.398081733190338e-14 :                  0.0\n",
      "         31 :                  0.0 :  2.6645352591003757e-15 :                  0.0\n",
      "         32 :                  0.0 : -2.6645352591003757e-15 :                  0.0\n",
      "         33 :                  0.0 :   4.440892098500626e-15 :                  0.0\n",
      "         34 :                  0.0 :   4.440892098500626e-15 :                  0.0\n",
      "         35 :                  0.0 :  1.3322676295501878e-14 :                  0.0\n",
      "         36 :                  0.0 :   3.552713678800501e-15 :                  0.0\n",
      "         37 :                  0.0 :   9.769962616701378e-15 :                  0.0\n",
      "         38 :                  0.0 :  -4.440892098500626e-15 :                  0.0\n",
      "    c_PF : Size=39\n",
      "        Key : Lower : Body                    : Upper\n",
      "          0 :   0.0 : -1.2514433933574765e-12 :   0.0\n",
      "          1 :   0.0 :  -3.637978807091713e-12 :   0.0\n",
      "          2 :   0.0 :  1.1639578190170141e-12 :   0.0\n",
      "          3 :   0.0 :  1.4551915228366852e-11 :   0.0\n",
      "          4 :   0.0 :  2.9103830456733704e-11 :   0.0\n",
      "          5 :   0.0 : -1.0913936421275139e-11 :   0.0\n",
      "          6 :   0.0 :  1.7927881401647028e-11 :   0.0\n",
      "          7 :   0.0 :  2.9833913117727207e-12 :   0.0\n",
      "          8 :   0.0 :  -1.309674590999066e-12 :   0.0\n",
      "          9 :   0.0 :  1.0913936421275139e-11 :   0.0\n",
      "         10 :   0.0 :  1.8189894035458565e-12 :   0.0\n",
      "         11 :   0.0 :   2.462197112862441e-12 :   0.0\n",
      "         12 :   0.0 :  -7.275957614183426e-12 :   0.0\n",
      "         13 :   0.0 :   3.637978807091713e-12 :   0.0\n",
      "         14 :   0.0 :  4.3653969328261155e-12 :   0.0\n",
      "         15 :   0.0 : -2.0954793455985055e-11 :   0.0\n",
      "         16 :   0.0 :   3.637978807091713e-12 :   0.0\n",
      "         17 :   0.0 :  1.7461587731304462e-12 :   0.0\n",
      "         18 :   0.0 :                     0.0 :   0.0\n",
      "         19 :   0.0 : -4.3653969328261155e-12 :   0.0\n",
      "         20 :   0.0 : -2.0374812947920873e-12 :   0.0\n",
      "         21 :   0.0 :   7.275957614183426e-12 :   0.0\n",
      "         22 :   0.0 :  -5.093259147770368e-12 :   0.0\n",
      "         23 :   0.0 :  -4.307665335545607e-12 :   0.0\n",
      "         24 :   0.0 :  -9.313438908975513e-12 :   0.0\n",
      "         25 :   0.0 :  2.1465051958102777e-12 :   0.0\n",
      "         26 :   0.0 : -2.3283597272438783e-12 :   0.0\n",
      "         27 :   0.0 :  1.3096190798478347e-12 :   0.0\n",
      "         28 :   0.0 :    2.76489942052649e-12 :   0.0\n",
      "         29 :   0.0 :   9.920952948050399e-13 :   0.0\n",
      "         30 :   0.0 :  3.3724134596013755e-12 :   0.0\n",
      "         31 :   0.0 :  -4.285460875053104e-12 :   0.0\n",
      "         32 :   0.0 :   1.305622276959184e-12 :   0.0\n",
      "         33 :   0.0 : -1.8243184740640572e-12 :   0.0\n",
      "         34 :   0.0 :  -5.269562564080843e-12 :   0.0\n",
      "         35 :   0.0 :  1.2656542480726785e-12 :   0.0\n",
      "         36 :   0.0 :  3.2178704145735537e-12 :   0.0\n",
      "         37 :   0.0 : -2.2861712523081223e-12 :   0.0\n",
      "         38 :   0.0 :   3.675282300719118e-12 :   0.0\n"
     ]
    }
   ],
   "source": [
    "model.display()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "dual : Direction=Suffix.IMPORT, Datatype=Suffix.FLOAT\n",
      "    Key        : Value\n",
      "     c_Lmax[0] :  -5.173499417519515e-09\n",
      "    c_Lmax[10] : -1.0418281750048098e-08\n",
      "    c_Lmax[11] : -1.2989281905691557e-08\n",
      "    c_Lmax[12] :  -5.897919730023587e-09\n",
      "    c_Lmax[13] :  -2.086925136876417e-09\n",
      "    c_Lmax[14] :   -8.13353602922655e-09\n",
      "    c_Lmax[15] :   -6.91429127541709e-09\n",
      "    c_Lmax[16] :  -6.853730421450688e-09\n",
      "    c_Lmax[17] : -2.3298548775740888e-08\n",
      "    c_Lmax[18] : -1.5844898008006302e-08\n",
      "    c_Lmax[19] : -3.2590877988859266e-09\n",
      "     c_Lmax[1] :  -7.119937510491267e-09\n",
      "    c_Lmax[20] : -1.0189377784184044e-08\n",
      "    c_Lmax[21] :  -9.987462167817988e-09\n",
      "    c_Lmax[22] : -1.5398013473608167e-08\n",
      "    c_Lmax[23] :  -8.505154514148243e-09\n",
      "    c_Lmax[24] :  -5.540473465502049e-09\n",
      "    c_Lmax[25] :  -1.430825971417092e-08\n",
      "    c_Lmax[26] :  -4.705398501407228e-09\n",
      "    c_Lmax[27] :  -5.342032374499292e-09\n",
      "    c_Lmax[28] :  -7.725990221684723e-09\n",
      "    c_Lmax[29] : -1.1627033938768815e-08\n",
      "     c_Lmax[2] : -1.0282603445728146e-07\n",
      "    c_Lmax[30] :  -9.819505280237523e-09\n",
      "    c_Lmax[31] :  -6.987019314066522e-09\n",
      "    c_Lmax[32] :  -3.275137771202952e-09\n",
      "    c_Lmax[33] : -3.6123827936719227e-09\n",
      "    c_Lmax[34] : -3.3326707655669507e-09\n",
      "    c_Lmax[35] :  -8.994775167491026e-09\n",
      "    c_Lmax[36] : -3.2590258275003474e-09\n",
      "    c_Lmax[37] :   -2.18319897673875e-08\n",
      "    c_Lmax[38] : -3.4366021424979328e-09\n",
      "    c_Lmax[39] : -1.1745614822270339e-08\n",
      "     c_Lmax[3] :  -7.557195816748373e-09\n",
      "    c_Lmax[40] : -3.4741728012132108e-09\n",
      "    c_Lmax[41] : -1.2686007165655086e-08\n",
      "    c_Lmax[42] :    -7.6992799306367e-09\n",
      "    c_Lmax[43] :  -7.159870327989912e-09\n",
      "    c_Lmax[44] :  -5.869096546513629e-09\n",
      "    c_Lmax[45] : -2.7335462671716487e-09\n",
      "     c_Lmax[4] :  -3.259113398363457e-09\n",
      "     c_Lmax[5] : -1.3861189920482072e-08\n",
      "     c_Lmax[6] : -1.0082701537587564e-08\n",
      "     c_Lmax[7] :  -7.274236803388102e-09\n",
      "     c_Lmax[8] :  -6.625627514343238e-09\n",
      "     c_Lmax[9] : -2.9729471781067515e-09\n",
      "     c_Lmin[0] :   2.346137324702204e-08\n",
      "    c_Lmin[10] :   3.877507987624023e-09\n",
      "    c_Lmin[11] :   3.611460294277785e-09\n",
      "    c_Lmin[12] :   5.214392868909402e-08\n",
      "    c_Lmin[13] :   4.373345844998491e-09\n",
      "    c_Lmin[14] :   4.330645367863154e-09\n",
      "    c_Lmin[15] :   4.778993818992638e-09\n",
      "    c_Lmin[16] :    4.80836212088666e-09\n",
      "    c_Lmin[17] :  5.1817892413705395e-09\n",
      "    c_Lmin[18] :    5.78731763760864e-09\n",
      "    c_Lmin[19] :  2.1250111556286493e-08\n",
      "     c_Lmin[1] :   3.956446073136702e-09\n",
      "    c_Lmin[20] :   1.015759118227316e-08\n",
      "    c_Lmin[21] :  1.0366515939218318e-08\n",
      "    c_Lmin[22] :   5.848993552819857e-09\n",
      "    c_Lmin[23] :   8.450636402580587e-09\n",
      "    c_Lmin[24] :  1.8043103993668917e-08\n",
      "    c_Lmin[25] :  6.0233821370890265e-09\n",
      "    c_Lmin[26] :   4.276162366210339e-08\n",
      "    c_Lmin[27] :   2.052527669939767e-08\n",
      "    c_Lmin[28] :   9.392200673934937e-09\n",
      "    c_Lmin[29] :   6.671031152192435e-09\n",
      "     c_Lmin[2] :   5.351300480688745e-09\n",
      "    c_Lmin[30] :   7.458672915419655e-09\n",
      "    c_Lmin[31] :   4.744852773842098e-09\n",
      "    c_Lmin[32] :  2.0594573890352012e-08\n",
      "    c_Lmin[33] :  1.2976793785884955e-08\n",
      "    c_Lmin[34] :  1.8576389886379966e-08\n",
      "    c_Lmin[35] :   8.017053566722778e-09\n",
      "    c_Lmin[36] :  2.1252707989736655e-08\n",
      "    c_Lmin[37] :   5.260078477517464e-09\n",
      "    c_Lmin[38] :  1.5897186456208342e-08\n",
      "    c_Lmin[39] :   6.632476503695531e-09\n",
      "     c_Lmin[3] :   1.556047973220606e-08\n",
      "    c_Lmin[40] :  1.5139979820140702e-08\n",
      "    c_Lmin[41] :   6.366001093304516e-09\n",
      "    c_Lmin[42] :    9.43146145561671e-09\n",
      "    c_Lmin[43] :  1.0390324848864654e-08\n",
      "    c_Lmin[44] :   1.526040911220049e-08\n",
      "    c_Lmin[45] :    9.43033839829192e-09\n",
      "     c_Lmin[4] :   2.124904169773555e-08\n",
      "     c_Lmin[5] :   8.035564782204822e-09\n",
      "     c_Lmin[6] :  1.0265900409841923e-08\n",
      "     c_Lmin[7] :   1.015855183233472e-08\n",
      "     c_Lmin[8] :  2.1899271181982663e-08\n",
      "     c_Lmin[9] :   7.383172259592871e-09\n",
      "       c_PF[0] :      1351.6919885944483\n",
      "      c_PF[10] :      1351.6919885903967\n",
      "      c_PF[11] :      1351.6919885901627\n",
      "      c_PF[12] :      1351.6919885895811\n",
      "      c_PF[13] :      1351.6919886075907\n",
      "      c_PF[14] :      1351.6919886109322\n",
      "      c_PF[15] :      1351.6919885998534\n",
      "      c_PF[16] :      1351.6919886094854\n",
      "      c_PF[17] :       1351.691988621691\n",
      "      c_PF[18] :      1351.6919885617972\n",
      "      c_PF[19] :        1351.69198856404\n",
      "       c_PF[1] :       1351.691988554998\n",
      "      c_PF[20] :      1351.6919885866103\n",
      "      c_PF[21] :      1351.6919885733798\n",
      "      c_PF[22] :       1351.691988575737\n",
      "      c_PF[23] :      1351.6919885973396\n",
      "      c_PF[24] :      1351.6919885532968\n",
      "      c_PF[25] :       1351.691988582079\n",
      "      c_PF[26] :        1351.69198859917\n",
      "      c_PF[27] :      1351.6919885833402\n",
      "      c_PF[28] :      1351.6919885749023\n",
      "      c_PF[29] :      1351.6919885370082\n",
      "       c_PF[2] :      1351.6919886336327\n",
      "      c_PF[30] :       1351.691988627487\n",
      "      c_PF[31] :       1351.691988557911\n",
      "      c_PF[32] :      1351.6919885444781\n",
      "      c_PF[33] :      1351.6919885546756\n",
      "      c_PF[34] :      1351.6919885553862\n",
      "      c_PF[35] :      1351.6919885632765\n",
      "      c_PF[36] :       1351.691988541631\n",
      "      c_PF[37] :       1351.691988568206\n",
      "      c_PF[38] :      1351.6919886104845\n",
      "       c_PF[3] :      1351.6919886317166\n",
      "       c_PF[4] :      1351.6919886329642\n",
      "       c_PF[5] :       1351.691988629773\n",
      "       c_PF[6] :      1351.6919886357598\n",
      "       c_PF[7] :      1351.6919886378669\n",
      "       c_PF[8] :      1351.6919886213116\n",
      "       c_PF[9] :      1351.6919885759016\n",
      "     c_Pinj[0] :     -1351.6919885944483\n",
      "    c_Pinj[10] :     -1351.6919885903967\n",
      "    c_Pinj[11] :     -1351.6919885901627\n",
      "    c_Pinj[12] :     -1351.6919885895811\n",
      "    c_Pinj[13] :     -1351.6919886075907\n",
      "    c_Pinj[14] :     -1351.6919886109322\n",
      "    c_Pinj[15] :     -1351.6919885998534\n",
      "    c_Pinj[16] :     -1351.6919886094856\n",
      "    c_Pinj[17] :     -1351.6919886216913\n",
      "    c_Pinj[18] :     -1351.6919885617972\n",
      "    c_Pinj[19] :       -1351.69198856404\n",
      "     c_Pinj[1] :      -1351.691988554998\n",
      "    c_Pinj[20] :     -1351.6919885866107\n",
      "    c_Pinj[21] :     -1351.6919885733798\n",
      "    c_Pinj[22] :     -1351.6919885757372\n",
      "    c_Pinj[23] :     -1351.6919885973396\n",
      "    c_Pinj[24] :     -1351.6919885532968\n",
      "    c_Pinj[25] :      -1351.691988582079\n",
      "    c_Pinj[26] :       -1351.69198859917\n",
      "    c_Pinj[27] :     -1351.6919885833402\n",
      "    c_Pinj[28] :     -1351.6919885749023\n",
      "    c_Pinj[29] :     -1351.6919885370082\n",
      "     c_Pinj[2] :      -1351.691988633633\n",
      "    c_Pinj[30] :      -1351.691988627487\n",
      "    c_Pinj[31] :      -1351.691988557911\n",
      "    c_Pinj[32] :     -1351.6919885444781\n",
      "    c_Pinj[33] :     -1351.6919885546756\n",
      "    c_Pinj[34] :     -1351.6919885553862\n",
      "    c_Pinj[35] :     -1351.6919885632765\n",
      "    c_Pinj[36] :      -1351.691988541631\n",
      "    c_Pinj[37] :      -1351.691988568206\n",
      "    c_Pinj[38] :     -1351.6919886104845\n",
      "     c_Pinj[3] :     -1351.6919886317169\n",
      "     c_Pinj[4] :      -1351.691988632964\n",
      "     c_Pinj[5] :     -1351.6919886297733\n",
      "     c_Pinj[6] :     -1351.6919886357598\n",
      "     c_Pinj[7] :     -1351.6919886378669\n",
      "     c_Pinj[8] :     -1351.6919886213116\n",
      "     c_Pinj[9] :     -1351.6919885759019\n",
      "     c_Pmax[0] : -1.3407703650077375e-08\n",
      "     c_Pmax[1] :      -29.69197608235188\n",
      "     c_Pmax[2] :  -7.901857901600595e-08\n",
      "     c_Pmax[3] :     -17.691976592213795\n",
      "     c_Pmax[4] :     -305.69197844804023\n",
      "     c_Pmax[5] : -1.9310568988243243e-07\n",
      "     c_Pmax[6] :      -161.6919770439529\n",
      "     c_Pmax[7] :      -193.6919773323394\n",
      "     c_Pmax[8] : -2.4887232934255322e-08\n",
      "     c_Pmax[9] : -1.1576259917656575e-08\n",
      "     c_Pmin[0] :   7.699982645048915e-09\n",
      "     c_Pmin[1] :    7.87218394929204e-09\n",
      "     c_Pmin[2] :   7.699839329273999e-09\n",
      "     c_Pmin[3] :    7.79007291903239e-09\n",
      "     c_Pmin[4] :  1.0013059877835816e-08\n",
      "     c_Pmin[5] :   7.699492552335724e-09\n",
      "     c_Pmin[6] :   8.769889155587413e-09\n",
      "     c_Pmin[7] :   9.018734570722497e-09\n",
      "     c_Pmin[8] :   7.700041035945482e-09\n",
      "     c_Pmin[9] :   7.700216908109779e-09\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 = 13.517\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",
      " [[6.60845994]\n",
      " [6.46000006]\n",
      " [6.60845994]\n",
      " [6.52000006]\n",
      " [5.08000005]\n",
      " [6.60845994]\n",
      " [5.80000006]\n",
      " [5.64000006]\n",
      " [6.60845994]\n",
      " [6.60845994]]\n",
      "Cost = \n",
      " [[4565.6280797 ]\n",
      " [4367.16008303]\n",
      " [4565.6280794 ]\n",
      " [4446.84007983]\n",
      " [2733.24005296]\n",
      " [4565.62807861]\n",
      " [3538.20006864]\n",
      " [3350.36006501]\n",
      " [4565.62807983]\n",
      " [4565.62808021]]\n",
      "MC = \n",
      " [[6.90845994]\n",
      " [6.76000006]\n",
      " [6.90845994]\n",
      " [6.82000006]\n",
      " [5.38000005]\n",
      " [6.90845994]\n",
      " [6.10000006]\n",
      " [5.94000006]\n",
      " [6.90845994]\n",
      " [6.90845994]]\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": [
      "41263.9407472222\n",
      "41263.9407472222\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": [
      "PYPOWER Version 5.1.4, 27-June-2018 -- DC Optimal Power Flow\n",
      "Python Interior Point Solver - PIPS, Version 1.0, 07-Feb-2011\n",
      "Converged!\n",
      "\n",
      "Converged in 0.06 seconds\n",
      "Objective Function Value = 41263.94 $/hr\n",
      "================================================================================\n",
      "|     System Summary                                                           |\n",
      "================================================================================\n",
      "\n",
      "How many?                How much?              P (MW)            Q (MVAr)\n",
      "---------------------    -------------------  -------------  -----------------\n",
      "Buses             39     Total Gen Capacity    7367.0           0.0 to 0.0\n",
      "Generators        10     On-line Capacity      7367.0           0.0 to 0.0\n",
      "Committed Gens    10     Generation (actual)   6254.2               0.0\n",
      "Loads             21     Load                  6254.2               0.0\n",
      "  Fixed           21       Fixed               6254.2               0.0\n",
      "  Dispatchable     0       Dispatchable           0.0 of 0.0        0.0\n",
      "Shunts             0     Shunt (inj)              0.0               0.0\n",
      "Branches          46     Losses (I^2 * Z)         0.00              0.00\n",
      "Transformers      12     Branch Charging (inj)     -                0.0\n",
      "Inter-ties         6     Total Inter-tie Flow   836.0               0.0\n",
      "Areas              3\n",
      "\n",
      "                          Minimum                      Maximum\n",
      "                 -------------------------  --------------------------------\n",
      "Voltage Magnitude   0.982 p.u. @ bus 31         1.064 p.u. @ bus 36  \n",
      "Voltage Angle     -18.84 deg   @ bus 39         9.01 deg   @ bus 36  \n",
      "Lambda P           13.52 $/MWh @ bus 30        13.52 $/MWh @ bus 8   \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.039  -14.747       -         -       97.60      0.00    13.517     -\n",
      "    2  1.048   -5.726       -         -         -         -      13.517     -\n",
      "    3  1.031   -9.626       -         -      322.00      0.00    13.517     -\n",
      "    4  1.004  -11.250       -         -      500.00      0.00    13.517     -\n",
      "    5  1.006  -10.521       -         -         -         -      13.517     -\n",
      "    6  1.008   -9.760       -         -         -         -      13.517     -\n",
      "    7  0.998  -12.439       -         -      233.80      0.00    13.517     -\n",
      "    8  0.998  -13.163       -         -      522.00      0.00    13.517     -\n",
      "    9  1.038  -16.577       -         -        6.50      0.00    13.517     -\n",
      "   10  1.018   -7.022       -         -         -         -      13.517     -\n",
      "   11  1.013   -7.963       -         -         -         -      13.517     -\n",
      "   12  1.001   -7.943       -         -        8.53      0.00    13.517     -\n",
      "   13  1.015   -7.710       -         -         -         -      13.517     -\n",
      "   14  1.012   -9.272       -         -         -         -      13.517     -\n",
      "   15  1.016   -9.299       -         -      320.00      0.00    13.517     -\n",
      "   16  1.033   -7.587       -         -      329.00      0.00    13.517     -\n",
      "   17  1.034   -8.831       -         -         -         -      13.517     -\n",
      "   18  1.032   -9.594       -         -      158.00      0.00    13.517     -\n",
      "   19  1.050   -2.224       -         -         -         -      13.517     -\n",
      "   20  0.991   -3.666       -         -      680.00      0.00    13.517     -\n",
      "   21  1.032   -4.862       -         -      274.00      0.00    13.517     -\n",
      "   22  1.050    0.161       -         -         -         -      13.517     -\n",
      "   23  1.045   -0.029       -         -      247.50      0.00    13.517     -\n",
      "   24  1.038   -7.390       -         -      308.60      0.00    13.517     -\n",
      "   25  1.058   -4.873       -         -      224.00      0.00    13.517     -\n",
      "   26  1.053   -7.963       -         -      139.00      0.00    13.517     -\n",
      "   27  1.038   -9.641       -         -      281.00      0.00    13.517     -\n",
      "   28  1.050   -6.312       -         -      206.00      0.00    13.517     -\n",
      "   29  1.050   -4.003       -         -      283.50      0.00    13.517     -\n",
      "   30  1.050    1.299    660.85      0.00       -         -      13.517     -\n",
      "   31  0.982    0.000*   646.00      0.00      9.20      0.00    13.517     -\n",
      "   32  0.984    1.080    660.85      0.00       -         -      13.517     -\n",
      "   33  0.997    3.452    652.00      0.00       -         -      13.517     -\n",
      "   34  1.012    1.620    508.00      0.00       -         -      13.517     -\n",
      "   35  1.049    5.711    660.85      0.00       -         -      13.517     -\n",
      "   36  1.064    9.010    580.00      0.00       -         -      13.517     -\n",
      "   37  1.028    2.812    564.00      0.00       -         -      13.517     -\n",
      "   38  1.026    2.051    660.85      0.00       -         -      13.517     -\n",
      "   39  1.030  -18.836    660.85      0.00   1104.00      0.00    13.517     -\n",
      "                        --------  --------  --------  --------\n",
      "               Total:   6254.23      0.00   6254.23      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      2   -383.08      0.00    383.08      0.00     0.000      0.00\n",
      "   1      1     39    285.48      0.00   -285.48      0.00     0.000      0.00\n",
      "   2      2      3    450.81      0.00   -450.81      0.00     0.000      0.00\n",
      "   3      2     25   -173.05      0.00    173.05      0.00     0.000      0.00\n",
      "   4      2     30   -660.85      0.00    660.85      0.00     0.000      0.00\n",
      "   5      3      4    133.05      0.00   -133.05      0.00     0.000      0.00\n",
      "   6      3     18     -4.24      0.00      4.24      0.00     0.000      0.00\n",
      "   7      4      5    -99.32      0.00     99.32      0.00     0.000      0.00\n",
      "   8      4     14   -267.62      0.00    267.62      0.00     0.000      0.00\n",
      "   9      5      6   -510.99      0.00    510.99      0.00     0.000      0.00\n",
      "  10      5      8    411.66      0.00   -411.66      0.00     0.000      0.00\n",
      "  11      6      7    508.31      0.00   -508.31      0.00     0.000      0.00\n",
      "  12      6     11   -382.50      0.00    382.50      0.00     0.000      0.00\n",
      "  13      6     31   -636.80      0.00    636.80      0.00     0.000      0.00\n",
      "  14      7      8    274.51      0.00   -274.51      0.00     0.000      0.00\n",
      "  15      8      9    164.17      0.00   -164.17      0.00     0.000      0.00\n",
      "  16      9     39    157.67      0.00   -157.67      0.00     0.000      0.00\n",
      "  17     10     11    381.72      0.00   -381.72      0.00     0.000      0.00\n",
      "  18     10     13    279.13      0.00   -279.13      0.00     0.000      0.00\n",
      "  19     10     32   -660.85      0.00    660.85      0.00     0.000      0.00\n",
      "  20     12     11      0.78      0.00     -0.78      0.00     0.000      0.00\n",
      "  21     12     13     -9.31      0.00      9.31      0.00     0.000      0.00\n",
      "  22     13     14    269.82      0.00   -269.82      0.00     0.000      0.00\n",
      "  23     14     15      2.20      0.00     -2.20      0.00     0.000      0.00\n",
      "  24     15     16   -317.80      0.00    317.80      0.00     0.000      0.00\n",
      "  25     16     17    243.94      0.00   -243.94      0.00     0.000      0.00\n",
      "  26     16     19   -480.00      0.00    480.00      0.00     0.000      0.00\n",
      "  27     16     21   -352.29      0.00    352.29      0.00     0.000      0.00\n",
      "  28     16     24    -58.46      0.00     58.46      0.00     0.000      0.00\n",
      "  29     17     18    162.24      0.00   -162.24      0.00     0.000      0.00\n",
      "  30     17     27     81.70      0.00    -81.70      0.00     0.000      0.00\n",
      "  31     19     20    172.00      0.00   -172.00      0.00     0.000      0.00\n",
      "  32     19     33   -652.00      0.00    652.00      0.00     0.000      0.00\n",
      "  33     20     34   -508.00      0.00    508.00      0.00     0.000      0.00\n",
      "  34     21     22   -626.29      0.00    626.29      0.00     0.000      0.00\n",
      "  35     22     23     34.56      0.00    -34.56      0.00     0.000      0.00\n",
      "  36     22     35   -660.85      0.00    660.85      0.00     0.000      0.00\n",
      "  37     23     24    367.06      0.00   -367.06      0.00     0.000      0.00\n",
      "  38     23     36   -580.00      0.00    580.00      0.00     0.000      0.00\n",
      "  39     25     26    166.95      0.00   -166.95      0.00     0.000      0.00\n",
      "  40     25     37   -564.00      0.00    564.00      0.00     0.000      0.00\n",
      "  41     26     27    199.30      0.00   -199.30      0.00     0.000      0.00\n",
      "  42     26     28    -60.79      0.00     60.79      0.00     0.000      0.00\n",
      "  43     26     29   -110.56      0.00    110.56      0.00     0.000      0.00\n",
      "  44     28     29   -266.79      0.00    266.79      0.00     0.000      0.00\n",
      "  45     29     38   -660.85      0.00    660.85      0.00     0.000      0.00\n",
      "                                                             --------  --------\n",
      "                                                    Total:     0.000      0.00\n",
      "\n",
      "================================================================================\n",
      "|     Generation Constraints                                                   |\n",
      "================================================================================\n",
      " Gen   Bus                Active Power Limits\n",
      "  #     #    Pmin mu    Pmin       Pg       Pmax    Pmax mu\n",
      "----  -----  -------  --------  --------  --------  -------\n",
      "   1    31      -        0.00    646.00    646.00    0.297\n",
      "   3    33      -        0.00    652.00    652.00    0.177\n",
      "   4    34      -        0.00    508.00    508.00    3.057\n",
      "   6    36      -        0.00    580.00    580.00    1.617\n",
      "   7    37      -        0.00    564.00    564.00    1.937\n"
     ]
    }
   ],
   "source": [
    "import pypower.api as pp\n",
    "ppc = case39()\n",
    "pp_dcopf = pp.rundcopf(ppc)\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": [
      "[[-9.76000000e+01 -5.24839748e-14 -3.22000000e+02 -5.00000000e+02\n",
      "   1.81675080e-12 -1.93406730e-12 -2.33800000e+02 -5.22000000e+02\n",
      "  -6.50000000e+00  1.49094301e-12 -5.73523412e-13 -8.53000000e+00\n",
      "  -1.26461458e-13 -8.94720261e-13 -3.20000000e+02 -3.29000000e+02\n",
      "   9.89100668e-13 -1.58000000e+02  2.15700133e-13 -6.80000000e+02\n",
      "  -2.74000000e+02 -1.15204693e-13 -2.47500000e+02 -3.08600000e+02\n",
      "  -2.24000000e+02 -1.39000000e+02 -2.81000000e+02 -2.06000000e+02\n",
      "  -2.83500000e+02  6.60846000e+02  6.36800000e+02  6.60846000e+02\n",
      "   6.52000000e+02  5.08000000e+02  6.60846000e+02  5.80000000e+02\n",
      "   5.64000000e+02  6.60846000e+02 -4.43154000e+02]]\n",
      "[ -97.6           0.         -322.         -500.            0.\n",
      "    0.         -233.8        -522.           -6.5           0.\n",
      "    0.           -8.53          0.            0.         -320.\n",
      " -329.            0.         -158.            0.         -680.\n",
      " -274.            0.         -247.5        -308.6        -224.\n",
      " -139.         -281.         -206.         -283.5         660.84600006\n",
      "  636.79999991  660.84600005  651.99999986  507.99999999  660.84600002\n",
      "  579.99999998  563.99999999  660.84600006 -443.15399992]\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": [],
   "source": [
    "pp_cost = np.sum(Cg0) + np.sum(Cg1 * pp_pg) + np.sum(Cg2 * pp_pg**2)"
   ]
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "a48edb53158dd8a459f3e7f13fbcd78db9c99f511fd1a201dbc67d4baa3f7a7a"
  },
  "kernelspec": {
   "display_name": "Python 3.9.7 64-bit ('power': conda)",
   "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
 },
 "nbformat": 4,
 "nbformat_minor": 2
}