04_Power_prediction.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Example notebook to post-process and combine WRF results for wind farm power outputs\n",
    "\n",
    "# The example shown is for power output in a wind farm in the small Denmark domain"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# import necessary packages\n",
    "\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import glob\n",
    "import pylab as plt\n",
    "import random"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# read in norm of initial wind velocity at future location of hypothetical wind farm\n",
    "df_vel_norm = pd.read_csv('vel_avg_norm.csv')\n",
    "\n",
    "# read in arbitrary file to get location of wind farm array\n",
    "farm_array = np.load(\"/rds/general/project/wrfwindpower/live/wrf_results_square/fixed_farm/200701_01/power.npy\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "# read in theoretical power curve of wind turbines used\n",
    "\n",
    "from scipy.interpolate import interp1d\n",
    "f = open('power_curve.txt').read()\n",
    "power_file = [line.split(' ') for line in f.split('\\n')]\n",
    "vel_txt = [0]\n",
    "for i in power_file[:-1]:\n",
    "    vel_txt.append(float(i[0]))\n",
    "for _ in range(10):\n",
    "    vel_txt.append(vel_txt[-1]+1)\n",
    "power_txt = [0]\n",
    "for i in power_file[:-1]:\n",
    "    power_txt.append(float(i[2]))\n",
    "for _ in range(10):\n",
    "    power_txt.append(power_txt[-1])\n",
    "power_curve = interp1d(vel_txt, [i*1000 for i in power_txt])\n",
    "\n",
    "power_curve_power = power_curve(df_vel_norm['norm_vel']).sum()*136"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# read in outputs from WRF simulations for whole of 2007 - this is what we use as the truth\n",
    "\n",
    "fname_fixed_list = []\n",
    "power_fixed_list = []\n",
    "\n",
    "path = \"/rds/general/project/wrfwindpower/live/wrf_results_square/fixed_farm/*/\"\n",
    "\n",
    "path_txt = str(path) + '/power.npy'\n",
    "for fname in glob.glob(path_txt):\n",
    "    fname_fixed_list.append(fname)\n",
    "\n",
    "fname_fixed_list.sort()\n",
    "\n",
    "for i in fname_fixed_list:\n",
    "    power_fixed_list.append(np.load(i))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/apps/jupyterhub/2019-04-29/miniconda/lib/python3.6/site-packages/ipykernel_launcher.py:20: RuntimeWarning: invalid value encountered in arccos\n",
      "/apps/jupyterhub/2019-04-29/miniconda/lib/python3.6/site-packages/ipykernel_launcher.py:22: RuntimeWarning: invalid value encountered in arccos\n",
      "/apps/jupyterhub/2019-04-29/miniconda/lib/python3.6/site-packages/ipykernel_launcher.py:24: RuntimeWarning: invalid value encountered in arccos\n"
     ]
    }
   ],
   "source": [
    "# load in WRF outputs for the closest datapoints to the cluster centres for cluster found by clustering\n",
    "# on wind velocity on Europe-wide domain\n",
    "\n",
    "def load_file(folder_init, file_init):\n",
    "    folder = '/rds/general/project/wrfwindpower/live/' + folder_init\n",
    "    cluster_0 = np.load(folder + '_0/' + file_init + '.npy')\n",
    "    cluster_1 = np.load(folder + '_1/' + file_init + '.npy')\n",
    "    cluster_2 = np.load(folder + '_2/' + file_init + '.npy')\n",
    "    cluster_3 = np.load(folder + '_3/' + file_init + '.npy')\n",
    "    cluster_4 = np.load(folder + '_4/' + file_init + '.npy')\n",
    "    cluster_5 = np.load(folder + '_5/' + file_init + '.npy')\n",
    "    \n",
    "    return cluster_0, cluster_1, cluster_2, cluster_3, cluster_4, cluster_5\n",
    "\n",
    "power_cluster_0_wind_d01, power_cluster_1_wind_d01, power_cluster_2_wind_d01,\\\n",
    "    power_cluster_3_wind_d01, power_cluster_4_wind_d01, power_cluster_5_wind_d01 = load_file('wrf_results_wind_d01/fixed_sq_farm', 'power')\n",
    "\n",
    "u1_cluster_0_wind_d01, u1_cluster_1_wind_d01, u1_cluster_2_wind_d01,\\\n",
    "    u1_cluster_3_wind_d01, u1_cluster_4_wind_d01, u1_cluster_5_wind_d01 = load_file('wrf_results_wind_d01/fixed_sq_farm', 'u')\n",
    "\n",
    "u2_cluster_0_wind_d01, u2_cluster_1_wind_d01, u2_cluster_2_wind_d01,\\\n",
    "    u2_cluster_3_wind_d01, u2_cluster_4_wind_d01, u2_cluster_5_wind_d01 = load_file('wrf_results_wind_d01/fixed_sq_farm', 'v')\n",
    "\n",
    "# estimate the norm of the velocity at the location of the wind farm\n",
    "norm_vel_0_wind_d01 = np.sqrt((u1_cluster_0_wind_d01[0, 11, :, :72])**2 + (u2_cluster_0_wind_d01[0, 11, :153, :])**2)[farm_array[-1]!=0].mean()\n",
    "norm_vel_1_wind_d01 = np.sqrt((u1_cluster_1_wind_d01[0, 11, :, :72])**2 + (u2_cluster_1_wind_d01[0, 11, :153, :])**2)[farm_array[-1]!=0].mean()\n",
    "norm_vel_2_wind_d01 = np.sqrt((u1_cluster_2_wind_d01[0, 11, :, :72])**2 + (u2_cluster_2_wind_d01[0, 11, :153, :])**2)[farm_array[-1]!=0].mean()\n",
    "norm_vel_3_wind_d01 = np.sqrt((u1_cluster_3_wind_d01[0, 11, :, :72])**2 + (u2_cluster_3_wind_d01[0, 11, :153, :])**2)[farm_array[-1]!=0].mean()\n",
    "norm_vel_4_wind_d01 = np.sqrt((u1_cluster_4_wind_d01[0, 11, :, :72])**2 + (u2_cluster_4_wind_d01[0, 11, :153, :])**2)[farm_array[-1]!=0].mean()\n",
    "norm_vel_5_wind_d01 = np.sqrt((u1_cluster_5_wind_d01[0, 11, :, :72])**2 + (u2_cluster_5_wind_d01[0, 11, :153, :])**2)[farm_array[-1]!=0].mean()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# split real power values for each day in 2007 by regime\n",
    "df_wind_d01 = pd.read_csv('labels_all_6_wind_d01.csv')\n",
    "df_2007_wind_d01 = df_wind_d01[df_wind_d01['Year'] == 2007].reset_index(drop = True)\n",
    "\n",
    "final_fixed_power = [i[-1].sum() for i in power_fixed_list]\n",
    "\n",
    "combined_wind_d01 = pd.concat([df_2007_wind_d01, df_vel_norm], join = 'inner', axis = 1)\n",
    "\n",
    "reg_fixed_power_wind_d01 = pd.concat([combined_wind_d01[['Regime', 'norm_vel']], pd.DataFrame(final_fixed_power)], axis = 1)\n",
    "\n",
    "reg_fixed_power_0_wind_d01 = reg_fixed_power_wind_d01[reg_fixed_power_wind_d01['Regime'] == 0].reset_index(drop = True)\n",
    "reg_fixed_power_1_wind_d01 = reg_fixed_power_wind_d01[reg_fixed_power_wind_d01['Regime'] == 1].reset_index(drop = True)\n",
    "reg_fixed_power_2_wind_d01 = reg_fixed_power_wind_d01[reg_fixed_power_wind_d01['Regime'] == 2].reset_index(drop = True)\n",
    "reg_fixed_power_3_wind_d01 = reg_fixed_power_wind_d01[reg_fixed_power_wind_d01['Regime'] == 3].reset_index(drop = True)\n",
    "reg_fixed_power_4_wind_d01 = reg_fixed_power_wind_d01[reg_fixed_power_wind_d01['Regime'] == 4].reset_index(drop = True)\n",
    "reg_fixed_power_5_wind_d01 = reg_fixed_power_wind_d01[reg_fixed_power_wind_d01['Regime'] == 5].reset_index(drop = True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Compare the power predictions for each different method assessed in the paper\n",
    "\n",
    "# WRF is used as the real value, and then we have three separate methods of combining a reduced number \n",
    "# of WRF runs\n",
    "\n",
    "# Here we show the power predictions in the Denmark domain for the clusters obtained by clustering wind velocity \n",
    "# on the Europe-wide domain but the code would be the same for the other domains\n",
    "\n",
    "# The plot in the paper is slightly different as it combines all the results from the different clusters for each domain"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "92995407243.39062\n",
      "83077128877.33334\n",
      "91764363366.70206\n",
      "102165465519.0\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAEKCAYAAAARnO4WAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3Xl8VeW1//HPIkADGgMoerVMwaoISThAwhQQFAQcsOBQsF4RrKJSrVMdrtP1Z6+tWurAlaJ4VaqliuJF/anVggVBAjJImBGKRgQnEBlkksD6/XFOzo9AhgNk55Ds7/v1you999lnP+s5JOs85zl7r23ujoiI1Hy1kh2AiIhUDSV8EZGQUMIXEQkJJXwRkZBQwhcRCQklfBGRkAg04ZtZAzObaGYrzGy5mXUJsj0RESlb7YCP/wTwrrtfbGZ1gfoBtyciImWwoC68MrNjgIVAS9fVXSIiSRfkCL8lsB543szaAvOBG9192747mdlwYDjAUUcd1aFVq1YBhiSSPIVbCgFocUyL6IYNq6L/HndKlcXw42efAVA3I6PK2pRgzZ8/f4O7N05k3yBH+DnAbCDP3T8ysyeALe5+b1nPycnJ8Xnz5gUSj0iyDXt3GADP93s+uuH582IPvF1lMXx++RAAmr/4QpW1KcEys/nunpPIvkF+absWWOvuH8XWJwLtA2xPRETKEVjCd/evgS/M7LTYpl7AsqDaExGR8gV9ls4NwPjYGTqfAsMCbk9ERMoQaMJ39wIgobmlsuzevZu1a9eyc+fOSopK5NCkpqbSpEkT6tSpk+xQRA5J0CP8w7Z27VrS0tJo0aIFZpbscCSk3J3vvvuOtWvXkqEzXKSaOuJLK+zcuZNjjz1WyV6Sysw49thj9UlTqrUjPuEDSvZyRNDvoVR31SLhi4jI4VPCr8DNN9/M448/Hl/v27cvV111VXz91ltv5dFHH6VevXpEIhFat27NkCFD2L17NwDTpk0jPT2dSCRCJBKhd+/eVd4HERFQwq9Q165dyc/PB2Dv3r1s2LCBpUuXxh/Pz88nLy+Pk08+mYKCAhYvXszatWt55ZVX4vt0796dgoICCgoKmDJlSpX3QUQElPArlJeXF0/4S5cuJTMzk7S0NL7//nt27drF8uXLadiwYXz/lJQUOnbsyLp165IVsohIqY740zJL+Pud8PXiyj3mv2XBOQ+V+fBJJ51E7dq1WbNmDfn5+XTp0oV169Yxa9Ys0tPTyc7Opm7duvH9d+7cyUcffcQTTzwR3zZjxgwikQgAl1xyCXfffXfl9kFEJAHVK+EnSfEoPz8/n1tuuYV169aRn59Peno6Xbt2BWD16tVEIhFWrVrFxRdfTHZ2dvz53bt356233kpW+CIiQHVL+OWMxINUPI+/ePFiMjMzadq0KX/605845phjuPLKKwHic/hfffUVPXv25M033+SCCy5ISrwiIqXRHH4C8vLyeOutt2jUqBEpKSk0atSITZs2MWvWLLp0KXnXxhNPPJGHHnqIP/zhD0mKVkSkdEr4CcjKymLDhg107ty5xLb09HSOO+64A/YfMGAA27dvZ8aMGVUZpohIuarXlE6SpKSksGXLlhLbxo0bF19u0aIFS5Ysia+bGQsXLoyv9+zZM+gQRUQqpBG+iEhIKOGLiISEEr6ISEgo4YuIhIQSvohISCjhi4iEhBJ+Ar7++msGDx7MySefTOvWrTn33HNZuXIlhYWFZGZmHtIxx40bx5dfflnJkSbmyy+/5OKLL65wv6OPPrrU7a+//jrLli076HbLOp6IVA0l/Aq4OwMHDqRnz56sXr2aZcuW8fvf/55vvvnmsI57KAm/qKjosNosdtJJJzFx4sRDfv6hJvyqtGfPnmSHIHLEUcKvwNSpU6lTpw7XXnttfFskEqF79+4l9hs3bhzXX399fP38889n2rRp7Nmzh6FDh5KZmUlWVhaPPfYYEydOZN68eVx22WVEIhF27NjB/Pnz6dGjBx06dKBv37589dVXQPSirbvuuosePXqUqMAJ0at9N23ahLtz7LHH8sILLwBw+eWXM2XKFPbs2cNtt91Gbm4u2dnZPP300wAlPpls376dX/ziF2RnZzNo0CA6derEvHnz4m3cfffdtG3bls6dO/PNN9+Qn5/Pm2++yW233UYkEmH16tWsXr2afv360aFDB7p3786KFSsA+Oyzz+jSpQu5ubnce++9Zb7GAwYMoEOHDrRp04axY8cCMGbMGG6//fYSr+8NN9wAwF//+lc6duxIJBLhmmuuiSf3o48+mvvuu49OnToxa9YsHnjgAXJzc8nMzGT48OG4OwBz584lOzubLl26cNttt8Vfi7JeL5GaolpdafvwnIdZsXFFpR6zVaNW3NHxjjIfX7JkCR06dDjk4xcUFLBu3br4lbibNm2iQYMGPPnkk4wcOZKcnBx2797NDTfcwBtvvEHjxo2ZMGECd999N88991z8OR988MEBx87Ly2PmzJk0b96cli1bMmPGDIYMGcLs2bMZM2YMzz77LOnp6cydO5ddu3aRl5dHnz59Styb9c9//jMNGzZk0aJFLFmyJF7GGWDbtm107tyZBx98kNtvv51nnnmGe+65hwsuuIDzzz8/Pi3Uq1cvnnrqKU455RQ++ugjRowYwT//+U9uvPFGrrvuOoYMGcLo0aPLfI2ee+45GjVqxI4dO8jNzeWiiy7i4osvpkuXLjzyyCMA8ddk+fLlTJgwgZkzZ1KnTh1GjBjB+PHjGTJkCNu2bSMzM5MHHngAgNatW3PfffcB0TfBt956i/79+zNs2DDGjh1L165dufPOO+NxlPV6ZWRkHNL/vciRplol/OqoZcuWfPrpp9xwww2cd9559OnT54B9PvnkE5YsWcLZZ58NREeaJ554YvzxQYMGlXrs7t27M336dJo3b851113H2LFjWbduHY0aNeLoo4/mH//4B4sWLYpP32zevJlVq1Zx6qmnxo/x4YcfcuONNwKQmZlZoqxz3bp1Of/88wHo0KEDkydPPiCGH374gfz8fC655JL4tl27dgEwc+ZMXnvtNSCacO+4o/Q31lGjRjFp0iQAvvjiC1atWkXnzp1p2bIls2fP5pRTTuGTTz4hLy+P0aNHM3/+fHJzcwHYsWMHxx9/PBAtgXHRRRfFjzt16lQeeeQRtm/fzsaNG2nTpg3du3dn69at8bLWv/zlL+Olq8t6vZTwpaaoVgm/vJF4UNq0aZPQfHft2rXZu3dvfH3nzp0ANGzYkIULF/Lee+8xevRoXnnllfjIvZi706ZNG2bNmlXqsY866qhSt59xxhmMHj2aNWvW8OCDDzJp0iQmTpwYn25yd/77v/+bvn37lnheYWFhibbLUqdOnfingZSUlFK/Q9i7dy8NGjSgoKCg1GPs+2miNNOmTWPKlCnMmjWL+vXr07Nnz/hrN2jQIF555RVatWrFwIEDMTPcnSuuuKLUaqSpqamkpKQA0dd/xIgRzJs3j6ZNm3L//fezc+fOcvtb1uslUlNoDr8CZ511Frt27eKZZ56Jb5s7d+4BUywtWrSgoKCAvXv38sUXXzBnzhwANmzYwN69e7nooov43e9+x8cffwxAWloaW7duBeC0005j/fr18YS/e/fuEvfNLUvTpk3ZsGEDq1atomXLlnTr1o2RI0fGE37fvn0ZM2ZM/IbqK1euZNu2bSWO0a1bt/j9d5ctW8bixRXfUWzf2I855hgyMjJ49dVXgWjSLC4cl5eXx8svvwzA+PHjSz3W5s2badiwIfXr12fFihXMnj07/tiFF17I66+/zksvvRT/lNOrVy8mTpzIt99+C8DGjRv5/PPPDzhu8ZvGcccdxw8//BB/027YsCFpaWnxdorjS/T1EqnOlPArYGZMmjSJyZMnc/LJJ9OmTRvuv/9+TjrppBL75eXlkZGRQVZWFr/97W9p3749AOvWraNnz55EIhGGDh0aH5kOHTqUa6+9lkgkwp49e5g4cSJ33HEHbdu2JRKJxO+jW5FOnTrFp2i6d+/OunXr6NatGwBXXXUVrVu3pn379mRmZnLNNdccMEofMWIE69evJzs7m4cffpjs7GzS09PLbXPw4MH88Y9/pF27dqxevZrx48fz7LPP0rZtW9q0acMbb7wBwBNPPMHo0aPJzc1l8+bNpR6rX79+FBUVkZ2dzb333luiBHXDhg1p3bo1n3/+OR07dgSi8/L/9V//RZ8+fcjOzubss8+Of8G9rwYNGnD11VeTlZXFgAED4lNAEJ2rHz58OF26dMHd4/1N5PUSqc6svI+4h31ws0JgK7AHKHL3nPL2z8nJ8X3PEAFYvnw5p59+emAxht2ePXvYvXs3qamprF69ml69erFy5coS9+mtaX744Yf4NQEPPfQQX3311QFnQJXlcH4fh707DIDn+z0f3fD8ebEH3j6k4x2Kzy8fAkDzF1+osjYlWGY2v6LcWqwq5vDPdPcNVdCOHILt27dz5plnsnv3btydMWPG1OhkD/D222/zhz/8gaKiIpo3b17i3gYiNVm1+tJWKl9aWhr7f6qq6QYNGlTmmU8iNVnQc/gO/MPM5pvZ8IDbEhGRcgQ9ws9z9y/N7HhgspmtcPfp++4QeyMYDtCsWbOAwxERCa9AR/ju/mXs32+BSUDHUvYZ6+457p7TuHHjIMMREQm1wBK+mR1lZmnFy0AfYEn5zxIRkaAEOcI/AfjQzBYCc4C33f3dANsLTHUqj5xo6eNE7F8QTkSqt8Dm8N39U6BtUMevKsXlka+44or4VZkFBQV88803NG3a9JCPO27cODIzMw+4gKs8RUVF1K5d/n/Z4ZY+PtLs2bMnXi5BRA6PrrStwJFcHvmDDz4gEokQiURo164dW7duLfGpY9y4cQwYMID+/fuTkZHBk08+yaOPPkq7du3o3LkzGzdujLdx00030bVrVzIzM+NlIfa1fv16LrroInJzc8nNzWXmzJkH7FNYWEj37t1p37497du3j18tPGjQIN555534fkOHDuW1114rsxzxtGnTOPPMM/nlL39JVlYWUHoJZYheNXvqqafSs2dPrr766vj/QSLxioRNtToP/+vf/55dyyu3PPJPTm/Fv911V5mPH8nlkUeOHMno0aPJy8vjhx9+IDU1tdT4FyxYwM6dO/nZz37Gww8/zIIFC7j55pt54YUXuOmmm4BoKeT8/HymT5/OlVdeGY+32I033sjNN99Mt27dWLNmDX379mX58uUl9jn++OOZPHkyqamprFq1iksvvZR58+YxePBgJkyYwLnnnsuPP/7I+++/X275ZoA5c+awZMmSeKXK0koo79q1K16fKC0tjbPOOou2bdsmHK9I2FSrhF8dBVkeOS8vj1tuuYXLLruMCy+8kCZNmhywz5lnnklaWhppaWmkp6fTv39/IHrzlEWLFsX3u/TSS4FoBc4tW7awadOmEseZMmVKibtcbdmyha1bt5KWlhbftnv3bq6//noKCgpISUlh5cqVAJxzzjn85je/YdeuXbz77rucccYZ1KtXr8xyxHXr1qVjx44lyhKXVkL566+/pkePHjRq1AiASy65JN5mIvGKhE21SvjljcSDciSXR77zzjs577zzeOedd+jcuTNTpkw5YJT/k5/8JL5cq1at+HqtWrVKFAbbv4zx/ut79+5l1qxZ1KtXr9RYAB577DFOOOEEFi5cyN69e+OxpKam0rNnT9577z0mTJgQf3MpqxzxtGnTSvS5rBLK5dWBSiRekbDRHH4FjuTyyKtXryYrK4s77riDnJyc+K0FD8WECROA6A1R0tPTD6iY2adPH5588sn4emn17zdv3syJJ55IrVq1ePHFF0vcV3bw4ME8//zzzJgxI57gEy1HXFYJ5Y4dO/LBBx/w/fffU1RUFL/ZSqLxioSNEn4FjuTyyI8//jiZmZm0bduWevXqcc455xxyPxs2bEjXrl259tprefbZZw94fNSoUcybN4/s7Gxat27NU089dcA+I0aM4C9/+QudO3dm5cqVJUbpffr0Yfr06fTu3TtenC3RcsRllVD+6U9/yl133UWnTp3o3bs3rVu3jr9RJRKvSNgEWh75YKk8cnL07Nkz/gVydVNc6rioqIiBAwdy5ZVXMnDgwMDaU3lkOdIcTHlkjfClWrv//vuJRCJkZmaSkZHBgAEDkh2SyBGrWn1pK8GYNm1askM4ZCNHjkx2CCLVhkb4IiIhoYQvIhISSvgiIiGhhC8iEhLV7kvbFndW7ilshQ+dV+E+KSkpZGVlUVRUREZGBi+++CINGjQ4/LYLCzn//PMPqFsjIhIEjfATUK9ePQoKCliyZAmNGjVi9OjRyQ5JROSgKeEfpC5durBu3TogetFPr169aN++PVlZWbzxxhtAdOR++umnc/XVV9OmTRv69OnDjh07AJg/fz5t27alS5cuJd44du7cybBhw8jKyqJdu3ZMnToVSLzEsYhIRZTwD8KePXt4//33ueCCC4BoUbBJkybx8ccfM3XqVG699dZ4Qa9Vq1bx61//mqVLl9KgQYN4nZdhw4YxatSoAwqlFSf/xYsX89JLL3HFFVfEC7AtWbKEv/3tb8yZM4e7776b+vXrs2DBArp06cILL+iKSRFJjBJ+Anbs2EEkEuHYY49l48aN8TLG7s5dd91FdnY2vXv3Zt26dXzzzTcAZGRkEIlEAOjQoQOFhYVs3ryZTZs20aNHDwAuv/zyeBsffvhhfL1Vq1Y0b948Xuq3uMRx48aNDyhxXFhYWCWvgYhUf0r4CSiew//888/58ccf46Px8ePHs379eubPn09BQQEnnHBCfFS+b1nilJQUioqKcPcDyg4XK6+mUaIljkVEyqOEfxDS09MZNWoUI0eOZPfu3WzevJnjjz+eOnXqMHXqVD7//PNyn9+gQQPS09P58MMPgegbRrEzzjgjvr5y5UrWrFnDaaedFlxnRCR0qt1pmYmcRhmkdu3a0bZtW15++WUuu+wy+vfvT05ODpFIhFatWlX4/Oeff54rr7yS+vXrl7jxx4gRI7j22mvJysqidu3ajBs3rsTIXkTkcKk8sshBUHlkOdKoPLKIiBxACV9EJCSqRcI/kqadJLz0eyjV3RGf8FNTU/nuu+/0xyZJ5e589913pKamJjsUkUN2xJ+l06RJE9auXcv69euTHYqEXGpqKk2aNEl2GCKH7IhP+HXq1CEjIyPZYYiIVHuBT+mYWYqZLTCzt4JuS0REylYVc/g3AsuroB0RESlHoAnfzJoA5wH/E2Q7IiJSsaDn8B8HbgfSytrBzIYDwwGaNWsWcDhSGSr7rmMHI1mlNSqjz/WafVfiWC/Xja4PLufYyS4lIjVLYCN8Mzsf+Nbd55e3n7uPdfccd89p3LhxUOGIiIRekFM6ecAFZlYIvAycZWZ/DbA9EREpR2AJ393/w92buHsLYDDwT3f/96DaExGR8h3xV9qKiEjlqJILr9x9GjCtKtoSEZHSaYQvIhISSvgiIiGhhC8iEhJK+CIiIaGELyISEkr4IiIhkVDCN7PMoAMREZFgJTrCf8rM5pjZCDNrEGhEIiISiIQSvrt3Ay4DmgLzzOxvZnZ2oJGJiEilSngO391XAfcAdwA9gFFmtsLMLgwqOBERqTyJzuFnm9ljRO9cdRbQ391Pjy0/FmB8IiJSSRKtpfMk8Axwl7vvKN7o7l+a2T2BRCYiIpUq0YR/LrDD3fcAmFktINXdt7v7i4FFJyIilSbROfwpQL191uvHtomISDWRaMJPdfcfildiy/WDCUlERIKQaMLfZmbti1fMrAOwo5z9RUTkCJPoHP5NwKtm9mVs/URgUDAhiYhIEBJK+O4+18xaAacBBqxw992BRiYiIpXqYG5xmAu0iD2nnZnh7i8EEpWIiFS6hBK+mb0InAwUAHtimx1QwhcRqSYSHeHnAK3d3YMMRkREgpPoWTpLgH8LMhAREQlWoiP844BlZjYH2FW80d0vCCQqERGpdIkm/PuDDEJERIKX6GmZH5hZc+AUd59iZvWBlGBDExGRypRoeeSrgYnA07FNPwVeDyooERGpfIl+aftrIA/YAvGboRwfVFAiIlL5Ek34u9z9x+IVM6tN9Dz8MplZauw+uAvNbKmZ/Z/DCVRERA5Pogn/AzO7C6gXu5ftq8D/reA5u4Cz3L0tEAH6mVnnQw9VREQOR6IJ/05gPbAYuAZ4h+j9bcvkUcUllevEfnThlohIkiR6ls5eorc4fOZgDm5mKcB84GfAaHf/qJR9hgPDAZo1a3YwhxcRkYOQ6Fk6n5nZp/v/VPQ8d9/j7hGgCdDRzDJL2Wesu+e4e07jxo0PvgciIpKQg6mlUywVuARolGgj7r7JzKYB/YiWaRARkSqW0Ajf3b/b52eduz8OnFXec8yssZk1iC3XA3oDKw47YhEROSSJlkduv89qLaIj/rQKnnYi8JfYPH4t4BV3f+uQohQRkcOW6JTOn/ZZLgIKgV+U9wR3XwS0O7SwRESksiV6ls6ZQQciIiLBSnRK55byHnf3RysnHBERCcrBnKWTC7wZW+8PTAe+CCIoERGpfAdzA5T27r4VwMzuB15196uCCkxERCpXoqUVmgE/7rP+I9Ci0qMREZHAJDrCfxGYY2aTiNbDGQi8EFhUIiJS6RI9S+dBM/s70D22aZi7LwguLBERqWyJTukA1Ae2uPsTwFozywgoJhERCUCixdP+E7gD+I/YpjrAX4MKSkREKl+iI/yBwAXANgB3/5KKSyuIiMgRJNGE/6O7O7EbmJjZUcGFJCIiQUg04b9iZk8DDczsamAKB3kzFBERSa5Ez9IZGbuX7RbgNOA+d58caGQiIlKpKkz4sfLG77l7b0BJXkSkmqpwSsfd9wDbzSy9CuIREZGAJHql7U5gsZlNJnamDoC7/yaQqEREpNIlmvDfjv2IiEg1VW7CN7Nm7r7G3f9SVQGJiEgwKprDf714wcxeCzgWEREJUEUJ3/ZZbhlkICIiEqyKEr6XsSwiItVMRV/atjWzLURH+vViy8TW3d2PCTQ6ERGpNOUmfHdPqapAREQkWAdTD19ERKoxJXwRkZBQwhcRCQklfBGRkAgs4ZtZUzObambLzWypmd0YVFsiIlKxRGvpHIoi4FZ3/9jM0oD5ZjbZ3ZcF2KaIiJQhsBG+u3/l7h/HlrcCy4GfBtWeiIiUL8gRfpyZtQDaAR+V8thwYDhAs2bNDrmNFncmp5hn4UPnJaVdSF6fkymMfQ6bZP4fJ/PvuSoE/qWtmR0NvAbc5O5b9n/c3ce6e4675zRu3DjocEREQivQhG9mdYgm+/Hu/r9BtiUiIuUL8iwdA54Flrv7o0G1IyIiiQlyhJ8HXA6cZWYFsZ9zA2xPRETKEdiXtu7+ISXr6YuISBLpSlsRkZBQwhcRCQklfBGRkFDCFxEJCSV8EZGQUMIXEQkJJXwRkZBQwhcRCQklfBGRkFDCFxEJCSV8EZGQUMIXEQkJJXwRkZBQwhcRCQklfBGRkFDCFxEJCSV8EZGQUMIXEQkJJXwRkZBQwhcRCQklfBGRkFDCFxEJCSV8EZGQUMIXEQkJJXwRkZBQwhcRCQklfBGRkAgs4ZvZc2b2rZktCaoNERFJXJAj/HFAvwCPLyIiByGwhO/u04GNQR1fREQOjubwRURConayAzCz4cBwgGbNmiU5moPX4s63kx2C1GCV/fv18KffAdBDv7elStbfc+FD51VJO0kf4bv7WHfPcfecxo0bJzscEZEaK+kJX0REqkaQp2W+BMwCTjOztWb2q6DaEhGRigU2h+/ulwZ1bBEROXia0hERCQklfBGRkFDCFxEJCSV8EZGQUMIXEQkJJXwRkZBQwhcRCQklfBGRkFDCFxEJCSV8EZGQUMIXEQkJJXwRkZBQwhcRCQklfBGRkFDCFxEJCSV8EZGQUMIXEQkJJXwRkZBQwhcRCQklfBGRkFDCFxEJCSV8EZGQUMIXEQkJJXwRkZBQwhcRCQklfBGRkFDCFxEJCSV8EZGQCDThm1k/M/vEzP5lZncG2ZaIiJQvsIRvZinAaOAcoDVwqZm1Dqo9EREpX5Aj/I7Av9z9U3f/EXgZ+HmA7YmISDnM3YM5sNnFQD93vyq2fjnQyd2v32+/4cDw2OppwCeVGMZxwIZKPF51E+b+h7nvoP6Hqf/N3b1xIjvWDjAIK2XbAe8u7j4WGBtIAGbz3D0niGNXB2Huf5j7Dup/2PtfliCndNYCTfdZbwJ8GWB7IiJSjiAT/lzgFDPLMLO6wGDgzQDbExGRcgQ2pePuRWZ2PfAekAI85+5Lg2qvDIFMFVUjYe5/mPsO6n/Y+1+qwL60FRGRI4uutBURCQklfBGRkKgRCT/REg5mdrGZuZnVmNO1Eum7mf3CzJaZ2VIz+1tVxxikivpvZs3MbKqZLTCzRWZ2bjLiDIKZPWdm35rZkjIeNzMbFXttFplZ+6qOMUgJ9P+yWL8XmVm+mbWt6hiPOO5erX+IfiG8GmgJ1AUWAq1L2S8NmA7MBnKSHXdV9R04BVgANIytH5/suKu4/2OB62LLrYHCZMddif0/A2gPLCnj8XOBvxO9JqYz8FGyY67i/nfd5/f+nJrW/0P5qQkj/ERLOPwOeATYWZXBBSyRvl8NjHb37wHc/dsqjjFIifTfgWNiy+nUoGtB3H06sLGcXX4OvOBRs4EGZnZi1UQXvIr67+75xb/3RAd6TaoksCNYTUj4PwW+2Gd9bWxbnJm1A5q6+1tVGVgVqLDvwKnAqWY208xmm1m/KosueIn0/37g381sLfAOcEPVhHZESOT1CYtfEf20E2pBllaoKuWWcDCzWsBjwNCqCqgKJVK+ojbRaZ2eREc4M8ws0903BRxbVUik/5cC49z9T2bWBXgx1v+9wYeXdAmVN6npzOxMogm/W7JjSbaaMMKvqIRDGpAJTDOzQqJzmW/WkC9uEylfsRZ4w913u/tnRIvTnVJF8QUtkf7/CngFwN1nAalEC2uFQejLm5hZNvA/wM/d/btkx5NsNSHhl1vCwd03u/tx7t7C3VsQncu7wN3nJSfcSpVI+YrXgTMBzOw4olM8n1ZplMFJpP9rgF4AZnY60YS/vkqjTJ43gSGxs3U6A5vd/atkB1VVzKyrecZoAAAEFklEQVQZ8L/A5e6+MtnxHAmq/ZSOl1HCwcweAOa5e42t35Ng398D+pjZMmAPcFtNGekk2P9bgWfM7Gai0xlDPXbaRnVnZi8Rnao7LvYdxX8CdQDc/Smi31mcC/wL2A4MS06kwUig//cBxwJ/NjOAIg95BU2VVhARCYmaMKUjIiIJUMIXEQkJJXwRkZBQwhcRCQklfBGRAFVU5G2/fc8ws4/NrMjMLt7vsXfNbJOZHXLFACV8qfbMbI+ZFZjZEjN71czqJymOm8xsiJm1NbOCfbZfambbzaxObD3LzBbFll82s5pyIZyUbhyQaEmTNUSrApRW1faPwOWHE4gSvtQEO9w94u6ZwI/AtUE3aGYp+63XBq4k+oe6GGhuZmmxh7sCK4B2+6zPjC2PAW4POl5JntKKvJnZybER+3wzm2FmrWL7Frr7IuCA0h/u/j6w9XBiUcKXmmYG8DMAM7slNupfYmY3xbbdbma/iS0/Zmb/jC33MrO/xpb7mNms2EfrV83s6Nj2QjO7z8w+BC7Zr92zgI/dvShWp2cu0Cn2WAdgNNFET+zf/H3i7R17w5DwGAvc4O4dgN8Cf66KRpXwpcaIJc1zgMVm1oHolaWdiNZPujpWNXU60D32lBzg6NhUSzeiheWOA+4Bert7e2AecMs+zex0927u/vJ+zecB8/dZzwe6mtlRREdr0yiZ8GcCxN4c/gXo5hwhERtAdAVejU39PQ1USdlqjSqkJqi3z5z5DOBZ4DpgkrtvAzCz/yWa6McAHWLTLbuAj4km/u7Ab4i+ObQGZsYux68LzNqnrQllxHAisHyf9ZlEyzrMAOa6+2oz+5mZNQaOdvd96xl9C5xEyTcMqblqAZvcPVLVDSvhS02wY/8/Hotl6/25++5Y1dRhREfhi4gWlzuZaMI+GZjs7peW0da2smIgWpit2Gwgl+gnh+I3jLVEC7zll3wqqbHnSwi4+xYz+8zMLnH3V2O/q9nuvjDotjWlIzXVdGCAmdWPTasMJDraLn7st7F/ZxD9krcgVlRtNpBnZsXfA9Q3s1MTaG85se8OANx9K9Gbjwzl/yf8WcBNHJjwTwWWHmwHpXqIFXmbBZxmZmvN7FfAZcCvzGwh0f/7n8f2zY0VgrsEeNrMlu5znBnAq0Cv2HH6HmwsGuFLjeTuH5vZOGBObNP/uPuC2PIM4G5glrtvM7OdsW24+3ozGwq8ZGY/ie1/D1BRed2/Ay/ut20m0TrsxXedmgX8nn0SvpmdQPQTSmjKFodNOZ8WDzhV093nUsatGN29e2nbD4aqZYpUEjObBNzu7qsO4jk3A1vc/dngIhOJ0pSOSOW5k4M/22IT8JcAYhE5gEb4IiIhoRG+iEhIKOGLiISEEr6ISEgo4YuIhIQSvohISPw/yG0G/m7MkusAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Calculate power prediction from selecting 6 datapoints from the dataset randomly\n",
    "# This is done 20 times and then averaged, and is used to show the results that could be obtained\n",
    "# if datapoints were picked randomly instead of found using a robust clustering approach\n",
    "\n",
    "\n",
    "avg_fixed_power = []\n",
    "for i in range(int(np.floor(20))):\n",
    "    avg_fixed_power.append(np.sum([final_fixed_power[i] for i in [random.randint(0, 364) for _ in range(6)]])*365/6)\n",
    "    \n",
    "plt.hist(avg_fixed_power, label = 'Random')\n",
    "\n",
    "# Calculate the real power output for every datapoint in 2007\n",
    "plt.plot([reg_fixed_power_wind_d01[0].sum() for _ in range(16)], [i for i in range(16)], label = 'WRF')\n",
    "print(reg_fixed_power[0].sum())\n",
    "print(np.mean(avg_fixed_power))\n",
    "\n",
    "# Calculate power prediction using the more complex weighted approach\n",
    "power_0_factor_wind_d01 = power_cluster_0_wind_d01[-1].sum()*reg_fixed_power_0_wind_d01['norm_vel']/norm_vel_0_wind_d01\n",
    "power_0_factor_wind_d01[power_0_factor_wind_d01 > reg_fixed_power[0].max()] = reg_fixed_power[0].max()\n",
    "power_1_factor_wind_d01 = power_cluster_1_wind_d01[-1].sum()*reg_fixed_power_1_wind_d01['norm_vel']/norm_vel_1_wind_d01\n",
    "power_1_factor_wind_d01[power_1_factor_wind_d01 > reg_fixed_power[0].max()] = reg_fixed_power[0].max()\n",
    "power_2_factor_wind_d01 = power_cluster_2_wind_d01[-1].sum()*reg_fixed_power_2_wind_d01['norm_vel']/norm_vel_2_wind_d01\n",
    "power_2_factor_wind_d01[power_2_factor_wind_d01 > reg_fixed_power[0].max()] = reg_fixed_power[0].max()\n",
    "power_3_factor_wind_d01 = power_cluster_3_wind_d01[-1].sum()*reg_fixed_power_3_wind_d01['norm_vel']/norm_vel_3_wind_d01\n",
    "power_3_factor_wind_d01[power_3_factor_wind_d01 > reg_fixed_power[0].max()] = reg_fixed_power[0].max()\n",
    "power_4_factor_wind_d01 = power_cluster_4_wind_d01[-1].sum()*reg_fixed_power_4_wind_d01['norm_vel']/norm_vel_4_wind_d01\n",
    "power_4_factor_wind_d01[power_4_factor_wind_d01 > reg_fixed_power[0].max()] = reg_fixed_power[0].max()\n",
    "power_5_factor_wind_d01 = power_cluster_5_wind_d01[-1].sum()*reg_fixed_power_5_wind_d01['norm_vel']/norm_vel_5_wind_d01\n",
    "power_5_factor_wind_d01[power_5_factor_wind_d01 > reg_fixed_power[0].max()] = reg_fixed_power[0].max()\n",
    "\n",
    "fix_power_sum_wind_d01 = 0\n",
    "fix_power_sum_wind_d01 += power_0_factor_wind_d01.sum()\n",
    "fix_power_sum_wind_d01 += power_1_factor_wind_d01.sum()\n",
    "fix_power_sum_wind_d01 += power_2_factor_wind_d01.sum()\n",
    "fix_power_sum_wind_d01 += power_3_factor_wind_d01.sum()\n",
    "fix_power_sum_wind_d01 += power_4_factor_wind_d01.sum()\n",
    "fix_power_sum_wind_d01 += power_5_factor_wind_d01.sum()\n",
    "\n",
    "print(fix_power_sum_wind_d01)\n",
    "\n",
    "plt.plot([fix_power_sum_wind_d01 for _ in range(16)], [i for i in range(16)], label = 'Cluster weighted average')\n",
    "\n",
    "# Calculate the power output from using the simple averaged approach\n",
    "fix_power_sum_wind_d01_simple = 0\n",
    "fix_power_sum_wind_d01_simple += power_cluster_0_wind_d01[-1].sum()*len(reg_fixed_power_0_wind_d01)\n",
    "fix_power_sum_wind_d01_simple += power_cluster_1_wind_d01[-1].sum()*len(reg_fixed_power_1_wind_d01)\n",
    "fix_power_sum_wind_d01_simple += power_cluster_2_wind_d01[-1].sum()*len(reg_fixed_power_2_wind_d01)\n",
    "fix_power_sum_wind_d01_simple += power_cluster_3_wind_d01[-1].sum()*len(reg_fixed_power_3_wind_d01)\n",
    "fix_power_sum_wind_d01_simple += power_cluster_4_wind_d01[-1].sum()*len(reg_fixed_power_4_wind_d01)\n",
    "fix_power_sum_wind_d01_simple += power_cluster_5_wind_d01[-1].sum()*len(reg_fixed_power_5_wind_d01)\n",
    "\n",
    "print(fix_power_sum_wind_d01_simple)\n",
    "\n",
    "plt.plot([fix_power_sum_wind_d01_simple for _ in range(16)], [i for i in range(16)], label = 'Cluster simple average')\n",
    "plt.ylim([0, 6])\n",
    "plt.xlabel('Power (W)')\n",
    "plt.ylabel('Frequency')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  }
 ],
 "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.6.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}