SOPF with DR_case39.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, rundcpf, ppoption, case39, rundcopf\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"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Converged in 0.05 seconds\n",
      "Objective Function Value = 41263.94 $/hr\n",
      "================================================================================\n",
      "|     System Summary                                                           |\n",
      "================================================================================\n",
      "\n",
      "How many?                How much?              P (MW)            Q (MVAr)\n",
      "---------------------    -------------------  -------------  -----------------\n",
      "Buses             39     Total Gen Capacity    7367.0           0.0 to 0.0\n",
      "Generators        10     On-line Capacity      7367.0           0.0 to 0.0\n",
      "Committed Gens    10     Generation (actual)   6254.2               0.0\n",
      "Loads             21     Load                  6254.2               0.0\n",
      "  Fixed           21       Fixed               6254.2               0.0\n",
      "  Dispatchable     0       Dispatchable           0.0 of 0.0        0.0\n",
      "Shunts             0     Shunt (inj)              0.0               0.0\n",
      "Branches          46     Losses (I^2 * Z)         0.00              0.00\n",
      "Transformers      12     Branch Charging (inj)     -                0.0\n",
      "Inter-ties         6     Total Inter-tie Flow   836.0               0.0\n",
      "Areas              3\n",
      "\n",
      "                          Minimum                      Maximum\n",
      "                 -------------------------  --------------------------------\n",
      "Voltage Magnitude   0.982 p.u. @ bus 31         1.064 p.u. @ bus 36  \n",
      "Voltage Angle     -18.84 deg   @ bus 39         9.01 deg   @ bus 36  \n",
      "Lambda P           13.52 $/MWh @ bus 30        13.52 $/MWh @ bus 8   \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.039  -14.747       -         -       97.60      0.00    13.517     -\n",
      "    2  1.048   -5.726       -         -         -         -      13.517     -\n",
      "    3  1.031   -9.626       -         -      322.00      0.00    13.517     -\n",
      "    4  1.004  -11.250       -         -      500.00      0.00    13.517     -\n",
      "    5  1.006  -10.521       -         -         -         -      13.517     -\n",
      "    6  1.008   -9.760       -         -         -         -      13.517     -\n",
      "    7  0.998  -12.439       -         -      233.80      0.00    13.517     -\n",
      "    8  0.998  -13.163       -         -      522.00      0.00    13.517     -\n",
      "    9  1.038  -16.577       -         -        6.50      0.00    13.517     -\n",
      "   10  1.018   -7.022       -         -         -         -      13.517     -\n",
      "   11  1.013   -7.963       -         -         -         -      13.517     -\n",
      "   12  1.001   -7.943       -         -        8.53      0.00    13.517     -\n",
      "   13  1.015   -7.710       -         -         -         -      13.517     -\n",
      "   14  1.012   -9.272       -         -         -         -      13.517     -\n",
      "   15  1.016   -9.299       -         -      320.00      0.00    13.517     -\n",
      "   16  1.033   -7.587       -         -      329.00      0.00    13.517     -\n",
      "   17  1.034   -8.831       -         -         -         -      13.517     -\n",
      "   18  1.032   -9.594       -         -      158.00      0.00    13.517     -\n",
      "   19  1.050   -2.224       -         -         -         -      13.517     -\n",
      "   20  0.991   -3.666       -         -      680.00      0.00    13.517     -\n",
      "   21  1.032   -4.862       -         -      274.00      0.00    13.517     -\n",
      "   22  1.050    0.161       -         -         -         -      13.517     -\n",
      "   23  1.045   -0.029       -         -      247.50      0.00    13.517     -\n",
      "   24  1.038   -7.390       -         -      308.60      0.00    13.517     -\n",
      "   25  1.058   -4.873       -         -      224.00      0.00    13.517     -\n",
      "   26  1.053   -7.963       -         -      139.00      0.00    13.517     -\n",
      "   27  1.038   -9.641       -         -      281.00      0.00    13.517     -\n",
      "   28  1.050   -6.312       -         -      206.00      0.00    13.517     -\n",
      "   29  1.050   -4.003       -         -      283.50      0.00    13.517     -\n",
      "   30  1.050    1.299    660.85      0.00       -         -      13.517     -\n",
      "   31  0.982    0.000*   646.00      0.00      9.20      0.00    13.517     -\n",
      "   32  0.984    1.080    660.85      0.00       -         -      13.517     -\n",
      "   33  0.997    3.452    652.00      0.00       -         -      13.517     -\n",
      "   34  1.012    1.620    508.00      0.00       -         -      13.517     -\n",
      "   35  1.049    5.711    660.85      0.00       -         -      13.517     -\n",
      "   36  1.064    9.010    580.00      0.00       -         -      13.517     -\n",
      "   37  1.028    2.812    564.00      0.00       -         -      13.517     -\n",
      "   38  1.026    2.051    660.85      0.00       -         -      13.517     -\n",
      "   39  1.030  -18.836    660.85      0.00   1104.00      0.00    13.517     -\n",
      "                        --------  --------  --------  --------\n",
      "               Total:   6254.23      0.00   6254.23      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      2   -383.08      0.00    383.08      0.00     0.000      0.00\n",
      "   1      1     39    285.48      0.00   -285.48      0.00     0.000      0.00\n",
      "   2      2      3    450.81      0.00   -450.81      0.00     0.000      0.00\n",
      "   3      2     25   -173.05      0.00    173.05      0.00     0.000      0.00\n",
      "   4      2     30   -660.85      0.00    660.85      0.00     0.000      0.00\n",
      "   5      3      4    133.05      0.00   -133.05      0.00     0.000      0.00\n",
      "   6      3     18     -4.24      0.00      4.24      0.00     0.000      0.00\n",
      "   7      4      5    -99.32      0.00     99.32      0.00     0.000      0.00\n",
      "   8      4     14   -267.62      0.00    267.62      0.00     0.000      0.00\n",
      "   9      5      6   -510.99      0.00    510.99      0.00     0.000      0.00\n",
      "  10      5      8    411.66      0.00   -411.66      0.00     0.000      0.00\n",
      "  11      6      7    508.31      0.00   -508.31      0.00     0.000      0.00\n",
      "  12      6     11   -382.50      0.00    382.50      0.00     0.000      0.00\n",
      "  13      6     31   -636.80      0.00    636.80      0.00     0.000      0.00\n",
      "  14      7      8    274.51      0.00   -274.51      0.00     0.000      0.00\n",
      "  15      8      9    164.17      0.00   -164.17      0.00     0.000      0.00\n",
      "  16      9     39    157.67      0.00   -157.67      0.00     0.000      0.00\n",
      "  17     10     11    381.72      0.00   -381.72      0.00     0.000      0.00\n",
      "  18     10     13    279.13      0.00   -279.13      0.00     0.000      0.00\n",
      "  19     10     32   -660.85      0.00    660.85      0.00     0.000      0.00\n",
      "  20     12     11      0.78      0.00     -0.78      0.00     0.000      0.00\n",
      "  21     12     13     -9.31      0.00      9.31      0.00     0.000      0.00\n",
      "  22     13     14    269.82      0.00   -269.82      0.00     0.000      0.00\n",
      "  23     14     15      2.20      0.00     -2.20      0.00     0.000      0.00\n",
      "  24     15     16   -317.80      0.00    317.80      0.00     0.000      0.00\n",
      "  25     16     17    243.94      0.00   -243.94      0.00     0.000      0.00\n",
      "  26     16     19   -480.00      0.00    480.00      0.00     0.000      0.00\n",
      "  27     16     21   -352.29      0.00    352.29      0.00     0.000      0.00\n",
      "  28     16     24    -58.46      0.00     58.46      0.00     0.000      0.00\n",
      "  29     17     18    162.24      0.00   -162.24      0.00     0.000      0.00\n",
      "  30     17     27     81.70      0.00    -81.70      0.00     0.000      0.00\n",
      "  31     19     20    172.00      0.00   -172.00      0.00     0.000      0.00\n",
      "  32     19     33   -652.00      0.00    652.00      0.00     0.000      0.00\n",
      "  33     20     34   -508.00      0.00    508.00      0.00     0.000      0.00\n",
      "  34     21     22   -626.29      0.00    626.29      0.00     0.000      0.00\n",
      "  35     22     23     34.56      0.00    -34.56      0.00     0.000      0.00\n",
      "  36     22     35   -660.85      0.00    660.85      0.00     0.000      0.00\n",
      "  37     23     24    367.06      0.00   -367.06      0.00     0.000      0.00\n",
      "  38     23     36   -580.00      0.00    580.00      0.00     0.000      0.00\n",
      "  39     25     26    166.95      0.00   -166.95      0.00     0.000      0.00\n",
      "  40     25     37   -564.00      0.00    564.00      0.00     0.000      0.00\n",
      "  41     26     27    199.30      0.00   -199.30      0.00     0.000      0.00\n",
      "  42     26     28    -60.79      0.00     60.79      0.00     0.000      0.00\n",
      "  43     26     29   -110.56      0.00    110.56      0.00     0.000      0.00\n",
      "  44     28     29   -266.79      0.00    266.79      0.00     0.000      0.00\n",
      "  45     29     38   -660.85      0.00    660.85      0.00     0.000      0.00\n",
      "                                                             --------  --------\n",
      "                                                    Total:     0.000      0.00\n",
      "\n",
      "================================================================================\n",
      "|     Generation Constraints                                                   |\n",
      "================================================================================\n",
      " Gen   Bus                Active Power Limits\n",
      "  #     #    Pmin mu    Pmin       Pg       Pmax    Pmax mu\n",
      "----  -----  -------  --------  --------  --------  -------\n",
      "   1    31      -        0.00    646.00    646.00    0.297\n",
      "   3    33      -        0.00    652.00    652.00    0.177\n",
      "   4    34      -        0.00    508.00    508.00    3.057\n",
      "   6    36      -        0.00    580.00    580.00    1.617\n",
      "   7    37      -        0.00    564.00    564.00    1.937\n"
     ]
    }
   ],
   "source": [
    "# Run PyPower case to get initial condition\n",
    "ppc = case39()\n",
    "ppopt = ppoption(VERBOSE=0)\n",
    "pp_dcpf = rundcopf(ppc, ppopt)\n",
    "pp_bus = pp_dcpf[\"bus\"]\n",
    "pp_branch = pp_dcpf[\"branch\"]\n",
    "pp_gen = pp_dcpf[\"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",
    "\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": "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*300/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": "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": "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(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",
    "            # 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[j],Tset[0:num_TCL[j],j],delta,R[0:num_TCL[j],j],C[0:num_TCL[j],j],P_el[0:num_TCL[j],j],COP[0:num_TCL[j],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[j],Tset[0:num_TCL[j],j],delta,R[0:num_TCL[j],j],C[0:num_TCL[j],j],P_el[0:num_TCL[j],j],COP[0:num_TCL[j],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": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Parameters related to the demand side\n",
    "num_TCL = np.zeros(NL)\n",
    "for k in range(NL):\n",
    "    num_TCL[k] = 50*np.ceil(PL[k]*baseMVA)\n",
    "\n",
    "num_TCL = num_TCL.astype(int)\n",
    "\n",
    "max_num_TCL = np.max(num_TCL).astype(int)\n",
    "\n",
    "Tset = np.zeros((max_num_TCL, NL))\n",
    "R = np.zeros((max_num_TCL, NL))\n",
    "C = np.zeros((max_num_TCL, NL))\n",
    "COP = np.zeros((max_num_TCL, NL))\n",
    "P_el = np.zeros((max_num_TCL, NL))\n",
    "Tset_set = [21,21.5,22,22.5,23]\n",
    "\n",
    "for i in range(max_num_TCL): \n",
    "    for j in range(NL):\n",
    "        s = random.sample(Tset_set,1)\n",
    "        Tset[i,j] = s[0] \n",
    "\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(0)\n",
    "\n",
    "for i in range(NL): \n",
    "    R[:,i] = np.random.uniform(R_low,R_high,max_num_TCL)\n",
    "    C[:,i] = np.random.uniform(C_low,C_high,max_num_TCL)\n",
    "    COP[:,i] = np.random.uniform(COP_low,COP_high,max_num_TCL)\n",
    "    for j in range(max_num_TCL):\n",
    "        P_el[j,i] = np.ceil((Tset[j,i] + delta/2 - Tamb_min)/(R[j,i]*COP[j,i])) + 3\n",
    "# Sampling\n",
    "eps = 0.1\n",
    "beta = 0.1\n",
    "num_decision_variables = 5*NG + 5*NL\n",
    "num_samples = 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",
    "# np.random.seed(1)\n",
    "# Tamb_sample = np.random.normal(mu_Tamb,sigma_Tamb,num_samples)\n",
    "# PW_sample = np.random.choice(PW_data[i,:], size=num_samples)\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[j],Tset[0:num_TCL[j],j],delta,R[0:num_TCL[j],j],C[0:num_TCL[j],j],P_el[0:num_TCL[j],j],COP[0:num_TCL[j],j],baseMVA)\n",
    "    S0[j] = (Sdown[j] + Sup[j])/2    \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1567\n",
      "901\n"
     ]
    }
   ],
   "source": [
    "# Sampling\n",
    "eps = 0.2\n",
    "beta = 0.2\n",
    "num_decision_variables = 5*NG + 5*NL\n",
    "num_samples = 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",
    "print(num_samples)\n",
    "print(num_samples2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "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": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "def cost(model,CG0,CG1,CG2,CRG_up,CRG_down,CRL_up,CRL_down,NG,NL):\n",
    "    r_pg = np.zeros(NG)\n",
    "    r_rg_up = np.zeros(NG)\n",
    "    r_rg_down = np.zeros(NG)\n",
    "    r_rl_up = np.zeros(NL)\n",
    "    r_rl_down = np.zeros(NL)\n",
    "    r_pc = np.zeros(NL)\n",
    "    gen_cost = np.zeros(NG)\n",
    "    gen_r_cost = np.zeros(NG)\n",
    "    load_r_cost = np.zeros(NL)\n",
    "\n",
    "\n",
    "    for i in range(NG):\n",
    "        r_pg[i] = model.PG[i]()\n",
    "        r_rg_up[i] = model.RG_up[i]()\n",
    "        r_rg_down[i] = model.RG_down[i]()\n",
    "        gen_cost[i] = CG0[i] + CG1[i]*r_pg[i] + CG2[i]*r_pg[i]**2\n",
    "        gen_r_cost[i] = CRG_up[i]*r_rg_up[i] + CRG_down[i]*r_rg_down[i]\n",
    "\n",
    "    for i in range(NL):\n",
    "        r_rl_up[i] = model.RL_up[i]()\n",
    "        r_rl_down[i] = model.RL_down[i]()\n",
    "        r_pc[i] = model.PC[i]()\n",
    "        load_r_cost[i] = CRL_up[i]*r_rl_up[i] + CRL_down[i]*r_rl_down[i]\n",
    "\n",
    "    total_cost = np.sum(gen_cost + gen_r_cost) + np.sum(load_r_cost)\n",
    "\n",
    "    return r_pg,r_rg_up,r_rg_down,r_rl_up,r_rl_down,r_pc,gen_cost,gen_r_cost,load_r_cost,total_cost"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Hour  0\n",
      "S0:  [ 2743.7359902   9001.64174456 13968.95079242  6539.46595792\n",
      " 14572.70970699   194.37525354   250.37288175  8954.00481698\n",
      "  9194.94944642  4415.58057107 19000.92300132  7664.33178879\n",
      "  6936.62063254  8636.44830484  6289.95395026  3882.68249753\n",
      "  7858.68936008  5761.76962214  7934.89635988   278.15895765\n",
      " 30845.83988752]\n",
      "S up 0: 5487.471980\n",
      "S up 1: 18003.283489\n",
      "S up 2: 27937.901585\n",
      "S up 3: 13078.931916\n",
      "S up 4: 29145.419414\n",
      "S up 5: 388.750507\n",
      "S up 6: 500.745764\n",
      "S up 7: 17908.009634\n",
      "S up 8: 18389.898893\n",
      "S up 9: 8831.161142\n",
      "S up 10: 38001.846003\n",
      "S up 11: 15328.663578\n",
      "S up 12: 13873.241265\n",
      "S up 13: 17272.896610\n",
      "S up 14: 12579.907901\n",
      "S up 15: 7765.364995\n",
      "S up 16: 15717.378720\n",
      "S up 17: 11523.539244\n",
      "S up 18: 15869.792720\n",
      "S up 19: 556.317915\n",
      "S up 20: 61691.679775\n",
      "Elapsed:  5429.622940778732\n",
      "Hour  1\n",
      "S0:  [5.24918020e+03 8.09673415e+03 1.34200951e+04 5.46016019e+03\n",
      " 1.40527105e+04 6.25798479e-01 8.00929511e-01 8.05177460e+03\n",
      " 8.29865057e+03 3.19517339e+03 1.88165876e+04 6.65947474e+03\n",
      " 5.88439765e+03 7.72528519e+03 1.17390083e+04 2.62153496e+03\n",
      " 6.87326403e+03 4.63950103e+03 6.95971967e+03 8.83495932e-01\n",
      " 3.45853091e+04]\n",
      "S up 0: 5445.567372\n",
      "S up 1: 17865.350752\n",
      "S up 2: 27723.930801\n",
      "S up 3: 12978.753898\n",
      "S up 4: 28922.412463\n",
      "S up 5: 385.734926\n",
      "S up 6: 496.900862\n",
      "S up 7: 17770.719157\n",
      "S up 8: 18249.082058\n",
      "S up 9: 8763.458500\n",
      "S up 10: 37710.792562\n",
      "S up 11: 15211.289775\n",
      "S up 12: 13766.980524\n",
      "S up 13: 17140.269680\n",
      "S up 14: 12483.809154\n",
      "S up 15: 7705.887558\n",
      "S up 16: 15596.939706\n",
      "S up 17: 11435.128216\n",
      "S up 18: 15748.114382\n",
      "S up 19: 552.109705\n",
      "S up 20: 61219.438448\n",
      "Elapsed:  4325.677249908447\n",
      "Hour  2\n",
      "S0:  [2.05789990e+03 8.30909521e+03 1.38312618e+04 5.57465684e+03\n",
      " 1.44874997e+04 3.79038104e-01 4.82637742e-01 8.26245255e+03\n",
      " 8.51854784e+03 3.22508389e+03 1.94292924e+04 6.81876253e+03\n",
      " 6.01473907e+03 7.92376988e+03 1.15506385e+04 2.62998578e+03\n",
      " 7.04052617e+03 4.72331598e+03 7.13018649e+03 5.04691517e-01\n",
      " 3.25887806e+04]\n",
      "S up 0: 5403.662765\n",
      "S up 1: 17727.418015\n",
      "S up 2: 27509.960016\n",
      "S up 3: 12878.575879\n",
      "S up 4: 28699.405513\n",
      "S up 5: 382.719344\n",
      "S up 6: 493.055960\n",
      "S up 7: 17633.428678\n",
      "S up 8: 18108.265223\n",
      "S up 9: 8695.755858\n",
      "S up 10: 37419.739122\n",
      "S up 11: 15093.915971\n",
      "S up 12: 13660.719782\n",
      "S up 13: 17007.642745\n",
      "S up 14: 12387.710411\n",
      "S up 15: 7646.410122\n",
      "S up 16: 15476.500691\n",
      "S up 17: 11346.717187\n",
      "S up 18: 15626.436043\n",
      "S up 19: 547.901496\n",
      "S up 20: 60747.197119\n",
      "Elapsed:  4301.436253786087\n",
      "Hour  3\n",
      "S0:  [1.27905341e+03 8.22981275e+03 1.37519883e+04 5.49511991e+03\n",
      " 1.44082275e+04 3.41062030e-01 4.32968731e-01 8.18317089e+03\n",
      " 8.43926637e+03 3.14554298e+03 1.93500279e+04 6.73922747e+03\n",
      " 5.93520276e+03 7.84448792e+03 1.14969882e+04 2.55044803e+03\n",
      " 6.96099281e+03 4.64378171e+03 7.05065599e+03 4.57108282e-01\n",
      " 3.25096130e+04]\n",
      "S up 0: 5361.758153\n",
      "S up 1: 17589.485256\n",
      "S up 2: 27295.989199\n",
      "S up 3: 12778.397845\n",
      "S up 4: 28476.398531\n",
      "S up 5: 379.703761\n",
      "S up 6: 489.211058\n",
      "S up 7: 17496.138176\n",
      "S up 8: 17967.448367\n",
      "S up 9: 8628.053204\n",
      "S up 10: 37128.685636\n",
      "S up 11: 14976.542149\n",
      "S up 12: 13554.459023\n",
      "S up 13: 16875.015786\n",
      "S up 14: 12291.611655\n",
      "S up 15: 7586.932677\n",
      "S up 16: 15356.061656\n",
      "S up 17: 11258.306142\n",
      "S up 18: 15504.757683\n",
      "S up 19: 543.693286\n",
      "S up 20: 60274.955718\n",
      "Elapsed:  4414.740390300751\n",
      "Hour  4\n",
      "S0:  [2.98188409e+03 8.05673986e+03 1.33580339e+04 5.42884719e+03\n",
      " 1.39880246e+04 5.20605799e-01 4.62629177e-01 8.01196482e+03\n",
      " 8.25781662e+03 3.17325068e+03 1.87321570e+04 6.62339340e+03\n",
      " 5.85132692e+03 7.68682906e+03 1.16057846e+04 2.60196599e+03\n",
      " 6.83709765e+03 4.61156826e+03 6.92350532e+03 4.98584395e-01\n",
      " 3.13658518e+04]\n",
      "S up 0: 5333.821739\n",
      "S up 1: 17497.530061\n",
      "S up 2: 27153.341953\n",
      "S up 3: 12711.612473\n",
      "S up 4: 28327.727174\n",
      "S up 5: 377.693371\n",
      "S up 6: 486.647789\n",
      "S up 7: 17404.611151\n",
      "S up 8: 17873.570439\n",
      "S up 9: 8582.918090\n",
      "S up 10: 36934.649932\n",
      "S up 11: 14898.292915\n",
      "S up 12: 13483.618499\n",
      "S up 13: 16786.597789\n",
      "S up 14: 12227.545804\n",
      "S up 15: 7547.281038\n",
      "S up 16: 15275.768946\n",
      "S up 17: 11199.365430\n",
      "S up 18: 15423.638756\n",
      "S up 19: 540.887813\n",
      "S up 20: 59960.128040\n",
      "Elapsed:  4950.610293388367\n",
      "Hour  5\n",
      "S0:  [1.37077032e+03 8.04726791e+03 1.33485680e+04 5.41818819e+03\n",
      " 1.39785590e+04 3.76875241e-01 4.75521755e-01 8.00249251e+03\n",
      " 8.24834440e+03 3.16258936e+03 1.87226974e+04 6.61391926e+03\n",
      " 5.84066858e+03 7.67735644e+03 1.16097980e+04 2.59130039e+03\n",
      " 6.82762261e+03 4.60090523e+03 6.91402817e+03 5.09072706e-01\n",
      " 3.13560772e+04]\n",
      "S up 0: 5319.853529\n",
      "S up 1: 17451.552455\n",
      "S up 2: 27082.018316\n",
      "S up 3: 12678.219780\n",
      "S up 4: 28253.391481\n",
      "S up 5: 376.688176\n",
      "S up 6: 485.366155\n",
      "S up 7: 17358.847629\n",
      "S up 8: 17826.631466\n",
      "S up 9: 8560.350528\n",
      "S up 10: 36837.632061\n",
      "S up 11: 14859.168290\n",
      "S up 12: 13448.198230\n",
      "S up 13: 16742.388781\n",
      "S up 14: 12195.512873\n",
      "S up 15: 7527.455214\n",
      "S up 16: 15235.622584\n",
      "S up 17: 11169.895068\n",
      "S up 18: 15383.079284\n",
      "S up 19: 539.485076\n",
      "S up 20: 59802.714170\n",
      "Elapsed:  4714.280195474625\n",
      "Hour  6\n",
      "S0:  [1.27961532e+03 7.95532786e+03 1.32566312e+04 5.32703976e+03\n",
      " 1.38866225e+04 5.16794483e-01 5.09428350e-01 7.91055244e+03\n",
      " 8.15640438e+03 3.07143961e+03 1.86307639e+04 6.52197824e+03\n",
      " 5.74952046e+03 7.58541622e+03 1.15867533e+04 2.50014952e+03\n",
      " 6.73568150e+03 4.50975570e+03 6.82208669e+03 5.41811648e-01\n",
      " 3.12640994e+04]\n",
      "S up 0: 5319.853529\n",
      "S up 1: 17451.552455\n",
      "S up 2: 27082.018316\n",
      "S up 3: 12678.219780\n",
      "S up 4: 28253.391481\n",
      "S up 5: 376.688176\n",
      "S up 6: 485.366155\n",
      "S up 7: 17358.847629\n",
      "S up 8: 17826.631466\n",
      "S up 9: 8560.350528\n",
      "S up 10: 36837.632061\n",
      "S up 11: 14859.168290\n",
      "S up 12: 13448.198230\n",
      "S up 13: 16742.388781\n",
      "S up 14: 12195.512873\n",
      "S up 15: 7527.455214\n",
      "S up 16: 15235.622584\n",
      "S up 17: 11169.895068\n",
      "S up 18: 15383.079284\n",
      "S up 19: 539.485076\n",
      "S up 20: 59802.714170\n",
      "Elapsed:  5108.238531589508\n",
      "Hour  7\n",
      "S0:  [1.20592011e+03 7.88019254e+03 1.31814958e+04 5.25334705e+03\n",
      " 1.38114870e+04 4.81373930e-01 6.05197809e-01 7.83541701e+03\n",
      " 8.08126895e+03 2.99774696e+03 1.85556285e+04 6.44769361e+03\n",
      " 5.67582777e+03 7.51028078e+03 1.15434509e+04 2.42645628e+03\n",
      " 6.66058961e+03 4.43606247e+03 6.74695069e+03 6.45141678e-01\n",
      " 3.11889274e+04]\n",
      "S up 0: 5333.821739\n",
      "S up 1: 17497.530061\n",
      "S up 2: 27153.341953\n",
      "S up 3: 12711.612473\n",
      "S up 4: 28327.727174\n",
      "S up 5: 377.693371\n",
      "S up 6: 486.647789\n",
      "S up 7: 17404.611151\n",
      "S up 8: 17873.570439\n",
      "S up 9: 8582.918090\n",
      "S up 10: 36934.649932\n",
      "S up 11: 14898.292915\n",
      "S up 12: 13483.618499\n",
      "S up 13: 16786.597789\n",
      "S up 14: 12227.545804\n",
      "S up 15: 7547.281038\n",
      "S up 16: 15275.768946\n",
      "S up 17: 11199.365430\n",
      "S up 18: 15423.638756\n",
      "S up 19: 540.887813\n",
      "S up 20: 59960.128040\n",
      "Elapsed:  4715.591188669205\n",
      "Hour  8\n",
      "S0:  [1.21098903e+03 7.88415735e+03 1.31854568e+04 5.25843000e+03\n",
      " 1.38154472e+04 1.85379250e+00 2.36297426e+00 7.83938105e+03\n",
      " 8.08523288e+03 3.00283177e+03 1.85595858e+04 6.45277593e+03\n",
      " 5.68091058e+03 7.51424492e+03 1.10279151e+04 2.43153697e+03\n",
      " 6.66567034e+03 4.44114181e+03 6.75174600e+03 2.58623664e+00\n",
      " 3.11926740e+04]\n",
      "S up 0: 5431.599170\n",
      "S up 1: 17819.373174\n",
      "S up 2: 27652.607207\n",
      "S up 3: 12945.361226\n",
      "S up 4: 28848.076815\n",
      "S up 5: 384.729732\n",
      "S up 6: 495.619228\n",
      "S up 7: 17724.955666\n",
      "S up 8: 18202.143114\n",
      "S up 9: 8740.890953\n",
      "S up 10: 37613.774751\n",
      "S up 11: 15172.165175\n",
      "S up 12: 13731.560278\n",
      "S up 13: 17096.060703\n",
      "S up 14: 12451.776240\n",
      "S up 15: 7686.061746\n",
      "S up 16: 15556.793368\n",
      "S up 17: 11405.657874\n",
      "S up 18: 15707.554937\n",
      "S up 19: 550.706969\n",
      "S up 20: 61062.024675\n",
      "Elapsed:  4495.958758592606\n",
      "Hour  9\n",
      "S0:  [1.25965549e+03 8.21175272e+03 1.37339214e+04 5.47576036e+03\n",
      " 1.43901588e+04 2.36902674e+00 3.01487783e+00 8.16510914e+03\n",
      " 8.42120441e+03 3.12618696e+03 1.93319544e+04 6.71986693e+03\n",
      " 5.91584307e+03 7.82642629e+03 5.18044664e+03 2.53108013e+03\n",
      " 6.94162808e+03 4.62441146e+03 7.03128221e+03 3.29102547e+00\n",
      " 3.24904818e+04]\n",
      "S up 0: 5543.344805\n",
      "S up 1: 18187.193842\n",
      "S up 2: 28223.196021\n",
      "S up 3: 13212.502633\n",
      "S up 4: 29442.762077\n",
      "S up 5: 392.771283\n",
      "S up 6: 505.872300\n",
      "S up 7: 18091.063636\n",
      "S up 8: 18577.654709\n",
      "S up 9: 8921.431348\n",
      "S up 10: 38389.917333\n",
      "S up 11: 15485.162011\n",
      "S up 12: 14014.922280\n",
      "S up 13: 17449.732545\n",
      "S up 14: 12708.039592\n",
      "S up 15: 7844.668261\n",
      "S up 16: 15877.964103\n",
      "S up 17: 11641.420636\n",
      "S up 18: 16032.030533\n",
      "S up 19: 561.928864\n",
      "S up 20: 62321.335003\n",
      "Elapsed:  4483.101762771606\n",
      "Hour  10\n",
      "S0:  [1.26227260e+03 7.95975063e+03 1.32823347e+04 5.32656149e+03\n",
      " 1.39149460e+04 2.67399753e+00 3.40889992e+00 7.91495638e+03\n",
      " 8.16089880e+03 3.06158104e+03 1.86788173e+04 6.52587513e+03\n",
      " 5.75079921e+03 7.58970025e+03 1.03193052e+04 2.48790676e+03\n",
      " 6.73965386e+03 4.50586951e+03 6.82608456e+03 3.73099886e+00\n",
      " 3.13641215e+04]\n",
      "S up 0: 5696.995267\n",
      "S up 1: 18692.947915\n",
      "S up 2: 29007.756659\n",
      "S up 3: 13579.822543\n",
      "S up 4: 30260.455389\n",
      "S up 5: 403.828424\n",
      "S up 6: 519.970293\n",
      "S up 7: 18594.462729\n",
      "S up 8: 19093.983819\n",
      "S up 9: 9169.674705\n",
      "S up 10: 39457.114768\n",
      "S up 11: 15915.533213\n",
      "S up 12: 14404.545528\n",
      "S up 13: 17936.031914\n",
      "S up 14: 13060.402190\n",
      "S up 15: 8062.752510\n",
      "S up 16: 16319.574429\n",
      "S up 17: 11965.594839\n",
      "S up 18: 16478.185039\n",
      "S up 19: 577.358994\n",
      "S up 20: 64052.888947\n",
      "Elapsed:  4498.758327722549\n",
      "Hour  11\n",
      "S0:  [1.25207222e+03 7.89395365e+03 1.31723527e+04 5.28262364e+03\n",
      " 1.37997128e+04 2.63621216e+00 3.34763018e+00 7.84953044e+03\n",
      " 8.09343172e+03 3.03644568e+03 1.85240364e+04 6.47198105e+03\n",
      " 5.70333938e+03 7.52697380e+03 1.19149902e+04 2.46753340e+03\n",
      " 6.68398486e+03 4.46874343e+03 6.76969722e+03 3.66341292e+00\n",
      " 3.11040215e+04]\n",
      "S up 0: 5836.677963\n",
      "S up 1: 19152.725808\n",
      "S up 2: 29720.995878\n",
      "S up 3: 13913.750796\n",
      "S up 4: 31003.815326\n",
      "S up 5: 413.880398\n",
      "S up 6: 532.786691\n",
      "S up 7: 19052.099708\n",
      "S up 8: 19563.375409\n",
      "S up 9: 9395.351195\n",
      "S up 10: 40427.297344\n",
      "S up 11: 16306.781000\n",
      "S up 12: 14758.749599\n",
      "S up 13: 18378.123616\n",
      "S up 14: 13380.732871\n",
      "S up 15: 8261.011557\n",
      "S up 16: 16721.039634\n",
      "S up 17: 12260.299584\n",
      "S up 18: 16883.781322\n",
      "S up 19: 591.386435\n",
      "S up 20: 65627.033905\n",
      "Elapsed:  4338.190977096558\n",
      "Hour  12\n",
      "S0:  [1.24882096e+03 7.89068297e+03 1.31638927e+04 5.27935503e+03\n",
      " 1.37912533e+04 2.10200605e+00 2.61224830e+00 7.84625983e+03\n",
      " 8.09016223e+03 3.03317768e+03 1.85155737e+04 6.46871155e+03\n",
      " 5.70007024e+03 7.52370288e+03 1.20358988e+04 2.46427012e+03\n",
      " 6.68071655e+03 4.46547857e+03 6.76643111e+03 2.84626926e+00\n",
      " 3.10960448e+04]\n",
      "S up 0: 5948.424604\n",
      "S up 1: 19520.549688\n",
      "S up 2: 30291.589687\n",
      "S up 3: 14180.894535\n",
      "S up 4: 31598.505816\n",
      "S up 5: 421.922009\n",
      "S up 6: 543.039854\n",
      "S up 7: 19418.210842\n",
      "S up 8: 19938.890278\n",
      "S up 9: 9575.893151\n",
      "S up 10: 41203.446711\n",
      "S up 11: 16619.780560\n",
      "S up 12: 15042.114054\n",
      "S up 13: 18731.798459\n",
      "S up 14: 13636.998525\n",
      "S up 15: 8419.619475\n",
      "S up 16: 17042.213161\n",
      "S up 17: 12496.064376\n",
      "S up 18: 17208.259724\n",
      "S up 19: 602.608438\n",
      "S up 20: 66886.355232\n",
      "Elapsed:  6419.565704584122\n",
      "Hour  13\n",
      "S0:  [5.07296982e+03 7.68091912e+03 1.27118336e+04 5.18977666e+03\n",
      " 1.33103212e+04 1.95535389e+00 2.48362614e+00 7.63854095e+03\n",
      " 7.87121831e+03 3.04697730e+03 1.78172076e+04 6.32439363e+03\n",
      " 5.59112842e+03 7.33082906e+03 1.23345182e+04 2.50425956e+03\n",
      " 6.52664313e+03 4.41336340e+03 6.60841626e+03 2.70977240e+00\n",
      " 2.98192443e+04]\n",
      "S up 0: 6032.234952\n",
      "S up 1: 19796.418790\n",
      "S up 2: 30719.536897\n",
      "S up 3: 14381.253205\n",
      "S up 4: 32044.525618\n",
      "S up 5: 427.953241\n",
      "S up 6: 550.729759\n",
      "S up 7: 19692.795378\n",
      "S up 8: 20220.527647\n",
      "S up 9: 9711.300201\n",
      "S up 10: 41785.561258\n",
      "S up 11: 16854.531243\n",
      "S up 12: 15254.638311\n",
      "S up 13: 18997.055725\n",
      "S up 14: 13829.198605\n",
      "S up 15: 8538.575931\n",
      "S up 16: 17283.094347\n",
      "S up 17: 12672.888731\n",
      "S up 18: 17451.619575\n",
      "S up 19: 611.024979\n",
      "S up 20: 67830.850314\n",
      "Elapsed:  5723.99827170372\n",
      "Hour  14\n",
      "S0:  [3.32816045e+03 7.73505047e+03 1.28108647e+04 5.22205343e+03\n",
      " 1.34146039e+04 1.79056020e+00 2.28177238e+00 7.69230125e+03\n",
      " 7.92702020e+03 3.06045374e+03 1.79610304e+04 6.36662458e+03\n",
      " 5.62692610e+03 7.38188947e+03 1.23645745e+04 2.51298231e+03\n",
      " 6.57065062e+03 4.43883522e+03 6.65314613e+03 2.49734762e+00\n",
      " 3.00726499e+04]\n",
      "S up 0: 6060.171815\n",
      "S up 1: 19888.375422\n",
      "S up 2: 30862.186377\n",
      "S up 3: 14448.039620\n",
      "S up 4: 32193.199313\n",
      "S up 5: 429.963658\n",
      "S up 6: 553.293068\n",
      "S up 7: 19784.323818\n",
      "S up 8: 20314.407039\n",
      "S up 9: 9756.436014\n",
      "S up 10: 41979.599997\n",
      "S up 11: 16932.781696\n",
      "S up 12: 15325.479933\n",
      "S up 13: 19085.475065\n",
      "S up 14: 13893.265485\n",
      "S up 15: 8578.228197\n",
      "S up 16: 17363.388306\n",
      "S up 17: 12731.830352\n",
      "S up 18: 17532.739758\n",
      "S up 19: 613.830501\n",
      "S up 20: 68145.682912\n",
      "Elapsed:  4775.451537132263\n",
      "Hour  15\n",
      "S0:  [1.27428879e+03 8.22636022e+03 1.37485032e+04 5.49038056e+03\n",
      " 1.44047391e+04 3.25032302e+00 4.08646188e+00 8.17971402e+03\n",
      " 8.43581049e+03 3.14081729e+03 1.93465097e+04 6.73448074e+03\n",
      " 5.93046042e+03 7.84103026e+03 5.19505810e+03 2.54571312e+03\n",
      " 6.95624035e+03 4.63903212e+03 7.04589163e+03 4.37941628e+00\n",
      " 3.25049728e+04]\n",
      "S up 0: 6018.266536\n",
      "S up 1: 19750.440526\n",
      "S up 2: 30648.212237\n",
      "S up 3: 14347.860035\n",
      "S up 4: 31970.188854\n",
      "S up 5: 426.948034\n",
      "S up 6: 549.448106\n",
      "S up 7: 19647.031209\n",
      "S up 8: 20173.588003\n",
      "S up 9: 9688.732320\n",
      "S up 10: 41688.541997\n",
      "S up 11: 16815.406061\n",
      "S up 12: 15219.217540\n",
      "S up 13: 18952.846104\n",
      "S up 14: 13797.165202\n",
      "S up 15: 8518.749820\n",
      "S up 16: 17242.947413\n",
      "S up 17: 12643.417954\n",
      "S up 18: 17411.059529\n",
      "S up 19: 609.622219\n",
      "S up 20: 67673.434192\n",
      "Elapsed:  4415.850698232651\n",
      "Hour  16\n",
      "S0:  [1.25411207e+03 7.89600776e+03 1.31744115e+04 5.28467383e+03\n",
      " 1.38017730e+04 2.00615763e+00 2.55727717e+00 7.85158329e+03\n",
      " 8.09548627e+03 3.03849252e+03 1.85261023e+04 6.47403266e+03\n",
      " 5.70538979e+03 7.52902503e+03 1.16998168e+04 2.46957871e+03\n",
      " 6.68603628e+03 4.47079037e+03 6.77174730e+03 2.78332705e+00\n",
      " 3.11061465e+04]\n",
      "S up 0: 5920.487895\n",
      "S up 1: 19428.593561\n",
      "S up 2: 30148.940991\n",
      "S up 3: 14114.108486\n",
      "S up 4: 31449.832939\n",
      "S up 5: 419.911603\n",
      "S up 6: 540.476559\n",
      "S up 7: 19326.682903\n",
      "S up 8: 19845.011401\n",
      "S up 9: 9530.757585\n",
      "S up 10: 41009.409038\n",
      "S up 11: 16541.530536\n",
      "S up 12: 14971.272820\n",
      "S up 13: 18643.379599\n",
      "S up 14: 13572.932000\n",
      "S up 15: 8379.967427\n",
      "S up 16: 16961.919642\n",
      "S up 17: 12437.123078\n",
      "S up 18: 17127.139985\n",
      "S up 19: 599.802932\n",
      "S up 20: 66571.524363\n",
      "Elapsed:  4514.826112508774\n",
      "Hour  17\n",
      "S0:  [5.01372572e+03 7.67906054e+03 1.27152459e+04 5.18655022e+03\n",
      " 1.33137334e+04 1.77557985e+00 2.26310121e+00 7.63652182e+03\n",
      " 7.87008127e+03 3.04375147e+03 1.78206182e+04 6.32116676e+03\n",
      " 5.58790182e+03 7.32764343e+03 1.23160788e+04 2.50103540e+03\n",
      " 6.52341671e+03 4.41013868e+03 6.60519077e+03 2.47644465e+00\n",
      " 2.98226467e+04]\n",
      "S up 0: 5850.646266\n",
      "S up 1: 19198.703706\n",
      "S up 2: 29792.319969\n",
      "S up 3: 13947.143701\n",
      "S up 4: 31078.151496\n",
      "S up 5: 414.885598\n",
      "S up 6: 534.068334\n",
      "S up 7: 19097.863514\n",
      "S up 8: 19610.314679\n",
      "S up 9: 9417.918897\n",
      "S up 10: 40524.315831\n",
      "S up 11: 16345.905872\n",
      "S up 12: 14794.170089\n",
      "S up 13: 18422.332889\n",
      "S up 14: 13412.766016\n",
      "S up 15: 8280.837509\n",
      "S up 16: 16761.186249\n",
      "S up 17: 12289.770128\n",
      "S up 18: 16924.341046\n",
      "S up 19: 592.789183\n",
      "S up 20: 65784.448774\n",
      "Elapsed:  4855.422092914581\n",
      "Hour  18\n",
      "S0:  [3.45649017e+03 7.71721346e+03 1.27975733e+04 5.19950304e+03\n",
      " 1.34013106e+04 6.41628719e-01 8.16932577e-01 7.67430230e+03\n",
      " 7.90991023e+03 3.03790887e+03 1.79477269e+04 6.34407168e+03\n",
      " 5.60437505e+03 7.36271487e+03 1.22535147e+04 2.49043581e+03\n",
      " 6.54840039e+03 4.41628474e+03 6.63120772e+03 8.72659900e-01\n",
      " 3.00580826e+04]\n",
      "S up 0: 5808.741376\n",
      "S up 1: 19060.770077\n",
      "S up 2: 29578.347798\n",
      "S up 3: 13846.965035\n",
      "S up 4: 30855.143091\n",
      "S up 5: 411.870001\n",
      "S up 6: 530.223407\n",
      "S up 7: 18960.572162\n",
      "S up 8: 19469.496936\n",
      "S up 9: 9350.215823\n",
      "S up 10: 40233.260507\n",
      "S up 11: 16228.531314\n",
      "S up 12: 14687.908668\n",
      "S up 13: 18289.705132\n",
      "S up 14: 13316.666627\n",
      "S up 15: 8221.359681\n",
      "S up 16: 16640.746460\n",
      "S up 17: 12201.358538\n",
      "S up 18: 16802.661932\n",
      "S up 19: 588.580942\n",
      "S up 20: 65312.204393\n",
      "Elapsed:  6904.716997623444\n",
      "Hour  19\n",
      "S0:  [2.54044177e+03 7.72833188e+03 1.28086905e+04 5.21032764e+03\n",
      " 1.34124269e+04 8.09685249e-01 8.76022861e-01 7.68541989e+03\n",
      " 7.92102762e+03 3.04873437e+03 1.79588431e+04 6.35489635e+03\n",
      " 5.61519983e+03 7.37383259e+03 1.22204672e+04 2.50125481e+03\n",
      " 6.55951586e+03 4.42710364e+03 6.64231852e+03 9.36793619e-01\n",
      " 3.00661041e+04]\n",
      "S up 0: 5738.900018\n",
      "S up 1: 18830.881097\n",
      "S up 2: 29221.728137\n",
      "S up 3: 13680.000885\n",
      "S up 4: 30483.463068\n",
      "S up 5: 406.844013\n",
      "S up 6: 523.815207\n",
      "S up 7: 18731.753639\n",
      "S up 8: 19234.801107\n",
      "S up 9: 9237.377562\n",
      "S up 10: 39748.169149\n",
      "S up 11: 16032.907392\n",
      "S up 12: 14510.806607\n",
      "S up 13: 18068.659250\n",
      "S up 14: 13156.501262\n",
      "S up 15: 8122.230143\n",
      "S up 16: 16440.013829\n",
      "S up 17: 12054.006145\n",
      "S up 18: 16599.863761\n",
      "S up 19: 581.567220\n",
      "S up 20: 64525.131799\n",
      "Elapsed:  4793.162998199463\n",
      "Hour  20\n",
      "S0:  [1.23730364e+03 7.88163391e+03 1.31607905e+04 5.26784290e+03\n",
      " 1.37881509e+04 1.84603488e+00 2.35312592e+00 7.83704230e+03\n",
      " 8.08186921e+03 3.02166535e+03 1.85124725e+04 6.45719976e+03\n",
      " 5.68855834e+03 7.51326258e+03 9.95892932e+03 2.45275650e+03\n",
      " 6.66920451e+03 4.45396568e+03 6.75491867e+03 2.57546205e+00\n",
      " 3.10928133e+04]\n",
      "S up 0: 5683.027026\n",
      "S up 1: 18646.970217\n",
      "S up 2: 28936.432879\n",
      "S up 3: 13546.429784\n",
      "S up 4: 30186.119543\n",
      "S up 5: 402.823229\n",
      "S up 6: 518.688655\n",
      "S up 7: 18548.699121\n",
      "S up 8: 19047.044753\n",
      "S up 9: 9147.107101\n",
      "S up 10: 39360.096703\n",
      "S up 11: 15876.408512\n",
      "S up 12: 14369.125191\n",
      "S up 13: 17891.822830\n",
      "S up 14: 13028.369186\n",
      "S up 15: 8042.926645\n",
      "S up 16: 16279.427988\n",
      "S up 17: 11936.124422\n",
      "S up 18: 16437.625490\n",
      "S up 19: 575.956253\n",
      "S up 20: 63895.474764\n",
      "Elapsed:  4310.100108385086\n",
      "Hour  21\n",
      "S0:  [1.15897184e+03 8.11103785e+03 1.36331725e+04 5.37505115e+03\n",
      " 1.42894067e+04 1.93062491e+00 2.46129934e+00 8.06439318e+03\n",
      " 8.32048814e+03 3.02549181e+03 1.92311708e+04 6.61914972e+03\n",
      " 5.81513080e+03 7.72571138e+03 5.07973925e+03 2.43038906e+03\n",
      " 6.84090940e+03 4.52370626e+03 6.93057045e+03 2.69383631e+00\n",
      " 3.23896166e+04]\n",
      "S up 0: 5613.185878\n",
      "S up 1: 18417.081912\n",
      "S up 2: 28579.814267\n",
      "S up 3: 13379.466123\n",
      "S up 4: 29814.440618\n",
      "S up 5: 397.797254\n",
      "S up 6: 512.280474\n",
      "S up 7: 18319.881264\n",
      "S up 8: 18812.349611\n",
      "S up 9: 9034.269167\n",
      "S up 10: 38875.006769\n",
      "S up 11: 15680.785162\n",
      "S up 12: 14192.023646\n",
      "S up 13: 17670.777579\n",
      "S up 14: 12868.204304\n",
      "S up 15: 7943.797401\n",
      "S up 16: 16078.695944\n",
      "S up 17: 11788.772455\n",
      "S up 18: 16234.827909\n",
      "S up 19: 568.942554\n",
      "S up 20: 63108.404480\n",
      "Elapsed:  4476.483024358749\n",
      "Hour  22\n",
      "S0:  [1.23894727e+03 7.88083945e+03 1.31592446e+04 5.26950637e+03\n",
      " 1.37866051e+04 8.01618604e-01 1.02134603e+00 7.83641651e+03\n",
      " 8.08031771e+03 3.02332596e+03 1.85109345e+04 6.45886528e+03\n",
      " 5.69022273e+03 7.51385981e+03 8.58686755e+03 2.45441213e+03\n",
      " 6.67086915e+03 4.45562484e+03 6.75658134e+03 9.33209982e-01\n",
      " 3.10909410e+04]\n",
      "S up 0: 5557.313015\n",
      "S up 1: 18233.171442\n",
      "S up 2: 28294.519648\n",
      "S up 3: 13245.895321\n",
      "S up 4: 29517.097762\n",
      "S up 5: 393.776477\n",
      "S up 6: 507.153935\n",
      "S up 7: 18136.827148\n",
      "S up 8: 18624.593675\n",
      "S up 9: 8943.998905\n",
      "S up 10: 38486.935191\n",
      "S up 11: 15524.286630\n",
      "S up 12: 14050.342542\n",
      "S up 13: 17493.941539\n",
      "S up 14: 12740.072524\n",
      "S up 15: 7864.494083\n",
      "S up 16: 15918.110459\n",
      "S up 17: 11670.890991\n",
      "S up 18: 16072.589996\n",
      "S up 19: 563.331601\n",
      "S up 20: 62478.748850\n",
      "Elapsed:  4509.725719213486\n",
      "Hour  23\n",
      "S0:  [1.05726886e+03 8.00967931e+03 1.35318270e+04 5.27335729e+03\n",
      " 1.41880624e+04 8.44064712e-01 1.07252226e+00 7.96303527e+03\n",
      " 8.21913025e+03 2.92379259e+03 1.91298385e+04 6.51745891e+03\n",
      " 5.71343813e+03 7.62435318e+03 4.97804616e+03 2.32868856e+03\n",
      " 6.73921855e+03 4.42201124e+03 6.82921012e+03 1.12419707e+00\n",
      " 3.22883159e+04]\n",
      "Elapsed:  4290.281541585922\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[j],Tset[0:num_TCL[j],j],delta,R[0:num_TCL[j],j],C[0:num_TCL[j],j],P_el[0:num_TCL[j],j],COP[0:num_TCL[j],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_sample1 = Tamb[h]*np.ones(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",
    "    # PW_sample = PW_f*np.ones(nsample)\n",
    "    # PW_sample[PW_sample > 0.15] = 0.1\n",
    "    # PW_f = np.median(PW_sample)\n",
    "    \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",
    "    # print(\"S up: \", identified_VS_all[4])\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",
    "    # Tamb_samples2.append(Tamb_sample2)\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[j],Tset[0:num_TCL[j],j],delta,R[0:num_TCL[j],j],C[0:num_TCL[j],j],P_el[0:num_TCL[j],j],COP[0:num_TCL[j],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": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pickle\n",
    "\n",
    "case_name = \"case39_with_DR\"\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_sample, 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()\n",
    "\n",
    "# Write P_B_test\n",
    "file = open(\"Results/P_B_test_%s_%d\" %(case_name,nsample), \"wb\")\n",
    "pickle.dump(P_B_test, file)\n",
    "file.close()\n",
    "\n",
    "# Write Pup_test\n",
    "file = open(\"Results/Pup_test_%s_%d\" %(case_name,nsample), \"wb\")\n",
    "pickle.dump(Pup_test, file)\n",
    "file.close()\n",
    "\n",
    "# Write Pdown_test\n",
    "file = open(\"Results/Pdown_test_%s_%d\" %(case_name,nsample), \"wb\")\n",
    "pickle.dump(Pdown_test, file)\n",
    "file.close()\n",
    "\n",
    "# Write Sup_test\n",
    "file = open(\"Results/Sup_test_%s_%d\" %(case_name,nsample), \"wb\")\n",
    "pickle.dump(Sup_test, file)\n",
    "file.close()\n",
    "\n",
    "# Write Sdown_test\n",
    "file = open(\"Results/Sdown_test_%s_%d\" %(case_name,nsample), \"wb\")\n",
    "pickle.dump(Sdown_test, 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
}