experiments.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "using MatpowerCases\n",
    "using JuMP "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "using Gadfly\n",
    "using Interact"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "networkData = loadcase(\"case96\");\n",
    "# compared with Yury's data:\n",
    "# Line capacities are larger here\n",
    "# Costs are missing\n",
    "# 99 generators here versus 96 in yury's model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "networkData[\"branch\"][:,6];"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "networkData[\"gen\"];"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "type Gen\n",
    "    bus::Int\n",
    "    Pg  # real power output (MW)\n",
    "    Qg  # reactive power output (MVAr)\n",
    "    Qmax # max reactive output (MVAr)\n",
    "    Qmin # (MVAr)\n",
    "    Vg\n",
    "    mBase\n",
    "    status # <= 0 means out of service\n",
    "    Pmax # max real power output (MW)\n",
    "    Pmin # min real power output (MW)\n",
    "end\n",
    "\n",
    "type Bus\n",
    "    id::Int\n",
    "    bustype::Int # 1- PQ 2- PV 3- Reference 4- Isolated\n",
    "    Pd # real power demand (MW)\n",
    "    Qd # reactive power demand (MVAr)\n",
    "    Gs # shunt conductance (MW demanded at V = 1.0 p.u.)\n",
    "    Bs # shunt susceptance (MVAr injected at V = 1.0 p.u.)\n",
    "    area\n",
    "    Vm # voltage magnitude\n",
    "    Va # voltage angle\n",
    "    baseKV # base voltage (kV)\n",
    "    zone # loss zone\n",
    "    maxVM # maximum voltage magnitude\n",
    "    minVM # minimum voltage magnitude\n",
    "end\n",
    "\n",
    "type Branch\n",
    "    f::Int # from\n",
    "    t::Int # to\n",
    "    r # resistance\n",
    "    x # reactance\n",
    "    b # total line charging susceptance\n",
    "    rateA # MVA rating A (long term rating)\n",
    "    rateB # MVA rating B (short term rating)\n",
    "    rateC # MVA rating C (emergency rating)\n",
    "    ratio # transformer off nominal turns ratio\n",
    "    angle # transformer phase shift angle (degrees), positive => delay\n",
    "    status # 1 - in service, 0 - out of service\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "function process(networkData)\n",
    "    bus = Dict{Int,Bus}()\n",
    "    for i in 1:size(networkData[\"bus\"],1)\n",
    "        b = Bus(networkData[\"bus\"][i,:]...)\n",
    "        bus[b.id] = b\n",
    "    end\n",
    "    gen = Gen[]\n",
    "    for i in 1:size(networkData[\"gen\"],1)\n",
    "        g = Gen(networkData[\"gen\"][i,:]...)\n",
    "        push!(gen,g)\n",
    "    end\n",
    "    branch = Branch[]\n",
    "    for i in 1:size(networkData[\"branch\"],1)\n",
    "        b = Branch(networkData[\"branch\"][i,:]...)\n",
    "        push!(branch,b)\n",
    "    end\n",
    "    return bus, gen, branch\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "bus,gen,branch = process(networkData);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "showall([g.Qmax for g in gen])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "length(branch)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "function solve_dcopf(bus,gen,branch)\n",
    "    #Define the optimal power flow (OPF) model\n",
    "\n",
    "    N_gen = length(gen)\n",
    "    N_line = length(branch)\n",
    "    bus_id = collect(keys(bus))\n",
    "    N_bus = length(bus_id)\n",
    "    opf=Model()\n",
    "    \n",
    "    # Define decision variables    \n",
    "    @defVar(opf, gen[i].Pmin <= g[i=1:N_gen] <= gen[i].Pmax ) # real power output of generators\n",
    "    #@defVar(opf, w[v=1:N_wind] >=0 ) ; # wind power injection\n",
    "    @defVar(opf, -branch[i].rateA <= f[i=1:N_line] <= branch[i].rateA) # real power flows \n",
    "    @defVar(opf, theta[bus_id]) # bus angle \n",
    "\n",
    "\n",
    "    # Define the objective function\n",
    "    @setObjective(opf, Min, sum{g[i],i=1:N_gen})\n",
    "\n",
    "    # Define the power balance constraint\n",
    "    @addConstraint(opf, power_balance[b=bus_id],\n",
    "        sum{ g[i], i=1:N_gen; gen[i].bus == b} +\n",
    "        sum{ f[l], l=1:N_line; branch[l].t == b} -\n",
    "        sum{ f[l], l=1:N_line; branch[l].f == b} - bus[b].Pd == 0);\n",
    "    \n",
    "    # Calculate f[l]\n",
    "    @addConstraint(opf, def_flow[l=1:N_line], branch[l].x * f[l] == theta[branch[l].t] - theta[branch[l].f])\n",
    "\n",
    "    # Slack bus\n",
    "    for b in values(bus)\n",
    "        if b.bustype == 3\n",
    "            @addConstraint(opf, theta[b.id] == 0)  # set the slack bus\n",
    "        end\n",
    "    end\n",
    "\n",
    "    # Solve statement\n",
    "    status = solve(opf)\n",
    "    @show status\n",
    "    @show getValue(g)\n",
    "    @show getValue(f)\n",
    "\n",
    "    return status, getObjectiveValue(opf), getValue(g)[:], getValue(f)[:]\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "status, objval, gval, fval = solve_dcopf(bus,gen,branch)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "length(bus)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Plot dispatch of every generator\n",
    "plot(x=1:length(gen),y=gval, Geom.bar,\n",
    "Guide.XLabel(\"Index of generators \"), Guide.YLabel(\"Dispatch, MW\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "function solve_decoupled_opf(bus,gen,branch)\n",
    "    #Define the optimal power flow (OPF) model\n",
    "\n",
    "    N_gen = length(gen)\n",
    "    N_line = length(branch)\n",
    "    bus_id = collect(keys(bus))\n",
    "    N_bus = length(bus_id)\n",
    "    opf=Model()\n",
    "    \n",
    "    # Define decision variables    \n",
    "    @defVar(opf, gen[i].Pmin <= p[i=1:N_gen] <= gen[i].Pmax ) # real power output of generators\n",
    "    #@defVar(opf, gen[i].Qmin <= q[i=1:N_gen] <= gen[i].Qmax ) # reactive power output of generators\n",
    "    @defVar(opf, q[i=1:N_gen] >= 0) # infeasible unless we remove the bounds\n",
    "    #@defVar(opf, w[v=1:N_wind] >=0 ) ; # wind power injection\n",
    "    @defVar(opf, f_p[i=1:N_line]) # real power flows\n",
    "    @defVar(opf, f_q[i=1:N_line]) # reactive power flows\n",
    "    @defVar(opf, theta[bus_id]) # bus angle \n",
    "    @defVar(opf, v[bus_id]) # voltage\n",
    "\n",
    "\n",
    "    # Define the objective function\n",
    "    @setObjective(opf, Min, sum{p[i]+q[i],i=1:N_gen}) # missing cost info\n",
    "\n",
    "    # Define the real power balance constraint\n",
    "    @addConstraint(opf, power_balance[b=bus_id],\n",
    "        sum{ p[i], i=1:N_gen; gen[i].bus == b} +\n",
    "        sum{ f_p[l], l=1:N_line; branch[l].t == b} -\n",
    "        sum{ f_p[l], l=1:N_line; branch[l].f == b} - bus[b].Pd == 0);\n",
    "    \n",
    "    @addConstraint(opf, power_balance_q[b=bus_id],\n",
    "        sum{ q[i], i=1:N_gen; gen[i].bus == b} +\n",
    "        sum{ f_q[l], l=1:N_line; branch[l].t == b} -\n",
    "        sum{ f_q[l], l=1:N_line; branch[l].f == b} - bus[b].Qd == 0);\n",
    "    \n",
    "    # Calculate f_p[l]\n",
    "    @addConstraint(opf, def_flow[l=1:N_line], branch[l].x * f_p[l] == theta[branch[l].t] - theta[branch[l].f])\n",
    "    # Calculate f_q[l]\n",
    "    @addConstraint(opf, def_flow_q[l=1:N_line], branch[l].x * f_q[l] == v[branch[l].t] - v[branch[l].f])\n",
    "    \n",
    "    # line limits\n",
    "    for l in 1:N_line\n",
    "        #@addConstraint(opf, f_p[l]^2 + f_q[l]^2 <= 10*branch[l].rateA) # infeasible unless we remove line limits\n",
    "    end\n",
    "\n",
    "    # Slack bus\n",
    "    for b in values(bus)\n",
    "        if b.bustype == 3\n",
    "            @addConstraint(opf, theta[b.id] == 0)  # set the slack bus\n",
    "            @addConstraint(opf, v[b.id] == 1)\n",
    "        elseif b.bustype == 1 # PQ bus\n",
    "            @addConstraint(opf, b.minVM <= v[b.id] <= b.maxVM)\n",
    "        elseif b.bustype == 2 # PV bus\n",
    "            @addConstraint(opf, v[b.id] == 1)\n",
    "        end\n",
    "    end\n",
    "    \n",
    "\n",
    "    # Solve statement\n",
    "    status = solve(opf)\n",
    "    @show status\n",
    "    #@show getValue(p)\n",
    "    #@show getValue(f_p)\n",
    "\n",
    "    return status, getObjectiveValue(opf), getValue(p)[:], getValue(q)[:], getValue(f_p)[:], getValue(f_q)[:]\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "status_2, objval_2, pval_2, qval_2, fpval_2, fqval_q = solve_decoupled_opf(bus,gen,branch)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "showall(fpval_2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 0.3.8-pre",
   "language": "julia",
   "name": "julia-0.3"
  },
  "language_info": {
   "name": "julia",
   "version": "0.3.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}