SOPF without DR v2.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Packages"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pyomo.environ import *\n",
    "import numpy as np\n",
    "import random\n",
    "import pandas as pd\n",
    "from numpy import flatnonzero as find\n",
    "from pypower.api import case9, ext2int, bustypes, makeBdc, rundcopf, ppoption\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",
    "import matplotlib.pyplot as plt\n",
    "import pickle"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Power network parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Converged in 0.04 seconds\n",
      "Objective Function Value = 5216.03 $/hr\n",
      "================================================================================\n",
      "|     System Summary                                                           |\n",
      "================================================================================\n",
      "\n",
      "How many?                How much?              P (MW)            Q (MVAr)\n",
      "---------------------    -------------------  -------------  -----------------\n",
      "Buses              9     Total Gen Capacity     820.0           0.0 to 0.0\n",
      "Generators         3     On-line Capacity       820.0           0.0 to 0.0\n",
      "Committed Gens     3     Generation (actual)    315.0               0.0\n",
      "Loads              3     Load                   315.0               0.0\n",
      "  Fixed            3       Fixed                315.0               0.0\n",
      "  Dispatchable     0       Dispatchable           0.0 of 0.0        0.0\n",
      "Shunts             0     Shunt (inj)              0.0               0.0\n",
      "Branches           9     Losses (I^2 * Z)         0.00              0.00\n",
      "Transformers       0     Branch Charging (inj)     -                0.0\n",
      "Inter-ties         0     Total Inter-tie Flow     0.0               0.0\n",
      "Areas              1\n",
      "\n",
      "                          Minimum                      Maximum\n",
      "                 -------------------------  --------------------------------\n",
      "Voltage Magnitude   1.000 p.u. @ bus 1          1.000 p.u. @ bus 1   \n",
      "Voltage Angle      -5.43 deg   @ bus 9          6.04 deg   @ bus 2   \n",
      "Lambda P           24.04 $/MWh @ bus 2         24.04 $/MWh @ bus 9   \n",
      "Lambda Q            0.00 $/MWh @ bus 1          0.00 $/MWh @ bus 1   \n",
      "\n",
      "================================================================================\n",
      "|     Bus Data                                                                 |\n",
      "================================================================================\n",
      " Bus      Voltage          Generation             Load          Lambda($/MVA-hr)\n",
      "  #   Mag(pu) Ang(deg)   P (MW)   Q (MVAr)   P (MW)   Q (MVAr)     P        Q   \n",
      "----- ------- --------  --------  --------  --------  --------  -------  -------\n",
      "    1  1.000    0.000*    86.56      0.00       -         -      24.044     -\n",
      "    2  1.000    6.040    134.38      0.00       -         -      24.044     -\n",
      "    3  1.000    4.003     94.06      0.00       -         -      24.044     -\n",
      "    4  1.000   -2.857       -         -         -         -      24.044     -\n",
      "    5  1.000   -4.635       -         -       90.00      0.00    24.044     -\n",
      "    6  1.000    0.845       -         -         -         -      24.044     -\n",
      "    7  1.000   -1.338       -         -      100.00      0.00    24.044     -\n",
      "    8  1.000    1.228       -         -         -         -      24.044     -\n",
      "    9  1.000   -5.430       -         -      125.00      0.00    24.044     -\n",
      "                        --------  --------  --------  --------\n",
      "               Total:    315.00      0.00    315.00      0.00\n",
      "\n",
      "================================================================================\n",
      "|     Branch Data                                                              |\n",
      "================================================================================\n",
      "Brnch   From   To    From Bus Injection   To Bus Injection     Loss (I^2 * Z)  \n",
      "  #     Bus    Bus    P (MW)   Q (MVAr)   P (MW)   Q (MVAr)   P (MW)   Q (MVAr)\n",
      "-----  -----  -----  --------  --------  --------  --------  --------  --------\n",
      "   0      1      4     86.56      0.00    -86.56      0.00     0.000      0.00\n",
      "   1      4      5     33.74      0.00    -33.74      0.00     0.000      0.00\n",
      "   2      5      6    -56.26      0.00     56.26      0.00     0.000      0.00\n",
      "   3      3      6     94.06      0.00    -94.06      0.00     0.000      0.00\n",
      "   4      6      7     37.80      0.00    -37.80      0.00     0.000      0.00\n",
      "   5      7      8    -62.20      0.00     62.20      0.00     0.000      0.00\n",
      "   6      8      2   -134.38      0.00    134.38      0.00     0.000      0.00\n",
      "   7      8      9     72.17      0.00    -72.17      0.00     0.000      0.00\n",
      "   8      9      4    -52.83      0.00     52.83      0.00     0.000      0.00\n",
      "                                                             --------  --------\n",
      "                                                    Total:     0.000      0.00\n"
     ]
    }
   ],
   "source": [
    "# Run PyPower case to get initial condition\n",
    "ppc = case9()\n",
    "ppopt = ppoption(VERBOSE=0)\n",
    "pp_dcopf = rundcopf(ppc, ppopt)\n",
    "pp_bus = pp_dcopf[\"bus\"]\n",
    "pp_branch = pp_dcopf[\"branch\"]\n",
    "pp_gen = pp_dcopf[\"gen\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "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 = np.count_nonzero(pp_bus[:, PD]) # Number of loads\n",
    "NW = 1 # Number of wind farms\n",
    "\n",
    "# Map generator to bus\n",
    "CG = np.zeros((NB,NG))\n",
    "CG[gbus,range(3)] = 1\n",
    "\n",
    "# Map load to bus\n",
    "CL = np.zeros((NB,NL))\n",
    "load_bus = np.nonzero(pp_bus[:, PD])[0]\n",
    "for i in range(len(load_bus)):\n",
    "    CL[load_bus[i], i] = 1\n",
    "\n",
    "# Map wind farm to bus\n",
    "CW = np.zeros((NB,1))\n",
    "CW[5,0] = 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",
    "# Cost coefficient for wind power\n",
    "CW1 = 0\n",
    "CW0 = 0\n",
    "\n",
    "# Cost coefficients for generation reserve\n",
    "CRG_up = 10*np.ones((NG,1))*baseMVA\n",
    "CRG_down = 10*np.ones((NG,1))*baseMVA\n",
    "\n",
    "# Cost coefficients for load reserve\n",
    "CRL_up = 9.8*np.ones((NL,1))*baseMVA\n",
    "CRL_down = 9.8*np.ones((NL,1))*baseMVA\n",
    "\n",
    "# Uncontrollable load\n",
    "load = pp_bus[:,PD]\n",
    "PL_total = load[load > 0]\n",
    "PL = PL_total/baseMVA\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": "markdown",
   "metadata": {},
   "source": [
    "## Read temperature and wind data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Read temperature data\n",
    "weather_data = pd.read_csv(\"Weather_Data/NYC_Hourly_Normal_Temp.csv\")\n",
    "# Read wind power data\n",
    "PW_data = pd.read_csv(\"spring_wind.csv\",header=None).to_numpy()\n",
    "PW_data = PW_data/500*20/baseMVA\n",
    "# PW_data = PW_data[:, 0:4000]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "Tamb = (weather_data.iloc[:,8] - 32)/1.8\n",
    "Tamb_min = min(Tamb)\n",
    "Tamb_10PCTL = (weather_data.iloc[:,6] - 32)/1.8\n",
    "Tamb_90PCTL = (weather_data.iloc[:,7] - 32)/1.8\n",
    "std_Tamb_10PCTL = (Tamb - Tamb_10PCTL)/1.282 \n",
    "std_Tamb_90PCTL = (Tamb_90PCTL - Tamb)/1.282 \n",
    "std_Tamb = (std_Tamb_10PCTL + std_Tamb_90PCTL)/10"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Virtual storage model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "def VS(dt,Tamb,num_TCL,Tset,delta,R,C,P_el,COP,baseMVA):\n",
    "    Tdown = np.zeros(num_TCL)\n",
    "    Tup = np.zeros(num_TCL)\n",
    "    T_ON = np.zeros(num_TCL)\n",
    "    T_OFF = np.zeros(num_TCL)\n",
    "    D = np.zeros(num_TCL)\n",
    "    B = 0\n",
    "    Pdown = 0\n",
    "    Pup = 0\n",
    "    Sdown = 0\n",
    "    Sup = 0\n",
    "    for i in range(0,num_TCL):\n",
    "        Tdown[i] = Tset[i] - delta/2\n",
    "        Tup[i] = Tset[i] + delta/2\n",
    "        T_ON[i] = C[i]*R[i]*np.log((Tdown[i] - Tamb - R[i]*P_el[i]*COP[i])/(Tup[i] - Tamb - R[i]*P_el[i]*COP[i]))\n",
    "        if ((Tdown[i] - Tamb - R[i]*P_el[i]*COP[i])/(Tup[i] - Tamb - R[i]*P_el[i]*COP[i]) < 0):\n",
    "            print(Tdown[i] - Tamb - R[i]*P_el[i]*COP[i])\n",
    "            print(Tup[i] - Tamb - R[i]*P_el[i]*COP[i])\n",
    "        T_OFF[i] = C[i]*R[i]*np.log((Tup[i] - Tamb)/(Tdown[i] - Tamb))\n",
    "        D[i] = T_ON[i]/(T_ON[i] + T_OFF[i])\n",
    "        if ( Tamb < Tdown[i] and Tamb + R[i]*P_el[i]*COP[i] > Tup[i] ):\n",
    "            B = B + P_el[i]*D[i]\n",
    "            Pup = Pup + P_el[i]\n",
    "            Sup = Sup + dt*P_el[i]*(1 - D[i])\n",
    "    return [B/1000/baseMVA,Pdown/1000/baseMVA,Pup/1000/baseMVA,Sdown,Sup]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Optimization problem"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 82,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Optimization: Problem formulation\n",
    "def optimizer(NG,NL,Tamb_sample1,Tamb_sample2,PW_sample,num_TCL,Tset,delta,R,C,P_el,COP,S0,PL,\\\n",
    "    Pmax,Pmin,Lmax,Lmin,CG,CW,CL,CG0,CG1,CG2,CW0,CW1,CRG_up,CRG_down,CRL_up,CRL_down,nsample,PW_f):\n",
    "    model = ConcreteModel()\n",
    "    # Defining the variables\n",
    "    model.PG = Var(range(NG),within = NonNegativeReals, initialize=0)\n",
    "    model.PC = Var(range(NL),within = NonNegativeReals, initialize=0)\n",
    "    model.RG_down = Var(range(NG),within = NonNegativeReals, initialize=0)\n",
    "    model.RG_up = Var(range(NG),within = NonNegativeReals, initialize=0)\n",
    "    model.RG = Var(range(NG),range(nsample),within = Reals, initialize=0)\n",
    "    # model.RL_down = Var(range(NL),within = NonNegativeReals, initialize=0)\n",
    "    # model.RL_up = Var(range(NL),within = NonNegativeReals, initialize=0)\n",
    "    # model.RL = Var(range(NL),range(nsample),within = Reals, initialize=0)\n",
    "    model.dG_down = Var(range(NG),within = NonNegativeReals, initialize=0)\n",
    "    model.dG_up = Var(range(NG),within = NonNegativeReals, initialize=0)\n",
    "    # model.dL_down = Var(range(NL),within = NonNegativeReals, initialize=0)\n",
    "    # model.dL_up = Var(range(NL),within = NonNegativeReals, initialize=0)\n",
    "    model.Va = Var(range(NB), within=Reals, initialize=0, bounds=(-1,1))\n",
    "    model.s1 = Var(range(nsample), within=NonNegativeReals, initialize=0)\n",
    "    model.s2 = Var(range(nsample), within=NonNegativeReals, initialize=0)\n",
    "    model.s3 = Var(range(nsample), within=NonNegativeReals, initialize=0)\n",
    "    model.s4 = Var(range(nsample), within=NonNegativeReals, initialize=0)\n",
    "    model.s5 = Var(range(nsample), within=NonNegativeReals, initialize=0)\n",
    "    model.s6 = Var(range(nsample), within=NonNegativeReals, initialize=0)\n",
    "    # model.s7 = Var(range(nsample), within=NonNegativeReals, initialize=0)\n",
    "    # model.s8 = Var(range(nsample), within=NonNegativeReals, initialize=0)\n",
    "    model.d = Var(range(nsample), within=NonNegativeReals, initialize=0)\n",
    "    model.e = Var(range(nsample), within=NonNegativeReals, initialize=0)\n",
    "    model.dual = Suffix(direction=Suffix.IMPORT)\n",
    "\n",
    "    ## Defining the deterministic constraints\n",
    "    model.c = ConstraintList()\n",
    "\n",
    "    # Power mismatch distribution\n",
    "    model.c.add(sum(model.dG_down[i] for i in range(NG)) == 1)\n",
    "    model.c.add(sum(model.dG_up[i] for i in range(NG)) == 1)\n",
    "    \n",
    "    # Line flow limit constraints\n",
    "    for i in range(NBr):\n",
    "        model.c.add(sum(Bf[i,j]*model.Va[j] for j in range(NB)) <= Lmax[i])\n",
    "        model.c.add(sum(Bf[i,j]*model.Va[j] for j in range(NB)) >= Lmin[i])\n",
    "    \n",
    "    PB_all = np.zeros((NL,nsample))\n",
    "    P_mis_wind_all = np.zeros(nsample)\n",
    "    identified_VS_all = np.zeros((NL,nsample,5))\n",
    "    \n",
    "    # Defining the CC constraints\n",
    "    for i in range(nsample):\n",
    "\n",
    "        for j in range(NG):\n",
    "        # Generation capacity limit\n",
    "            model.c.add(model.PG[j] + model.RG[j,i] >= Pmin[j])\n",
    "            model.c.add(model.PG[j] + model.RG[j,i] <= Pmax[j])\n",
    "            model.c.add(model.PG[j] + model.RG_up[j] <= Pmax[j])\n",
    "            model.c.add(model.PG[j] - model.RG_down[j] >= Pmin[j])\n",
    "            # Reserve capacity limit\n",
    "            model.c.add(model.RG[j,i] >= -model.RG_down[j])\n",
    "            model.c.add(model.RG[j,i] <= model.RG_up[j])\n",
    "    \n",
    "        # DC Power flow constraint\n",
    "        for j in range(NB):\n",
    "            model.c.add(sum(CG[j,k]*model.PG[k] + CG[j,k]*model.RG[k,i] for k in range(NG)) + CW[j]*PW_sample[i]\\\n",
    "            - sum(CL[j,k]*PL[k] + CL[j,k]*model.PC[k] for k in range(NL)) >= sum(B[j,k]*model.Va[k] for k in range(NB)) - model.s1[i])\n",
    "            model.c.add(sum(CG[j,k]*model.PG[k] + CG[j,k]*model.RG[k,i] for k in range(NG)) + CW[j]*PW_sample[i]\\\n",
    "            - sum(CL[j,k]*PL[k] + CL[j,k]*model.PC[k] for k in range(NL)) <= sum(B[j,k]*model.Va[k] for k in range(NB)) + model.s2[i])\n",
    "            # model.c.add(sum(CG[j,k]*model.PG[k] + CG[j,k]*model.RG[k,i] for k in range(NG)) + CW[j]*PW_sample[i]\\\n",
    "            # - sum(CL[j,k]*PL[k] + CL[j,k]*model.PC[k] for k in range(NL)) >= sum(B[j,k]*model.Va[k] for k in range(NB)) - model.s1[i])\n",
    "            # model.c.add(sum(CG[j,k]*model.PG[k] + CG[j,k]*model.RG[k,i] for k in range(NG)) + CW[j]*PW_sample[i]\\\n",
    "            # - sum(CL[j,k]*PL[k] + CL[j,k]*model.PC[k] for k in range(NL)) <= sum(B[j,k]*model.Va[k] for k in range(NB)) + model.s2[i])\n",
    "        \n",
    "        P_mis_wind = PW_sample[i] - PW_f\n",
    "        P_mis_wind_all[i] = P_mis_wind\n",
    "\n",
    "        PB = np.zeros(NL)\n",
    "\n",
    "        # print(i)\n",
    "\n",
    "        # P_mis_load = 0\n",
    "        for j in range(NL):\n",
    "            identified_VS = VS(0.25,Tamb_sample1[i],num_TCL,Tset[:,j],delta,R[:,j],C[:,j],P_el[:,j],COP[:,j],baseMVA)\n",
    "            # print(\"VS0\", j, identified_VS[0])\n",
    "            identified_VS_all[j,i,:] = identified_VS\n",
    "            PB[j] = identified_VS[0]\n",
    "        \n",
    "        PB_all[:,i] = PB\n",
    "        \n",
    "        model.c.add(model.e[i] >= P_mis_wind - sum(model.PC[j] for j in range(NL)) + sum(PB[j] for j in range(NL)))\n",
    "        model.c.add(model.d[i] >= - P_mis_wind + sum(model.PC[j] for j in range(NL)) - sum(PB[j] for j in range(NL)))\n",
    "        \n",
    "        for j in range(NG):\n",
    "            model.c.add(model.RG[j,i] == model.dG_up[j]*model.d[i] - model.dG_down[j]*model.e[i])\n",
    "        \n",
    "\n",
    "        for j in range(NL):\n",
    "            # model.c.add(model.RL[j,i] == model.dL_up[j]*model.e[i] - model.dL_down[j]*model.d[i])\n",
    "            identified_VS1 = VS(0.25,Tamb_sample1[i],num_TCL,Tset[:,j],delta,R[:,j],C[:,j],P_el[:,j],COP[:,j],baseMVA)\n",
    "            # identified_VS2 = VS(0.25,Tamb_sample2[i],num_TCL,Tset[:,j],delta,R[:,j],C[:,j],P_el[:,j],COP[:,j],baseMVA)\n",
    "\n",
    "            model.c.add(model.PC[j] >= identified_VS1[1])\n",
    "            model.c.add(model.PC[j] <= identified_VS1[2])\n",
    "            # model.c.add(model.RL[j,i] >= -model.RL_down[j])\n",
    "            # model.c.add(model.RL[j,i] <= model.RL_up[j])\n",
    "            model.c.add(S0[j] + 0.25*(model.PC[j] - identified_VS1[0])*1000*baseMVA >= identified_VS1[3] - model.s3[i])\n",
    "            model.c.add(S0[j] + 0.25*(model.PC[j] - identified_VS1[0])*1000*baseMVA  <= identified_VS1[4] + model.s4[i])\n",
    "            model.c.add(S0[j] + 0.75*(model.PC[j] - identified_VS1[0])*1000*baseMVA  + 0.25*(model.PC[j]\\\n",
    "                 - identified_VS1[0])*1000*baseMVA  >= identified_VS1[3] - model.s5[i])\n",
    "            model.c.add(S0[j] + 0.75*(model.PC[j] - identified_VS1[0])*1000*baseMVA  + 0.25*(model.PC[j]\\\n",
    "                 - identified_VS1[0])*1000*baseMVA  <= identified_VS1[4] + model.s6[i])\n",
    "\n",
    "            # model.c.add(S0[j] + 0.75*(model.PC[j] - identified_VS2[0])*1000*baseMVA  + 0.25*(model.PC[j] + model.RL[j,i]\\\n",
    "            #      - identified_VS2[0])*1000*baseMVA  >= identified_VS2[3] - model.s7[i])\n",
    "            # model.c.add(S0[j] + 0.75*(model.PC[j] - identified_VS2[0])*1000*baseMVA  + 0.25*(model.PC[j] + model.RL[j,i]\\\n",
    "            #      - identified_VS2[0])*1000*baseMVA  <= identified_VS2[4] + model.s8[i])\n",
    "    \n",
    "    def cost(model, CG0, CG1, CG2, CRG_down, CRG_up, nsample, CRL_down, CRL_up):\n",
    "\n",
    "        cost = sum(CG0[i] for i in range(NG))\\\n",
    "                + sum(CG1[i]*model.PG[i] for i in range(NG))\\\n",
    "                + sum(CG2[i]*model.PG[i]**2 for i in range(NG))\\\n",
    "                + sum(CRG_up[i]*model.RG_up[i] for i in range(NG))\\\n",
    "                + sum(CRG_down[i]*model.RG_down[i] for i in range(NG))\\\n",
    "                + 1000*sum(model.s1[i] + model.s2[i] for i in range(nsample))\\\n",
    "                + 1e10*sum(model.s3[i] + model.s4[i] + model.s5[i] + model.s6[i] for i in range(nsample))\\\n",
    "                + 1000*sum(model.d[i] + model.e[i] for i in range(nsample))\n",
    "                # + sum(CRL_up[i]*model.RL_up[i] for i in range(NL))\\\n",
    "                # + sum(CRL_down[i]*model.RL_down[i] for i in range(NL))\\\n",
    "                # + 1e10*sum(model.s7[i] + model.s8[i] for i in range(nsample))\\\n",
    "                \n",
    "                \n",
    "        return cost\n",
    "\n",
    "    model.obj = Objective(expr=cost(model, CG0, CG1, CG2, CRG_down, CRG_up, nsample, CRL_down, CRL_up), sense=minimize)\n",
    "\n",
    "    \n",
    "    return model, PB_all, P_mis_wind_all, identified_VS_all"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 84,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of samples 1:  647\n",
      "Number of samples 2:  453\n"
     ]
    }
   ],
   "source": [
    "# Parameters related to the demand side\n",
    "num_TCL = 5000\n",
    "Tset = np.zeros((num_TCL,NL))\n",
    "R = np.zeros((num_TCL,NL))\n",
    "C = np.zeros((num_TCL,NL))\n",
    "COP = np.zeros((num_TCL,NL))\n",
    "P_el = np.zeros((num_TCL,NL))\n",
    "Tset_set = [21,21.5,22,22.5,23]\n",
    "for i in range(num_TCL): \n",
    "    for j in range(NL):\n",
    "        s = random.sample(Tset_set,1)\n",
    "        Tset[i,j] = s[0] \n",
    "delta = 1\n",
    "R_nom = 2\n",
    "R_low = 0.8*R_nom\n",
    "R_high = 1.2*R_nom\n",
    "C_nom = 2\n",
    "C_low = 0.8*C_nom\n",
    "C_high = 1.2*C_nom\n",
    "COP_nom = 2.5\n",
    "COP_low = 0.8*COP_nom\n",
    "COP_high = 1.2*COP_nom\n",
    "\n",
    "np.random.seed(1)\n",
    "\n",
    "for i in range(NL): \n",
    "    R[:,i] = np.random.uniform(R_low,R_high,num_TCL)\n",
    "    C[:,i] = np.random.uniform(C_low,C_high,num_TCL)\n",
    "    COP[:,i] = np.random.uniform(COP_low,COP_high,num_TCL)\n",
    "    for j in range(num_TCL):\n",
    "        P_el[j,i] = np.ceil((Tset[j,i] + delta/2 - Tamb_min)/(R[j,i]*COP[j,i])) + 3\n",
    "\n",
    "# Sampling\n",
    "eps = 0.1\n",
    "beta = 0.1\n",
    "num_decision_variables = 5*NG + 5*NL\n",
    "num_samples1 = np.ceil((2/eps)*(np.log(1/beta) + num_decision_variables)).astype(int)\n",
    "num_samples2 = np.ceil((num_decision_variables + 1 + np.log(1/beta) + np.sqrt(2*(num_decision_variables+1)*np.log(1/beta)))/eps).astype(int)\n",
    "\n",
    "print(\"Number of samples 1: \", num_samples1)\n",
    "print(\"Number of samples 2: \", num_samples2)\n",
    "\n",
    "P_B = np.zeros((NL,1))\n",
    "Pdown = np.zeros((NL,1))\n",
    "Pup = np.zeros((NL,1))\n",
    "Sdown = np.zeros((NL,1))\n",
    "Sup = np.zeros((NL,1))\n",
    "S0 = np.zeros((NL,1))\n",
    "\n",
    "for j in range(NL):\n",
    "    [P_B[j],Pdown[j],Pup[j],Sdown[j],Sup[j]] = \\\n",
    "        VS(0.25,Tamb[0],num_TCL,Tset[:,j],delta,R[:,j],C[:,j],P_el[:,j],COP[:,j],baseMVA)\n",
    "    S0[j] = (Sdown[j] + Sup[j])/2    \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Single period problem"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 85,
   "metadata": {},
   "outputs": [],
   "source": [
    "# # Single period\n",
    "\n",
    "# # Define the model\n",
    "# nsample= 50\n",
    "# np.random.seed(1)\n",
    "# Tamb_sample = np.random.normal(Tamb[0],std_Tamb[0],nsample)\n",
    "# PW_sample = np.random.choice(PW_data[0,:], size=nsample)\n",
    "# PW_f = np.median(PW_data[0,:])\n",
    "# # PW_f = np.median(PW_data)\n",
    "# model, PB_all, P_mis_wind, identified_VS_all = optimizer(NG,NL,\n",
    "# Tamb_sample,PW_sample,num_TCL,Tset,delta,R,C,P_el,COP,S0,PL,\\\n",
    "#     Pmax,Pmin,Lmax,Lmin,CG,CW,CL,CG0,CG1,CG2,CW0,CW1,CRG_up,CRG_down,CRL_up,CRL_down,nsample,PW_f)\n",
    "\n",
    "# # Solve the problem\n",
    "# solver = SolverFactory(\"bonmin\")\n",
    "# solver.solve(model)\n",
    "\n",
    "# model.display()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Multi-period problem"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 86,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Read S0 data from OPF with DR\n",
    "# file = open(\"Results/S0all_453\",\"rb\")\n",
    "# S0_all_DR = pickle.load(file)\n",
    "# file.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 107,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Hour  0\n",
      "S0:  [2794.48257109 2798.6107005  2790.34622915]\n",
      "S up 0: 5588.965142\n",
      "S up 1: 5597.221401\n",
      "S up 2: 5580.692458\n",
      "Elapsed:  119.62601232528687\n",
      "Hour  1\n",
      "S0:  [2993.72054089 2999.19363934 3000.11664782]\n",
      "S up 0: 5546.203297\n",
      "S up 1: 5554.381374\n",
      "S up 2: 5537.839264\n",
      "Elapsed:  158.44917559623718\n",
      "Hour  2\n",
      "S0:  [3094.78089704 3113.02737514 3100.88079848]\n",
      "S up 0: 5503.441452\n",
      "S up 1: 5511.541347\n",
      "S up 2: 5494.986069\n",
      "Elapsed:  135.3354456424713\n",
      "Hour  3\n",
      "S0:  [3123.59019831 3126.96765805 3110.50478846]\n",
      "S up 0: 5460.679600\n",
      "S up 1: 5468.701315\n",
      "S up 2: 5452.132865\n",
      "Elapsed:  119.95629453659058\n",
      "Hour  4\n",
      "S0:  [3034.36016865 3037.94551439 3020.63056531]\n",
      "S up 0: 5432.171692\n",
      "S up 1: 5440.141287\n",
      "S up 2: 5423.564056\n",
      "Elapsed:  119.33100700378418\n",
      "Hour  5\n",
      "S0:  [3005.85400579 3009.38726516 2992.06347762]\n",
      "S up 0: 5417.917735\n",
      "S up 1: 5425.861270\n",
      "S up 2: 5409.279648\n",
      "Elapsed:  187.31128478050232\n",
      "Hour  6\n",
      "S0:  [3001.5076641  3005.01488006 2987.68690256]\n",
      "S up 0: 5417.917735\n",
      "S up 1: 5425.861270\n",
      "S up 2: 5409.279648\n",
      "Elapsed:  239.21683406829834\n",
      "Hour  7\n",
      "S0:  [3001.50781059 3005.01502655 2987.68693819]\n",
      "S up 0: 5432.171692\n",
      "S up 1: 5440.141287\n",
      "S up 2: 5423.564056\n",
      "Elapsed:  123.4524757862091\n",
      "Hour  8\n",
      "S0:  [3005.85350955 3009.38676893 2992.06298149]\n",
      "S up 0: 5531.949349\n",
      "S up 1: 5540.101365\n",
      "S up 2: 5523.554866\n",
      "Elapsed:  134.9678089618683\n",
      "Hour  9\n",
      "S0:  [3093.79678525 3099.45304994 3100.40665211]\n",
      "S up 0: 5645.980947\n",
      "S up 1: 5654.341450\n",
      "S up 2: 5637.830060\n",
      "Elapsed:  126.4932770729065\n",
      "Hour  10\n",
      "S0:  [3118.96907449 3123.78479662 3114.19931909]\n",
      "S up 0: 5802.774597\n",
      "S up 1: 5811.421777\n",
      "S up 2: 5794.958645\n",
      "Elapsed:  122.85708546638489\n",
      "Hour  11\n",
      "S0:  [3164.34779044 3169.4880973  3161.28203668]\n",
      "S up 0: 5945.314732\n",
      "S up 1: 5954.222537\n",
      "S up 2: 5937.803257\n",
      "Elapsed:  119.2367651462555\n",
      "Hour  12\n",
      "S0:  [3237.69962966 3242.96446971 3234.04884854]\n",
      "S up 0: 6059.347325\n",
      "S up 1: 6068.463636\n",
      "S up 2: 6052.079428\n",
      "Elapsed:  118.83356618881226\n",
      "Hour  13\n",
      "S0:  [3331.64479982 3337.13818276 3327.12712187]\n",
      "S up 0: 6144.872138\n",
      "S up 1: 6154.144835\n",
      "S up 2: 6137.786923\n",
      "Elapsed:  118.51557302474976\n",
      "Hour  14\n",
      "S0:  [3464.83696828 3470.604979   3459.80672945]\n",
      "S up 0: 6173.380491\n",
      "S up 1: 6182.705318\n",
      "S up 2: 6166.356170\n",
      "Elapsed:  121.53012943267822\n",
      "Hour  15\n",
      "S0:  [3518.72216676 3524.36219452 3512.18098494]\n",
      "S up 0: 6130.617978\n",
      "S up 1: 6139.864610\n",
      "S up 2: 6123.502316\n",
      "Elapsed:  121.09373307228088\n",
      "Hour  16\n",
      "S0:  [3686.14749312 3691.16359914 3674.6614379 ]\n",
      "S up 0: 6030.839128\n",
      "S up 1: 6039.903312\n",
      "S up 2: 6023.510337\n",
      "Elapsed:  119.60745429992676\n",
      "Hour  17\n",
      "S0:  [3725.76889994 3730.61772962 3713.51651089]\n",
      "S up 0: 5959.568779\n",
      "S up 1: 5968.502647\n",
      "S up 2: 5952.087752\n",
      "Elapsed:  119.1927638053894\n",
      "Hour  18\n",
      "S0:  [3634.2908851  3638.97254411 3621.8430859 ]\n",
      "S up 0: 5916.806658\n",
      "S up 1: 5925.662337\n",
      "S up 2: 5909.234288\n",
      "Elapsed:  119.03307437896729\n",
      "Hour  19\n",
      "S0:  [3591.53576119 3596.13930541 3578.9966253 ]\n",
      "S up 0: 5845.536580\n",
      "S up 1: 5854.261946\n",
      "S up 2: 5837.811972\n",
      "Elapsed:  120.73143196105957\n",
      "Hour  20\n",
      "S0:  [3429.28533987 3433.59227703 3416.39961192]\n",
      "S up 0: 5788.520612\n",
      "S up 1: 5797.141729\n",
      "S up 2: 5780.674211\n",
      "Elapsed:  119.57761859893799\n",
      "Hour  21\n",
      "S0:  [3261.06648821 3265.06588088 3247.82141957]\n",
      "S up 0: 5717.250743\n",
      "S up 1: 5725.741552\n",
      "S up 2: 5709.252101\n",
      "Elapsed:  119.36332607269287\n",
      "Hour  22\n",
      "S0:  [3301.01204546 3305.08466899 3287.85242238]\n",
      "S up 0: 5660.234902\n",
      "S up 1: 5668.621465\n",
      "S up 2: 5652.114464\n",
      "Elapsed:  123.56386160850525\n",
      "Hour  23\n",
      "S0:  [3132.7980007  3136.56304566 3119.27886521]\n",
      "Elapsed:  129.16775727272034\n"
     ]
    }
   ],
   "source": [
    "# Multi-period \n",
    "import time\n",
    "\n",
    "# S0_all = S0_all_DR\n",
    "S0_all = np.zeros((24,NL))\n",
    "R_PC = np.zeros((24,NL))\n",
    "models = []\n",
    "Tamb_samples1 = []\n",
    "Tamb_samples2 = []\n",
    "PW_samples = []\n",
    "PW_fs = []\n",
    "PB_alls = []\n",
    "P_mis_wind_alls = []\n",
    "identified_VS_alls = []\n",
    "P_B_test = []\n",
    "Pdown_test = []\n",
    "Pup_test = []\n",
    "Sdown_test = []\n",
    "Sup_test = []\n",
    "\n",
    "# Initialize at hour 0\n",
    "for j in range(NL):\n",
    "    [P_B[j],Pdown[j],Pup[j],Sdown[j],Sup[j]] = \\\n",
    "        VS(0.25,Tamb[0],num_TCL,Tset[:,j],delta,R[:,j],C[:,j],P_el[:,j],COP[:,j],baseMVA)\n",
    "    S0_all[0,j] = (Sdown[j] + Sup[j])/2  \n",
    "\n",
    "nsample= 453\n",
    "\n",
    "for h in range(24):\n",
    "\n",
    "    print(\"Hour \", h)\n",
    "    t_start = time.time()\n",
    "\n",
    "    print(\"S0: \", S0_all[h])\n",
    "\n",
    "    np.random.seed(1)\n",
    "    Tamb_sample1 = np.random.normal(Tamb[h],std_Tamb[h],nsample)\n",
    "    Tamb_sample2 = np.random.normal(Tamb[h+1],std_Tamb[h+1],nsample)\n",
    "    PW_sample = np.random.choice(PW_data[h,:], size=nsample)\n",
    "    PW_f = np.median(PW_data[h, :])\n",
    "\n",
    "    # Define the problem\n",
    "    model, PB_all, P_mis_wind_all, identified_VS_all = optimizer(NG,NL,Tamb_sample1,Tamb_sample2,PW_sample,num_TCL,Tset,delta,R,C,P_el,COP,S0_all[h,:],PL,\\\n",
    "        Pmax,Pmin,Lmax,Lmin,CG,CW,CL,CG0,CG1,CG2,CW0,CW1,CRG_up,CRG_down,CRL_up,CRL_down,nsample,PW_f)\n",
    "\n",
    "    # Solve the problem\n",
    "    solver = SolverFactory(\"bonmin\")\n",
    "    solver.solve(model)\n",
    "    r_pc = model.PC[:]()\n",
    "    R_PC[h,:] = r_pc\n",
    "\n",
    "    models.append(model)\n",
    "    Tamb_samples1.append(Tamb_sample1)\n",
    "    PW_samples.append(PW_sample)\n",
    "    PW_fs.append(PW_f)\n",
    "    PB_alls.append(PB_all)\n",
    "    P_mis_wind_alls.append(P_mis_wind_all)\n",
    "    identified_VS_alls.append(identified_VS_all)\n",
    "\n",
    "    if h != 23:\n",
    "        for j in range(NL):\n",
    "            [P_B[j],Pdown[j],Pup[j],Sdown[j],Sup[j]] = \\\n",
    "                VS(0.25,Tamb[h],num_TCL,Tset[:,j],delta,R[:,j],C[:,j],P_el[:,j],COP[:,j],baseMVA)\n",
    "            S0_all[h+1,j] = S0_all[h,j] + r_pc[j]*1000*baseMVA - P_B[j]*1000*baseMVA\n",
    "        \n",
    "            print(\"S up %d: %f\" %(j, Sup[j]))\n",
    "\n",
    "    t_end = time.time()\n",
    "    print(\"Elapsed: \", t_end - t_start)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 112,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "PC : Size=3, Index=PC_index\n",
      "    Key : Lower : Value               : Upper : Fixed : Stale : Domain\n",
      "      0 :     0 : 0.21552930057504951 :  None : False : False : NonNegativeReals\n",
      "      1 :     0 :  0.2156447374589551 :  None : False : False : NonNegativeReals\n",
      "      2 :     0 :  0.2158650084265998 :  None : False : False : NonNegativeReals\n"
     ]
    }
   ],
   "source": [
    "model.PC.pprint()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 113,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "PG : Size=3, Index=PG_index\n",
      "    Key : Lower : Value              : Upper : Fixed : Stale : Domain\n",
      "      0 :     0 : 1.0714170284240991 :  None : False : False : NonNegativeReals\n",
      "      1 :     0 : 1.5467279062987633 :  None : False : False : NonNegativeReals\n",
      "      2 :     0 : 1.1262174757106598 :  None : False : False : NonNegativeReals\n"
     ]
    }
   ],
   "source": [
    "models[2].PG.pprint()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 114,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "RG_up : Size=3, Index=RG_up_index\n",
      "    Key : Lower : Value                : Upper : Fixed : Stale : Domain\n",
      "      0 :     0 : 0.034175305559122404 :  None : False : False : NonNegativeReals\n",
      "      1 :     0 :  0.03417530672466326 :  None : False : False : NonNegativeReals\n",
      "      2 :     0 :  0.03417530507559341 :  None : False : False : NonNegativeReals\n",
      "RG_down : Size=3, Index=RG_down_index\n",
      "    Key : Lower : Value                : Upper : Fixed : Stale : Domain\n",
      "      0 :     0 :  0.05764701072062478 :  None : False : False : NonNegativeReals\n",
      "      1 :     0 : 0.057647004903107125 :  None : False : False : NonNegativeReals\n",
      "      2 :     0 : 0.057647010348583425 :  None : False : False : NonNegativeReals\n"
     ]
    }
   ],
   "source": [
    "model.RG_up.pprint()\n",
    "model.RG_down.pprint()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Save the results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 111,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pickle\n",
    "\n",
    "case_name = \"case9_without_DR_v6\"\n",
    "\n",
    "# Write models\n",
    "file = open(\"Results/models_%s_%d\" %(case_name,nsample), \"wb\")\n",
    "pickle.dump(models, file)\n",
    "file.close()\n",
    "\n",
    "# Write S0_all\n",
    "file = open(\"Results/S0all_%s_%d\" %(case_name,nsample),\"wb\")\n",
    "pickle.dump(S0_all, file)\n",
    "file.close()\n",
    "\n",
    "# Write R_PC\n",
    "file = open(\"Results/R_PC_%s_%d\" %(case_name,nsample), \"wb\")\n",
    "pickle.dump(R_PC, file)\n",
    "file.close()\n",
    "\n",
    "# Write Tamb_samples1\n",
    "file = open(\"Results/Tamb_samples1_%s_%d\" %(case_name,nsample), \"wb\")\n",
    "pickle.dump(Tamb_samples1, file)\n",
    "file.close()\n",
    "\n",
    "# Write Tamb_samples2\n",
    "file = open(\"Results/Tamb_samples2_%s_%d\" %(case_name,nsample), \"wb\")\n",
    "pickle.dump(Tamb_samples2, file)\n",
    "file.close()\n",
    "\n",
    "# Write PW_samples\n",
    "file = open(\"Results/PW_samples_%s_%d\" %(case_name,nsample), \"wb\")\n",
    "pickle.dump(PW_samples, file)\n",
    "file.close()\n",
    "\n",
    "# Write PW_fs\n",
    "file = open(\"Results/PW_fs_%s_%d\" %(case_name,nsample), \"wb\")\n",
    "pickle.dump(PW_fs, file)\n",
    "file.close()\n",
    "\n",
    "# Write PB\n",
    "file = open(\"Results/PB_alls_%s_%d\" %(case_name,nsample), \"wb\")\n",
    "pickle.dump(PB_alls, file)\n",
    "file.close()\n",
    "\n",
    "# Write identified_VS_alls\n",
    "file = open(\"Results/VS_alls_%s_%d\" %(case_name,nsample), \"wb\")\n",
    "pickle.dump(identified_VS_alls, file)\n",
    "file.close()\n",
    "\n",
    "# Write Pmis_winds\n",
    "file = open(\"Results/P_mis_winds_%s_%d\" %(case_name,nsample), \"wb\")\n",
    "pickle.dump(P_mis_wind_alls, file)\n",
    "file.close()"
   ]
  }
 ],
 "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.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}