Market 3 bus Iterative Paper (Vectorized).ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "66eb3eff",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<style>.container { width:100% !important; }</style>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from IPython.display import display, HTML\n",
    "display(HTML(\"<style>.container { width:100% !important; }</style>\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "327b2320",
   "metadata": {},
   "outputs": [],
   "source": [
    "CONGESTION = 1\n",
    "\n",
    "SUPPLY_BID = 0\n",
    "LOAD_ATTACK = 0\n",
    "\n",
    "LOAD1 = 0\n",
    "DA_LOAD2 = 75       # in MWs\n",
    "\n",
    "RT_LOAD2 = 90 + LOAD_ATTACK # in MWs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "33239253",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "impedance_pu: 0.0001890359168241966, sus_pu: 5290.0\n",
      "Day Ahead Market(DAM) Data:\n",
      "\n",
      "Updated DAM Bus Dataframe:\n",
      "   bus_i  type  pd     qd  gs  bs  area  vm  va  basekv  zone  vmax  vmin\n",
      "0      1     2   0   0.00   0   0     1   1   0     230     1   1.1   0.9\n",
      "1      2     2  75   0.00   0   0     1   1   0     230     1   1.1   0.9\n",
      "2      3     1   0  98.61   0   0     1   1   0     230     1   1.1   0.9\n",
      "\n",
      "Updated DAM_gen:\n",
      "   bus   pg  qg   qmax   qmin  vg  mbase  status  pmax  pmin  ...  qc1min  \\\n",
      "0    1   40   0   30.0  -30.0   1    100       1  1000     0  ...       0   \n",
      "1    3  170   0  127.5 -127.5   1    100       1  1000     0  ...       0   \n",
      "\n",
      "   qc1max  qc2min  qc2max  ramp_agc  ramp_10  ramp_30  ramp_q  apf  id  \n",
      "0       0       0       0         0        0        0       0    0   1  \n",
      "1       0       0       0         0        0        0       0    0   2  \n",
      "\n",
      "[2 rows x 22 columns]\n",
      "\n",
      "Updated DAM_gencost:\n",
      "   model  startup  shutdown  n   x1  y1  id\n",
      "0      2        0         0  2  0.3   3   1\n",
      "1      2        0         0  2  0.8   8   2\n",
      "\n",
      "Updated Branch:\n",
      "   fbus  tbus        r         x        b  ratea  rateb  ratec  ratio  angle  \\\n",
      "0     1     3  0.00281  0.000189  0.00712    500    500    500      0      0   \n",
      "1     1     2  0.00281  0.000189  0.00712    500    500    500      0      0   \n",
      "2     2     3  0.00281  0.000189  0.00712    500    500    500      0      0   \n",
      "\n",
      "   status  angmin  angmax  id     sus  \n",
      "0       1    -360     360   1  5290.0  \n",
      "1       1    -360     360   2  5290.0  \n",
      "2       1    -360     360   3  5290.0  \n"
     ]
    }
   ],
   "source": [
    "import os\n",
    "import copy\n",
    "import pandas as pd\n",
    "\n",
    "def input_data(datadir):\n",
    "    \"\"\"\n",
    "    Reading 3-bus data from the CSV files in <datadir> \n",
    "    and processing them as dataframes.\n",
    "    \n",
    "    return (dataframes): bus, gen, gencost, branch\n",
    "    \"\"\"\n",
    "    \n",
    "    # Read the CSV files into dataframes\n",
    "    gen = pd.read_csv(os.path.join(datadir, 'gen.csv'))\n",
    "    gencost = pd.read_csv(os.path.join(datadir, 'gencost.csv'))\n",
    "    branch = pd.read_csv(os.path.join(datadir, 'branch.csv'))\n",
    "    bus = pd.read_csv(os.path.join(datadir, 'bus.csv'))\n",
    "\n",
    "    # Rename all columns to lowercase\n",
    "    gen.columns = gen.columns.str.lower()\n",
    "    gencost.columns = gencost.columns.str.lower()\n",
    "    branch.columns = branch.columns.str.lower()\n",
    "    bus.columns = bus.columns.str.lower()\n",
    "\n",
    "    # Create generator ids\n",
    "    gen['id'] = range(1, len(gen) + 1)\n",
    "    gencost['id'] = [1, 2]\n",
    "\n",
    "    # Create line ids\n",
    "    branch['id'] = range(1, len(branch) + 1)\n",
    "\n",
    "    # Calculate the susceptance of each line\n",
    "    branch['sus'] = 1 / branch['x']\n",
    "    \n",
    "    return bus, gen, gencost, branch\n",
    "\n",
    "\n",
    "# conversion to per-unit from Ohms\n",
    "baseKV = 230 * 1e3  # 230 kV\n",
    "mbase  = 100 * 1e6  # 100 MVA\n",
    "zbase  = ( baseKV * baseKV ) / mbase\n",
    "x = 0.1\n",
    "xpu = x/zbase\n",
    "sus = 1/xpu\n",
    "print(f\"impedance_pu: {xpu}, sus_pu: {sus}\")\n",
    "\n",
    "print(\"Day Ahead Market(DAM) Data:\\n\")\n",
    "DAM_bus, DAM_gen, DAM_gencost, branch = input_data(datadir = 'opf_data')\n",
    "\n",
    "# modify the data as per the requirements of the current system\n",
    "print(\"Updated DAM Bus Dataframe:\")\n",
    "DAM_bus[\"pd\"] = [LOAD1, DA_LOAD2, 0]      # demand on each bus (in MWs)\n",
    "print(DAM_bus)\n",
    "print(\"\")\n",
    "\n",
    "print(\"Updated DAM_gen:\")\n",
    "DAM_gen.loc[1, \"bus\"] = 3       # add generator to bus 3 \n",
    "print(DAM_gen)\n",
    "print(\"\")\n",
    "\n",
    "print(\"Updated DAM_gencost:\")\n",
    "DAM_gencost[\"x1\"] = [0.3, 0.8]  # generator cost curves\n",
    "DAM_gencost[\"y1\"] = [  3,   8]\n",
    "print(DAM_gencost)\n",
    "print(\"\")\n",
    "\n",
    "# print(\"Updated DAM_branch:\")\n",
    "branch['x'] = xpu               # in pu via per unit calculations\n",
    "branch['sus'] = 1 / branch['x']\n",
    "\n",
    "print(\"Updated Branch:\")\n",
    "print(branch)\n",
    "\n",
    "if CONGESTION:\n",
    "    branch.loc[((branch['fbus'] == 1) & (branch['tbus'] == 2)), 'ratea'] = 47  # creating congestion\n",
    "#    branch.loc[((branch['fbus'] == 2) & (branch['tbus'] == 3)), 'ratea'] = 47  # creating congestion"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fe258a43",
   "metadata": {},
   "source": [
    "### Construct all the relevant matrices from the above dataframes for optimization problem formulation (Overall approach to Construction):\n",
    "\n",
    "#### Optimization variable:\n",
    "- Ensure the dimension of `p = [x, y, v]` is fixed.\n",
    "- Both the below decisions to easily compute `p_{-i}`\n",
    "    -  Make `x` dynamically sized - `len(G)`\n",
    "    -  Make `y` and `v` constant sized - `len(N)` each\n",
    "    \n",
    "#### Objective Function:\n",
    "- (Done) `A`: Coefficient matrix for the constraints\n",
    "- (Done) `b`: Right-hand side vector for the constraints\n",
    "\n",
    "#### Modeling Constraints:\n",
    "- (Done) `S`: System matrix, possibly including incidence matrix `Inci` and matrix `B`\n",
    "- (Done) `c`: Capacity bounds associated with the constraint modeled as `lcb` and `ucb`\n",
    "- `\\Phi`: Hardcode based on the modeling choice of variable `p`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "d91739aa",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[0.3 0.  0.  0.  0.  0.  0.  0. ]\n",
      " [0.  0.8 0.  0.  0.  0.  0.  0. ]\n",
      " [0.  0.  0.  0.  0.  0.  0.  0. ]\n",
      " [0.  0.  0.  0.  0.  0.  0.  0. ]\n",
      " [0.  0.  0.  0.  0.  0.  0.  0. ]\n",
      " [0.  0.  0.  0.  0.  0.  0.  0. ]\n",
      " [0.  0.  0.  0.  0.  0.  0.  0. ]\n",
      " [0.  0.  0.  0.  0.  0.  0.  0. ]]\n",
      "[3. 8. 0. 0. 0. 0. 0. 0.]\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "\n",
    "# Modeling Objective Function constants\n",
    "G_DAM = DAM_gen['id'].values      # 1, 2\n",
    "N = DAM_bus['bus_i'].values   # all buses in the network: 1, 2, 3\n",
    "\n",
    "dim_p = len(G_DAM) + len(N) + len(N)\n",
    "\n",
    "# Modeling Objective Function Constants (A, b):\n",
    "A = np.zeros( ( dim_p, dim_p ) ) # initialize A as zeros # make the genreators dynamic\n",
    "b = np.zeros( dim_p )\n",
    "\n",
    "for i in G_DAM:  # go through each generator\n",
    "    alpha_i  = DAM_gencost.loc[DAM_gencost['id']==i, 'x1'].values[0]   # alpha_i associated \n",
    "    beta_i   = DAM_gencost.loc[DAM_gencost['id']==i, 'y1'].values[0]   #  beta_i associated\n",
    "    \n",
    "    bus_index = DAM_gen.loc[DAM_gen['id']==i, 'bus'].values[0]  # bus index associated with curr generator\n",
    "    \n",
    "    A[    i-1,     i-1] = alpha_i   # for physical supply (0, 1, 2)\n",
    "    #A[len(G_DAM) + len(N) + bus_index-1, len(G_DAM) + len(N) + bus_index-1] = alpha_i   # for physical supply (6, 7, 8) - hard coded for the 3-bus\n",
    "    \n",
    "    b[    i-1] = beta_i\n",
    "    #b[len(G_DAM) + len(N) + bus_index-1] = beta_i   # for physical supply (6, 7, 8) - hard coded for the 3-bus\n",
    "\n",
    "print( A )\n",
    "print( b )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "2b162cab",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Susceptance Matrix:\n",
      "[[10580. -5290. -5290.]\n",
      " [-5290. 10580. -5290.]\n",
      " [-5290. -5290. 10580.]]\n",
      "Pseudo-inverse of B:\n",
      "[[ 4.20079815e-05 -2.10039908e-05 -2.10039908e-05]\n",
      " [-2.10039908e-05  4.20079815e-05 -2.10039908e-05]\n",
      " [-2.10039908e-05 -2.10039908e-05  4.20079815e-05]]\n",
      "Incidence Matrix incidence:\n",
      "[[ 5290.     0. -5290.]\n",
      " [ 5290. -5290.     0.]\n",
      " [    0.  5290. -5290.]]\n",
      "Capacity Bounds:\n",
      "[[500]\n",
      " [ 47]\n",
      " [500]] [[-500]\n",
      " [ -47]\n",
      " [-500]]\n",
      "Shift Factor Matrix S:\n",
      "[[ 3.33333333e-01  7.80218988e-17 -3.33333333e-01]\n",
      " [ 3.33333333e-01 -3.33333333e-01  2.42522473e-17]\n",
      " [ 3.63885354e-17  3.33333333e-01 -3.33333333e-01]] (3, 3)\n",
      "[[1 0 1 0 0 1 0 0]\n",
      " [0 0 0 1 0 0 1 0]\n",
      " [0 1 0 0 1 0 0 1]] (3, 8)\n"
     ]
    }
   ],
   "source": [
    "# Modeling Contraints (S, c, \\Phi):\n",
    "\n",
    "# construct B in N x N \n",
    "B = np.zeros( ( len(N), len(N) ) ) # initialize B as zeors\n",
    "for i in N:  # go through each bus:\n",
    "    #print(i)\n",
    "    tbus_values = branch.loc[branch['fbus'] == i, 'tbus'].values  # find all the other buses connected to this bus and their sus value\n",
    "    #print(tbus_values)\n",
    "    for curr_to_bus in tbus_values:\n",
    "        curr_sus = branch.loc[(branch['fbus'] == i) & (branch['tbus'] == curr_to_bus), 'sus'].values[0]\n",
    "        #print(curr_sus)\n",
    "        B[i-1, curr_to_bus-1] = 0 - curr_sus\n",
    "        B[curr_to_bus-1, i-1] = 0 - curr_sus\n",
    "        \n",
    "    curr_index = i-1\n",
    "    curr_vector = B[curr_index]\n",
    "    curr_total = sum(curr_vector[:curr_index]) + sum(curr_vector[curr_index+1:]) # sum of all except curr_index\n",
    "    B[i-1, i-1] =  0 - curr_total # to be the sum of the remaining elements in   \n",
    "    #print(B) \n",
    "    #print(\"\")\n",
    "B_plus = np.linalg.pinv(B)   # Compute B^+ for computing downstream power-flows\n",
    "print(\"Susceptance Matrix:\")\n",
    "print(B) # in N x N \n",
    "print(\"Pseudo-inverse of B:\")\n",
    "print(B_plus)\n",
    "\n",
    "\n",
    "# construct Incidence in L x N\n",
    "ucb = []  \n",
    "lcb = []\n",
    "incidence = []\n",
    "for l in branch.index:\n",
    "    from_bus = branch.loc[l, 'fbus']\n",
    "    to_bus = branch.loc[l, 'tbus']\n",
    "    curr_sus = branch.loc[(branch['fbus'] == from_bus) & (branch['tbus'] == to_bus), 'sus'].values[0]\n",
    "    \n",
    "    vector = np.zeros( len(N) ) # Create vector for this branch\n",
    "    vector[from_bus - 1] =     curr_sus\n",
    "    vector[to_bus - 1]   = 0 - curr_sus\n",
    "    \n",
    "    incidence.append(vector)     \n",
    "    ucb.append(     branch.loc[l, 'ratea'] )\n",
    "    lcb.append( 0 - branch.loc[l, 'ratea'] )\n",
    "    \n",
    "incidence = np.array(incidence)        # convert to numpy\n",
    "ucb = np.array(ucb).reshape(-1, 1)  \n",
    "lcb = np.array(lcb).reshape(-1, 1)  \n",
    "\n",
    "print(\"Incidence Matrix incidence:\")\n",
    "print(incidence)\n",
    "print(\"Capacity Bounds:\")\n",
    "print(ucb, lcb)\n",
    "\n",
    "# S in L x N = incidence * B_plus (as per the paper)\n",
    "S = np.dot(incidence, B_plus)    \n",
    "print(\"Shift Factor Matrix S:\")\n",
    "print(S, S.shape)\n",
    "\n",
    "# construct \\Phi (hardcoded for now)\n",
    "#  -------------- x ---- y -------- v -------\n",
    "PHI = np.array([ [1, 0,  1,  0,  0, 1, 0, 0],\n",
    "                 [0, 0,  0,  1,  0, 0, 1, 0],\n",
    "                 [0, 1,  0,  0,  1, 0, 0, 1]  ])\n",
    "\n",
    "print(PHI, PHI.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e321051c",
   "metadata": {},
   "source": [
    "## DAM Optimization (Vectorized format):\n",
    "\n",
    "Objective function:\n",
    "$$\n",
    "\\min \\quad 0.5 p^T A p + b^T p\n",
    "$$\n",
    "\n",
    "Power balance constraint:\n",
    "$$\n",
    "\\quad \\quad \\quad \\quad \\quad 1^T p = 0 \\quad : \\lambda \\in \\mathbb{R}\n",
    "$$ \n",
    "\n",
    "Line constraints:\n",
    "$$\n",
    "\\quad \\quad \\quad \\quad \\quad -c \\leq S \\Phi p \\leq c \\quad : \\mu^-, \\mu^+ \\in \\mathbb{R}^{L}\n",
    "$$\n",
    "\n",
    "Generation limits:\n",
    "$$\n",
    "\\quad \\quad \\quad \\quad \\quad p^{min} \\leq p \\leq p^{max} \\in \\mathbb{R}^{G_{DAM}}\n",
    "$$\n",
    "\n",
    "the optimization variable are cleared energy of physical, defined as:\n",
    "$$\n",
    "p \\triangleq [x, y, v] ^ T\n",
    "$$\n",
    "\n",
    "and the dimensions are:\n",
    "$$\n",
    "A \\in \\mathbb{R}^{p_{dim}\\times p_{dim}},\\;b \\in \\mathbb{R}^{p_{dim}} \\quad p_{dim} = G_{DAM} + N + N\n",
    "$$\n",
    "$$x\\in \\mathbb{R}^{G_{DAM}} \\;- \\text{Physical Supply Variable} $$ \n",
    "$$y\\in \\mathbb{R}^{N} \\;- \\text{Physical Load Constant} $$ \n",
    "$$v\\in \\mathbb{R}^{N} \\;- \\text{Virtual Supply bid Constant} $$\n",
    "$$S \\in \\mathbb{R}^{L \\times N} \\;-\\quad \\text{Shift Factor Matrix}$$\n",
    "$$\\Phi \\in \\mathbb{R}^{N \\times p_{dim}} $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "12059226",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/opt/homebrew/Caskroom/miniforge/base/envs/PGOC/lib/python3.9/site-packages/scipy/__init__.py:146: UserWarning: A NumPy version >=1.16.5 and <1.23.0 is required for this version of SciPy (detected version 1.26.4\n",
      "  warnings.warn(f\"A NumPy version >={np_minversion} and <{np_maxversion}\"\n"
     ]
    }
   ],
   "source": [
    "import cvxpy as cp\n",
    "\n",
    "def DAM_dcopf_vectorized(gen, branch, gencost, bus, S, ucb, lcb):\n",
    "    \n",
    "    # Define base MVA (for conversion between per unit and MWs)\n",
    "    baseMVA = gen['mbase'].iloc[0]\n",
    "    \n",
    "    x = cp.Variable( len(G_DAM) )      # dynamic based on generator locations\n",
    "    y = 0 - bus[\"pd\"].to_numpy()   # dynamically read load data from dataframe - negative for load values\n",
    "    v = np.array([0, SUPPLY_BID, 0])        # manually change to add virtual bids\n",
    "        \n",
    "    # decision variable in MWs\n",
    "    p = cp.Variable( dim_p )\n",
    "\n",
    "    # Define constraints on the concatenated variable - dim_p\n",
    "    constraints = [ p[0:len(G_DAM)] == x,\n",
    "                    p[len(G_DAM):len(G_DAM)+len(y)] == y,\n",
    "                    p[len(G_DAM)+len(y) :] == v       ]\n",
    "    \n",
    "    objective = cp.Minimize( 0.5 * cp.quad_form(p, A) + cp.matmul(b, p) ) # 1a\n",
    "    \n",
    "    # constraint - 1 (equality constraints)\n",
    "    ones_vector = np.ones( dim_p )  # A vector of ones with the same length as p\n",
    "    constraints += [ cp.matmul(ones_vector, p) == 0 ]  # Adding the balance constraint\n",
    "    \n",
    "    def find_flows(p):\n",
    "        p_normalized = p / baseMVA               # in per unit\n",
    "        pp_vector = np.matmul(PHI, p_normalized)\n",
    "        flows = np.matmul(S, pp_vector)           # power-flow vector of size L - still in per unit\n",
    "        flows_scaled = baseMVA * flows             # power-flow vector of size L - in MWs\n",
    "        flows_scaled = np.reshape(flows_scaled, (-1, 1))\n",
    "        return flows_scaled\n",
    "    \n",
    "    # constraint - 2 (inequality capacity constraints)\n",
    "    flows_scaled = find_flows(p)\n",
    "    constraints.append(\n",
    "        flows_scaled <= ucb\n",
    "    )\n",
    "    constraints.append(\n",
    "        flows_scaled >= lcb\n",
    "    )\n",
    "    \n",
    "    # constraint - 3 (inequality Max, Min generation constraints)\n",
    "    for g in G_DAM: # 1, 2\n",
    "        constraints.append( p[g-1] <= gen.loc[g-1, 'pmax'] ) # G1-g[0] AND G3-g[2] only \n",
    "    for g in G_DAM: # 1, 2\n",
    "        constraints.append( p[g-1] >= gen.loc[g-1, 'pmin'] ) # G1-g[0] AND G3-g[2] only\n",
    "        \n",
    "    # Define the problem and solve\n",
    "    problem = cp.Problem(objective, constraints)\n",
    "    problem.solve(solver=cp.OSQP)     # is an open-source C library for solving convex quadratic programs\n",
    "    #problem.solve(solver=cp.XPRESS ) \n",
    "    \n",
    "    # Access the optimized values\n",
    "    optimized_p = p.value\n",
    "    #print(optimized_p)\n",
    "\n",
    "    # Prepare generation output\n",
    "    optimized_x = optimized_p[0:len(G_DAM)]\n",
    "    generation = pd.DataFrame({\n",
    "        'id':   gen['id'],\n",
    "        'node': gen['bus'],\n",
    "        'gen':  optimized_x\n",
    "    })\n",
    "    \n",
    "    # prepare power-flow output\n",
    "    optimized_flows = find_flows(optimized_p)\n",
    "    flows = pd.DataFrame({\n",
    "        'fbus': branch['fbus'],\n",
    "        'tbus': branch['tbus'],\n",
    "        'flow': np.squeeze(optimized_flows)\n",
    "    })\n",
    "    \n",
    "    # prepare prices\n",
    "#     for constraint in constraints[3 : 5+1]:\n",
    "#         print(constraint.dual_value) \n",
    "        \n",
    "    llambda  = abs( constraints[3].dual_value ) # not sure if abs is correct or flipping?\n",
    "    mu_plus  = constraints[4].dual_value\n",
    "    mu_minus = constraints[5].dual_value\n",
    "    mu = mu_plus - mu_minus\n",
    "    ones_vector = np.ones( len(N) ).reshape(-1, 1)\n",
    "    pi = llambda * ones_vector - np.matmul(S.T, mu)\n",
    "    # Extract the prices (dual values of the balance constraints)\n",
    "    prices = pd.DataFrame({\n",
    "        'node': bus['bus_i'],\n",
    "        'value': np.squeeze(pi),\n",
    "        '  mu_plus': np.squeeze(mu_plus),\n",
    "        '  mu_minus': np.squeeze(mu_minus),\n",
    "        '  mu(dual)': np.squeeze(mu)\n",
    "    })\n",
    "            \n",
    "    # Return the solution and objective value\n",
    "    return {\n",
    "        'generation': generation.round(7),\n",
    "        'flows': flows.round(5),\n",
    "        'prices': prices,\n",
    "        'p': optimized_p,\n",
    "        'cost': '{:.5f}'.format(problem.value),\n",
    "        'status': problem.status\n",
    "    }"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "6d6d854c",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Day Ahead Market(DAM) Results:\n",
      "\n",
      "generation:\n",
      "   id  node        gen\n",
      "0   1     1  59.090909\n",
      "1   2     3  15.909091\n",
      "\n",
      "flows:\n",
      "   fbus  tbus      flow\n",
      "0     1     3  14.39394\n",
      "1     1     2  44.69697\n",
      "2     2     3 -30.30303\n",
      "\n",
      "prices:\n",
      "   node      value    mu_plus    mu_minus    mu(dual)\n",
      "0     1  20.727273        0.0         0.0         0.0\n",
      "1     2  20.727273        0.0         0.0         0.0\n",
      "2     3  20.727273        0.0         0.0         0.0\n",
      "\n",
      "p:\n",
      "[ 5.90909091e+01  1.59090909e+01 -1.60405221e-20 -7.50000000e+01\n",
      " -1.60405221e-20 -1.60405221e-20 -1.60405221e-20 -1.60405221e-20]\n",
      "\n",
      "cost:\n",
      "929.54545\n",
      "\n",
      "status:\n",
      "optimal\n",
      "\n",
      "optimized_x = [59.0909091 15.9090909]\n"
     ]
    }
   ],
   "source": [
    "print(\"Day Ahead Market(DAM) Results:\\n\")\n",
    "\n",
    "DAM = DAM_dcopf_vectorized(DAM_gen, branch, DAM_gencost, DAM_bus, S, ucb, lcb)\n",
    "for key, value in DAM.items():\n",
    "    if key==\"generation\":\n",
    "        optimized_x = value[\"gen\"].to_numpy()\n",
    "    if key==\"p\":\n",
    "        p = value\n",
    "        #p = np.where(np.abs(p) < 1e-2, 0, p)\n",
    "    if key==\"prices\":\n",
    "        mu = value[\"  mu(dual)\"].to_numpy()\n",
    "    print(f\"{key}:\") \n",
    "    print(value)\n",
    "    print()\n",
    "    \n",
    "print(f\"optimized_x = {optimized_x}\") # optimized vector to be used in DAM"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "51375eec",
   "metadata": {},
   "source": [
    "### Post Optimality Analysis:\n",
    "\n",
    "$$p_{-i} \\triangleq [x, y, v_{-i}] ^ T \\quad \\quad p_{-i_{dim}}$$\n",
    "$$\\bar D \\in \\mathbb{R}^{\\bar L \\times L} \\;$$\n",
    "$$S \\in \\mathbb{R}^{L \\times N}$$\n",
    "$$\\Phi_{-i} \\in \\mathbb{R}^{N \\times {p_{-i}}_{dim}}$$\n",
    "$$A_{-i} \\in \\mathbb{R}^{{p_{-i}}_{dim} \\times {p_{-i}}_{dim}} $$\n",
    " \n",
    "$$ \\bar D S \\in \\mathbb{R}^{\\bar L \\times N}$$\n",
    "$$X = \\bar D S \\Phi_{-i}\\in \\mathbb{R}^{\\bar L \\times {p_{-i}}_{dim}} \\;$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "ce6ea542",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "p_minusi:\n",
      "[ 5.90909091e+01  1.59090909e+01 -1.60405221e-20 -7.50000000e+01\n",
      " -1.60405221e-20 -1.60405221e-20 -1.60405221e-20]\n",
      "PHI_minusi:\n",
      "[[1 0 1 0 0 1 0]\n",
      " [0 0 0 1 0 0 0]\n",
      " [0 1 0 0 1 0 1]]\n",
      "[ 59.09090909 -75.          15.90909091]\n",
      "flows:\n",
      "[ 14.39393939  44.6969697  -30.3030303 ]\n",
      "Dbar: [] with shape (0, 3)\n",
      "X: [] of Shape (0, 7)\n",
      "(1, 7) (7, 1)\n",
      "H: [[4.58333333]] of Shape (1, 1)\n",
      "H_inverse: [[0.21818182]] of Shape (1, 1)\n",
      "[[-1.]\n",
      " [-1.]\n",
      " [-1.]] (3, 1)\n",
      "[[1.]] (1, 1)\n",
      "\n",
      "DAM_sesitivity:\n",
      " [[-0.21818182]\n",
      " [-0.21818182]\n",
      " [-0.21818182]]\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "<ipython-input-8-319e7463f7e4>:24: RuntimeWarning: divide by zero encountered in divide\n",
      "  LAMBDA = np.where(A_minusi != 0, 1 / A_minusi, 0)\n"
     ]
    }
   ],
   "source": [
    "p_minusi   = np.concatenate((p[:-2], p[-1:]))\n",
    "p_minusi_dim = len(p_minusi)\n",
    "PHI_minusi = np.delete(PHI, -2, axis=1)\n",
    "A_minusi = np.delete(np.delete(A, -2, axis=0), -2, axis=1) # delete row and column\n",
    "\n",
    "ppvector_minusi = np.matmul(PHI_minusi, p_minusi)\n",
    "flows = np.matmul(S, ppvector_minusi)\n",
    "\n",
    "print(f\"p_minusi:\\n{p_minusi}\")\n",
    "print(f\"PHI_minusi:\\n{PHI_minusi}\")\n",
    "print(ppvector_minusi)\n",
    "print(f\"flows:\\n{flows}\")\n",
    "\n",
    "Lbar = np.count_nonzero(mu)  # how many non-zero entries in mu? - Lbar\n",
    "nonzero_indices_mu = np.nonzero(mu) # find the index of non-zero in mu? - using which construct \n",
    "L = incidence.shape[0]\n",
    "Dbar = np.zeros((Lbar, L))\n",
    "Dbar[np.arange(Lbar), nonzero_indices_mu[0]] = 1\n",
    "print(f\"Dbar: {Dbar} with shape {Dbar.shape}\")\n",
    "\n",
    "X = np.matmul(Dbar, np.matmul(S, PHI_minusi))\n",
    "print(f\"X: {X} of Shape {X.shape}\")\n",
    "\n",
    "LAMBDA = np.where(A_minusi != 0, 1 / A_minusi, 0)\n",
    "#print(LAMBDA, LAMBDA.shape)\n",
    "\n",
    "ones_vector_p = np.ones( p_minusi_dim ).reshape(-1, 1) # column vector\n",
    "H_left  = np.vstack( (ones_vector_p.T,  X) )\n",
    "H_right = np.vstack( (ones_vector_p.T, -X) ).T\n",
    "print(H_left.shape, H_right.shape)\n",
    "H = np.matmul(H_left, np.matmul(LAMBDA, H_right))\n",
    "print(f\"H: {H} of Shape {H.shape}\")\n",
    "\n",
    "H_inverse = np.linalg.inv(H)\n",
    "print(f\"H_inverse: {H_inverse} of Shape {H_inverse.shape}\")\n",
    "\n",
    "ones_vector_N = np.ones( len(N) ).reshape(-1, 1)        # column vector\n",
    "temp3 = np.vstack( (-ones_vector_N.T, np.matmul(Dbar, S)) ).T\n",
    "print(temp3, temp3.shape)\n",
    "\n",
    "zeros_vector_Lbar = np.zeros( Lbar ).reshape(-1, 1)      # column vector\n",
    "DAM_sesitivity_right = np.vstack( (1, zeros_vector_Lbar) )\n",
    "print(DAM_sesitivity_right, DAM_sesitivity_right.shape)\n",
    "\n",
    "DAM_sesitivity = np.matmul(temp3, np.matmul(H_inverse, DAM_sesitivity_right))\n",
    "\n",
    "print(f\"\\nDAM_sesitivity:\\n {DAM_sesitivity}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "061c437e",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "c517baf0",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Real Time Market(RTM) Data:\n",
      "\n",
      "Updated Bus Dataframe:\n",
      "   bus_i  type  pd     qd  gs  bs  area  vm  va  basekv  zone  vmax  vmin\n",
      "0      1     2   0   0.00   0   0     1   1   0     230     1   1.1   0.9\n",
      "1      2     2  90   0.00   0   0     1   1   0     230     1   1.1   0.9\n",
      "2      3     1   0  98.61   0   0     1   1   0     230     1   1.1   0.9\n",
      "\n",
      "Updated RTM_gen:\n",
      "   bus   pg  qg   qmax   qmin  vg  mbase  status  pmax  pmin  ...  qc1min  \\\n",
      "0    1   40   0   30.0  -30.0   1    100       1  1000     0  ...       0   \n",
      "1    2  170   0  127.5 -127.5   1    100       1  1000     0  ...       0   \n",
      "2    3  100   0    0.0    0.0   0    100       1  1000     0  ...       0   \n",
      "\n",
      "   qc1max  qc2min  qc2max  ramp_agc  ramp_10  ramp_30  ramp_q  apf  id  \n",
      "0       0       0       0         0        0        0       0    0   1  \n",
      "1       0       0       0         0        0        0       0    0   2  \n",
      "2       0       0       0         0        0        0       0    0   3  \n",
      "\n",
      "[3 rows x 22 columns]\n",
      "\n",
      "Updated RTM_gencost:\n",
      "   model  startup  shutdown  n   x1  y1  id\n",
      "0      2        0         0  2  1.8  10   1\n",
      "1      2        0         0  2  1.7   5   2\n",
      "2      2        0         0  2  0.1  14   3\n",
      "\n",
      "Updated RTM_branch:\n",
      "   fbus  tbus        r         x        b  ratea  rateb  ratec  ratio  angle  \\\n",
      "0     1     3  0.00281  0.000189  0.00712    500    500    500      0      0   \n",
      "1     1     2  0.00281  0.000189  0.00712     47    500    500      0      0   \n",
      "2     2     3  0.00281  0.000189  0.00712    500    500    500      0      0   \n",
      "\n",
      "   status  angmin  angmax  id     sus  \n",
      "0       1    -360     360   1  5290.0  \n",
      "1       1    -360     360   2  5290.0  \n",
      "2       1    -360     360   3  5290.0  \n"
     ]
    }
   ],
   "source": [
    "print(\"Real Time Market(RTM) Data:\\n\")\n",
    "\n",
    "datadir = 'opf_data'\n",
    "RTM_bus, RTM_gen, RTM_gencost, _ = input_data(datadir)\n",
    "\n",
    "print(\"Updated Bus Dataframe:\")\n",
    "RTM_bus[\"pd\"] = [LOAD1, RT_LOAD2, 0]  # remember this is positive. Keep in mind the sign when using later.\n",
    "print(RTM_bus)\n",
    "print(\"\")\n",
    "\n",
    "print(\"Updated RTM_gen:\")\n",
    "# Create a DataFrame with the new row data\n",
    "new_row = pd.DataFrame({'bus': [3], 'pg': [100], 'qg': [0], 'qmax': [0], 'qmin': [0], 'vg': [0], 'mbase': [100], 'status': [1], 'pmax': [1000], 'pmin': [0], 'qc1max': [0], 'qc1min': [0], 'qc2max': [0], 'qc2min': [0], 'ramp_agc': [0], 'ramp_10': [0], 'ramp_30': [0], 'ramp_q': [0], 'apf': [0], 'id': [3]})\n",
    "RTM_gen = pd.concat([RTM_gen, new_row], ignore_index=True) # Append the new row to RTM_gen\n",
    "print(RTM_gen)\n",
    "print(\"\")\n",
    "\n",
    "print(\"Updated RTM_gencost:\")\n",
    "last_row = RTM_gencost.iloc[-1] # Get the last row of the DataFrame\n",
    "RTM_gencost = pd.concat([RTM_gencost, pd.DataFrame([last_row])], ignore_index=True) # Concatenate the original DataFrame with the duplicated row\n",
    "RTM_gencost[\"x1\"] = [1.8, 1.7,  0.1]\n",
    "RTM_gencost[\"y1\"] = [ 10,   5,   14]\n",
    "RTM_gencost[\"id\"] = [  1,   2,    3]\n",
    "print(RTM_gencost)\n",
    "print(\"\")\n",
    "\n",
    "print(\"Updated RTM_branch:\")\n",
    "print(branch)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "88879061",
   "metadata": {},
   "source": [
    "### Construct all the relevant matrices from the above dataframes for optimization problem formulation (Overall approach to Construction):\n",
    "\n",
    "#### Optimization variable:\n",
    "- Ensure the dimension of `z` is fixed.\n",
    "    \n",
    "#### Objective Function:\n",
    "- (Done) `C`: Coefficient matrix for the constraints\n",
    "- (Done) `d`: Right-hand side vector for the constraints\n",
    "\n",
    "#### Modeling Constraints:\n",
    "- `l` for \n",
    "- `\\Psi` for `x` from day ahead\n",
    "- `\\Theta` for `z` from day ahead\n",
    "- `\\Omega` for `l` from day ahead"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "74ed05ed",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[1.8 0.  0. ]\n",
      " [0.  1.7 0. ]\n",
      " [0.  0.  0.1]]\n",
      "[10.  5. 14.]\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "\n",
    "# Modeling Objective Function constants\n",
    "G_RTM = RTM_gen['id'].values      # 1, 2, 3\n",
    "\n",
    "dim_x = len(optimized_x)\n",
    "dim_z = len(G_RTM)\n",
    "\n",
    "# Modeling Objective Function Constants (A, b):\n",
    "C = np.zeros( ( dim_z, dim_z ) ) # initialize A as zeros # make the genreators dynamic\n",
    "d = np.zeros( dim_z )\n",
    "\n",
    "for i in G_RTM:  # go through each generator\n",
    "    alpha_i  = RTM_gencost.loc[RTM_gencost['id']==i, 'x1'].values[0]   # alpha_i associated \n",
    "    beta_i   = RTM_gencost.loc[RTM_gencost['id']==i, 'y1'].values[0]   #  beta_i associated\n",
    "    bus_index = RTM_gen.loc[RTM_gen['id']==i, 'bus'].values[0]  # bus index associated with curr generator\n",
    "    \n",
    "    C[    i-1,     i-1] = alpha_i   # for physical supply (0, 1, 2)\n",
    "    d[    i-1] = beta_i\n",
    "\n",
    "print( C )\n",
    "print( d )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "bda91642",
   "metadata": {},
   "outputs": [],
   "source": [
    "# for x\n",
    "PSI = np.array([ [1, 0],\n",
    "                 [0, 0],\n",
    "                 [0, 1]  ])\n",
    "# for z\n",
    "THETA = np.array([  [1, 0, 0],\n",
    "                    [0, 1, 0],\n",
    "                    [0, 0, 1]  ])\n",
    "# for l\n",
    "OMEGA = np.array([  [1, 0, 0],\n",
    "                    [0, 1, 0],\n",
    "                    [0, 0, 1]  ])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "49e44e2d",
   "metadata": {},
   "source": [
    "## Real-Time Market Optimization (Vectorized format):\n",
    "\n",
    "$$\n",
    "\\min \\quad 0.5 z^T C z + d^T z\n",
    "$$\n",
    "\n",
    "Power balance constraint:\n",
    "$$\n",
    "1^T z + 1^T x = 1^T l \\quad : \\delta\n",
    "$$\n",
    "\n",
    "Line constraints:\n",
    "$$\n",
    "-c \\leq S(\\Psi x + \\Theta z - \\Omega l) \\leq c \\quad : \\eta^-, \\eta^+\n",
    "$$\n",
    "\n",
    "Generation limits:\n",
    "$$\n",
    "z^{min} \\leq z \\leq z^{max}\n",
    "$$\n",
    "\n",
    "the optimization variable are adjustment to physical energy:\n",
    "$$\n",
    "z \n",
    "$$\n",
    "\n",
    "and the dimensions are:\n",
    "$$\n",
    "C \\in \\mathbb{R}^{z_{dim}\\times z_{dim}},\\;d \\in \\mathbb{R}^{z_{dim}} \\quad z_{dim} = G_{RTM}\n",
    "$$\n",
    "$$z\\in \\mathbb{R}^{G_{RTM}} \\;- \\text{Physical Supply Adjustment Variable} $$ \n",
    "$$\\Theta \\in \\mathbb{R}^{N \\times G_{RTM}} $$\n",
    "\n",
    "$$x\\in \\mathbb{R}^{G_{DAM}} \\;- \\text{Physical Supply Variable} $$ \n",
    "$$\\Psi \\in \\mathbb{R}^{N \\times G_{DAM}} $$\n",
    "\n",
    "$$l\\in \\mathbb{R}^{N} \\;- \\text{Physical Supply Variable} $$ \n",
    "$$\\Omega \\in \\mathbb{R}^{N \\times N} $$\n",
    "\n",
    "$$S \\in \\mathbb{R}^{L \\times N} \\;-\\quad \\text{Shift Factor Matrix}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "7c9c22b3",
   "metadata": {},
   "outputs": [],
   "source": [
    "import cvxpy as cp\n",
    "\n",
    "def RTM_dcopf_vectorized(gen, branch, gencost, bus, S, ucb, lcb, optimized_x):\n",
    "    \n",
    "    # Define base MVA (for conversion between per unit and MWs)\n",
    "    baseMVA = gen['mbase'].iloc[0]\n",
    "    \n",
    "    # decision variable in MWs\n",
    "    z = cp.Variable( len(G_RTM) )      # dynamic based on generator locations\n",
    "    load_RTM = bus[\"pd\"].to_numpy()    # dynamically read load data from dataframe - negative for load values\n",
    "    dim_load = len(load_RTM)\n",
    "            \n",
    "    objective = cp.Minimize( 0.5 * cp.quad_form(z, C) + cp.matmul(d, z) ) # 6a\n",
    "    \n",
    "    # constraint - 1 (equality constraints)\n",
    "    ones_vector_z = np.ones( dim_z )\n",
    "    ones_vector_x = np.ones( dim_x )\n",
    "    ones_load_RTM = np.ones( dim_load )\n",
    "    constraints   = [ cp.matmul(ones_vector_z, z) + cp.matmul(ones_vector_x, optimized_x) == cp.matmul(ones_load_RTM, load_RTM) ]  # Adding the balance constraint\n",
    "    \n",
    "    def find_flows(z, optimized_x, load_RTM):\n",
    "        z_normalized = z / baseMVA               # in per unit\n",
    "        x_normalized = optimized_x / baseMVA               # in per unit\n",
    "        load_normalized = load_RTM / baseMVA               # in per unit\n",
    "        \n",
    "        pp_vector  = np.matmul(THETA,  z_normalized) \n",
    "        pp_vector += np.matmul(PSI,    x_normalized)\n",
    "        pp_vector  -= np.matmul(OMEGA, load_normalized)\n",
    "        flows = np.matmul(S, pp_vector)           # power-flow vector of size L - still in per unit\n",
    "        flows_scaled = baseMVA * flows            # power-flow vector of size L - in MWs\n",
    "        flows_scaled = np.reshape(flows_scaled, (-1, 1))\n",
    "        return flows_scaled\n",
    "    \n",
    "    # constraint - 2 (inequality capacity constraints)\n",
    "    flows_scaled = find_flows(z, optimized_x, load_RTM)\n",
    "    constraints.append(\n",
    "        flows_scaled <= ucb\n",
    "    )\n",
    "    constraints.append(\n",
    "        flows_scaled >= lcb\n",
    "    )\n",
    "    \n",
    "    # constraint - 3 (inequality Max, Min generation constraints)\n",
    "    for g in G_RTM: # 1, 2\n",
    "        constraints.append( z[g-1] <= gen.loc[g-1, 'pmax'] ) # G1-g[0] AND G3-g[2] only\n",
    "        \n",
    "    # Define the problem and solve\n",
    "    problem = cp.Problem(objective, constraints)\n",
    "    problem.solve(solver=cp.OSQP)     # is an open-source C library for solving convex quadratic programs\n",
    "    #problem.solve(solver=cp.XPRESS ) \n",
    "    \n",
    "    # Access the optimized values\n",
    "    optimized_z = z.value\n",
    "    #print(optimized_z)\n",
    "\n",
    "    # Prepare generation output\n",
    "    generation = pd.DataFrame({\n",
    "        'id':   gen['id'],\n",
    "        'node': gen['bus'],\n",
    "        'gen':  optimized_z\n",
    "    })\n",
    "    \n",
    "    # prepare power-flow output\n",
    "    optimized_flows = find_flows(optimized_z, optimized_x, load_RTM)\n",
    "    flows = pd.DataFrame({\n",
    "        'fbus': branch['fbus'],\n",
    "        'tbus': branch['tbus'],\n",
    "        'flow': np.squeeze(optimized_flows)\n",
    "    })\n",
    "    \n",
    "    # prepare prices\n",
    "#     for constraint in constraints[3 : 5+1]:\n",
    "#         print(constraint.dual_value) \n",
    "        \n",
    "    delta  = abs( constraints[0].dual_value )\n",
    "    eta_plus  = constraints[1].dual_value\n",
    "    eta_minus = constraints[2].dual_value\n",
    "    eta = eta_plus - eta_minus\n",
    "    ones_vector = np.ones( len(N) ).reshape(-1, 1)\n",
    "    sigma = delta * ones_vector - np.matmul(S.T, eta)\n",
    "    # Extract the prices (dual values of the balance constraints)\n",
    "    prices = pd.DataFrame({\n",
    "        'node': bus['bus_i'],\n",
    "        'value': np.squeeze(sigma),\n",
    "        '  eta_plus': np.squeeze(eta_plus),\n",
    "        '  eta_minus': np.squeeze(eta_minus),\n",
    "        '  eta(dual)': np.squeeze(eta)\n",
    "        \n",
    "    })\n",
    "            \n",
    "    # Return the solution and objective value\n",
    "    return {\n",
    "        'generation': generation.round(7),\n",
    "        'flows': flows.round(5),\n",
    "        'prices': prices,\n",
    "        'z': optimized_z,\n",
    "        'cost': '{:.5f}'.format(problem.value),\n",
    "        'status': problem.status\n",
    "    }"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "3b109e4d",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Real Time Market(RTM):\n",
      "\n",
      "generation:\n",
      "   id  node       gen\n",
      "0   1     1  0.160839\n",
      "1   2     2  8.251748\n",
      "2   3     3  6.587413\n",
      "\n",
      "flows:\n",
      "   fbus  tbus      flow\n",
      "0     1     3  12.25175\n",
      "1     1     2  47.00000\n",
      "2     2     3 -34.74825\n",
      "\n",
      "prices:\n",
      "   node      value    eta_plus    eta_minus    eta(dual)\n",
      "0     1  10.289510    0.000000          0.0     0.000000\n",
      "1     2  19.027972   13.107692          0.0    13.107692\n",
      "2     3  14.658741    0.000000          0.0     0.000000\n",
      "\n",
      "z:\n",
      "[0.16083916 8.25174826 6.58741259]\n",
      "\n",
      "cost:\n",
      "195.16154\n",
      "\n",
      "status:\n",
      "optimal\n",
      "\n"
     ]
    }
   ],
   "source": [
    "print(\"Real Time Market(RTM):\\n\")\n",
    "\n",
    "RTM = RTM_dcopf_vectorized(RTM_gen, branch, RTM_gencost, RTM_bus, S, ucb, lcb, optimized_x)\n",
    "\n",
    "for key, value in RTM.items():\n",
    "    if key==\"prices\":\n",
    "        eta = value[\"  eta(dual)\"].to_numpy()\n",
    "    print(f\"{key}:\") \n",
    "    print(value)\n",
    "    print()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0f221c6a",
   "metadata": {},
   "source": [
    "### Post Optimality Analysis (VB Sensitivity):\n",
    "\n",
    "$$\\bar R \\in \\mathbb{R}^{\\bar L \\times L} \\;$$\n",
    "$$S \\in \\mathbb{R}^{L \\times N}$$\n",
    "$$\\hat \\Theta \\in \\mathbb{R}^{N \\times G_{RTM}} $$\n",
    " \n",
    "$$ \\bar R S \\in \\mathbb{R}^{\\bar L \\times N}$$\n",
    "$$Y = \\bar R S \\hat \\Theta\\in \\mathbb{R}^{\\bar L \\times G_{RTM}} \\;$$\n",
    "\n",
    "$$ \\bar R S \\in \\mathbb{R}^{\\bar L \\times N}$$\n",
    "$$\\Psi \\in \\mathbb{R}^{N \\times G_{DAM}} $$\n",
    "$$ \\bar R S \\Psi \\in \\mathbb{R}^{\\bar L \\times G_{DAM}}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "d2e189f8",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Rbar: [[0. 1. 0.]] of shape (1, 3)\n",
      "Y: [[ 3.33333333e-01 -3.33333333e-01  2.42522473e-17]] of Shape (1, 3)\n",
      "[[ 0.55555556  0.          0.        ]\n",
      " [ 0.          0.58823529  0.        ]\n",
      " [ 0.          0.         10.        ]] (3, 3)\n",
      "(2, 3) (3, 2)\n",
      "E: [[ 1.11437908e+01  1.08932462e-02]\n",
      " [-1.08932462e-02 -1.27087872e-01]] of Shape (2, 2)\n",
      "E_inverse: [[ 8.97435897e-02  7.69230769e-03]\n",
      " [-7.69230769e-03 -7.86923077e+00]] of Shape (2, 2)\n",
      "[[ 1.00000000e+00 -3.33333333e-01]\n",
      " [ 1.00000000e+00  3.33333333e-01]\n",
      " [ 1.00000000e+00 -2.42522473e-17]] (3, 2)\n",
      "[[1.00000000e+00 1.00000000e+00]\n",
      " [3.33333333e-01 2.42522473e-17]] (2, 2)\n",
      "(3, 2)\n",
      "(7, 1)\n",
      "\n",
      "RTM_sesitivity:\n",
      " [[ 0.73006993]\n",
      " [-0.54685315]\n",
      " [ 0.09160839]]\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "<ipython-input-14-8387614e4752>:10: RuntimeWarning: divide by zero encountered in divide\n",
      "  GAMMA = np.where(C != 0, 1 / C, 0)\n"
     ]
    }
   ],
   "source": [
    "Lbar = np.count_nonzero(eta)  # how many non-zero entries in mu? - Lbar\n",
    "nonzero_indices_eta = np.nonzero(eta) # find the index of non-zero in mu? - using which construct \n",
    "Rbar = np.zeros((Lbar, L))\n",
    "Rbar[np.arange(Lbar), nonzero_indices_eta[0]] = 1\n",
    "print(f\"Rbar: {Rbar} of shape {Rbar.shape}\")\n",
    "\n",
    "Y = np.matmul(Rbar, np.matmul(S, THETA))\n",
    "print(f\"Y: {Y} of Shape {Y.shape}\")\n",
    "\n",
    "GAMMA = np.where(C != 0, 1 / C, 0)\n",
    "print(GAMMA, GAMMA.shape)\n",
    "\n",
    "ones_vector_z = np.ones( dim_z ).reshape(-1, 1)\n",
    "temp1 = np.vstack((ones_vector_z.T, Y))\n",
    "temp2 = np.vstack((ones_vector_z.T, -Y)).T\n",
    "print(temp1.shape, temp2.shape)\n",
    "E = np.matmul(temp1, np.matmul(GAMMA, temp2))\n",
    "print(f\"E: {E} of Shape {E.shape}\")\n",
    "\n",
    "E_inverse = np.linalg.inv(E)\n",
    "print(f\"E_inverse: {E_inverse} of Shape {E_inverse.shape}\")\n",
    "\n",
    "\n",
    "ones_vector_N = np.ones( len(N) ).reshape(-1, 1)\n",
    "temp3 = np.vstack( (ones_vector_N.T, -np.matmul(Rbar, S)) ).T\n",
    "print(temp3, temp3.shape)\n",
    "\n",
    "ones_vector_DAM = np.ones( len(G_DAM) ).reshape(-1, 1)\n",
    "temp4 = np.vstack( (ones_vector_DAM.T, np.matmul(np.matmul(Rbar, S), PSI)) )\n",
    "print(temp4, temp4.shape)\n",
    "\n",
    "firstPart =  np.matmul(np.matmul(temp3, E_inverse), temp4)\n",
    "print(firstPart.shape)\n",
    "secondPart = np.matmul(np.matmul(LAMBDA, H_right), np.matmul(H_inverse, DAM_sesitivity_right))\n",
    "print(secondPart.shape)\n",
    "\n",
    "Khat = np.array([[1, 0, 0, 0, 0, 0, 0],\n",
    "                 [0, 1, 0, 0, 0, 0, 0]])  # Example K matrix\n",
    "\n",
    "RTM_sesitivity = np.matmul(firstPart, np.matmul(Khat, secondPart))\n",
    "print(f\"\\nRTM_sesitivity:\\n {RTM_sesitivity}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "587a2111",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-0.94825175],\n",
       "       [ 0.32867133],\n",
       "       [-0.30979021]])"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sesitivity_VB = DAM_sesitivity - RTM_sesitivity\n",
    "sesitivity_VB"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "8392343a",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Bus 1: 20.727272727272727, 10.28951048153846, 10.43776\n",
      "Bus 2: 20.727272727272727, 19.027972035897502, 1.69930\n",
      "Bus 3: 20.727272727272727, 14.658741258717981, 6.06853\n"
     ]
    }
   ],
   "source": [
    "RTM['prices']['value'] = pd.to_numeric(RTM['prices']['value'])\n",
    "DAM['prices']['value'] = pd.to_numeric(DAM['prices']['value'])\n",
    "\n",
    "count = 1\n",
    "for dam_price, rtm_price in zip(DAM['prices']['value'], RTM['prices']['value']):\n",
    "    print(f\"Bus {count}: {dam_price}, {rtm_price}, {dam_price - rtm_price:.5f}\")\n",
    "    if count == 2:\n",
    "        if SUPPLY_BID>0:\n",
    "            print(f\"Bus-2 Profit: {(dam_price - rtm_price)*SUPPLY_BID - (LOAD_ATTACK)*rtm_price}\")\n",
    "    count += 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7528473a",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.18"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}