SVR.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "51f84407",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "from datetime import datetime\n",
    "from multiprocessing import Lock, Process, freeze_support\n",
    "\n",
    "from sklearn.svm import SVR\n",
    "from statsmodels.tools.eval_measures import rmse\n",
    "from sklearn.pipeline import make_pipeline\n",
    "from sklearn.preprocessing import StandardScaler\n",
    "import os\n",
    "import time\n",
    "\n",
    "# objective function for PSO algorithm\n",
    "# returns RMSE error for PSO particles\n",
    "def SVRfx(particles, Xtrain, ytrain, Xtest, ytest):\n",
    "    rmses = []\n",
    "    for params in particles:\n",
    "        model = make_pipeline(\n",
    "            StandardScaler(),\n",
    "            SVR(kernel='rbf', epsilon=params[0], gamma=params[1], C=params[2])\n",
    "        )\n",
    "        model.fit(Xtrain,ytrain)\n",
    "        predictions = model.predict(Xtest)\n",
    "        rmses.append(rmse(ytest,predictions))\n",
    "    return rmses\n",
    "\n",
    "PSOn = 7 # run PSO optimization of SVR parameters each 7 days, note: the smaller the number the more computationally expensive calculation\n",
    "datasets = 'datasets' # folder name with datasets for individual periods in data folder\n",
    "FILE_NAME = 'SVR_PSO.csv' # name of file with results\n",
    "\n",
    "# PSO parameters\n",
    "max_bound = np.array([1,1,1000]) # epsilon, gamma, C\n",
    "min_bound = np.array([0.0001,0.0001,0.1])\n",
    "bounds = (min_bound, max_bound) # lower and upper bounds for epsilon, gamma and C\n",
    "options = {'c1': 0.5, 'c2': 0.5, 'w':0.2} # acceleration coefficients c1 and c2, inertia weight w\n",
    "\n",
    "# division of hours into 8 groups for paralelization\n",
    "hours_l = list(range(24))\n",
    "hour_groups = np.array_split(hours_l, 8)\n",
    "\n",
    "import pyswarms as ps\n",
    "from sklearn.model_selection import train_test_split\n",
    "\n",
    "write_lock = Lock()\n",
    "# SVR calculation for given hours - each hourly period is forecasted by a separate SVR model trained \n",
    "# only on data from that period\n",
    "def SVR_calc(hours):\n",
    "    global PSOn, datasets, FILE_NAME\n",
    "    global bounds, options\n",
    "\n",
    "    for hour in hours:\n",
    "        print(f'{str(hour).zfill(2)}:00')\n",
    "        # load data from given period\n",
    "        data = pd.read_csv(os.getcwd() + f'/data/SVR_dataset/{str(hour).zfill(2)}_00.csv',parse_dates=['datetime'],index_col='datetime').dropna()\n",
    "        #Dropping ewma only for curve features purpose\n",
    "        data.drop(['ewma'], axis = 1, inplace = True)\n",
    "\n",
    "        # create list of dates to forecast\n",
    "        dates = data.resample('D').count()[['grid1-loss']]\n",
    "        dates = dates[dates['grid1-loss']>0]\n",
    "        DATES = dates.loc['12-2019':].reset_index()['datetime'].apply(datetime.date).to_frame()['datetime']\n",
    "        dates = dates.reset_index()['datetime'].apply(datetime.date).to_frame()\n",
    "\n",
    "        i = 0\n",
    "        for date in DATES:\n",
    "            print(date)\n",
    "            test_date = date\n",
    "            date_index = dates.index[dates['datetime']==test_date][0]\n",
    "\n",
    "            # create train set for PSO parameter optimization from 90% of all dates preceding the current test date\n",
    "            # validation set is the last 10% of dates preceding the current test date\n",
    "            X_valtrain, X_val, y_valtrain, y_val = train_test_split(data[:str(data.iloc[date_index-1].name.date())].iloc[:,1:],\n",
    "                                                                    np.ravel(data[:str(data.iloc[date_index-1].name.date())].iloc[:,0]),\n",
    "                                                                    test_size=.1, shuffle=False, stratify=None)\n",
    "\n",
    "            # create train set from all dates preceding the current test date\n",
    "            # test set is the 1 current test date\n",
    "            X_train, X_test, y_train, y_test = train_test_split(data[:str(test_date)].iloc[:, 1:],\n",
    "                                                                np.ravel(data[:str(test_date)].iloc[:, 0]),\n",
    "                                                                test_size=1, shuffle=False, stratify=None)\n",
    "\n",
    "            # if the current iteration is the PSOnth date, optimize SVR parameters with PSO\n",
    "            if i % PSOn == 0:\n",
    "                optimizer = ps.single.GlobalBestPSO(n_particles=30, dimensions=3, options=options, bounds=bounds)\n",
    "                cost, pos = optimizer.optimize(SVRfx, iters=30, Xtrain=X_valtrain,ytrain=y_valtrain,Xtest=X_val,ytest=y_val,verbose=False)\n",
    "            i = i + 1\n",
    "\n",
    "            # create SVR model with feature normalization to (0,1) and forecast the next 15-minute period\n",
    "            model = make_pipeline(\n",
    "                StandardScaler(),\n",
    "                SVR(kernel='rbf', epsilon=pos[0], gamma=pos[1], C=pos[2])\n",
    "            )\n",
    "            model.fit(X_train,y_train)\n",
    "            predictions = model.predict(X_test)\n",
    "\n",
    "            # create forecasts dataframe and append it to file with results\n",
    "            forecasts = data[:str(test_date)].loc[str(test_date)][['grid1-loss']].reset_index() if hour > 0 else data[:str(test_date)].loc[[str(test_date)]][['grid1-loss']].reset_index()\n",
    "            forecasts.rename({'grid1-loss':'real'},axis=1,inplace=True)\n",
    "            forecasts['predicted'] = predictions\n",
    "\n",
    "            # In your SVR_calc function, when writing to the file:\n",
    "            with write_lock:\n",
    "                forecasts.to_csv(f'{FILE_NAME}', mode='a', index=False, header=False)\n",
    "\n",
    "if __name__ == \"__main__\":\n",
    "    freeze_support() \n",
    "    \n",
    "    # create empty file for results with header\n",
    "    f = open(f'{FILE_NAME}', \"w\")\n",
    "    f.write(\"datetime,real,predicted\\n\")\n",
    "    f.close()\n",
    "    \n",
    "    processes = list()\n",
    "    # run 8 paralel procecess \n",
    "    for hours in hour_groups:\n",
    "        p = Process(target = SVR_calc, args = (hours,))\n",
    "        processes.append(p)\n",
    "        p.start()\n",
    "    \n",
    "    # Wait for all processes to finish\n",
    "    for p in processes:\n",
    "        p.join()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8d8c4bc3",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.10.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}