imbalance_reserves-end_to_end-withARB.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "5e9d1bed",
   "metadata": {},
   "source": [
    "Load packages"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "a1593cd1",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "import pyomo.environ as pyo\n",
    "from datetime import date\n",
    "import pandas as pd\n",
    "import ast\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7a2319d8",
   "metadata": {},
   "source": [
    "Define DA model with Imbalance Reserves"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "8d6954ca",
   "metadata": {},
   "outputs": [],
   "source": [
    "DAIRmodel = pyo.AbstractModel()\n",
    "\n",
    "#Sets\n",
    "DAIRmodel.G = pyo.RangeSet(1, 5)\n",
    "DAIRmodel.S = pyo.RangeSet(1, 5)\n",
    "DAIRmodel.R = pyo.RangeSet(1, 4)\n",
    "\n",
    "#Parameters\n",
    "DAIRmodel.VC = pyo.Param(DAIRmodel.G) #Variable Costs of generators\n",
    "DAIRmodel.CAP = pyo.Param(DAIRmodel.G) #Capacity of generators\n",
    "DAIRmodel.RR= pyo.Param(DAIRmodel.G)  #Ramp rate\n",
    "DAIRmodel.DEMAND=pyo.Param(within=pyo.NonNegativeIntegers) #Demand\n",
    "DAIRmodel.D1=pyo.Param(within=pyo.NonNegativeIntegers) #linear cost coefficient for demand slack variable\n",
    "DAIRmodel.D2=pyo.Param(within=pyo.NonNegativeIntegers) #quadratic cost coefficient for demand slack variable\n",
    "DAIRmodel.REDA=pyo.Param(within=pyo.NonNegativeIntegers) #Renewable day-ahead\n",
    "DAIRmodel.virtualcost=pyo.Param()\n",
    "DAIRmodel.VUB=pyo.Param()\n",
    "DAIRmodel.VLB=pyo.Param()\n",
    "\n",
    "#Parameters related to IR\n",
    "DAIRmodel.avup= pyo.Param(DAIRmodel.G) #Capacity bid for IR up\n",
    "DAIRmodel.avdn= pyo.Param(DAIRmodel.G) #Capacity bid for IR down\n",
    "DAIRmodel.PEN=pyo.Param(DAIRmodel.R,within=pyo.NonNegativeIntegers) #Scarcity price for IR up\n",
    "DAIRmodel.PEND=pyo.Param(DAIRmodel.R,within=pyo.NonNegativeReals) #Scarcity price for IR down\n",
    "DAIRmodel.probTU=pyo.Param(DAIRmodel.R) #probability of scarcity if IR up at level R is unmet\n",
    "DAIRmodel.probTD=pyo.Param(DAIRmodel.R) #probability of surplus if IR down at level R is unmet\n",
    "\n",
    "#Parameters related to IR requirement estimation - typically exogenous to model\n",
    "DAIRmodel.prob=pyo.Param(DAIRmodel.S)\n",
    "DAIRmodel.RE = pyo.Param(DAIRmodel.S)\n",
    "\n",
    "# Decision variables\n",
    "DAIRmodel.xDA=pyo.Var(DAIRmodel.G, domain=pyo.NonNegativeReals)\n",
    "DAIRmodel.rgDA=pyo.Var(domain=pyo.NonNegativeReals)\n",
    "DAIRmodel.d=pyo.Var(domain=pyo.NonNegativeReals)\n",
    "DAIRmodel.arb=pyo.Var()\n",
    "#Decision variables for IR\n",
    "DAIRmodel.isu  = pyo.Var(DAIRmodel.G, domain=pyo.NonNegativeReals) #IR up provided by supplier g\n",
    "DAIRmodel.isd  = pyo.Var(DAIRmodel.G, domain=pyo.NonNegativeReals) #IR down provided by supplier g\n",
    "DAIRmodel.siu  = pyo.Var(DAIRmodel.R, domain=pyo.NonNegativeReals) #Shortfall of IR up\n",
    "DAIRmodel.sid  = pyo.Var(DAIRmodel.R, domain=pyo.NonNegativeReals) #Shortfall of IR down\n",
    "\n",
    "#Decision variables for IR supply by renewables\n",
    "DAIRmodel.iupre =pyo.Var(domain=pyo.NonNegativeReals)\n",
    "DAIRmodel.idnre =pyo.Var(domain=pyo.NonNegativeReals)\n",
    "\n",
    "def obj_expression(m):\n",
    "    return sum(m.VC[g]*m.xDA[g] for g in m.G)\\\n",
    "+sum(m.probTU[r]*m.PEN[r]*m.siu[r] for r in m.R)\\\n",
    "+sum(m.probTD[r]*m.PEND[r]*m.sid[r] for r in m.R)\\\n",
    "+sum(m.avup[g]*m.isu[g] for g in m.G)\\\n",
    "+sum(m.avdn[g]*m.isd[g] for g in m.G)+m.D1 *m.d+ m.D2 *m.d*m.d +m.virtualcost*m.arb\n",
    "DAIRmodel.OBJ = pyo.Objective(rule=obj_expression)\n",
    "\n",
    "def DA_energy_balance(model):\n",
    "    return sum(model.xDA[g] for g in model.G)+ model.rgDA +model.arb +model.d == model.DEMAND\n",
    "DAIRmodel.Con3 = pyo.Constraint(expr=DA_energy_balance) \n",
    "\n",
    "def DA_IRup_balance(model):\n",
    "    return sum(model.isu[g] for g in model.G) +sum(model.siu[t] for t in model.R) +model.iupre >= np.mean(np.array([model.RE[s] for s in model.S]))-min(model.RE[s] for s in model.S)\n",
    "DAIRmodel.Con4UP= pyo.Constraint(rule=DA_IRup_balance)\n",
    "\n",
    "def shortage_IR_up(model,r):\n",
    "    return model.siu[r]<=max(0, min(np.mean(np.array([model.RE[s] for s in model.S])),model.RE[r+1])-model.RE[r])\n",
    "DAIRmodel.Con4UPstep= pyo.Constraint(DAIRmodel.R,rule=shortage_IR_up)\n",
    "\n",
    "def DA_IRdn_balance(model):\n",
    "    return sum(model.isd[g] for g in model.G)+sum(model.sid[t] for t in model.R) +model.idnre>= max(model.RE[s] for s in model.S)-np.mean(np.array([model.RE[s] for s in model.S]))\n",
    "DAIRmodel.Con4DN = pyo.Constraint(rule=DA_IRdn_balance)\n",
    "\n",
    "def shortage_IR_dn(model,r):\n",
    "    return model.sid[r]<=max(0, model.RE[r+1]-max(np.mean(np.array([model.RE[s] for s in model.S])),model.RE[r]))\n",
    "DAIRmodel.Con4DNstep= pyo.Constraint(DAIRmodel.R,rule=shortage_IR_dn)\n",
    "\n",
    "\n",
    "\n",
    "def RE_dn(model):\n",
    "    return model.idnre<= model.rgDA\n",
    "DAIRmodel.Con5 = pyo.Constraint(expr=RE_dn)\n",
    "\n",
    "def RE_up(model):\n",
    "    return model.iupre +model.rgDA<= np.mean(np.array([model.RE[s] for s in model.S]))\n",
    "DAIRmodel.Con6 = pyo.Constraint(expr=RE_up)                                                                                                            \n",
    "                                                                                                            \n",
    "                                                                                               \n",
    "def RRUP(model,g):\n",
    "    return model.isu[g]<=model.RR[g]\n",
    "DAIRmodel.Con10UP = pyo.Constraint(DAIRmodel.G, rule=RRUP)  \n",
    "\n",
    "def RRDN(model,g):\n",
    "    return model.isd[g] <=model.RR[g]\n",
    "DAIRmodel.Con10DN = pyo.Constraint(DAIRmodel.G, rule=RRDN) \n",
    "\n",
    "def DA_capacity_cons(model, g):\n",
    "    return model.xDA[g] +model.isu[g]<= model.CAP[g]\n",
    "DAIRmodel.Con11 = pyo.Constraint(DAIRmodel.G, rule=DA_capacity_cons)\n",
    "\n",
    "def DA_down_cons(model, g):\n",
    "    return  model.isd[g] <= model.xDA[g]\n",
    "DAIRmodel.Con12 = pyo.Constraint(DAIRmodel.G, rule=DA_down_cons)\n",
    "\n",
    "def VB_UB(model):\n",
    "    return model.arb<=model.VUB\n",
    "DAIRmodel.Con13 = pyo.Constraint(rule=VB_UB)\n",
    "\n",
    "def VB_LB(model):\n",
    "    return model.arb>=model.VLB\n",
    "DAIRmodel.Con14 = pyo.Constraint(rule=VB_LB)\n",
    "\n",
    "\n",
    "DAIRmodel.dual = pyo.Suffix(direction=pyo.Suffix.IMPORT)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "45d423cd",
   "metadata": {},
   "source": [
    "Stochastic Program that is used to simulate the impact of DA decisions on RT"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "a71efc99",
   "metadata": {},
   "outputs": [],
   "source": [
    "RTsim = pyo.AbstractModel()\n",
    "#Sets\n",
    "RTsim.S = pyo.RangeSet(1, 5)\n",
    "RTsim.G = pyo.RangeSet(1, 5)\n",
    "\n",
    "#Parameters\n",
    "RTsim.VC = pyo.Param(RTsim.G)\n",
    "RTsim.CAP = pyo.Param(RTsim.G)\n",
    "RTsim.RR= pyo.Param(RTsim.G)\n",
    "\n",
    "RTsim.prob=pyo.Param(RTsim.S)\n",
    "RTsim.RE = pyo.Param(RTsim.S)\n",
    "RTsim.DEMAND=pyo.Param(within=pyo.NonNegativeIntegers)\n",
    "RTsim.D1=pyo.Param(within=pyo.NonNegativeIntegers)\n",
    "RTsim.D2=pyo.Param(within=pyo.NonNegativeIntegers)\n",
    "RTsim.xDA=pyo.Param(RTsim.G, domain=pyo.NonNegativeReals)\n",
    "RTsim.REDA=pyo.Param( domain=pyo.NonNegativeReals)\n",
    "RTsim.PEN=pyo.Param(within=pyo.NonNegativeIntegers)\n",
    "RTsim.PENDN=pyo.Param()\n",
    "RTsim.ARB=pyo.Param()\n",
    "RTsim.DAdr=pyo.Param()\n",
    "\n",
    "# Variables\n",
    "RTsim.xup = pyo.Var(RTsim.S, RTsim.G, domain=pyo.NonNegativeReals)\n",
    "RTsim.xdn = pyo.Var(RTsim.S, RTsim.G, domain=pyo.NonNegativeReals)\n",
    "RTsim.d= pyo.Var(RTsim.S, domain=pyo.NonNegativeReals)\n",
    "RTsim.rgup = pyo.Var(RTsim.S, domain=pyo.NonNegativeReals)\n",
    "RTsim.rgdn = pyo.Var(RTsim.S, domain=pyo.NonNegativeReals)\n",
    "RTsim.sdup = pyo.Var(RTsim.S, domain=pyo.NonNegativeReals)\n",
    "RTsim.sddn = pyo.Var(RTsim.S, domain=pyo.NonNegativeReals)\n",
    "\n",
    "#Objective\n",
    "def obj_expression(m):\n",
    "    return sum(m.prob[s]*m.VC[g]*(m.xup[s,g]-m.xdn[s,g]) for g in m.G for s in m.S)\\\n",
    "-m.PENDN*sum(m.prob[s]*m.sdup[s] for s in m.S)+m.PEN*sum(m.prob[s]*m.sddn[s] for s in m.S)\\\n",
    "+sum(m.prob[s]*(m.D1*(m.DAdr+m.d[s])+m.D2*(m.DAdr+m.d[s])*(m.DAdr+m.d[s])) for s in m.S)\\\n",
    "-sum(m.prob[s]*(m.D1*(m.DAdr)+m.D2*(m.DAdr)*(m.DAdr)) for s in m.S)\n",
    "RTsim.OBJ = pyo.Objective(rule=obj_expression)\n",
    "\n",
    "#Constraints\n",
    "def RT_energy_balance(model,s):\n",
    "    return sum(model.xup[s,g]-model.xdn[s,g] for g in model.G)+ model.rgup[s]-model.rgdn[s]+model.d[s] == model.ARB\n",
    "\n",
    "def RT_capacity_cons(model,s, g):\n",
    "    return model.xDA[g]+model.xup[s,g] <= model.CAP[g]\n",
    "\n",
    "def RT_capacity_min(model,s, g):\n",
    "    return model.xDA[g]-model.xdn[s,g] >=0\n",
    "\n",
    "def RT_RE_availability(model,s):\n",
    "    return model.rgup[s]-model.rgdn[s]+model.sdup[s]-model.sddn[s] ==model.RE[s]-model.REDA \n",
    "\n",
    "def RT_ramp_up(model,s,g):\n",
    "    return model.xup[s,g] <=model.RR[g]\n",
    "\n",
    "def RT_ramp_dn(model,s,g):\n",
    "    return model.xdn[s,g] <= model.RR[g]\n",
    "\n",
    "# Numbering of constraints follows the RT problem\n",
    "RTsim.Con3 = pyo.Constraint(RTsim.S, rule=RT_energy_balance)\n",
    "RTsim.Con4 = pyo.Constraint(RTsim.S, rule= RT_RE_availability)\n",
    "RTsim.Con5up = pyo.Constraint(RTsim.S,RTsim.G, rule=RT_ramp_up)\n",
    "RTsim.Con5dn = pyo.Constraint(RTsim.S, RTsim.G,rule=RT_ramp_dn)\n",
    "RTsim.Con6 = pyo.Constraint(RTsim.S,RTsim.G, rule=RT_capacity_cons)\n",
    "RTsim.Con7 = pyo.Constraint(RTsim.S,RTsim.G, rule=RT_capacity_min)\n",
    "\n",
    "#Record duals\n",
    "RTsim.dual = pyo.Suffix(direction=pyo.Suffix.IMPORT)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e82f21b2",
   "metadata": {},
   "source": [
    "Use CPLEX to solve the problem [quadratic objective function due to elastic demand]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "f69773f5",
   "metadata": {},
   "outputs": [],
   "source": [
    "opt = pyo.SolverFactory('cplex',executable='C:/Program Files/IBM/ILOG/CPLEX_Studio2211/cplex/bin/x64_win64/cplex')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "f403857a",
   "metadata": {},
   "outputs": [],
   "source": [
    "input_data=pd.read_csv('Input_files/IR_input_sectionV.csv',index_col=0)\n",
    "Summary_cost=pd.DataFrame()\n",
    "Summary_price_convergence=pd.DataFrame()\n",
    "Summary_margin_intermittency=pd.DataFrame()\n",
    "Summary_costD=pd.DataFrame()\n",
    "Summary_costCURT=pd.DataFrame()\n",
    "for l in input_data.columns:\n",
    "        dataDA=input_data[l].to_dict()\n",
    "        for k in range(0,len(input_data[l])):\n",
    "            dataDA[input_data.index[k]]=ast.literal_eval(input_data[l][input_data.index[k]])\n",
    "        dataDA = {None: dataDA}\n",
    "        i = DAIRmodel.create_instance(dataDA)\n",
    "        opt.solve(i)\n",
    "        #Values to pass from DA to RT module\n",
    "        Energy=pd.DataFrame()\n",
    "        E_grid_data = {(G): v.value for (G), v in i.xDA.items()}\n",
    "        Energy['en']= pd.DataFrame.from_dict(E_grid_data, orient=\"index\", columns=[\"en\"])\n",
    "        Energy['en']=round(Energy['en'],2)\n",
    "        temp=i.rgDA.value\n",
    "        Total=pd.DataFrame()\n",
    "        Total.at['cost','DA']=sum(np.multiply(np.array([i.VC[k] for k in i.G]),np.array([i.xDA[g].value for g in i.G])))\n",
    "        Total.at['price','DA']=i.dual[i.Con3]\n",
    "        E_grid_data = {(G): v.value for (G), v in i.isu.items()}\n",
    "        df = pd.DataFrame.from_dict(E_grid_data, orient=\"index\", columns=[\"isu\"])\n",
    "        E_grid_data = {(G): v.value for (G), v in i.isd.items()}\n",
    "        df ['isd']= pd.DataFrame.from_dict(E_grid_data, orient=\"index\", columns=[\"isd\"])\n",
    "        Prices=pd.DataFrame()\n",
    "        Prices.at[0,'up']=i.dual[i.Con4UP]\n",
    "        Prices.at[0,'down']=i.dual[i.Con4DN]\n",
    "        Gross_margins=pd.DataFrame()\n",
    "        Gross_margins[\"en\"]=np.multiply(i.dual[i.Con3]-np.array([i.VC[k] for k in i.G]), Energy['en'])\n",
    "        Gross_margins.at[0,\"en\"]=i.rgDA.value*i.dual[i.Con3]\n",
    "        Gross_margins[\"up\"]=i.dual[i.Con4UP]*df['isu']\n",
    "        Gross_margins[\"down\"]=i.dual[i.Con4DN]*df['isd']\n",
    "            \n",
    "        dataRT = {None:{\n",
    "            'RE': dataDA[None][\"RE\"],\n",
    "            'CAP': dataDA[None][\"CAP\"],\n",
    "            'VC': dataDA[None][\"VC\"],\n",
    "            'DEMAND': dataDA[None][\"DEMAND\"],\n",
    "            'D1': dataDA[None][\"D1\"],\n",
    "            'D2': dataDA[None][\"D2\"],\n",
    "            'prob':{1: 0.2, 2:0.2, 3:0.2, 4:0.2, 5:0.2},\n",
    "              'RR': dataDA[None][\"RR\"],\n",
    "            'xDA':Energy['en'],\n",
    "             'PEN':{None: dataDA[None][\"PEN\"][1]},\n",
    "             'PENDN':{None: dataDA[None][\"PEND\"][1]},\n",
    "             'REDA':{None: temp},\n",
    "             'ARB':{None: i.arb.value},\n",
    "             'DAdr':{None: i.d.value}\n",
    "            }}\n",
    "        iRT = RTsim.create_instance(dataRT)\n",
    "        opt.solve(iRT)\n",
    "        RTmargins=pd.DataFrame()\n",
    "        for s in iRT.S:  \n",
    "            RTmargins[s]=(np.multiply(iRT.dual[iRT.Con3[s]]-iRT.prob[s]*np.array([iRT.VC[k] for k in iRT.G]),np.array([iRT.xup[s,k].value-iRT.xdn[s,k].value  for k in i.G])))\n",
    "            Total.at['cost',s]=iRT.prob[s]*sum(np.multiply(np.array([i.VC[k] for k in i.G]),np.array([iRT.xup[s,g].value-iRT.xdn[s,g].value for g in iRT.G])))\n",
    "            Total.at['price',s]=iRT.dual[iRT.Con3[s]]\n",
    "            Total.at['unmet_demand',s]=iRT.prob[s]*(iRT.D1*(i.d.value+iRT.d[s].value)+iRT.D2*(iRT.d[s].value+i.d.value)*(i.d.value+iRT.d[s].value))\n",
    "            Total.at['curtail cost',s]=iRT.prob[s]*(iRT.PENDN*iRT.sdup[s].value)\n",
    "        price_convergence=Total.at['price','DA']-Total[Total.columns[1:]].loc['price'].sum()\n",
    "        Total_margin=pd.DataFrame()\n",
    "        for k in range(1,6):\n",
    "             Total_margin[k]=iRT.prob[k]*Gross_margins.sum(axis=1)[0:5].values+RTmargins[k]\n",
    "             Total_margin.at['RE',k]=(iRT.dual[iRT.Con3[k]]*(iRT.rgup[k].value-iRT.rgdn[k].value)) +iRT.prob[k]*temp*i.dual[i.Con3]\n",
    "        writer = pd.ExcelWriter('Test_output_files/IR_'+l+\"_v\"+date.today().strftime(\"%Y_%m_%d\")+'revision.xlsx')      \n",
    "        round(Total,2).to_excel(writer,'System_metrics')\n",
    "        round(RTmargins,2).to_excel(writer,'RTmargins')\n",
    "        round(Gross_margins,2).to_excel(writer,'Grossmargins')\n",
    "        round(df,2).to_excel(writer,'IR_supply_AWARDS')\n",
    "        round(Energy,2).to_excel(writer,'DA_energy')\n",
    "        round(Total_margin,2).to_excel(writer,'Totalmargins')\n",
    "        round(Prices,2).to_excel(writer,'Prices')\n",
    "        pd.DataFrame(dataDA).to_excel(writer,'input_DA')\n",
    "        pd.DataFrame(dataRT).to_excel(writer,'input_RT')\n",
    "        writer.close()\n",
    "        Summary_cost[l]=Total.loc['cost'].transpose()\n",
    "        Summary_costD[l]=Total.loc['unmet_demand'].transpose()\n",
    "        Summary_price_convergence[l]=Total.loc['price'].transpose()\n",
    "        Summary_margin_intermittency[l]=Total_margin.apply(lambda row: sum(row>0.01),axis=1)\n",
    "        Summary_costCURT[l]=Total.loc['curtail cost'].transpose()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "2c51fdd7",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>up</th>\n",
       "      <th>down</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>30.000001</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "          up  down\n",
       "0  30.000001   0.0"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Prices"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4252be05",
   "metadata": {},
   "outputs": [],
   "source": [
    "i.dual[i.Con4UP]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "186aba99",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.8.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}