Market 3 bus Iterative (Reference Bus).ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "7d32466e",
   "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": "f244f201",
   "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": "markdown",
   "id": "c0f55eb1",
   "metadata": {},
   "source": [
    "## DAM Optimization:\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 \n",
    "$$ \n",
    "\n",
    "Line constraints:\n",
    "$$\n",
    "\\quad \\quad \\quad \\quad \\quad -c \\leq S \\Phi p \\leq c \\quad : \\mu^-, \\mu^+\n",
    "$$\n",
    "\n",
    "Generation limits:\n",
    "$$\n",
    "\\quad \\quad \\quad \\quad \\quad p^{min} \\leq p \\leq p^{max}\n",
    "$$\n",
    "\n",
    "where the optimization variables are cleared energy of physical and virtual supplies and loads in the system, defined as:\n",
    "\n",
    "$$\n",
    "p \\triangleq [x \\quad y \\quad v \\quad w]^T\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "ec005504",
   "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 os\n",
    "import copy\n",
    "import cvxpy as cp\n",
    "import pandas as pd\n",
    "\n",
    "def dcopf(gen, branch, gencost, bus, supply_vb):\n",
    "    \n",
    "    # Define sets\n",
    "    G = gen['id'].values\n",
    "    N = bus['bus_i'].values\n",
    "\n",
    "    # Define base MVA\n",
    "    baseMVA = gen['mbase'].iloc[0]\n",
    "\n",
    "    # Decision variables\n",
    "    GEN = cp.Variable(   len(G) ) # Generation\n",
    "    THETA = cp.Variable( len(N) )             # Voltage phase angle of bus\n",
    "    FLOW = cp.Variable( (len(N), len(N)) )    # Flows between all pairs of nodes\n",
    "    \n",
    "    # Objective function\n",
    "    objective = cp.Minimize( 0.5 * cp.sum(cp.multiply(gencost['x1'], cp.square(GEN))) + cp.sum(cp.multiply(gencost['y1'], GEN)) )\n",
    "    #objective = cp.Minimize( cp.sum( cp.multiply(gencost['x1'], GEN) ) )\n",
    "\n",
    "    # Define the constraints\n",
    "    constraints = []\n",
    "\n",
    "    # Slack bus reference angle constraint\n",
    "    constraints.append(THETA[0] == 0)\n",
    "    \n",
    "    # Supply-demand balance constraints\n",
    "    for i in N:\n",
    "        gen_indices = gen.loc[gen['bus'] == i, 'id'].values\n",
    "        \n",
    "        print( f\"i: {i}, gen_indices: {gen_indices}\")\n",
    "        print( f\"gen_indices-1: {gen_indices-1}\" ) \n",
    "        print( f\"GEN[gen_indices-1]: {GEN[gen_indices-1]}\" ) \n",
    "        print(\" \")\n",
    "        \n",
    "        pd_value = bus.loc[bus['bus_i'] == i, 'pd'].values[0]           # demand in MW\n",
    "        fbus_values = branch.loc[branch['fbus'] == i, 'tbus'].values\n",
    "        fbus_values_list = list(fbus_values) # Convert generator expression to a list\n",
    "        flow_vars = [ FLOW[i-1, j-1] for j in fbus_values_list ]        # List of CVXPY flow variables\n",
    "        \n",
    "        if i==2: # Virtual bidding at Bus-2\n",
    "            constraints.append(\n",
    "                cp.sum(GEN[gen_indices-1]) + supply_vb - (pd_value) == cp.sum(flow_vars)\n",
    "            )\n",
    "        else:\n",
    "            constraints.append(\n",
    "                cp.sum(GEN[gen_indices-1]) - pd_value == cp.sum(flow_vars)\n",
    "            )\n",
    "    \n",
    "    # Max generation constraints\n",
    "    for g in G: # 1, 2\n",
    "        constraints.append( GEN[g-1] <= gen.loc[g-1, 'pmax'] ) # also in MWs -- G1-g[0] AND G3-g[2] only \n",
    "    \n",
    "    # Flow constraints on each branch\n",
    "    for l in branch.index:\n",
    "        from_bus = branch.loc[l, 'fbus']\n",
    "        to_bus = branch.loc[l, 'tbus']\n",
    "        susceptance = branch.loc[l, 'sus']\n",
    "        constraints.append( FLOW[from_bus-1, to_bus-1] == baseMVA * susceptance * (THETA[from_bus-1] - THETA[to_bus-1]) )\n",
    "    \n",
    "    # Max line flow constraints\n",
    "    for l in branch.index:\n",
    "        from_bus = branch.loc[l, 'fbus']\n",
    "        to_bus = branch.loc[l, 'tbus']\n",
    "        constraints.append(FLOW[from_bus-1, to_bus-1] <= branch.loc[l, 'ratea'])\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.ECOS)  # is an open-source C library for solving convex second-order and exponential cone programs. \n",
    "    problem.solve(solver=cp.XPRESS ) \n",
    "    \n",
    "    # Prepare the output data\n",
    "    generation = pd.DataFrame({\n",
    "        'id': gen['id'],\n",
    "        'node': gen['bus'],\n",
    "        'gen': GEN.value\n",
    "    })\n",
    "    \n",
    "    angles = THETA.value\n",
    "    \n",
    "    flows = pd.DataFrame({\n",
    "        'fbus': branch['fbus'],\n",
    "        'tbus': branch['tbus'],\n",
    "        'flow': baseMVA * branch['sus'] * (angles[branch['fbus']-1] - angles[branch['tbus']-1])\n",
    "    })\n",
    "    \n",
    "    # Extract the prices (dual values of the balance constraints)\n",
    "    prices = pd.DataFrame({\n",
    "        'node': bus['bus_i'],\n",
    "        'value': [constraint.dual_value for constraint in constraints[1:len(N)+1]]\n",
    "    })\n",
    "    \n",
    "    # Return the solution and objective value\n",
    "    return {\n",
    "        'generation': generation.round(7),\n",
    "        'angles': [round(angle, 7) for angle in angles],\n",
    "        'flows': flows.round(5),\n",
    "        'prices': prices.applymap(lambda x: '{:.5f}'.format(x)),\n",
    "        'cost': '{:.5f}'.format(problem.value),\n",
    "        'status': problem.status\n",
    "    }"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "fa362483",
   "metadata": {},
   "outputs": [],
   "source": [
    "def input_data(datadir):\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",
    "    #gen['id'] = [1, 3]\n",
    "    gencost['id'] = [1, 2]\n",
    "\n",
    "    # Create line ids\n",
    "    branch['id'] = range(1, len(branch) + 1)\n",
    "\n",
    "    # Add set of rows for reverse direction with same parameters\n",
    "    branch2 = copy.deepcopy(branch)\n",
    "    branch2['fbus'], branch2['tbus'] = branch2['tbus'], branch2['fbus']\n",
    "    branch2 = branch2[branch.columns]  # Reorder columns to match original branch dataframe\n",
    "    branch = pd.concat([branch, branch2], ignore_index=True)\n",
    "\n",
    "    # Calculate the susceptance of each line\n",
    "    # Assuming reactance is much greater than resistance (x >> 0)\n",
    "    # Treat susceptance as the reciprocal of reactance (x)\n",
    "    branch['sus'] = 1 / branch['x']\n",
    "    \n",
    "    return bus, gen, gencost, branch"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "46f916aa",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "impedance_pu: 0.0001890359168241966, sus_pu: 5290.0\n"
     ]
    }
   ],
   "source": [
    "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}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5971a4c5",
   "metadata": {},
   "source": [
    "$ Z_{\\text{base}} = \\frac{{V_{\\text{base}}^2}}{{S_{\\text{base}}}}, X_{\\text{pu}} = \\frac{{X}}{{Z_{\\text{base}}}} $"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "8d5a63ae",
   "metadata": {},
   "outputs": [],
   "source": [
    "# import numpy as np\n",
    "\n",
    "# datadir = 'opf_data'\n",
    "# DAM_bus, DAM_gen, DAM_gencost, branch = input_data(datadir)\n",
    "\n",
    "# G = DAM_gen['id'].values\n",
    "# #G = np.array([1, 3])\n",
    "# N = DAM_bus['bus_i'].values\n",
    "# G"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "8de8ee76",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "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 DAM_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",
      "3     3     1  0.00281  0.000189  0.00712    500    500    500      0      0   \n",
      "4     2     1  0.00281  0.000189  0.00712     47    500    500      0      0   \n",
      "5     3     2  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",
      "3       1    -360     360   1  5290.0  \n",
      "4       1    -360     360   2  5290.0  \n",
      "5       1    -360     360   3  5290.0  \n"
     ]
    }
   ],
   "source": [
    "print(\"Day Ahead Market(DAM) Data:\\n\")\n",
    "\n",
    "datadir = 'opf_data'\n",
    "DAM_bus, DAM_gen, DAM_gencost, branch = input_data(datadir)\n",
    "\n",
    "print(\"Updated DAM Bus Dataframe:\")\n",
    "DAM_bus[\"pd\"] = [LOAD1, DA_LOAD2, 0]\n",
    "print(DAM_bus)\n",
    "print(\"\")\n",
    "\n",
    "print(\"Updated DAM_gen:\")\n",
    "DAM_gen.loc[1, \"bus\"] = 3 # generator to bus 3 \n",
    "print(DAM_gen)\n",
    "print(\"\")\n",
    "\n",
    "print(\"Updated DAM_gencost:\")\n",
    "DAM_gencost[\"x1\"] = [0.3, 0.8]\n",
    "DAM_gencost[\"y1\"] = [  3,   8]\n",
    "print(DAM_gencost)\n",
    "print(\"\")\n",
    "\n",
    "print(\"Updated DAM_branch:\")\n",
    "branch['x'] = xpu   # pu via per unit calculations\n",
    "branch['sus'] = 1 / branch['x']\n",
    "\n",
    "if congestion:\n",
    "    branch.loc[((branch['fbus'] == 1) & (branch['tbus'] == 2)) | ((branch['fbus'] == 2) & (branch['tbus'] == 1)), 'ratea'] = 47  # creating congestion\n",
    "\n",
    "print(branch)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "c2396a0d",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Day Ahead Market(DAM) Results:\n",
      "\n",
      "i: 1, gen_indices: [1]\n",
      "gen_indices-1: [0]\n",
      "GEN[gen_indices-1]: var1[0]\n",
      " \n",
      "i: 2, gen_indices: []\n",
      "gen_indices-1: []\n",
      "GEN[gen_indices-1]: var1[]\n",
      " \n",
      "i: 3, gen_indices: [2]\n",
      "gen_indices-1: [1]\n",
      "GEN[gen_indices-1]: var1[1]\n",
      " \n",
      "Using the Community license in this session. If you have a full Xpress license, pass the full path to your license file to xpress.init(). If you want to use the FICO Community license and no longer want to see this message, use the following code before using the xpress module:\n",
      "  xpress.init('/opt/homebrew/Caskroom/miniforge/base/envs/PGOC/lib/python3.9/site-packages/xpress/license/community-xpauth.xpr')\n",
      "generation:\n",
      "   id  node        gen\n",
      "0   1     1  59.090909\n",
      "1   2     3  15.909091\n",
      "\n",
      "angles:\n",
      "[0.0, -8.45e-05, -2.72e-05]\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",
      "3     3     1 -14.39394\n",
      "4     2     1 -44.69697\n",
      "5     3     2  30.30303\n",
      "\n",
      "prices:\n",
      "      node      value\n",
      "0  1.00000  -20.72727\n",
      "1  2.00000  -20.72727\n",
      "2  3.00000  -20.72727\n",
      "\n",
      "cost:\n",
      "929.54545\n",
      "\n",
      "status:\n",
      "optimal\n",
      "\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "<ipython-input-3-cd910b6aaf57>:102: FutureWarning: DataFrame.applymap has been deprecated. Use DataFrame.map instead.\n",
      "  'prices': prices.applymap(lambda x: '{:.5f}'.format(x)),\n"
     ]
    }
   ],
   "source": [
    "print(\"Day Ahead Market(DAM) Results:\\n\")\n",
    "\n",
    "DAM = dcopf(DAM_gen, branch, DAM_gencost, DAM_bus, SUPPLY_BID)\n",
    "\n",
    "for key, value in DAM.items():\n",
    "    print(f\"{key}:\") \n",
    "    print(value)\n",
    "    print()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ca6eb18b",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "58af140e",
   "metadata": {},
   "source": [
    "## Real-Time Market Optimization\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",
    "LMPs:\n",
    "$$\n",
    "\\mathbf{\\sigma}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "c3860784",
   "metadata": {},
   "outputs": [],
   "source": [
    "def RTM_dcopf(gen, branch, gencost, bus, **kwargs):\n",
    "    \n",
    "    # Define sets\n",
    "    G = gen['id'].values\n",
    "    N = bus['bus_i'].values\n",
    "\n",
    "    # Define base MVA\n",
    "    baseMVA = gen['mbase'].iloc[0]\n",
    "\n",
    "    # Decision variables\n",
    "    GEN = cp.Variable(   len(G) )    # Generation\n",
    "    THETA = cp.Variable( len(N) )             # Voltage phase angle of bus\n",
    "    FLOW = cp.Variable( (len(N), len(N)) )    # Flows between all pairs of nodes\n",
    "\n",
    "    # Objective function\n",
    "    objective = cp.Minimize( 0.5 * cp.sum(cp.multiply(gencost['x1'], cp.square(GEN))) + cp.sum(cp.multiply(gencost['y1'], GEN)) )\n",
    "    #objective = cp.Minimize( cp.sum( cp.multiply(gencost['x1'], GEN) ) )\n",
    "\n",
    "    # Define the constraints\n",
    "    constraints = []\n",
    "\n",
    "    # Slack bus reference angle constraint\n",
    "    constraints.append(THETA[0] == 0)\n",
    "\n",
    "    # Supply-demand balance constraints\n",
    "    for i in N:\n",
    "        gen_indices = gen.loc[gen['bus'] == i, 'id'].values\n",
    "        pd_value = bus.loc[bus['bus_i'] == i, 'pd'].values[0]\n",
    "        fbus_values = branch.loc[branch['fbus'] == i, 'tbus'].values\n",
    "        fbus_values_list = list(fbus_values)  # Convert generator expression to a list\n",
    "        flow_vars = [ FLOW[i-1, j-1] for j in fbus_values_list ]  # List of CVXPY flow variables\n",
    "        \n",
    "        if i==1:\n",
    "            constraints.append(\n",
    "                cp.sum(GEN[gen_indices-1]) + kwargs.get('arg1') - pd_value == cp.sum(flow_vars)\n",
    "        )\n",
    "        elif i==3:\n",
    "            constraints.append(\n",
    "                cp.sum(GEN[gen_indices-1]) + kwargs.get('arg2') - pd_value == cp.sum(flow_vars)\n",
    "            )\n",
    "        else:\n",
    "            constraints.append(\n",
    "                cp.sum(GEN[gen_indices-1]) - pd_value == cp.sum(flow_vars)\n",
    "            )     \n",
    "            \n",
    "    # Max generation constraints\n",
    "    for g in G:\n",
    "        constraints.append(GEN[g-1] <= gen.loc[g-1, 'pmax'])\n",
    "        \n",
    "    #constraints.append(GEN[0] >= kwargs.get('arg1') )\n",
    "    #constraints.append(GEN[2] >= kwargs.get('arg2') )\n",
    "\n",
    "    # Flow constraints on each branch\n",
    "    for l in branch.index:\n",
    "        from_bus = branch.loc[l, 'fbus']\n",
    "        to_bus = branch.loc[l, 'tbus']\n",
    "        susceptance = branch.loc[l, 'sus']\n",
    "        constraints.append( FLOW[from_bus-1, to_bus-1] == baseMVA * susceptance * (THETA[from_bus-1] - THETA[to_bus-1]) )\n",
    "    \n",
    "    # Max line flow constraints\n",
    "    for l in branch.index:\n",
    "        from_bus = branch.loc[l, 'fbus']\n",
    "        to_bus = branch.loc[l, 'tbus']\n",
    "        constraints.append( FLOW[from_bus-1, to_bus-1] <= branch.loc[l, 'ratea'] )\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.ECOS)  # is an open-source C library for solving convex second-order and exponential cone programs. \n",
    "    problem.solve(solver=cp.XPRESS ) \n",
    "     \n",
    "    \n",
    "    # Prepare the output data\n",
    "    generation = pd.DataFrame({\n",
    "        'id': gen['id'],\n",
    "        'node': gen['bus'],\n",
    "        'gen': GEN.value\n",
    "    })\n",
    "    \n",
    "    angles = THETA.value\n",
    "    \n",
    "    flows = pd.DataFrame({\n",
    "        'fbus': branch['fbus'],\n",
    "        'tbus': branch['tbus'],\n",
    "        'flow': baseMVA * branch['sus'] * (angles[branch['fbus']-1] - angles[branch['tbus']-1])\n",
    "    })\n",
    "    \n",
    "    # Extract the prices (dual values of the balance constraints)\n",
    "    prices = pd.DataFrame({\n",
    "        'node': bus['bus_i'],\n",
    "        'value': [constraint.dual_value for constraint in constraints[1:len(N)+1]]\n",
    "    })\n",
    "    \n",
    "    # Return the solution and objective value\n",
    "    return {\n",
    "        'generation': generation.round(7),\n",
    "        'angles': [round(angle, 7) for angle in angles],\n",
    "        'flows': flows.round(5),\n",
    "        'prices': prices.applymap(lambda x: '{:.5f}'.format(x)),\n",
    "        'cost': '{:.5f}'.format(problem.value),\n",
    "        'status': problem.status\n",
    "    }"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "03c5d832",
   "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",
      "3     3     1  0.00281  0.000189  0.00712    500    500    500      0      0   \n",
      "4     2     1  0.00281  0.000189  0.00712     47    500    500      0      0   \n",
      "5     3     2  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",
      "3       1    -360     360   1  5290.0  \n",
      "4       1    -360     360   2  5290.0  \n",
      "5       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]\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",
    "#RTM_branch['x'] = 0.00189\n",
    "#if congestion:\n",
    "#    RTM_branch.loc[((RTM_branch['fbus'] == 1) & (RTM_branch['tbus'] == 2)) | ((RTM_branch['fbus'] == 2) & (RTM_branch['tbus'] == 1)), 'ratea'] = 47  # creating congestion\n",
    "print(branch)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "b2fea23b",
   "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",
      "angles:\n",
      "[0.0, -8.88e-05, -2.32e-05]\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",
      "3     3     1 -12.25175\n",
      "4     2     1 -47.00000\n",
      "5     3     2  34.74825\n",
      "\n",
      "prices:\n",
      "      node      value\n",
      "0  1.00000  -10.28951\n",
      "1  2.00000  -19.02797\n",
      "2  3.00000  -14.65874\n",
      "\n",
      "cost:\n",
      "195.16154\n",
      "\n",
      "status:\n",
      "optimal\n",
      "\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "<ipython-input-9-8e0f91e04fa5>:99: FutureWarning: DataFrame.applymap has been deprecated. Use DataFrame.map instead.\n",
      "  'prices': prices.applymap(lambda x: '{:.5f}'.format(x)),\n"
     ]
    }
   ],
   "source": [
    "print(\"Real Time Market(RTM):\\n\")\n",
    "\n",
    "RTM = RTM_dcopf(RTM_gen, branch, RTM_gencost, RTM_bus, arg1=DAM['generation'][\"gen\"][0], arg2=DAM['generation'][\"gen\"][1])\n",
    "\n",
    "for key, value in RTM.items():\n",
    "    print(f\"{key}:\") \n",
    "    print(value)\n",
    "    print()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "7169705a",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Bus 1: 20.72727, 10.28951, 10.43776\n",
      "Bus 2: 20.72727, 19.02797, 1.69930\n",
      "Bus-2 Profit: 0.0\n",
      "Bus 3: 20.72727, 14.65874, 6.06853\n"
     ]
    }
   ],
   "source": [
    "RTM['prices']['value'] = pd.to_numeric(RTM['prices']['value'])\n",
    "RTM['prices']['value'] = RTM['prices']['value'].apply(lambda x: -x)\n",
    "\n",
    "DAM['prices']['value'] = pd.to_numeric(DAM['prices']['value'])\n",
    "DAM['prices']['value'] = DAM['prices']['value'].apply(lambda x: -x)\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",
    "        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": "c5ac6e78",
   "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
}