gurobi.ipynb

{
 "metadata": {
  "name": "",
  "signature": "sha256:5945faa8716449b3fe9ee52acca6466f05ce7745e373493c0192001b3d2e9fc2"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from gurobipy import Model, GRB, quicksum\n",
      "from jabr import build_constraint_matrix, convert_matrix_to_mosek\n",
      "from loadcase import load_case\n",
      "from numpy import real, imag, array, sqrt, pi"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 18
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "#case = load_case('cases/case5_renumber_tree.m')\n",
      "case = load_case('cases/case14_tree.m')\n",
      "#case = load_case('cases/case15_v2.m')\n",
      "#case = load_case('cases/case85_v2.m')\n",
      "#case = load_case('cases/case118_v2.m')\n",
      "#case = load_case('cases/case_ieee30_tree.m')\n",
      "#case = load_case('cases/case_ieee30_tree_v2.m')\n",
      "G, B = case.G, case.B\n",
      "P = real(case.demands)\n",
      "Q = imag(case.demands)\n",
      "branches = case.branch_list\n",
      "n = len(case.demands)\n",
      "vhat = case.vhat\n",
      "s2 = 2**.5\n",
      "gens = {}"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 23
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "#external_gens = {2:1.045, 5:1.01, 8:1.01, 11:1.082, 13:1.071}\n",
      "external_gens = {2:1.045, 3:1.01, 6:1.07, 8:1.09}\n",
      "#external_gens = {1: 0.955, 4: 0.998, 6: 0.990, 8: 1.015, 10: 1.050, 12: 0.990, 15: 0.970, 18: 0.973, 19: 0.962, 24: 0.992, 25: 1.050, 26: 1.015, 27: 0.968, 31: 0.967, 32: 0.963, 34: 0.984, 36: 0.980, 40: 0.970, 42: 0.985, 46: 1.005, 49: 1.025, 54: 0.955, 55: 0.952, 56: 0.954, 59: 0.985, 61: 0.995, 62: 0.998, 65: 1.005, 66: 1.050, 69: 1.035, 70: 0.984, 72: 0.980, 73: 0.991, 74: 0.958, 76: 0.943, 77: 1.006, 80: 1.040, 85: 0.985, 87: 1.015, 89: 1.005, 90: 0.985, 91: 0.980, 92: 0.990, 99: 1.010, 100: 1.017, 103: 1.010, 104: 0.971, 105: 0.965, 107: 0.952, 110: 0.973, 111: 0.980, 112: 0.975, 113: 0.993, 116: 1.005}\n",
      "gens = {case.i2e.index(g):v for g, v in external_gens.iteritems()}\n",
      "for gen, v in sorted(gens.iteritems()):\n",
      "    print('%3d %5.3f' % (gen, v))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "  1 1.045\n",
        "  2 1.010\n",
        "  5 1.070\n",
        "  7 1.090\n"
       ]
      }
     ],
     "prompt_number": 24
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "m = Model(\"jabr\")\n",
      "u = [m.addVar(name='u_%d'%i) for i in range(n)]\n",
      "R = {(i, j): m.addVar(name='R_%d_%d' % (i, j)) for i, j in branches}\n",
      "I = {(i, j): m.addVar(lb=-GRB.INFINITY, name='I_%d_%d' % (i, j)) for i, j in branches}\n",
      "for i, j in branches:\n",
      "    R[j, i] = R[i, j]\n",
      "    I[j, i] = I[i, j]\n",
      "m.update()\n",
      "m.addConstr(u[0] == vhat*vhat/s2, 'u0')\n",
      "for gen, v in gens.iteritems():\n",
      "    m.addConstr(u[gen] == v*v/s2, 'u%d' % gen)\n",
      "for i, j in branches:\n",
      "    m.addQConstr(2*u[i]*u[j] >= R[i,j]*R[i,j] + I[i,j]*I[i,j], 'cone_%d_%d' % (i, j))\n",
      "k = lambda i: (j for j in B[i, :].nonzero()[1])\n",
      "s = lambda i, j: 1 if i < j else -1\n",
      "for i in range(1, n):\n",
      "    m.addConstr(-s2*u[i]*G[i, :].sum() + quicksum(G[i,j]*R[i,j] + B[i,j]*s(i,j)*I[i,j] for j in k(i)) == P[i],\n",
      "                'real_flow_%d_%d' % (i, j))\n",
      "    if i in gens:\n",
      "        continue\n",
      "    m.addConstr(s2*u[i]*B[i, :].sum() + quicksum(-B[i,j]*R[i,j] + G[i,j]*s(i,j)*I[i,j] for j in k(i)) == Q[i],\n",
      "                'reac_flow_%d_%d' % (i, j))\n",
      "m.setObjective(quicksum(R[i,j] for i, j in branches), sense=GRB.MAXIMIZE)\n",
      "m.params.outputFlag = 0\n",
      "#m.setParam('BarConvTol', 1e-1)\n",
      "#m.setParam('DualReductions', 0)\n",
      "#m.setParam('BarHomogeneous', 1)\n",
      "m.optimize()\n",
      "print(m.status)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "2\n"
       ]
      }
     ],
     "prompt_number": 25
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "for i, x in enumerate(u):\n",
      "    y = x.getAttr('x')\n",
      "    print('%2d %f %.3f %d' % (i, (s2*y)**.5, y, 1 if i in gens else 0))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        " 0 1.060000 0.795 0\n",
        " 1 1.045000 0.772 1\n",
        " 2 1.010000 0.721 1\n",
        " 3 0.987370 0.689 0\n",
        " 4 0.997318 0.703 0\n",
        " 5 1.070000 0.810 1\n",
        " 6 1.009382 0.720 0\n",
        " 7 1.090000 0.840 1\n",
        " 8 0.975697 0.673 0\n",
        " 9 0.967655 0.662 0\n",
        "10 1.063496 0.800 0\n",
        "11 1.058978 0.793 0\n",
        "12 1.054283 0.786 0\n",
        "13 0.940518 0.625 0\n"
       ]
      }
     ],
     "prompt_number": 26
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "ans = {}\n",
      "for i, x in enumerate(u):\n",
      "    y = x.getAttr('x')\n",
      "    ans[case.i2e[i]] = y\n",
      "for i in sorted(ans.keys()):\n",
      "    y = ans[i]\n",
      "    print('%f' % (s2*y)**.5, end=' ')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "0.955000 0.980472 0.961572 0.998000 1.000238 0.990000 0.988690 1.015000 1.023043 1.050000 0.980390 0.990000 0.960078 0.986182 0.970000 0.948927 0.980796 0.973000 0.962000 0.919712 0.926202 0.945818 0.992908 0.992000 1.050000 1.015000 0.968000 0.948590 0.938225 0.987268 0.967000 0.963000 0.958661 0.984000 0.973965 0.980000 0.982412 0.978610 0.960911 0.970000 0.953767 0.985000 0.887454 0.920242 0.944878 1.005000 1.002062 1.015978 1.025000 1.011785 0.977020 0.970174 0.933770 0.955000 0.952000 0.954000 1.001990 0.971617 0.985000 0.992461 0.995000 0.998000 0.985913 0.992243 1.005000 1.050000 1.036778 1.010163 1.035000 0.984000 0.986261 0.980000 0.991000 0.958000 0.954734 0.943000 1.006000 1.000077 1.011476 1.040000 1.020536 0.963628 0.962438 0.976994 0.985000 0.983921 1.015000 0.990531 1.005000 0.985000 0.980000 0.990000 0.980754 0.967956 0.947796 0.997081 1.029095 1.023089 1.010000 1.017000 0.952092 0.975669 1.010000 0.971000 0.965000 0.946465 0.952000 0.959310 0.957528 0.973000 0.980000 0.975000 0.993000 0.959951 0.958697 1.005000 0.971351 0.941883 "
       ]
      }
     ],
     "prompt_number": 31
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "for p, q in zip(P, Q):\n",
      "    print('%7.3f %7.3f' % (p, q))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "  0.000   0.000\n",
        "-18.300  -5.184\n",
        "  2.400   1.200\n",
        "  7.600   1.600\n",
        " 94.200 -13.590\n",
        "  0.000   0.000\n",
        " 22.800  10.900\n",
        " 30.000 -68.944\n",
        "  0.000   0.000\n",
        "  5.800   2.000\n",
        "  0.000 -26.019\n",
        " 11.200   7.500\n",
        "  0.000 -17.640\n",
        "  6.200   1.600\n",
        "  8.200   2.500\n",
        "  3.500   1.800\n",
        "  9.000   5.800\n",
        "  3.200   0.900\n",
        "  9.500   3.400\n",
        "  2.200   0.700\n",
        " 17.500  11.200\n",
        "  0.000   0.000\n",
        "  3.200   1.600\n",
        "  8.700   6.700\n",
        "  0.000   0.000\n",
        "  3.500   2.300\n",
        "  0.000   0.000\n",
        "  0.000   0.000\n",
        "  2.400   0.900\n",
        " 10.600   1.900\n"
       ]
      }
     ],
     "prompt_number": 16
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from numpy.testing import assert_allclose"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 28
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%load_ext autoreload\n",
      "%autoreload 2"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 1
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import numpy as np\n",
      "import jabr\n",
      "import loadcase\n",
      "from collections import defaultdict, deque\n",
      "from numpy import pi, array, empty, set_printoptions"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 166
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "case = loadcase.load_case('cases/case118_v2.m')\n",
      "u, R, I = jabr.build_gurobi_model(case)\n",
      "V, theta_bus = jabr.recover_original_variables(u, I)\n",
      "matpower_equivalent = empty((len(V), 2))\n",
      "for i, (theta, v) in enumerate(zip(theta_bus, V)):\n",
      "    matpower_equivalent[case.i2e[i]-1] = [v, theta*180/pi]\n",
      "    #print('%2d %7.3f %7.3f' % (case.i2e[i], v, theta*180/pi))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 127
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "theta_bus"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 173,
       "text": [
        "[0,\n",
        " -0.4347134729659961,\n",
        " -0.4161570164716648,\n",
        " -0.41097743895066197,\n",
        " -0.31459509702899663,\n",
        " -0.3118184235941396,\n",
        " -0.3503439217766781,\n",
        " -0.35428773039166506,\n",
        " -0.2256784182421615,\n",
        " -0.09504751359704322,\n",
        " 0.03828248064638265,\n",
        " -0.3986213963446618,\n",
        " -0.4051983522664112,\n",
        " -0.4212464254468908,\n",
        " -0.4151163655351874,\n",
        " -0.4025149044443071,\n",
        " -0.40396086979980184,\n",
        " -0.36044776588705496,\n",
        " -0.41795099627348753,\n",
        " -0.4403661817086273,\n",
        " -0.37174548210956093,\n",
        " -0.35444902137695194,\n",
        " -0.3216342887717713,\n",
        " -0.25630355758593215,\n",
        " -0.2630414521109072,\n",
        " -0.12378427317677462,\n",
        " -0.08688369156919024,\n",
        " -0.3388425541760361,\n",
        " -0.37475141063384754,\n",
        " -0.39911814958692854,\n",
        " -0.266454286061434,\n",
        " -0.4364190445135751,\n",
        " -0.39330349591433766,\n",
        " -0.42412351287731387,\n",
        " -0.40326974975383656,\n",
        " -0.40951313251924687,\n",
        " -0.4114850696120179,\n",
        " -0.3934062019788969,\n",
        " -0.29844977213533974,\n",
        " -0.41998637518094817,\n",
        " -0.5231287770018425,\n",
        " -0.8568768256222524,\n",
        " -0.8080527953677621,\n",
        " -0.5962505639324057,\n",
        " -0.5473546451309743,\n",
        " -0.5159836261277151,\n",
        " -0.5488893074954835,\n",
        " -0.09119653452270865,\n",
        " -0.3318225970917832,\n",
        " -0.3240145854042134,\n",
        " -0.343241692791405,\n",
        " -0.38057422952202213,\n",
        " -0.3906695938630217,\n",
        " -0.6246634989265959,\n",
        " -0.5964377181744267,\n",
        " -0.437804096201168,\n",
        " -0.5292735166635318,\n",
        " -0.3577007072350632,\n",
        " -0.3888573803046458,\n",
        " -0.29026546978291023,\n",
        " -0.1732607194108599,\n",
        " -0.16267742307624466,\n",
        " -0.27266610034421723,\n",
        " -0.211344852928258,\n",
        " -0.17096164713015555,\n",
        " -0.11113874035250763,\n",
        " -0.11206311339043132,\n",
        " -0.1367317751317863,\n",
        " -0.07019436236147913,\n",
        " -0.10047460310341132,\n",
        " -0.10807532087608622,\n",
        " -0.13028514792875298,\n",
        " -0.11188167492283022,\n",
        " -0.219327572251852,\n",
        " -0.18515246756645756,\n",
        " -0.28512568182124776,\n",
        " -0.18666649002704383,\n",
        " -0.19444569574059417,\n",
        " -0.10007415910972703,\n",
        " -0.07871759318886007,\n",
        " -0.0738442820461552,\n",
        " -0.19342996122745731,\n",
        " -0.17908397917402136,\n",
        " -0.10105769319639624,\n",
        " -0.09592746319452701,\n",
        " -0.11910344791168981,\n",
        " -0.11495030090553282,\n",
        " 0.06060359020809605,\n",
        " 0.09354429359040091,\n",
        " -0.07966168111483253,\n",
        " -0.06089386317686637,\n",
        " -0.049623192258329354,\n",
        " -0.05824375398056989,\n",
        " -0.31960752338091875,\n",
        " -0.3350165056865919,\n",
        " -0.14030155034306369,\n",
        " -0.09026926315944876,\n",
        " -0.11144486314494761,\n",
        " -0.414590229541627,\n",
        " -0.3811115001717124,\n",
        " -0.08551173404538984,\n",
        " -0.06295689526329507,\n",
        " -0.41152375779648515,\n",
        " -0.6691022061526621,\n",
        " -0.7076067713182306,\n",
        " -0.4734869110640965,\n",
        " -0.8134781216446494,\n",
        " -0.7140729276261113,\n",
        " -0.7162382656992378,\n",
        " -0.5464526781667565,\n",
        " -0.5177406000143351,\n",
        " -0.6005442288130756,\n",
        " -0.4167328345981923,\n",
        " -0.3981616291283095,\n",
        " -0.3551726919563401,\n",
        " -0.07715949548807974,\n",
        " -0.43158141420573076,\n",
        " -0.20038579789230043]"
       ]
      }
     ],
     "prompt_number": 173
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "for i in range(118):\n",
      "    x = matpower_equivalent[i, :]\n",
      "    print('%f %f' % tuple(x))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "0.955000 -24.907247\n",
        "0.980472 -23.844041\n",
        "0.961572 -23.547273\n",
        "0.998000 -18.024971\n",
        "1.000238 -17.865880\n",
        "0.990000 -20.073228\n",
        "0.988689 -20.299192\n",
        "1.015000 -12.930421\n",
        "1.023043 -5.445821\n",
        "1.050000 2.193425\n",
        "0.980390 -22.839324\n",
        "0.990000 -23.216155\n",
        "0.960077 -24.135642\n",
        "0.986182 -23.784416\n",
        "0.970000 -23.062405\n",
        "0.948927 -23.145253\n",
        "0.980795 -20.652136\n",
        "0.973000 -23.946828\n",
        "0.962000 -25.231124\n",
        "0.919710 -21.299447\n",
        "0.926200 -20.308433\n",
        "0.945816 -18.428287\n",
        "0.992908 -14.685112\n",
        "0.992000 -15.071165\n",
        "1.050000 -7.092316\n",
        "1.015000 -4.978069\n",
        "0.968000 -19.414248\n",
        "0.948589 -21.471674\n",
        "0.938225 -22.867785\n",
        "0.987268 -15.266706\n",
        "0.967000 -25.004969\n",
        "0.963000 -22.534630\n",
        "0.958661 -24.300487\n",
        "0.984000 -23.105655\n",
        "0.973965 -23.463374\n",
        "0.980000 -23.576358\n",
        "0.982412 -22.540515\n",
        "0.978609 -17.099912\n",
        "0.960911 -24.063447\n",
        "0.970000 -29.973071\n",
        "0.953766 -49.095426\n",
        "0.985000 -46.298015\n",
        "0.887453 -34.162641\n",
        "0.920242 -31.361111\n",
        "0.944877 -29.563684\n",
        "1.005000 -31.449041\n",
        "1.002062 -5.225177\n",
        "1.015977 -19.012034\n",
        "1.025000 -18.564668\n",
        "1.011784 -19.666300\n",
        "0.977019 -21.805297\n",
        "0.970172 -22.383719\n",
        "0.933770 -35.790582\n",
        "0.955000 -34.173364\n",
        "0.952000 -25.084327\n",
        "0.954000 -30.325139\n",
        "1.001990 -20.494741\n",
        "0.971616 -22.279887\n",
        "0.985000 -16.630986\n",
        "0.992461 -9.927108\n",
        "0.995000 -9.320730\n",
        "0.998000 -15.622617\n",
        "0.985912 -12.109168\n",
        "0.992243 -9.795381\n",
        "1.005000 -6.367781\n",
        "1.050000 -6.420743\n",
        "1.036778 -7.834154\n",
        "1.010162 -4.021841\n",
        "1.035000 0.000000\n",
        "0.984000 -5.756771\n",
        "0.986261 -6.192260\n",
        "0.980000 -7.464789\n",
        "0.991000 -6.410348\n",
        "0.958000 -12.566544\n",
        "0.954734 -10.608455\n",
        "0.943000 -16.336498\n",
        "1.006000 -10.695202\n",
        "1.000076 -11.140918\n",
        "1.011475 -5.733827\n",
        "1.040000 -4.510186\n",
        "1.020536 -4.230966\n",
        "0.963627 -11.082720\n",
        "0.962437 -10.260756\n",
        "0.976994 -5.790179\n",
        "0.985000 -5.496239\n",
        "0.983920 -6.824125\n",
        "1.015000 -6.586167\n",
        "0.990531 3.472330\n",
        "1.005000 5.359693\n",
        "0.985000 -4.564278\n",
        "0.980000 -3.488961\n",
        "0.990000 -2.843199\n",
        "0.980754 -3.337121\n",
        "0.967955 -18.312162\n",
        "0.947796 -19.195032\n",
        "0.997081 -8.038687\n",
        "1.029095 -5.172048\n",
        "1.023089 -6.385320\n",
        "1.010000 -23.754270\n",
        "1.017000 -21.836080\n",
        "0.952091 -4.899461\n",
        "0.975669 -3.607164\n",
        "1.010000 -23.578574\n",
        "0.971000 -38.336732\n",
        "0.965000 -40.542882\n",
        "0.946465 -27.128802\n",
        "0.952000 -46.608863\n",
        "0.959309 -40.913365\n",
        "0.957527 -41.037430\n",
        "0.973000 -31.309432\n",
        "0.980000 -29.664351\n",
        "0.975000 -34.408650\n",
        "0.993000 -23.877033\n",
        "0.959951 -22.812981\n",
        "0.958697 -20.349896\n",
        "1.005000 -4.420913\n",
        "0.971351 -24.727794\n",
        "0.941883 -11.481260\n"
       ]
      }
     ],
     "prompt_number": 158
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "x = [(0.955000, -24.907252),(0.980472, -23.844046),(0.961572, -23.547277),(0.998000, -18.024976),(1.000238, -17.865884),(0.990000, -20.073233),(0.988689, -20.299196),(1.015000, -12.930425),(1.023043, -5.445826),(1.050000, 2.193420),(0.980390, -22.839328),(0.990000, -23.216160),(0.960077, -24.135647),(0.986182, -23.784421),(0.970000, -23.062410),(0.948927, -23.145257),(0.980795, -20.652140),(0.973000, -23.946833),(0.962000, -25.231128),(0.919710, -21.299451),(0.926200, -20.308437),(0.945816, -18.428291),(0.992908, -14.685116),(0.992000, -15.071169),(1.050000, -7.092321),(1.015000, -4.978073),(0.968000, -19.414252),(0.948589, -21.471678),(0.938225, -22.867790),(0.987268, -15.266710),(0.967000, -25.004973),(0.963000, -22.534634),(0.958661, -24.300492),(0.984000, -23.105659),(0.973965, -23.463378),(0.980000, -23.576362),(0.982412, -22.540519),(0.978609, -17.099916),(0.960911, -24.063451),(0.970000, -29.973075),(0.953766, -49.095428),(0.985000, -46.298017),(0.887453, -34.162643),(0.920242, -31.361113),(0.944877, -29.563686),(1.005000, -31.449043),(1.002062, -5.225177),(1.015977, -19.012037),(1.025000, -18.564671),(1.011784, -19.666303),(0.977019, -21.805299),(0.970172, -22.383721),(0.933770, -35.790584),(0.955000, -34.173366),(0.952000, -25.084330),(0.954000, -30.325142),(1.001990, -20.494743),(0.971616, -22.279889),(0.985000, -16.630989),(0.992461, -9.927111),(0.995000, -9.320732),(0.998000, -15.622620),(0.985912, -12.109171),(0.992243, -9.795384),(1.005000, -6.367784),(1.050000, -6.420746),(1.036778, -7.834157),(1.010162, -4.021841),(1.035000, 0.000000),(0.984000, -5.756771),(0.986261, -6.192260),(0.980000, -7.464789),(0.991000, -6.410348),(0.958000, -12.566544),(0.954734, -10.608455),(0.943000, -16.336499),(1.006000, -10.695203),(1.000076, -11.140918),(1.011475, -5.733827),(1.040000, -4.510186),(1.020536, -4.230966),(0.963627, -11.082721),(0.962437, -10.260757),(0.976994, -5.790181),(0.985000, -5.496240),(0.983920, -6.824126),(1.015000, -6.586168),(0.990531, 3.472328),(1.005000, 5.359691),(0.985000, -4.564280),(0.980000, -3.488963),(0.990000, -2.843202),(0.980754, -3.337123),(0.967955, -18.312165),(0.947796, -19.195035),(0.997081, -8.038687),(1.029095, -5.172048),(1.023089, -6.385321),(1.010000, -23.754273),(1.017000, -21.836083),(0.952091, -4.899464),(0.975669, -3.607166),(1.010000, -23.578577),(0.971000, -38.336736),(0.965000, -40.542885),(0.946465, -27.128805),(0.952000, -46.608866),(0.959309, -40.913368),(0.957527, -41.037433),(0.973000, -31.309435),(0.980000, -29.664354),(0.975000, -34.408653),(0.993000, -23.877037),(0.959951, -22.812985),(0.958697, -20.349900),(1.005000, -4.420914),(0.971351, -24.727798),(0.941883, -11.481261)]\n",
      "matpower_answer = array(x)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 170
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "set_printoptions(precision=5)\n",
      "print(abs(matpower_answer-matpower_equivalent).max(axis=0))\n",
      "x = matpower_answer-matpower_equivalent\n",
      "for i in range(118):\n",
      "    y = x[i, :]\n",
      "    if abs(y[0]) > .01:\n",
      "        print(i+1, y)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "[  4.66161e-07   5.34142e-06]\n"
       ]
      }
     ],
     "prompt_number": 172
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def recover_original_variables(u, I):\n",
      "    \"\"\" given Jabr variables u and I, return bus voltages and angles \"\"\"\n",
      "    V = sqrt(sqrt(2) * array(u))\n",
      "    theta_branch = {(i, j): asin(I[(i, j)]/(V[i]*V[j])) for (i, j) in I}\n",
      "    theta_bus = recover_bus_angles(theta_branch)\n",
      "    return V, theta_bus"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 161
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "np.real(case.demands)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 167,
       "text": [
        "array([   0.,   51.,   20.,   39.,   39.,    0.,   52.,   19.,   28.,    0., -450.,   70.,  -38.,\n",
        "         34.,   14.,   90.,   25.,   11.,   60.,   45.,   18.,   14.,   10.,    7.,   13., -220.,\n",
        "       -314.,   71.,   17.,   24.,    0.,   36.,   59.,   23.,   59.,   33.,   31.,    0.,    0.,\n",
        "         27.,   66.,   37.,   96.,   18.,   16.,   53.,    9.,   34.,   20., -117.,   17.,   17.,\n",
        "         18.,   23.,   65.,   63.,   84.,   12.,   12.,  122.,   78., -160.,   77.,    0.,    0.,\n",
        "       -391., -353.,   28.,    0.,   66.,    0.,   12.,    6.,   68.,   47.,   68.,   61.,   71.,\n",
        "         39., -347.,    0.,   54.,   20.,   11.,   24.,   21.,   -4.,   48., -607.,  163.,   10.,\n",
        "         65.,   12.,   30.,   42.,   38.,   15.,   34.,   42., -215.,   22.,    5.,  -17.,   38.,\n",
        "         31.,   43.,   50.,    2.,    8.,   39.,  -36.,   68.,    6.,    8.,   22.,  184.,   20.,\n",
        "         33.])"
       ]
      }
     ],
     "prompt_number": 167
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from jabr import solve"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 180
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "answer = solve('cases/case15_v2.m')\n",
      "for bus in sorted(answer.keys()):\n",
      "    v, t = answer[bus]\n",
      "    print('%3d %7.3f %7.3f' % (bus, v, 180/pi*t))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "  1   1.000   0.000\n",
        "  2   0.971   0.032\n",
        "  3   0.957   0.049\n",
        "  4   0.951   0.056\n",
        "  5   0.950   0.068\n",
        "  6   0.958   0.189\n",
        "  7   0.956   0.216\n",
        "  8   0.957   0.205\n",
        "  9   0.968   0.072\n",
        " 10   0.967   0.085\n",
        " 11   0.950   0.131\n",
        " 12   0.946   0.182\n",
        " 13   0.945   0.198\n",
        " 14   0.949   0.085\n",
        " 15   0.948   0.087\n"
       ]
      }
     ],
     "prompt_number": 181
    }
   ],
   "metadata": {}
  }
 ]
}