SOPF without DR_case39.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Packages"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "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 case39, ext2int, bustypes, makeBdc, rundcpf, ppoption, rundcopf, case30, case9, case118, case300\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": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Converged in 0.04 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, OPF_ALG_DC=200)\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": 6,
   "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 = 0\n",
    "CW0 = 0\n",
    "\n",
    "# Cost coefficients for generation reserve\n",
    "CRG_up = np.ones((NG,1))*10*baseMVA\n",
    "CRG_down = np.ones((NG,1))*10*baseMVA\n",
    "\n",
    "# Cost coefficients for load reserve\n",
    "CRL_up = 1000*np.ones((NL,1))*baseMVA\n",
    "CRL_down = 1000*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": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([10.4 ,  6.46,  7.25,  6.52,  5.08,  6.87,  5.8 ,  5.64,  8.65,\n",
       "       11.  ])"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Pmax"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "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"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "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)/2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "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": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Optimization: Problem formulation\n",
    "def optimizer(NG,NL,Tamb_sample,PW_sample,PW_data,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)\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)\n",
    "    # model.RL_up = Var(range(NL),within = NonNegativeReals)\n",
    "    # model.RL = Var(range(NL),within = NonNegativeReals)\n",
    "    # model.dG_down = Var(range(NG),within = NonNegativeReals)\n",
    "    # model.dG_up = Var(range(NG),within = NonNegativeReals)\n",
    "    # model.dL_down = Var(range(NL),within = NonNegativeReals)\n",
    "    # model.dL_up = Var(range(NL),within = NonNegativeReals)\n",
    "    # model.P_mis = Var(range(nsample), within=Reals)\n",
    "    model.Va = Var(range(NB), within=Reals, initialize=0, bounds=(-1,1))\n",
    "    # model.y = Var(range(20), within=Binary)\n",
    "    # model.a = Var(range(20), within=NonNegativeReals, initialize=0)\n",
    "    # model.s1 = Var(range(nsample), within=NonNegativeReals, initialize=0)\n",
    "    # model.s2 = Var(range(nsample), within=NonNegativeReals, initialize=0)\n",
    "\n",
    "    model.dual = Suffix(direction=Suffix.IMPORT)\n",
    "\n",
    "    ## Defining the deterministic constraints\n",
    "    model.c = ConstraintList()\n",
    "\n",
    "    \n",
    "\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",
    "    # model.c.add(sum(model.dG_down[i] for i in range(NG)) == 1)\n",
    "    # model.c.add(sum(model.dG_up[i] for i in range(NG)) == 1)\n",
    "    \n",
    "    # Line flow limit constraints\n",
    "    for i in range(NBr):\n",
    "        model.c.add(sum(Bf[i,j]*model.Va[j] for j in range(NB)) <= Lmax[i])\n",
    "        model.c.add(sum(Bf[i,j]*model.Va[j] for j in range(NB)) >= Lmin[i])\n",
    "    \n",
    "    # i = 100\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] >= Pmin[j])\n",
    "            model.c.add(model.PG[j] <= 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",
    "            ### Add PC RL\n",
    "            model.c.add(sum(CG[j,k]*model.PG[k] for k in range(NG))\\\n",
    "                - sum(CL[j,k]*PL[k] for k in range(NL)) == sum(B[j,k]*model.Va[k] for k 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] 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] for k in range(NL)) <= sum(B[j,k]*model.Va[k] for k in range(NB)) + model.s2[i])\n",
    "\n",
    "        \n",
    "        # P_mis_wind = PW_sample[i] - PW_f\n",
    "\n",
    "        # P_mis_load = 0\n",
    "        # for j in range(NL):\n",
    "        #     identified_VS = VS(0.25,Tamb_sample[j],num_TCL,Tset[:,j],delta,R[:,j],C[:,j],P_el[:,j],COP[:,j],baseMVA)\n",
    "        #     P_mis_load = P_mis_load + (model.PC[j] - identified_VS[0])\n",
    "        # P_mis = P_mis_wind - P_mis_load\n",
    "        \n",
    "        # P_mis = P_mis_wind\n",
    "        # d = max(-P_mis, 0)\n",
    "        # e = max(P_mis, 0)\n",
    "        \n",
    "        # model.c.add(model.RG[i] == model.dG_up[i]*max(-model.P_mis,0) - model.dG_down[i]*max(model.P_mis,0))\n",
    "        # for j in range(NG):\n",
    "            # model.c.add(model.RG[j] == model.dG_up[j]*model.a[2*i] - model.dG_down[j]*model.a[2*i+1])\n",
    "            # model.c.add(model.RG[j,i] == model.dG_up[j]*d - model.dG_down[j]*e)\n",
    "            # model.c.add(model.a[2*i+1] <= model.y[2*i+1]*P_mis[i])\n",
    "            # model.c.add(model.a[2*i] <= -model.y[2*i]*P_mis[i])\n",
    "        \n",
    "\n",
    "    # for j in range(NL):\n",
    "    #     # model.c.add(model.RL[j] == model.dL_up[j]*max(model.P_mis,0) - model.dL_down[j]*max(-model.P_mis,0))\n",
    "    #     model.c.add(model.RL[j] == model.dL_up[j]*P_mis - model.dL_down[j]*(-P_mis))\n",
    "    #     identified_VS = VS(0.25,Tamb_sample[j],num_TCL,Tset[:,j],delta,R[:,j],C[:,j],P_el[:,j],COP[:,j],baseMVA)\n",
    "    #     model.c.add(model.PC[i] + model.RL[i] >= identified_VS[1])\n",
    "    #     model.c.add(model.PC[i] + model.RL[i] <= identified_VS[2])\n",
    "    #     model.c.add(model.RL[i] >= -model.RL_down[i])\n",
    "    #     model.c.add(model.RL[i] <= model.RL_up[i])\n",
    "    #     model.c.add(S0[j] + 0.25*(model.PC[i] + model.RL[i] - identified_VS[0])*1000*baseMVA >= identified_VS[3])\n",
    "    #     model.c.add(S0[j] + 0.25*(model.PC[i] + model.RL[i] - identified_VS[0])*1000*baseMVA  <= identified_VS[4])\n",
    "    #     model.c.add(S0[j] + 0.75*(model.PC[i] - identified_VS[0])*1000*baseMVA  + 0.25*(model.PC[i] + model.RL[i] - identified_VS[0])*1000*baseMVA  >= identified_VS[3])\n",
    "    #     model.c.add(S0[j] + 0.75*(model.PC[i] - identified_VS[0])*1000*baseMVA  + 0.25*(model.PC[i] + model.RL[i] - identified_VS[0])*1000*baseMVA  <= identified_VS[4])         \n",
    "    \n",
    "    def cost(model, CG0, CG1, CG2, CRG_down, CRG_up, nsample, CRL_down, CRL_up):\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",
    "        #     + 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",
    "        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",
    "                \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",
    "        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     "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Parameters related to the demand side\n",
    "num_TCL = 3000\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",
    "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])) + 1\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",
    "i = 0\n",
    "mu_Tamb = Tamb[i]\n",
    "sigma_Tamb = std_Tamb[i] \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[i],num_TCL,Tset[:,j],delta,R[:,j],C[:,j],P_el[:,j],COP[:,j],baseMVA)\n",
    "    S0[j] = (Sdown[j] + Sup[j])/2    \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define the model\n",
    "nsample= 1\n",
    "# Pmax = 10*Pmax\n",
    "PW_f = np.median(PW_data)\n",
    "model = optimizer(NG,NL,Tamb_sample,PW_sample,PW_data,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)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Model unknown\n",
      "\n",
      "  Variables:\n",
      "    PG : Size=10, Index=PG_index\n",
      "        Key : Lower : Value              : Upper : Fixed : Stale : Domain\n",
      "          0 :     0 :  6.608459941000233 :  None : False : False : NonNegativeReals\n",
      "          1 :     0 :  6.460000064599663 :  None : False : False : NonNegativeReals\n",
      "          2 :     0 :  6.608459941000192 :  None : False : False : NonNegativeReals\n",
      "          3 :     0 :  6.520000065199435 :  None : False : False : NonNegativeReals\n",
      "          4 :     0 :  5.080000050799967 :  None : False : False : NonNegativeReals\n",
      "          5 :     0 :  6.608459941000068 :  None : False : False : NonNegativeReals\n",
      "          6 :     0 :  5.800000057999938 :  None : False : False : NonNegativeReals\n",
      "          7 :     0 :  5.640000056399948 :  None : False : False : NonNegativeReals\n",
      "          8 :     0 : 6.6084599410002465 :  None : False : False : NonNegativeReals\n",
      "          9 :     0 :   6.60845994100031 :  None : False : False : NonNegativeReals\n",
      "    Va : Size=39, Index=Va_index\n",
      "        Key : Lower : Value                 : Upper : Fixed : Stale : Domain\n",
      "          0 :    -1 :   -0.1327754642806673 :     1 : False : False :  Reals\n",
      "          1 :    -1 :  0.024671596272499402 :     1 : False : False :  Reals\n",
      "          2 :    -1 :  -0.04340108180809893 :     1 : False : False :  Reals\n",
      "          3 :    -1 :  -0.07174155580124031 :     1 : False : False :  Reals\n",
      "          4 :    -1 :  -0.05902809960026357 :     1 : False : False :  Reals\n",
      "          5 :    -1 :  -0.04574245893094485 :     1 : False : False :  Reals\n",
      "          6 :    -1 :  -0.09250687133591312 :     1 : False : False :  Reals\n",
      "          7 :    -1 :  -0.10513427753839724 :     1 : False : False :  Reals\n",
      "          8 :    -1 :  -0.16472839900174166 :     1 : False : False :  Reals\n",
      "          9 :    -1 :  0.002036072116914492 :     1 : False : False :  Reals\n",
      "         10 :    -1 :  -0.01437786863426882 :     1 : False : False :  Reals\n",
      "         11 :    -1 : -0.014038523406235917 :     1 : False : False :  Reals\n",
      "         12 :    -1 : -0.009966364878203016 :     1 : False : False :  Reals\n",
      "         13 :    -1 : -0.037218284916886925 :     1 : False : False :  Reals\n",
      "         14 :    -1 :  -0.03769541053439441 :     1 : False : False :  Reals\n",
      "         15 :    -1 : -0.007822091677462159 :     1 : False : False :  Reals\n",
      "         16 :    -1 : -0.029533174166201867 :     1 : False : False :  Reals\n",
      "         17 :    -1 : -0.042836971499390517 :     1 : False : False :  Reals\n",
      "         18 :    -1 :   0.08577791058452612 :     1 : False : False :  Reals\n",
      "         19 :    -1 :   0.06061775132762805 :     1 : False : False :  Reals\n",
      "         20 :    -1 :   0.03973638840973609 :     1 : False : False :  Reals\n",
      "         21 :    -1 :    0.1274162936853491 :     1 : False : False :  Reals\n",
      "         22 :    -1 :   0.12409844186931163 :     1 : False : False :  Reals\n",
      "         23 :    -1 : -0.004372894906656162 :     1 : False : False :  Reals\n",
      "         24 :    -1 :  0.039553837440928497 :     1 : False : False :  Reals\n",
      "         25 :    -1 : -0.014371235341224774 :     1 : False : False :  Reals\n",
      "         26 :    -1 :  -0.04366784786394861 :     1 : False : False :  Reals\n",
      "         27 :    -1 :  0.014442370060481075 :     1 : False : False :  Reals\n",
      "         28 :    -1 :   0.05472738781503294 :     1 : False : False :  Reals\n",
      "         29 :    -1 :   0.14727504932790625 :     1 : False : False :  Reals\n",
      "         30 :    -1 :   0.12460154279709615 :     1 : False : False :  Reals\n",
      "         31 :    -1 :    0.1434571148543186 :     1 : False : False :  Reals\n",
      "         32 :    -1 :   0.18484279157516634 :     1 : False : False :  Reals\n",
      "         33 :    -1 :     0.152880712250257 :     1 : False : False :  Reals\n",
      "         34 :    -1 :   0.22427979527055758 :     1 : False : False :  Reals\n",
      "         35 :    -1 :    0.2818584434469099 :     1 : False : False :  Reals\n",
      "         36 :    -1 :   0.17367303878211923 :     1 : False : False :  Reals\n",
      "         37 :    -1 :   0.16039666227162688 :     1 : False : False :  Reals\n",
      "         38 :    -1 :  -0.20414618237870058 :     1 : False : False :  Reals\n",
      "\n",
      "  Objectives:\n",
      "    obj : Size=1, Index=None, Active=True\n",
      "        Key  : Active : Value\n",
      "        None :   True : 41263.940746896915\n",
      "\n",
      "  Constraints:\n",
      "    c : Size=151\n",
      "        Key : Lower                : Body                    : Upper\n",
      "          1 :                 None :     -3.8308287239213312 :                  6.0\n",
      "          2 :                 -6.0 :     -3.8308287239213312 :                 None\n",
      "          3 :                 None :      2.8548287239213304 :                 10.0\n",
      "          4 :                -10.0 :      2.8548287239213304 :                 None\n",
      "          5 :                 None :       4.508124376198564 :                  5.0\n",
      "          6 :                 -5.0 :       4.508124376198564 :                 None\n",
      "          7 :                 None :      -1.730493159119662 :                  5.0\n",
      "          8 :                 -5.0 :      -1.730493159119662 :                 None\n",
      "          9 :                 None :      -6.608459941000233 :                  9.0\n",
      "         10 :                 -9.0 :      -6.608459941000233 :                 None\n",
      "         11 :                 None :      1.3305386851240084 :                  5.0\n",
      "         12 :                 -5.0 :      1.3305386851240084 :                 None\n",
      "         13 :                 None :    -0.04241430892544473 :                  5.0\n",
      "         14 :                 -5.0 :    -0.04241430892544473 :                 None\n",
      "         15 :                 None :     -0.9932387657013075 :                  6.0\n",
      "         16 :                 -6.0 :     -0.9932387657013075 :                 None\n",
      "         17 :                 None :     -2.6762225491746805 :                  5.0\n",
      "         18 :                 -5.0 :     -2.6762225491746805 :                 None\n",
      "         19 :                 None :      -5.109861795891817 :                 12.0\n",
      "         20 :                -12.0 :      -5.109861795891817 :                 None\n",
      "         21 :                 None :       4.116623030190507 :                  9.0\n",
      "         22 :                 -9.0 :       4.116623030190507 :                 None\n",
      "         23 :                 None :       5.083088304887856 :                  9.0\n",
      "         24 :                 -9.0 :       5.083088304887856 :                 None\n",
      "         25 :                 None :     -3.8249500361800033 :                  4.8\n",
      "         26 :                 -4.8 :     -3.8249500361800033 :                 None\n",
      "         27 :                 None :      -6.368000064599663 :                 18.0\n",
      "         28 :                -18.0 :      -6.368000064599663 :                 None\n",
      "         29 :                 None :       2.745088304887851 :                  9.0\n",
      "         30 :                 -9.0 :       2.745088304887851 :                 None\n",
      "         31 :                 None :      1.6417113350783583 :                  9.0\n",
      "         32 :                 -9.0 :      1.6417113350783583 :                 None\n",
      "         33 :                 None :      1.5767113350783566 :                  9.0\n",
      "         34 :                 -9.0 :      1.5767113350783566 :                 None\n",
      "         35 :                 None :      3.8171955235310024 :                  6.0\n",
      "         36 :                 -6.0 :      3.8171955235310024 :                 None\n",
      "         37 :                 None :      2.7912644174691876 :                  6.0\n",
      "         38 :                 -6.0 :      2.7912644174691876 :                 None\n",
      "         39 :                 None :      -6.608459941000192 :                  9.0\n",
      "         40 :                 -9.0 :      -6.608459941000192 :                 None\n",
      "         41 :                 None :    0.007754512649000289 :                  5.0\n",
      "         42 :                 -5.0 :    0.007754512649000289 :                 None\n",
      "         43 :                 None :     -0.0930545126490003 :                  5.0\n",
      "         44 :                 -5.0 :     -0.0930545126490003 :                 None\n",
      "         45 :                 None :      2.6982099048201893 :                  6.0\n",
      "         46 :                 -6.0 :      2.6982099048201893 :                 None\n",
      "         47 :                 None :    0.021987355645506135 :                  6.0\n",
      "         48 :                 -6.0 :    0.021987355645506135 :                 None\n",
      "         49 :                 None :     -3.1780126443544945 :                  6.0\n",
      "         50 :                 -6.0 :     -3.1780126443544945 :                 None\n",
      "         51 :                 None :      2.4394474706449114 :                  6.0\n",
      "         52 :                 -6.0 :      2.4394474706449114 :                 None\n",
      "         53 :                 None :      -4.800000115999399 :                  6.0\n",
      "         54 :                 -6.0 :      -4.800000115999399 :                 None\n",
      "         55 :                 None :        -3.5228503768295 :                  6.0\n",
      "         56 :                 -6.0 :        -3.5228503768295 :                 None\n",
      "         57 :                 None :      -0.584609622170508 :                  6.0\n",
      "         58 :                 -6.0 :      -0.584609622170508 :                 None\n",
      "         59 :                 None :      1.6224143089254448 :                  6.0\n",
      "         60 :                 -6.0 :      1.6224143089254448 :                 None\n",
      "         61 :                 None :       0.817033161719465 :                  6.0\n",
      "         62 :                 -6.0 :       0.817033161719465 :                 None\n",
      "         63 :                 None :       1.719999949200032 :                  9.0\n",
      "         64 :                 -9.0 :       1.719999949200032 :                 None\n",
      "         65 :                 None :      -6.520000065199434 :                  9.0\n",
      "         66 :                 -9.0 :      -6.520000065199434 :                 None\n",
      "         67 :                 None :      -5.080000050799967 :                  9.0\n",
      "         68 :                 -9.0 :      -5.080000050799967 :                 None\n",
      "         69 :                 None :        -6.2628503768295 :                  9.0\n",
      "         70 :                 -9.0 :        -6.2628503768295 :                 None\n",
      "         71 :                 None :      0.3456095641705694 :                  6.0\n",
      "         72 :                 -6.0 :      0.3456095641705694 :                 None\n",
      "         73 :                 None :      -6.608459941000069 :                  9.0\n",
      "         74 :                 -9.0 :      -6.608459941000069 :                 None\n",
      "         75 :                 None :      3.6706096221705082 :                  6.0\n",
      "         76 :                 -6.0 :      3.6706096221705082 :                 None\n",
      "         77 :                 None :      -5.800000057999937 :                  9.0\n",
      "         78 :                 -9.0 :      -5.800000057999937 :                 None\n",
      "         79 :                 None :       1.669506897280287 :                  6.0\n",
      "         80 :                 -6.0 :       1.669506897280287 :                 None\n",
      "         81 :                 None :     -5.6400000563999475 :                  9.0\n",
      "         82 :                 -9.0 :     -5.6400000563999475 :                 None\n",
      "         83 :                 None :      1.9929668382805334 :                  6.0\n",
      "         84 :                 -6.0 :      1.9929668382805334 :                 None\n",
      "         85 :                 None :     -0.6078819705001234 :                  6.0\n",
      "         86 :                 -6.0 :     -0.6078819705001234 :                 None\n",
      "         87 :                 None :     -1.1055779705001234 :                  6.0\n",
      "         88 :                 -6.0 :     -1.1055779705001234 :                 None\n",
      "         89 :                 None :     -2.6678819705001233 :                  6.0\n",
      "         90 :                 -6.0 :     -2.6678819705001233 :                 None\n",
      "         91 :                 None :     -6.6084599410002465 :                 12.0\n",
      "         92 :                -12.0 :     -6.6084599410002465 :                 None\n",
      "         93 :                  0.0 :       6.608459941000233 :                 None\n",
      "         94 :                 None :       6.608459941000233 :                 10.4\n",
      "         95 :                  0.0 :       6.460000064599663 :                 None\n",
      "         96 :                 None :       6.460000064599663 :                 6.46\n",
      "         97 :                  0.0 :       6.608459941000192 :                 None\n",
      "         98 :                 None :       6.608459941000192 :                 7.25\n",
      "         99 :                  0.0 :       6.520000065199435 :                 None\n",
      "        100 :                 None :       6.520000065199435 :                 6.52\n",
      "        101 :                  0.0 :       5.080000050799967 :                 None\n",
      "        102 :                 None :       5.080000050799967 :                 5.08\n",
      "        103 :                  0.0 :       6.608459941000068 :                 None\n",
      "        104 :                 None :       6.608459941000068 :                 6.87\n",
      "        105 :                  0.0 :       5.800000057999938 :                 None\n",
      "        106 :                 None :       5.800000057999938 :                  5.8\n",
      "        107 :                  0.0 :       5.640000056399948 :                 None\n",
      "        108 :                 None :       5.640000056399948 :                 5.64\n",
      "        109 :                  0.0 :      6.6084599410002465 :                 None\n",
      "        110 :                 None :      6.6084599410002465 :                 8.65\n",
      "        111 :                  0.0 :        6.60845994100031 :                 None\n",
      "        112 :                 None :        6.60845994100031 :                 11.0\n",
      "        113 :               -0.976 :     -0.9760000000000009 :               -0.976\n",
      "        114 :                  0.0 :  -8.881784197001252e-16 :                  0.0\n",
      "        115 :                -3.22 :     -3.2199999999999998 :                -3.22\n",
      "        116 :                 -5.0 :     -4.9999999999999964 :                 -5.0\n",
      "        117 :                  0.0 :  -7.105427357601002e-15 :                  0.0\n",
      "        118 :                  0.0 :   4.440892098500626e-15 :                  0.0\n",
      "        119 :               -2.338 :      -2.338000000000001 :               -2.338\n",
      "        120 :                -5.22 :      -5.220000000000001 :                -5.22\n",
      "        121 :               -0.065 :    -0.06500000000000128 :               -0.065\n",
      "        122 :                  0.0 : -1.7763568394002505e-15 :                  0.0\n",
      "        123 :                  0.0 :   5.551115123125783e-16 :                  0.0\n",
      "        124 : -0.08529999999999999 :    -0.08530000000000001 : -0.08529999999999999\n",
      "        125 :                  0.0 :   2.220446049250313e-15 :                  0.0\n",
      "        126 :                  0.0 : -2.4424906541753444e-15 :                  0.0\n",
      "        127 :                 -3.2 :     -3.2000000000000015 :                 -3.2\n",
      "        128 :                -3.29 :     -3.2900000000000014 :                -3.29\n",
      "        129 :                  0.0 :  -8.881784197001252e-16 :                  0.0\n",
      "        130 :                -1.58 :                   -1.58 :                -1.58\n",
      "        131 :                  0.0 :                     0.0 :                  0.0\n",
      "        132 :                 -6.8 :                    -6.8 :                 -6.8\n",
      "        133 :                -2.74 :     -2.7399999999999993 :                -2.74\n",
      "        134 :                  0.0 :                     0.0 :                  0.0\n",
      "        135 :               -2.475 :     -2.4749999999999996 :               -2.475\n",
      "        136 :  -3.0860000000000003 :                  -3.086 :  -3.0860000000000003\n",
      "        137 :                -2.24 :     -2.2399999999999993 :                -2.24\n",
      "        138 :                -1.39 :     -1.3900000000000001 :                -1.39\n",
      "        139 :                -2.81 :     -2.8099999999999987 :                -2.81\n",
      "        140 :                -2.06 :                   -2.06 :                -2.06\n",
      "        141 :               -2.835 :      -2.835000000000001 :               -2.835\n",
      "        142 :                  0.0 :                     0.0 :                  0.0\n",
      "        143 :                  0.0 :   8.881784197001252e-16 :                  0.0\n",
      "        144 :                  0.0 :                     0.0 :                  0.0\n",
      "        145 :                  0.0 :   8.881784197001252e-16 :                  0.0\n",
      "        146 :                  0.0 :                     0.0 :                  0.0\n",
      "        147 :                  0.0 :  -8.881784197001252e-16 :                  0.0\n",
      "        148 :                  0.0 :   8.881784197001252e-16 :                  0.0\n",
      "        149 :                  0.0 :   8.881784197001252e-16 :                  0.0\n",
      "        150 :                  0.0 :                     0.0 :                  0.0\n",
      "        151 :                  0.0 : -1.7763568394002505e-15 :                  0.0\n"
     ]
    }
   ],
   "source": [
    "# Solve the problem\n",
    "solver = SolverFactory(\"bonmin\")\n",
    "solver.solve(model)\n",
    "\n",
    "model.display()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "41263.940746896915"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model.obj()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[660.8459941 ]\n",
      " [646.00000646]\n",
      " [660.8459941 ]\n",
      " [652.00000652]\n",
      " [508.00000508]\n",
      " [660.8459941 ]\n",
      " [580.0000058 ]\n",
      " [564.00000564]\n",
      " [660.8459941 ]\n",
      " [660.8459941 ]]\n",
      "[[4565.62807741]\n",
      " [4367.1600854 ]\n",
      " [4565.62807741]\n",
      " [4446.84008698]\n",
      " [2733.24005314]\n",
      " [4565.62807741]\n",
      " [3538.20006902]\n",
      " [3350.36006531]\n",
      " [4565.62807741]\n",
      " [4565.62807741]]\n",
      "[[690.8459941 ]\n",
      " [676.00000646]\n",
      " [690.8459941 ]\n",
      " [682.00000652]\n",
      " [538.00000508]\n",
      " [690.8459941 ]\n",
      " [610.0000058 ]\n",
      " [594.00000564]\n",
      " [690.8459941 ]\n",
      " [690.8459941 ]]\n"
     ]
    }
   ],
   "source": [
    "# Print generation, cost and marginal cost\n",
    "r_pg = np.zeros((NG,1))\n",
    "r_cost = np.zeros((NG,1))\n",
    "r_mc = np.zeros((NG,1))\n",
    "for i in range(NG):\n",
    "    r_pg[i] = model.PG[i]()\n",
    "    r_cost[i] = CG0[i] + CG1[i]*r_pg[i] + CG2[i]*r_pg[i]**2\n",
    "    r_mc[i] = CG1[i] + CG2[i]*r_pg[i]\n",
    "\n",
    "print(r_pg*baseMVA)\n",
    "print(r_cost)\n",
    "print(r_mc)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([10.4 ,  6.46,  7.25,  6.52,  5.08,  6.87,  5.8 ,  5.64,  8.65,\n",
       "       11.  ])"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Pmax"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "PG : Size=10, Index=PG_index\n",
      "    Key : Lower : Value              : Upper : Fixed : Stale : Domain\n",
      "      0 :     0 :  6.608459941000233 :  None : False : False : NonNegativeReals\n",
      "      1 :     0 :  6.460000064599663 :  None : False : False : NonNegativeReals\n",
      "      2 :     0 :  6.608459941000192 :  None : False : False : NonNegativeReals\n",
      "      3 :     0 :  6.520000065199435 :  None : False : False : NonNegativeReals\n",
      "      4 :     0 :  5.080000050799967 :  None : False : False : NonNegativeReals\n",
      "      5 :     0 :  6.608459941000068 :  None : False : False : NonNegativeReals\n",
      "      6 :     0 :  5.800000057999938 :  None : False : False : NonNegativeReals\n",
      "      7 :     0 :  5.640000056399948 :  None : False : False : NonNegativeReals\n",
      "      8 :     0 : 6.6084599410002465 :  None : False : False : NonNegativeReals\n",
      "      9 :     0 :   6.60845994100031 :  None : False : False : NonNegativeReals\n"
     ]
    }
   ],
   "source": [
    "model.PG.pprint()"
   ]
  }
 ],
 "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
}