SOPF with DR.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.17 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": 3,
   "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(NG)] = 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 = 100\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",
    "CRL_up = np.array([20, 5, 9.8])*baseMVA\n",
    "CRL_down = np.array([20, 5, 9.8])*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": 4,
   "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": 5,
   "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": 6,
   "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": 7,
   "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)) + sum(model.dL_up[j] for j in range(NL)) == 1)\n",
    "    model.c.add(sum(model.dG_up[i] for i in range(NG)) + sum(model.dL_down[j] for j in range(NL)) == 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(9)) <= Lmax[i])\n",
    "        model.c.add(sum(Bf[i,j]*model.Va[j] for j in range(9)) >= 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] + CL[j,k]*model.RL[k,i] 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] + CL[j,k]*model.RL[k,i] 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] + model.RL[j,i] >= identified_VS1[1])\n",
    "            model.c.add(model.PC[j] + model.RL[j,i] <= 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] + model.RL[j,i] - identified_VS1[0])*1000*baseMVA >= identified_VS1[3] - model.s3[i])\n",
    "            model.c.add(S0[j] + 0.25*(model.PC[j] + model.RL[j,i] - 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] + model.RL[j,i]\\\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] + model.RL[j,i]\\\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": 8,
   "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    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Single period problem"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "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],num_samples)\n",
    "# PW_sample = np.random.choice(PW_data[0,:], size=num_samples)\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": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Hour  0\n",
      "S0:  [5590.96647353 5589.30897543 5594.15000414]\n",
      "S up 0: 11181.932947\n",
      "S up 1: 11178.617951\n",
      "S up 2: 11188.300008\n",
      "Elapsed:  468.2229142189026\n",
      "Hour  1\n",
      "S0:  [5646.71791292 5645.40562951 5649.49903462]\n",
      "S up 0: 11096.139237\n",
      "S up 1: 11092.839451\n",
      "S up 2: 11102.482446\n",
      "Elapsed:  458.6143114566803\n",
      "Hour  2\n",
      "S0:  [5678.89478962 5677.66627196 5681.12596915]\n",
      "S up 0: 11010.345524\n",
      "S up 1: 11007.060950\n",
      "S up 2: 11016.664882\n",
      "Elapsed:  426.34666299819946\n",
      "Hour  3\n",
      "S0:  [5746.60121099 5745.41817516 5748.43761086]\n",
      "S up 0: 10924.551796\n",
      "S up 1: 10921.282434\n",
      "S up 2: 10930.847303\n",
      "Elapsed:  430.1574065685272\n",
      "Hour  4\n",
      "S0:  [5819.67278075 5818.06227097 5823.25010438]\n",
      "S up 0: 10867.355961\n",
      "S up 1: 10864.096743\n",
      "S up 2: 10873.635570\n",
      "Elapsed:  424.1109175682068\n",
      "Hour  5\n",
      "S0:  [5943.93289105 5942.21692084 5947.6332074 ]\n",
      "S up 0: 10838.758038\n",
      "S up 1: 10835.503891\n",
      "S up 2: 10845.029697\n",
      "Elapsed:  429.9609866142273\n",
      "Hour  6\n",
      "S0:  [6119.98350757 6118.11394116 6124.00909461]\n",
      "S up 0: 10838.758038\n",
      "S up 1: 10835.503891\n",
      "S up 2: 10845.029697\n",
      "Elapsed:  481.5819787979126\n",
      "Hour  7\n",
      "S0:  [6431.38657233 6581.53852016 6099.10531958]\n",
      "S up 0: 10867.355961\n",
      "S up 1: 10864.096743\n",
      "S up 2: 10873.635570\n",
      "Elapsed:  505.3042585849762\n",
      "Hour  8\n",
      "S0:  [6549.11354183 6536.3287039  6578.98842335]\n",
      "S up 0: 11067.541334\n",
      "S up 1: 11064.246618\n",
      "S up 2: 11073.876592\n",
      "Elapsed:  488.23130464553833\n",
      "Hour  9\n",
      "S0:  [6744.17344534 6742.01207476 6748.57649428]\n",
      "S up 0: 11296.324580\n",
      "S up 1: 11292.989305\n",
      "S up 2: 11302.723445\n",
      "Elapsed:  504.7632791996002\n",
      "Hour  10\n",
      "S0:  [6841.5102284  6837.80046377 6844.18786521]\n",
      "S up 0: 11610.901923\n",
      "S up 1: 11607.510897\n",
      "S up 2: 11617.388262\n",
      "Elapsed:  463.9126718044281\n",
      "Hour  11\n",
      "S0:  [6968.98111825 6966.86213395 6973.89272335]\n",
      "S up 0: 11896.882210\n",
      "S up 1: 11893.440518\n",
      "S up 2: 11903.448085\n",
      "Elapsed:  467.99730825424194\n",
      "Hour  12\n",
      "S0:  [7153.9040838  7151.70950033 7158.73536781]\n",
      "S up 0: 12125.667396\n",
      "S up 1: 12122.185182\n",
      "S up 2: 12132.296910\n",
      "Elapsed:  455.97529196739197\n",
      "Hour  13\n",
      "S0:  [7310.13232883 7307.95186437 7315.09214726]\n",
      "S up 0: 12297.257017\n",
      "S up 1: 12293.744419\n",
      "S up 2: 12303.934268\n",
      "Elapsed:  495.4715700149536\n",
      "Hour  14\n",
      "S0:  [7558.23594562 7555.94920704 7563.35695947]\n",
      "S up 0: 12354.453720\n",
      "S up 1: 12350.930994\n",
      "S up 2: 12361.146885\n",
      "Elapsed:  497.49036145210266\n",
      "Hour  15\n",
      "S0:  [7733.3952608  7731.00984555 7738.65052349]\n",
      "S up 0: 12268.658698\n",
      "S up 1: 12265.151163\n",
      "S up 2: 12275.327992\n",
      "Elapsed:  520.5958242416382\n",
      "Hour  16\n",
      "S0:  [7973.76160594 7971.22495528 7979.19268898]\n",
      "S up 0: 12068.471003\n",
      "S up 1: 12064.998918\n",
      "S up 2: 12075.084607\n",
      "Elapsed:  457.4467580318451\n",
      "Hour  17\n",
      "S0:  [8200.86170631 8198.20832244 8206.1886169 ]\n",
      "S up 0: 11925.480305\n",
      "S up 1: 11922.033547\n",
      "S up 2: 11932.054135\n",
      "Elapsed:  574.4669396877289\n",
      "Hour  18\n",
      "S0:  [8338.54675852 8335.92637012 8343.82319342]\n",
      "S up 0: 11839.686060\n",
      "S up 1: 11836.254500\n",
      "S up 2: 11846.236027\n",
      "Elapsed:  469.3959469795227\n",
      "Hour  19\n",
      "S0:  [8408.52555898 8405.92126757 8413.77889459]\n",
      "S up 0: 11696.695896\n",
      "S up 1: 11693.289669\n",
      "S up 2: 11703.206095\n",
      "Elapsed:  479.4063038825989\n",
      "Hour  20\n",
      "S0:  [8434.56364026 8431.88926579 8439.72642578]\n",
      "S up 0: 11582.303950\n",
      "S up 1: 11578.917992\n",
      "S up 2: 11588.782336\n",
      "Elapsed:  468.32554721832275\n",
      "Hour  21\n",
      "S0:  [8361.65683417 8359.00887749 8366.73334308]\n",
      "S up 0: 11439.314195\n",
      "S up 1: 11435.953576\n",
      "S up 2: 11445.752819\n",
      "Elapsed:  477.5917112827301\n",
      "Hour  22\n",
      "S0:  [8267.99735226 8265.34212319 8273.09202666]\n",
      "S up 0: 11324.922494\n",
      "S up 1: 11321.582151\n",
      "S up 2: 11331.329312\n",
      "Elapsed:  449.4817180633545\n",
      "Hour  23\n",
      "S0:  [8063.66620395 8061.06313853 8068.67103818]\n",
      "Elapsed:  483.6703941822052\n"
     ]
    }
   ],
   "source": [
    "# Multi-period \n",
    "import time\n",
    "\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= num_samples2\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",
    "\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",
    "        P_B_test.append(P_B)\n",
    "        Pdown_test.append(Pdown)\n",
    "        Pup_test.append(Pup)\n",
    "        Sdown_test.append(Sdown)\n",
    "        Sup_test.append(Sup)\n",
    "\n",
    "    t_end = time.time()\n",
    "    print(\"Elapsed: \", t_end - t_start)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Save the results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "case_name = \"case9_with_DR_double_load\"\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
}