problems example elaborated. its noted however that make_notebooks is not...

problems example elaborated. its noted however that make_notebooks is not working due to a bug in the drivers.ipynb

.gitignore

@@ -112,3 +112,6 @@ venv.bak/
 # mypy
 .mypy_cache/
 /.githooks
+*.DS_Store
+__pycache__
+*.egg-info

_notebooks/elements/problems.ipynb

@@ -4,7 +4,8 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "# Problems"
+    "# Problems\n",
+    "In Topfarm the problem is the collection/container of components and drivers. Optimizing a problem is executing a workflow - finding the best feasible solution under a series of constraints using a specific driver, and initial conditions. To the problem is also attached a plotting routine so one can follow the optimization path as well as a recorder so that intermediate variable values are accessible after the optimization is finished. The Topfarm Problem inherits its fundamental nature from the OpenMDAO problem, and is being adapted so it can connect the given workflow with the driver. For example if the user specifies a boundary constraint and this is not supported by the driver, the Problem is equipped with a penalty component that will deter the driver giving unfeasible solutions. Or if your workflow does not have gradients for all components and a gradient based driver is specified, finite differencing is applied to obtain the gradients. "
    ]
   },
   {
@@ -18,6 +19,13 @@
     "Make sure you first run the code below in order to set up and initialize needed variables."
    ]
   },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "**First we install topfarm and import supporting libraries in Python numpy and matplotlib**"
+   ]
+  },
   {
    "cell_type": "code",
    "execution_count": 1,
@@ -28,33 +36,52 @@
     "try:  # install Topfarm if needed\n",
     "    import topfarm\n",
     "except ModuleNotFoundError:\n",
-    "    !pip install topfarm"
+    "    !pip install topfarm\n",
+    "import numpy as np"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 2,
+   "execution_count": 3,
    "metadata": {},
    "outputs": [],
    "source": [
-    "import matplotlib.pyplot as plt \n",
-    "\n",
-    "from topfarm import TopFarmProblem\n",
-    "from topfarm.easy_drivers import EasyScipyOptimizeDriver\n",
-    "from topfarm.examples.iea37 import get_iea37_initial, get_iea37_constraints, get_iea37_cost\n",
-    "from topfarm.plotting import NoPlot, XYPlotComp"
+    "# non-updating, inline plots\n",
+    "#%matplotlib inline\n",
+    "# ...or updating plots in new window\n",
+    "%matplotlib qt"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 3,
+   "execution_count": 4,
    "metadata": {},
    "outputs": [],
    "source": [
-    "# non-updating, inline plots\n",
-    "%matplotlib inline\n",
-    "# ...or updating plots in new window\n",
-    "# %matplotlib qt"
+    "import matplotlib.pyplot as plt "
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "**Next we import and initialize several functions and classes from Topfarm to set up the problem including:**\n",
+    "* **TopFarmProblem - overall topfarm problem class to which the objectives, design variables, and constraints are added**\n",
+    "* **EasyScipyOptimizeDriver - a subclass of ScipyOptimizeDriver which is configured for the given workflow**\n",
+    "* **get_iea37_initial, get_iea37_constraints, get_iea37_cost - functions to get the initial layout, the constraints and the cost function for the IEA task 37 benchmark example**\n",
+    "* **NoPlot, XYPlotComp - plotting components to visualize the results**"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 5,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from topfarm import TopFarmProblem\n",
+    "from topfarm.easy_drivers import EasyScipyOptimizeDriver\n",
+    "from topfarm.examples.iea37 import get_iea37_initial, get_iea37_constraints, get_iea37_cost\n",
+    "from topfarm.plotting import NoPlot, XYPlotComp"
    ]
   },
   {
@@ -70,42 +97,73 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 4,
+   "execution_count": 12,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from topfarm.examples.iea37 import *\n",
+    "def get_iea37_cost(n_wt=9, n_wd=16):\n",
+    "    \"\"\"Cost component that wraps the IEA 37 AEP calculator\"\"\"\n",
+    "    wd = np.linspace(0., 360., n_wd, endpoint=False)\n",
+    "    site = IEA37Site(n_wt)\n",
+    "    wind_turbines = IEA37_WindTurbines()\n",
+    "    wake_model = IEA37SimpleBastankhahGaussian(site, wind_turbines)\n",
+    "    aep_calc = PyWakeAEP(wake_model)\n",
+    "    return aep_calc.get_TopFarm_cost_component(n_wt, wd=wd)\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 26,
    "metadata": {},
    "outputs": [],
    "source": [
-    "def optimize_iea37_locs(n_wt, driver):\n",
+    "def optimize_iea37_locs(n_wt, n_wd, driver, state=None):\n",
+    "    \"\"\"\n",
+    "    Parameters\n",
+    "    ----------\n",
+    "    - n_wt: int\n",
+    "        Number of wind turbines\n",
+    "\n",
+    "    - n_wd: int\n",
+    "        Number of wind directions to consider for the AEP\n",
+    "        \n",
+    "    - driver: TopfarmDriver instance\n",
+    "        The optimization algorithm to use\n",
+    "        \n",
+    "    - state: dict(x=[], y=[]) [default=None]\n",
+    "        State to start from the optimization\n",
+    "    \n",
+    "    Returns\n",
+    "    -------\n",
+    "    - state: The state after the optimization\n",
+    "\n",
+    "    \"\"\"\n",
+    "    \n",
     "    initial = get_iea37_initial(n_wt)\n",
+    "    design_vars = dict(zip('xy', (initial[:, :2]).T))\n",
+    "    \n",
     "    tf = TopFarmProblem(\n",
-    "            dict(zip('xy', (initial[:, :2]).T)),\n",
-    "            get_iea37_cost(n_wt),\n",
+    "            design_vars,\n",
+    "            get_iea37_cost(n_wt, n_wd=n_wd),\n",
     "            constraints=get_iea37_constraints(n_wt),\n",
     "            driver=driver,\n",
     "            plot_comp=XYPlotComp())\n",
-    "    tf.optimize()"
+    "    \n",
+    "    if not state:\n",
+    "        _, state = tf.evaluate()\n",
+    "        \n",
+    "    _, state, _ = tf.optimize(state)\n",
+    "    return state, "
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 5,
+   "execution_count": 27,
    "metadata": {},
-   "outputs": [
-    {
-     "ename": "TypeError",
-     "evalue": "__init__() missing 1 required positional argument: 'windTurbines'",
-     "output_type": "error",
-     "traceback": [
-      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
-      "\u001b[1;31mTypeError\u001b[0m                                 Traceback (most recent call last)",
-      "\u001b[1;32m<ipython-input-5-df00a1eefbe1>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m\u001b[0m\n\u001b[1;32m----> 1\u001b[1;33m \u001b[0moptimize_iea37_locs\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;36m9\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mEasyScipyOptimizeDriver\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mdisp\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;32mFalse\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m",
-      "\u001b[1;32m<ipython-input-4-c9615be881ef>\u001b[0m in \u001b[0;36moptimize_iea37_locs\u001b[1;34m(n_wt, driver)\u001b[0m\n\u001b[0;32m      3\u001b[0m     tf = TopFarmProblem(\n\u001b[0;32m      4\u001b[0m             \u001b[0mdict\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mzip\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;34m'xy'\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;33m(\u001b[0m\u001b[0minitial\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;33m:\u001b[0m\u001b[1;36m2\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mT\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m----> 5\u001b[1;33m             \u001b[0mget_iea37_cost\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mn_wt\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m      6\u001b[0m             \u001b[0mconstraints\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mget_iea37_constraints\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mn_wt\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      7\u001b[0m             \u001b[0mdriver\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mdriver\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
-      "\u001b[1;32mc:\\mmpe\\programming\\python\\topfarm\\topfarm2\\topfarm\\examples\\iea37\\_iea37.py\u001b[0m in \u001b[0;36mget_iea37_cost\u001b[1;34m(n_wt)\u001b[0m\n\u001b[0;32m     50\u001b[0m     \u001b[0msite\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mIEA37Site\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mn_wt\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     51\u001b[0m     \u001b[0mwind_turbines\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mIEA37_WindTurbines\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m---> 52\u001b[1;33m     \u001b[0mwake_model\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mIEA37SimpleBastankhahGaussian\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mwind_turbines\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m     53\u001b[0m     \u001b[0maep_calc\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mPyWakeAEP\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0msite\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mwind_turbines\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mwake_model\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     54\u001b[0m     \u001b[1;32mreturn\u001b[0m \u001b[0maep_calc\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mget_TopFarm_cost_component\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mn_wt\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mwd\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mwd\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
-      "\u001b[1;31mTypeError\u001b[0m: __init__() missing 1 required positional argument: 'windTurbines'"
-     ]
-    }
-   ],
+   "outputs": [],
    "source": [
-    "optimize_iea37_locs(9, EasyScipyOptimizeDriver(disp=False))"
+    "state = optimize_iea37_locs(9, n_wd=16, EasyScipyOptimizeDriver(disp=False))"
    ]
   },
   {
@@ -114,13 +172,38 @@
    "source": [
     "The final optimized output is much lower than the value reported in Baker et al. (2019), which is 257.790 GWh. Moreover, the layout does not match the figures given in Appendix A in the same reference. This is due to the fact that the SLSQP optimizer was attracted to a local minimum. To find the global optimum, more advanced optimization procedures should be used. This benchmark is discussed in more detail in the validation report linked in TOPFARM's documentation."
    ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 15,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "{'x': array([-680.80198103, -854.32314245, -695.66777104,   89.8040452 ,\n",
+       "          11.54426432,  -64.22919445,  708.93660969,  860.87425685,\n",
+       "         677.89454027]),\n",
+       " 'y': array([-588.64986431,  283.07590553,  571.00468684, -895.50836595,\n",
+       "           9.60717992,  897.70519161, -554.44466226, -262.4795496 ,\n",
+       "         591.9957705 ])}"
+      ]
+     },
+     "execution_count": 15,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "state"
+   ]
   }
  ],
  "metadata": {
   "kernelspec": {
-   "display_name": "Python 3",
+   "display_name": "Python [conda env:topfarm_clean]",
    "language": "python",
-   "name": "python3"
+   "name": "conda-env-topfarm_clean-py"
   },
   "language_info": {
    "codemirror_mode": {
@@ -132,9 +215,9 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.6.7"
+   "version": "3.7.3"
   }
  },
  "nbformat": 4,
- "nbformat_minor": 2
+ "nbformat_minor": 4
 }

_notebooks/elements/roads_and_cables.ipynb

+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "%%capture\n",
+    "try:\n",
+    "    !git clone https://gitlab.windenergy.dtu.dk/TOPFARM/workshop-december-2019.git\n",
+    "except:\n",
+    "    pass\n",
+    "# !cd workshop-december-2019\n",
+    "!pip install -e ./workshop-december-2019\n",
+    "# !cd wind2loads\n",
+    "!pip install ./workshop-december-2019/wind2loads"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "In colab, use the \"inline\" backend"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# non-updating, inline plots\n",
+    "%matplotlib inline\n",
+    "# ...or updating plots in new window\n",
+    "#%matplotlib qt"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Let's import a few classes"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import numpy as np\n",
+    "# from openmdao.api import view_model\n",
+    "# from py_wake.examples.data.iea37._iea37 import IEA37_WindTurbines, IEA37Site\n",
+    "from topfarm.cost_models.cost_model_wrappers import CostModelComponent\n",
+    "from topfarm import TopFarmGroup, TopFarmProblem\n",
+    "from topfarm.easy_drivers import EasyRandomSearchDriver, EasyScipyOptimizeDriver\n",
+    "from topfarm.drivers.random_search_driver import RandomizeTurbinePosition_Circle\n",
+    "# from topfarm.constraint_components.boundary import CircleBoundaryConstraint\n",
+    "from topfarm.constraint_components.spacing import SpacingConstraint\n",
+    "from topfarm.constraint_components.boundary import XYBoundaryConstraint\n",
+    "from topfarm.cost_models.electrical.simple_msp import ElNetLength, ElNetCost \n",
+    "\n",
+    "\n",
+    "import matplotlib.pylab as plt\n",
+    "\n",
+    "# Convenient plot component for plotting the cables from simple_msp in jupyter notebooks\n",
+    "from workshop.cabling import XYCablePlotComp, get_site"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Setting up the site to optimize\n",
+    "We will use the IEA-37 site, using the DTU 10MW reference turbine"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from workshop.cabling import get_site\n",
+    "from py_wake.examples.data.DTU10MW_RWT import DTU10MW\n",
+    "\n",
+    "\n",
+    "\n",
+    "site = get_site()\n",
+    "n_wt = len(site.initial_position)\n",
+    "windTurbines = DTU10MW()\n",
+    "Drotor_vector = [windTurbines.diameter()] * n_wt \n",
+    "power_rated_vector = [float(windTurbines.power(20)/1000)] * n_wt \n",
+    "hub_height_vector = [windTurbines.hub_height()] * n_wt \n",
+    "rated_rpm_array = 12. * np.ones([n_wt])\n",
+    "\n",
+    "print('Number of turbines:', n_wt)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Setting up the AEP calculator\n",
+    "- Using the Gaussian wake model from Bastankhah & Porté Agel\n",
+    "- Based on 16 wind direction to speed things up (not critical here because we will be using the RandomSearch algorithm)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from py_wake.wake_models.gaussian import IEA37SimpleBastankhahGaussian\n",
+    "from py_wake.aep_calculator import AEPCalculator\n",
+    "\n",
+    "## We use the Gaussian wake model\n",
+    "wake_model = IEA37SimpleBastankhahGaussian(site, windTurbines)\n",
+    "AEPCalc = AEPCalculator(wake_model)\n",
+    "\n",
+    "## The AEP is calculated using n_wd wind directions\n",
+    "n_wd = 16\n",
+    "wind_directions = np.linspace(0., 360., n_wd, endpoint=False)\n",
+    "\n",
+    "def aep_func(x, y, **kwargs):\n",
+    "    \"\"\"A simple function that takes as input the x,y position of the turbines and return the AEP\"\"\"\n",
+    "    return AEPCalc.calculate_AEP(x_i=x, y_i=y, wd=wind_directions).sum(-1).sum(-1)*10**6"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Setting up the NREL IRR cost model\n",
+    "Based on the 2006 NREL report"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from topfarm.cost_models.economic_models.turbine_cost import economic_evaluation as EE_NREL\n",
+    "\n",
+    "def irr_nrel(aep, electrical_connection_cost, **kwargs):\n",
+    "    return EE_NREL(Drotor_vector, power_rated_vector, hub_height_vector, aep, electrical_connection_cost).calculate_irr()\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Setting up the DTU IRR cost model\n",
+    "Based on Witold's recent work"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from topfarm.cost_models.economic_models.dtu_wind_cm_main import economic_evaluation as EE_DTU\n",
+    "\n",
+    "distance_from_shore = 10.0 # [km]\n",
+    "energy_price = 0.2 / 7.4 # [DKK/kWh] / [DKK/EUR] -> [EUR/kWh]\n",
+    "project_duration = 20 # [years]\n",
+    "water_depth_array = 20 * np.ones([n_wt])\n",
+    "Power_rated_array = np.array(power_rated_vector)/1.0E3 # [MW]\n",
+    "\n",
+    "ee_dtu = EE_DTU(distance_from_shore, energy_price, project_duration)\n",
+    "\n",
+    "\n",
+    "def irr_dtu(aep, electrical_connection_cost, **kwargs):\n",
+    "    ee_dtu.calculate_irr(\n",
+    "                    rated_rpm_array, \n",
+    "                    Drotor_vector, \n",
+    "                    Power_rated_array,\n",
+    "                    hub_height_vector, \n",
+    "                    water_depth_array, \n",
+    "                    aep, \n",
+    "                    electrical_connection_cost)\n",
+    "    return ee_dtu.IRR"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Setting up the Topfarm problem"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "## Some user options\n",
+    "#@markdown Which IRR Cost model to use\n",
+    "IRR_COST = 'DTU' #@param [\"DTU\", \"NREL\"]\n",
+    "\n",
+    "#@markdown Minimum spacing between the turbines\n",
+    "min_spacing = 2 #@param {type:\"slider\", min:2, max:10, step:1}\n",
+    "\n",
+    "#@markdown Minimum spacing between the turbines\n",
+    "cable_cost_per_meter = 750. # {type:\"slider\", min:0, max:10000, step:1}\n",
+    "\n",
+    "## Electrical grid cable components (Minimum spanning tree from Topfarm report 2010)\n",
+    "elnetlength = ElNetLength(n_wt=n_wt)\n",
+    "elnetcost = ElNetCost(n_wt=n_wt, output_key='electrical_connection_cost', cost_per_meter=cable_cost_per_meter)\n",
+    "\n",
+    "# The Topfarm IRR cost model components\n",
+    "irr_dtu_comp = CostModelComponent(input_keys=['aep', ('electrical_connection_cost', 0.0)],   n_wt=n_wt, \n",
+    "                                  cost_function=irr_dtu, output_key=\"irr\", output_unit=\"%\",   objective=True, income_model=True)\n",
+    "irr_nrel_comp = CostModelComponent(input_keys=['aep', ('electrical_connection_cost', 0.0)],   n_wt=n_wt, \n",
+    "                                   cost_function=irr_nrel, output_key=\"irr\", output_unit=\"%\",   objective=True, income_model=True)\n",
+    "irr_cost_models = {'DTU': irr_dtu_comp, 'NREL': irr_nrel_comp}\n",
+    "\n",
+    "\n",
+    "## The Topfarm AEP component\n",
+    "aep_comp = CostModelComponent(input_keys=['x','y'], n_wt=n_wt, cost_function=aep_func, \n",
+    "                              output_key=\"aep\", output_unit=\"GWh\", objective=False, output_val=np.zeros(n_wt))\n",
+    "\n",
+    "## Plotting component\n",
+    "plot_comp = XYCablePlotComp(memory=0, plot_improvements_only=False, plot_initial=False)\n",
+    "\n",
+    "\n",
+    "## The group containing all the components\n",
+    "group = TopFarmGroup([aep_comp, elnetlength, elnetcost, irr_cost_models[IRR_COST]])\n",
+    "\n",
+    "problem = TopFarmProblem(\n",
+    "        design_vars=dict(zip('xy', site.initial_position.T)),\n",
+    "        cost_comp=group,\n",
+    "        driver=EasyRandomSearchDriver(randomize_func=RandomizeTurbinePosition_Circle(), max_iter=1000),\n",
+    "        constraints=[SpacingConstraint(min_spacing * windTurbines.diameter(0)),\n",
+    "                     XYBoundaryConstraint(site.boundary)],\n",
+    "        expected_cost=1.0,\n",
+    "        plot_comp=plot_comp)\n",
+    "\n",
+    "cost, state, recorder = problem.optimize()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Exercises\n",
+    "- Try to see what is the effect of increasing or decreasing the cost of the cable\n",
+    "- Change between IRR cost model. Ask Witold about the difference between DTU and NREL models"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "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.7.3"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 4
+}

docs/source/api_reference/dtucost.rst

@@ -5,7 +5,7 @@ DTU Cost Model
 
 
 
-.. autoclass:: topfarm.cost_models.dtu_wind_cm.dtu_wind_cm_main.economic_evaluation
+.. autoclass:: topfarm.cost_models.economic_models.dtu_wind_cm_main.economic_evaluation
     
     .. autosummary::
         __init__

docs/source/user_guide.rst

@@ -61,7 +61,8 @@ DTU Cost Model
 	called individually. It generally relies on curve fitting using input parameters such as 
 	rated power or water depth, and was tuned using data obtained from the 
 	industry. It supports three types of drivetrains: high-speed, medium-speed and direct-drive.
-	It also supports two types of offshore foundations: monopile and jacket. 
+	It also supports two types of offshore foundations: monopile and jacket. The model was originally created in Excel by Kenneth Thomse and
+	Søren Oemann Lind and Rasmus (publication forthcoming). The monetary units are 2017 Euro. 
 
 	
 Drivers

examples/docs/example_8_irr_opt_with_capacityconst_turbine_type.py

@@ -52,7 +52,7 @@ def main(obj=False, max_con_on=True):
     
     # ---------------- DEFINE SITE & SELECT WAKE MODEL -------------------
 #        site = UniformWeibullSite(p_wd=[50, 50], a=[9, 9], k=[2.3, 2.3], ti=.1, alpha=0, h_ref=100)
-    site = UniformWeibullSite(p_wd=[100], a=[9], k=[2.3], ti=.1, alpha=0, h_ref=100)
+    site = UniformWeibullSite(p_wd=[100], a=[9], k=[2.3], ti=.1)
     site.default_ws = [9]  # reduce the number of calculations
     site.default_wd = [0]  # reduce the number of calculations
 

examples/docs/example_9_irr_opt_on_weibull_site.py

+import numpy as np
+from openmdao.api import view_model
+from py_wake.examples.data.iea37._iea37 import IEA37_WindTurbines, IEA37Site
+from py_wake.wake_models.gaussian import IEA37SimpleBastankhahGaussian
+from py_wake.aep_calculator import AEPCalculator
+from topfarm.cost_models.economic_models.turbine_cost import economic_evaluation
+from topfarm.cost_models.cost_model_wrappers import CostModelComponent
+from topfarm import TopFarmGroup, TopFarmProblem
+from topfarm.easy_drivers import EasyRandomSearchDriver
+from topfarm.drivers.random_search_driver import RandomizeTurbinePosition_Circle
+from topfarm.constraint_components.boundary import CircleBoundaryConstraint
+from topfarm.plotting import XYPlotComp, NoPlot
+from topfarm.constraint_components.spacing import SpacingConstraint
+
+from py_wake.examples.data.DTU10MW_RWT import DTU10MW
+from py_wake.site import WaspGridSiteBase, UniformWeibullSite
+from py_wake.site.shear import PowerShear
+from topfarm.constraint_components.boundary import XYBoundaryConstraint
+
+
+def get_site():
+    f = [0.035972, 0.039487, 0.051674, 0.070002, 0.083645, 0.064348,
+         0.086432, 0.117705, 0.151576, 0.147379, 0.10012, 0.05166]
+    A = [9.176929, 9.782334, 9.531809, 9.909545, 10.04269, 9.593921,
+     9.584007, 10.51499, 11.39895, 11.68746, 11.63732, 10.08803]
+    k = [2.392578, 2.447266, 2.412109, 2.591797, 2.755859, 2.595703,
+     2.583984, 2.548828, 2.470703, 2.607422, 2.626953, 2.326172]
+    ti = 0.001
+    h_ref = 100
+    alpha = .1
+    site = UniformWeibullSite(f, A, k, ti, shear=PowerShear(h_ref=h_ref, alpha=alpha))
+    spacing = 2000
+    N = 5
+    theta = 76 # deg
+    dx = np.tan(np.radians(theta))
+    x = np.array([np.linspace(0,(N-1)*spacing,N)+i*spacing/dx for i in range(N)])
+    y = np.array(np.array([N*[i*spacing] for i in range(N)]))
+    initial_positions = np.column_stack((x.ravel(),y.ravel()))
+    eps = 2000
+    delta = 5
+    site.boundary = np.array([(0-delta, 0-delta),
+                     ((N-1)*spacing+eps, 0-delta),
+                     ((N-1)*spacing*(1+1/dx)+eps*(1+np.cos(np.radians(theta))), (N-1)*spacing+eps*np.sin(np.radians(theta))-delta),
+                     ((N-1)*spacing/dx+eps*np.cos(np.radians(theta)), (N-1)*spacing+eps*np.sin(np.radians(theta)))])
+    site.initial_position = initial_positions
+    return site
+
+def main():
+    if __name__ == '__main__':
+        plot_comp = XYPlotComp()
+        site = get_site()
+        n_wt = len(site.initial_position)
+        windTurbines = DTU10MW()
+        min_spacing = 2 * windTurbines.diameter(0)
+        wake_model = IEA37SimpleBastankhahGaussian(site, windTurbines)
+        Drotor_vector = [windTurbines.diameter()] * n_wt 
+        power_rated_vector = [float(windTurbines.power(20)/1000)] * n_wt 
+        hub_height_vector = [windTurbines.hub_height()] * n_wt 
+        AEPCalc = AEPCalculator(wake_model)         
+
+        def aep_func(x, y, **kwargs):
+            return AEPCalc.calculate_AEP(x_i=x, y_i=y).sum(-1).sum(-1)*10**6
+        
+        def irr_func(aep, **kwargs):
+            return economic_evaluation(Drotor_vector, power_rated_vector, hub_height_vector, aep).calculate_irr()
+        
+        aep_comp = CostModelComponent(input_keys=['x','y'], n_wt=n_wt, cost_function=aep_func, output_key="aep", output_unit="GWh", objective=False, output_val=np.zeros(n_wt))
+        irr_comp = CostModelComponent(input_keys=['aep'],   n_wt=n_wt, cost_function=irr_func, output_key="irr", output_unit="%",   objective=True, income_model=True)
+        group = TopFarmGroup([aep_comp, irr_comp])
+        problem = TopFarmProblem(
+                design_vars=dict(zip('xy', site.initial_position.T)),
+                cost_comp=group,
+                driver=EasyRandomSearchDriver(randomize_func=RandomizeTurbinePosition_Circle(), max_iter=50),
+            constraints=[SpacingConstraint(min_spacing),
+                         XYBoundaryConstraint(site.boundary),],
+                plot_comp=plot_comp)
+        cost, state, recorder = problem.optimize()
+
+main()

setup.py

@@ -29,7 +29,7 @@ setup(name='topfarm',
         'openmdao==2.6',  # for optimization
         'pytest',  # for testing
         'pytest-cov',  # for calculating coverage
-		'py_wake',  # for calculating AEP
+		'py_wake>=1.0.11',  # for calculating AEP
         'scipy',  # constraints
         'sphinx',  # generating documentation
         'sphinx_rtd_theme',  # docs theme

topfarm/cost_models/economic_models/dtu_wind_cm_main.py

@@ -15,14 +15,13 @@ class economic_evaluation():
         # Initialize dictionaries
         self.project_costs = {}
         self.project_costs_sums = {}
-
         self.discount_rate = discount_rate              # [-]
         self.distance_from_shore = distance_from_shore  # [km]
         self.energy_price = energy_price                # [Euro/kWh]
         self.project_duration = project_duration        # [years]
 
     def calculate_npv(self, rated_rpm_array, D_rotor_array, Power_rated_array,
-                      hub_height_array, water_depth_array, aep_array):
+                      hub_height_array, water_depth_array, aep_array, cabling_cost=None):
         '''
         Calculate Net Present Value [Euro]
 
@@ -41,7 +40,8 @@ class economic_evaluation():
             Power_rated_array,
             hub_height_array,
             aep_array,
-            water_depth_array)
+            water_depth_array,
+            cabling_cost)
 
         self.calculate_cash_flow()
 
@@ -50,7 +50,7 @@ class economic_evaluation():
         return self.NPV
 
     def calculate_irr(self, rated_rpm_array, D_rotor_array, Power_rated_array,
-                      hub_height_array, water_depth_array, aep_array):
+                      hub_height_array, water_depth_array, aep_array, cabling_cost=None):
         '''
         Calculate Internal Rate of Return [%]
 
@@ -69,7 +69,8 @@ class economic_evaluation():
             Power_rated_array,
             hub_height_array,
             aep_array,
-            water_depth_array)
+            water_depth_array,
+            cabling_cost)
 
         self.calculate_cash_flow()
 
@@ -97,7 +98,7 @@ class economic_evaluation():
                                     self.project_costs_sums["OPEX"]))
 
     def calculate_expenditures(self, rated_rpm, rotor_diameter, rated_power,
-                               hub_height, aep_array, water_depth):
+                               hub_height, aep_array, water_depth, cabling_cost=None):
         '''Calculate DEVEX, CAPEX, OPEX and ABEX [Euro]'''
 
         # Ensure that the variables are numpy arrays
@@ -107,7 +108,7 @@ class economic_evaluation():
         self.calculate_devex(rated_power)
 
         self.calculate_capex(rated_rpm, rotor_diameter, rated_power,
-                             hub_height, water_depth)
+                             hub_height, water_depth, cabling_cost)
 
         self.calculate_opex(rated_power)
 
@@ -132,7 +133,7 @@ class economic_evaluation():
             self.project_costs[myName].values())
 
     def calculate_capex(self, rated_rpm, rotor_diameter, rated_power,
-                        hub_height, water_depth):
+                        hub_height, water_depth, cabling_cost=None):
         '''Calculate CAPEX [Euro]'''
 
         rated_rpm = np.array(rated_rpm)
@@ -166,6 +167,9 @@ class economic_evaluation():
             50.0 *
             np.sum(rated_power) ** 2.0}
 
+        if cabling_cost:
+            self.project_costs[myName]["array_of_cables"] = cabling_cost
+
         self.project_costs_sums[myName] = sum(
             self.project_costs[myName].values())