wake_steering_and_loads.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Load Constrained Wake Steering Optimization"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Try this yourself](https://colab.research.google.com/github/DTUWindEnergy/TopFarm2/blob/master/docs/notebooks/wake_steering_and_loads.ipynb) (requires google account)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this example, we perform an active yaw (wake steering) optimization while considering the fatigue loading as constraint. Here, the design variables are changed to the yaw angles of the turbines, which are dependent on the number of turbines and the wind speed and wind direction cases selected. In addition, we need to specify the wake deflection model to represent the wake behind the turbine as accurate as possible.\n",
    "\n",
    "Load calculations of the Damage Equivalent Loads (DEL) and Lifetime Damage Equivalent Loads (LDEL) are computed via PyWake for the turbines and flow cases selected. Then, the load constraints can be included either with the `AEPMaxLoadCostModel` or as **post_constraints** in the `TopFarmProblem`. Both types of set ups are included in this example."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "Nv2ihk66SNgi"
   },
   "source": [
    "**Install TopFarm and PyWake if needed**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "executionInfo": {
     "elapsed": 67702,
     "status": "ok",
     "timestamp": 1623612898460,
     "user": {
      "displayName": "Mikkel Friis-Møller",
      "photoUrl": "https://lh3.googleusercontent.com/a-/AOh14GiVrnxS0oNYcvEfwYdFjWYAU5G0YxLXELnknXMi=s64",
      "userId": "10444369613733539918"
     },
     "user_tz": -120
    },
    "id": "HP8XVx4URYcr"
   },
   "outputs": [],
   "source": [
    "# Install TopFarm if needed\n",
    "import importlib\n",
    "if not importlib.util.find_spec(\"topfarm\"):\n",
    "    !pip install git+https://gitlab.windenergy.dtu.dk/TOPFARM/TopFarm2.git"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Note:** The `surrogates_interface` package is required to import `IEA34_130_1WT_Surrogate` from `py_wake`. If you encounter a `ModuleNotFoundError`, install `surrogates_interface` by running:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from py_wake.examples.data.iea34_130rwt import IEA34_130_1WT_Surrogate \n",
    "if not importlib.util.find_spec(\"surrogates_interface\"):\n",
    "    !pip install surrogates_interface==2.2.1 --index-url https://gitlab.windenergy.dtu.dk/api/v4/projects/2552/packages/pypi/simple --ignore-requires-python"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "TbtAki7QSZG4"
   },
   "source": [
    "**First we import basic Python elements and some Topfarm classes**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "executionInfo": {
     "elapsed": 4290,
     "status": "ok",
     "timestamp": 1623612902741,
     "user": {
      "displayName": "Mikkel Friis-Møller",
      "photoUrl": "https://lh3.googleusercontent.com/a-/AOh14GiVrnxS0oNYcvEfwYdFjWYAU5G0YxLXELnknXMi=s64",
      "userId": "10444369613733539918"
     },
     "user_tz": -120
    },
    "id": "5YqUNim5R3JG"
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from numpy import newaxis as na\n",
    "import time\n",
    "import matplotlib.pyplot as plt\n",
    "import warnings\n",
    "warnings.filterwarnings(\"ignore\")\n",
    "\n",
    "from topfarm.cost_models.cost_model_wrappers import AEPMaxLoadCostModelComponent\n",
    "from topfarm.easy_drivers import EasyScipyOptimizeDriver\n",
    "from topfarm import TopFarmProblem\n",
    "from topfarm.plotting import NoPlot\n",
    "\n",
    "from py_wake.examples.data.lillgrund import LillgrundSite\n",
    "from py_wake.deficit_models.gaussian import IEA37SimpleBastankhahGaussian, NiayifarGaussian\n",
    "from py_wake.turbulence_models.stf import STF2017TurbulenceModel\n",
    "from py_wake.examples.data.iea34_130rwt import IEA34_130_1WT_Surrogate \n",
    "from py_wake.deflection_models.jimenez import JimenezWakeDeflection\n",
    "from py_wake.superposition_models import MaxSum\n",
    "from py_wake.wind_turbines.power_ct_functions import SimpleYawModel"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "WjnaXixxSdZL"
   },
   "source": [
    "**Next, we select the site and wind turbine to use**.\n",
    "\n",
    "Usually, adding the loads as constraints into TOPFARM's problem requires an accurate calculation of the fatigue loading, which here is done by a surrogate of the IEA3 3.4MW turbine. In addition, it is necessary to specify a turbulence model (`STF2017TurbulenceModel`) that is adequate enough to represent the turbulence intensity in the site, which the surrogate model for the turbine will need for the load calculation. For the wake deflection, we will use the `JimenezWakeDeflection` model.\n",
    "\n",
    "For more information about wake models, please see: https://topfarm.pages.windenergy.dtu.dk/PyWake/notebooks/EngineeringWindFarmModels.html"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "executionInfo": {
     "elapsed": 1395,
     "status": "ok",
     "timestamp": 1623612904131,
     "user": {
      "displayName": "Mikkel Friis-Møller",
      "photoUrl": "https://lh3.googleusercontent.com/a-/AOh14GiVrnxS0oNYcvEfwYdFjWYAU5G0YxLXELnknXMi=s64",
      "userId": "10444369613733539918"
     },
     "user_tz": -120
    },
    "id": "E-TFCSUSSt5e"
   },
   "outputs": [],
   "source": [
    "site = LillgrundSite()\n",
    "windTurbines = IEA34_130_1WT_Surrogate()\n",
    "wfm = IEA37SimpleBastankhahGaussian(site, windTurbines,deflectionModel=JimenezWakeDeflection(), turbulenceModel=STF2017TurbulenceModel(addedTurbulenceSuperpositionModel=MaxSum()))\n",
    "\n",
    "#choosing the flow cases for the optimization - this will determine the speed and accuracy of the simulation\n",
    "wsp = np.asarray([10, 15])\n",
    "wdir = np.asarray([90])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Next, we set up the load constraint**\n",
    "\n",
    "In this example, we will calculate nominal loads and use this as a basis for the load constraint. The loads are represented by the Lifetime Damage Equivalent Loads (LDEL)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "executionInfo": {
     "elapsed": 1981,
     "status": "ok",
     "timestamp": 1623612906110,
     "user": {
      "displayName": "Mikkel Friis-Møller",
      "photoUrl": "https://lh3.googleusercontent.com/a-/AOh14GiVrnxS0oNYcvEfwYdFjWYAU5G0YxLXELnknXMi=s64",
      "userId": "10444369613733539918"
     },
     "user_tz": -120
    },
    "id": "hjG7faI1iZBh",
    "outputId": "618fbcaa-ac5b-431e-81b2-00005db43238"
   },
   "outputs": [],
   "source": [
    "%%capture\n",
    "x, y = site.initial_position.T\n",
    "\n",
    "#keeping only every second turbine as lillegrund turbines are approx. half the size of the iea 3.4MW\n",
    "x = x[::2]\n",
    "y = y[::2]\n",
    "n_wt = x.size\n",
    "\n",
    "#setting up the size of the yaw angle to represent the number of turbines, wind speeds and wind directions used\n",
    "i = n_wt\n",
    "k = wsp.size\n",
    "l = wdir.size\n",
    "yaw_zero = np.zeros((i, l, k))\n",
    "\n",
    "#choosing a load ratio for the constraint\n",
    "load_fact = 1.02\n",
    "simulationResult = wfm(x,y,wd=wdir, ws=wsp, yaw=yaw_zero, tilt=0)\n",
    "nom_loads = simulationResult.loads('OneWT')['LDEL'].values\n",
    "max_loads = nom_loads * load_fact"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "0v7KeHZvinpG"
   },
   "source": [
    "**Setting up optimization problem**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "executionInfo": {
     "elapsed": 13,
     "status": "ok",
     "timestamp": 1623612906112,
     "user": {
      "displayName": "Mikkel Friis-Møller",
      "photoUrl": "https://lh3.googleusercontent.com/a-/AOh14GiVrnxS0oNYcvEfwYdFjWYAU5G0YxLXELnknXMi=s64",
      "userId": "10444369613733539918"
     },
     "user_tz": -120
    },
    "id": "VXXj960rjJt4"
   },
   "outputs": [],
   "source": [
    "#parameters needed for the optimization\n",
    "maxiter = 5\n",
    "driver = EasyScipyOptimizeDriver(optimizer='SLSQP', maxiter=maxiter)\n",
    "yaw_min, yaw_max =  - 40, 40\n",
    "yaw_init = np.zeros((i, l, k))\n",
    "tol = 1e-8\n",
    "ec = 1e-4\n",
    "step = 1e-2\n",
    "\n",
    "#setting up cost function - aep and nominal loads calculation\n",
    "def aep_load_func(yaw_ilk):\n",
    "    simres = wfm(x, y, wd=wdir, ws=wsp, yaw=yaw_ilk, tilt=0)\n",
    "    aep = simres.aep().sum()\n",
    "    loads = simres.loads('OneWT')['LDEL'].values\n",
    "    return [aep, loads]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "S3Eq4hVNjWai"
   },
   "source": [
    "**Now, we can set up the gradient function used in the optimization**\n",
    "\n",
    "For some problems it is sufficient to rely on the automatic finite difference calculated by OpenMDAO or you can specify the explicit gradients from your model. In this case we don't have explicit gradients but the automatic finite difference is also inefficient, so we do a manual population of the Jacobian."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "executionInfo": {
     "elapsed": 11,
     "status": "ok",
     "timestamp": 1623612906114,
     "user": {
      "displayName": "Mikkel Friis-Møller",
      "photoUrl": "https://lh3.googleusercontent.com/a-/AOh14GiVrnxS0oNYcvEfwYdFjWYAU5G0YxLXELnknXMi=s64",
      "userId": "10444369613733539918"
     },
     "user_tz": -120
    },
    "id": "Ry6CRAZLj4VF"
   },
   "outputs": [],
   "source": [
    "s = nom_loads.shape[0]\n",
    "P_ilk = np.broadcast_to(simulationResult.P.values[na], (i, l, k))\n",
    "lifetime = float(60 * 60 * 24 * 365 * 20)\n",
    "f1zh = 10.0 ** 7.0\n",
    "lifetime_on_f1zh = lifetime / f1zh\n",
    "indices = np.arange(i * l * k).reshape((i, l, k))\n",
    "\n",
    "def aep_load_gradient(yaw_ilk):\n",
    "    simres0 = wfm(x, y, wd=wdir, ws=wsp, yaw=yaw_ilk, tilt=0)\n",
    "    aep0 = simres0.aep()\n",
    "    DEL0 = simulationResult.loads('OneWT')['DEL'].values\n",
    "    LDEL0 = simulationResult.loads('OneWT')['LDEL'].values\n",
    "    d_aep_d_yaw = np.zeros(i*l*k)\n",
    "    d_load_d_yaw = np.zeros((s * i, i * l * k))\n",
    "    for n in range(n_wt):\n",
    "        yaw_step = yaw_ilk.copy()\n",
    "        yaw_step = yaw_step.reshape(i, l, k)\n",
    "        yaw_step[n, :, :] += step\n",
    "        simres_fd = wfm(x, y, wd=wdir, ws=wsp, yaw=yaw_step, tilt=0)\n",
    "        aep_fd = simres_fd.aep()\n",
    "        d_aep_d_yaw[n * l * k : (n + 1) * l * k] = (((aep_fd.values - aep0.values) / step).sum((0))).ravel()\n",
    "        \n",
    "        DEL_fd = simres_fd.loads('OneWT')['DEL'].values\n",
    "        for _ls in range(s):\n",
    "            m = simulationResult.loads('OneWT').m.values[_ls]\n",
    "            for _wd in range(l):\n",
    "                for _ws in range(k):\n",
    "                    DEL_fd_fc = DEL0.copy()\n",
    "                    DEL_fd_fc[:, :, _wd, _ws] = DEL_fd[:, :, _wd, _ws]\n",
    "                    DEL_fd_fc = DEL_fd_fc[_ls, :, :, :]\n",
    "                    f = DEL_fd_fc.mean()\n",
    "                    LDEL_fd = (((P_ilk * (DEL_fd_fc/f) ** m).sum((1, 2)) * lifetime_on_f1zh) ** (1/m))*f\n",
    "                    d_load_d_yaw[n_wt * _ls : n_wt * (_ls + 1), indices[n, _wd, _ws]] = (LDEL_fd - LDEL0[_ls]) / step\n",
    "\n",
    "    return d_aep_d_yaw, d_load_d_yaw"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "TfbCb-z2j6nd"
   },
   "source": [
    "**Specifying the cost component and topfarm problem**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "executionInfo": {
     "elapsed": 12,
     "status": "ok",
     "timestamp": 1623612906115,
     "user": {
      "displayName": "Mikkel Friis-Møller",
      "photoUrl": "https://lh3.googleusercontent.com/a-/AOh14GiVrnxS0oNYcvEfwYdFjWYAU5G0YxLXELnknXMi=s64",
      "userId": "10444369613733539918"
     },
     "user_tz": -120
    },
    "id": "5nqGrkbukCfd"
   },
   "outputs": [],
   "source": [
    "cost_comp = AEPMaxLoadCostModelComponent(input_keys=[('yaw_ilk', np.zeros((i, l, k)))],\n",
    "                                          n_wt = n_wt,\n",
    "                                          aep_load_function = aep_load_func,\n",
    "                                          aep_load_gradient = aep_load_gradient,\n",
    "                                          max_loads = max_loads, \n",
    "                                          objective=True,\n",
    "                                          income_model=True,\n",
    "                                          output_keys=[('AEP', 0), ('loads', np.zeros((s, i)))]\n",
    "                                          )\n",
    "\n",
    "problem = TopFarmProblem(design_vars={'yaw_ilk': (yaw_init, yaw_min, yaw_max)},\n",
    "                          n_wt=n_wt,\n",
    "                          cost_comp=cost_comp,\n",
    "                          driver=EasyScipyOptimizeDriver(optimizer='SLSQP', maxiter=maxiter, tol=tol),\n",
    "                          plot_comp=NoPlot(),\n",
    "                          expected_cost=ec)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "Jt_GXqoHkV_W"
   },
   "source": [
    "**Now we run the optimization**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "executionInfo": {
     "elapsed": 538141,
     "status": "ok",
     "timestamp": 1623613843204,
     "user": {
      "displayName": "Mikkel Friis-Møller",
      "photoUrl": "https://lh3.googleusercontent.com/a-/AOh14GiVrnxS0oNYcvEfwYdFjWYAU5G0YxLXELnknXMi=s64",
      "userId": "10444369613733539918"
     },
     "user_tz": -120
    },
    "id": "kOQhnRLUSGIz",
    "outputId": "6e174217-d355-472b-e3ef-9be7dee337da"
   },
   "outputs": [],
   "source": [
    "%%capture\n",
    "cost, state, recorder = problem.optimize()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**We can also do some plotting to visualize the final yaw angles found by the optimizer for a given flow case**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 293
    },
    "executionInfo": {
     "elapsed": 10125,
     "status": "ok",
     "timestamp": 1623613892575,
     "user": {
      "displayName": "Mikkel Friis-Møller",
      "photoUrl": "https://lh3.googleusercontent.com/a-/AOh14GiVrnxS0oNYcvEfwYdFjWYAU5G0YxLXELnknXMi=s64",
      "userId": "10444369613733539918"
     },
     "user_tz": -120
    },
    "id": "miJIu_aMtvwy",
    "outputId": "980c94bb-a157-45aa-de6a-6a2322a9dc58"
   },
   "outputs": [],
   "source": [
    "%%capture\n",
    "simulationResult = wfm(x,y,wd=wdir[0], ws=wsp[0], yaw=state['yaw_ilk'][:,0,0], tilt=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "simulationResult.flow_map().plot_wake_map()\n",
    "plt.xlabel('x [m]')\n",
    "plt.ylabel('y [m]')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Note:\n",
    "\n",
    "Due to the computational expensive nature of this problem, the maximum number of iterations were restricted to 5, for educational purposes. To see a convergence in loading and AEP, the number of iterations must be increased."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Setting up the load constraints as post_constraints"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Alternatively to using the `AEPMaxLoadCostModelComponent`, the common `CostModelComponent` can also be used any load constrained problem. The main difference lies in the way the load constraints are specified within the `TopFarmProblem`. You should set up the load constraint as a **post_constraints** element. Below is an example on how to set up this type of problem."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from topfarm.cost_models.cost_model_wrappers import CostModelComponent\n",
    "\n",
    "#setting up the cost model component and topfarm problem\n",
    "\n",
    "cost_comp = CostModelComponent([('yaw_ilk', np.zeros((i, l, k)))],\n",
    "                            n_wt=n_wt,\n",
    "                            cost_function=aep_load_func,\n",
    "                            cost_gradient_function=aep_load_gradient, \n",
    "                            objective=True,\n",
    "                            output_keys=[('AEP', 0.0), ('loads', np.zeros((s, i)))])\n",
    "\n",
    "#parameters for optimization\n",
    "maxiter = 5\n",
    "tol = 1e-4\n",
    "ec = 1\n",
    "\n",
    "problem = TopFarmProblem(design_vars={'yaw_ilk': (yaw_init, yaw_min, yaw_max)},\n",
    "                          n_wt = n_wt,\n",
    "                          post_constraints=[('loads', {'upper': max_loads})],\n",
    "                          cost_comp=cost_comp,\n",
    "                          driver=EasyScipyOptimizeDriver(optimizer='SLSQP', maxiter=maxiter, tol=tol),\n",
    "                          plot_comp=NoPlot(),\n",
    "                          expected_cost=ec)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can evaluate the topfarm problem to make sure it is set up properly without having to run the optimization again."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%capture\n",
    "problem.evaluate()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Evaluating the topfarm problem yields the initial conditions in terms of AEP and position of the turbines:\n",
    "\n",
    "    (72.12110945807837,\n",
    "     {'yaw_ilk': array([[[0., 0.]],\n",
    "  \n",
    "         [[0., 0.]],\n",
    "  \n",
    "         [[0., 0.]],\n",
    "  \n",
    "         [[0., 0.]],\n",
    "  \n",
    "         [[0., 0.]],\n",
    "  \n",
    "         [[0., 0.]],\n",
    "  \n",
    "         [[0., 0.]],\n",
    "  \n",
    "         [[0., 0.]],\n",
    "  \n",
    "         [[0., 0.]],\n",
    "  \n",
    "         [[0., 0.]],\n",
    "  \n",
    "         [[0., 0.]],\n",
    "  \n",
    "         [[0., 0.]],\n",
    "  \n",
    "         [[0., 0.]],\n",
    "  \n",
    "         [[0., 0.]],\n",
    "  \n",
    "         [[0., 0.]],\n",
    "  \n",
    "         [[0., 0.]],\n",
    "  \n",
    "         [[0., 0.]],\n",
    "  \n",
    "         [[0., 0.]],\n",
    "  \n",
    "         [[0., 0.]],\n",
    "  \n",
    "         [[0., 0.]],\n",
    "  \n",
    "         [[0., 0.]],\n",
    "  \n",
    "         [[0., 0.]],\n",
    "  \n",
    "         [[0., 0.]],\n",
    "  \n",
    "         [[0., 0.]]])})"
   ]
  }
 ],
 "metadata": {
  "colab": {
   "authorship_tag": "ABX9TyOvPUxm8Sna0KLEzObl3Sp5",
   "name": "WakeSteeringOptimization.ipynb",
   "provenance": []
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}