Cables

Try this yourself (requires google account)

TOPFARM can use the Electrical Network Design package EDWIN to optimize the carray cabels as well as the substation position at each iteration of the layout optimization

Install TOPFARM if needed

# Install TopFarm if needed
import importlib
if not importlib.util.find_spec("topfarm"):
    %pip install git+https://gitlab.windenergy.dtu.dk/TOPFARM/TopFarm2.git

Import

import numpy as np

from topfarm.constraint_components.spacing import SpacingConstraint
from topfarm.constraint_components.boundary import XYBoundaryConstraint
from topfarm.easy_drivers import EasyScipyOptimizeDriver
from topfarm.cost_models.cost_model_wrappers import CostModelComponent
from topfarm._topfarm import TopFarmProblem, TopFarmGroup
from topfarm.cost_models.py_wake_wrapper import PyWakeAEPCostModelComponent
from topfarm.plotting import XYPlotComp, NoPlot
from topfarm.utils import plot_list_recorder
from topfarm.cost_models.electrical.optiwindnet_wrapper import WFNComponent

from py_wake.examples.data.iea37 import IEA37_WindTurbines
from py_wake import BastankhahGaussian
from py_wake.examples.data.hornsrev1 import Hornsrev1Site


from costmodels.finance import Depreciation, Technology, Product
from costmodels.project import Project
from costmodels.models import DTUOffshoreCostModel
from topfarm.examples.bathymetry_ficticio import (
    gaussian_surface,
    get_bathymetry_func_rect,
    plot_bathymetry_rect,
)

Inputs

plot = True
n_wt = 30
x_min = 0
x_max = 6000
y_min = -10000
y_max = 0
maxiter = 100
sigma = 3000.0
# Here you set up cables [<number of turbines can be connected>, <price in € per meter>]
cables = np.array([(2, 2000), (5, 2200)])
LIFETIME = 25  # years
el_price = 50  # fixed ppa price Euro per MWh
driver = EasyScipyOptimizeDriver(maxiter=maxiter)

Geomertry

rng1 = np.random.default_rng(1)
rng2 = np.random.default_rng(2)
initial = np.asarray([rng1.random(n_wt)*(x_max-x_min)+x_min, rng2.random(n_wt)*(y_max-y_min)+y_min]).T
x_init = initial[:,0]
y_init = initial[:,1]
x_ss_init = x_init.mean()
y_ss_init = y_init.mean()
turbines_pos=np.asarray([x_init, y_init]).T
substations_pos = np.asarray([[x_ss_init], [y_ss_init]]).T
boundary = np.array([(x_min, y_min), (x_max, y_min), (x_max, y_max), (x_min, y_max)])  # turbine boundaries

Site

windTurbines = IEA37_WindTurbines()
site = Hornsrev1Site()
wfm = BastankhahGaussian(site, windTurbines)

/builds/TOPFARM/TopFarm2/.pixi/envs/default/lib/python3.12/site-packages/py_wake/deficit_models/gaussian.py:124: UserWarning: The BastankhahGaussian model is not representative of the setup used in the literature. For this, use py_wake.literature.gaussian_models.Bastankhah_PorteAgel_2014 instead
  DeprecatedModel.__init__(self, 'py_wake.literature.gaussian_models.Bastankhah_PorteAgel_2014')

Bathymetry

g = gaussian_surface(sigma, x_min=x_min, x_max=x_max, y_min=y_min, y_max=y_max)
if plot:
    plot_bathymetry_rect(g, x_min, x_max, y_min, y_max)
bathymetry_interpolator = get_bathymetry_func_rect(g, x_min, x_max, y_min, y_max)
def water_depth_func(x, y, **kwargs):
    xnew, ynew = np.meshgrid(x, y)
    points = np.array([xnew.flatten(), ynew.flatten()]).T
    return - np.diag(bathymetry_interpolator(points).reshape(n_wt, n_wt).T)

Economy

Drotor_vector = [windTurbines.diameter()] * n_wt
power_rated_vector = [float(windTurbines.power(20)) * 1e-6] * n_wt
hub_height_vector = [windTurbines.hub_height()] * n_wt
simres = wfm(x_init, y_init)
aep_ref = simres.aep().values.sum()
RP_MW = windTurbines.power(20) * 1e-6
CF_ref = aep_ref * 1e3 / (RP_MW * 24 * 365 * n_wt)

cost_model = DTUOffshoreCostModel(
    rated_power=windTurbines.power(20) / 1e6,
    rotor_speed=10.0,
    rotor_diameter=windTurbines.diameter(),
    hub_height=windTurbines.hub_height(),
    lifetime=LIFETIME,
    capacity_factor=CF_ref,
    nwt=n_wt,
    profit=0,
)

wind_plant = Technology(
    name="wind",
    lifetime=LIFETIME,
    product=Product.SPOT_ELECTRICITY,
    opex=12600 * n_wt * RP_MW + 1.35 * aep_ref * 1000,  # Euro
    wacc=0.06,
    cost_model=cost_model,
)

project = Project(
    technologies=[wind_plant],
    product_prices={Product.SPOT_ELECTRICITY: el_price},
    depreciation=Depreciation(rate=(0, 1), year=(0, LIFETIME)),
)


def economic_func(aep, water_depth, cabling_cost, **kwargs):
    aep_scaled = aep * 10**3
    npv, aux = project.npv(
        productions={wind_plant.name: aep_scaled},
        cost_model_args={
            wind_plant.name: {"water_depth": water_depth, "aep": aep_scaled}
        },
        finance_args={"shared_capex": cabling_cost},
        return_aux=True,
    )
    return npv, {
        "LCOE": aux["LCOE"][0],
        "IRR": aux["IRR"],
        "CAPEX": aux["CAPEX"],
        "OPEX": np.mean(aux["OPEX"]),
    }


def economic_func_grad(aep, water_depth, cabling_cost, **kwargs):
    aep_scaled = aep * 10**3
    grad = project.npv_grad(
        productions={wind_plant.name: aep_scaled},
        cost_model_args={
            wind_plant.name: {"water_depth": water_depth, "aep": aep_scaled}
        },
        finance_args={"shared_capex": cabling_cost},
    )
    return (
        grad[0][wind_plant.name] * 1e3,  # dNPV/dAEP | because of AEP scaling at the top
        grad[1][wind_plant.name]["water_depth"],  # dNPV/dWaterDepth
        grad[2]["shared_capex"],  # dNPV/dCablingCost
    )

Components

aep_comp = PyWakeAEPCostModelComponent(wfm, n_wt, objective=False, output_key="aep")

water_depth_component = CostModelComponent(
    input_keys=[("x", x_init), ("y", y_init)],
    n_wt=n_wt,
    cost_function=water_depth_func,
    objective=False,
    output_keys=[("water_depth", np.zeros(n_wt))],
)

cable_component = WFNComponent(turbines_pos, substations_pos, cables)

npv_comp = CostModelComponent(
    input_keys=[("aep", 0), ("water_depth", np.zeros(n_wt)), ("cabling_cost", 0)],
    n_wt=n_wt,
    cost_function=economic_func,
    cost_gradient_function=economic_func_grad,
    output_keys=[("NPV", 0)],
    additional_output=[
        ("LCOE", 0),
        ("IRR", 0),
        ("CAPEX", 0),
        ("OPEX", 0),
        # ('LCOE_ref', 0)
    ],
    maximize=True,
    objective=True,
)

cost_comp = TopFarmGroup([aep_comp, water_depth_component, cable_component, npv_comp])

Problem Assembly

if plot:
    plot_comp = XYPlotComp()
else:
    plot_comp = NoPlot()

tf = TopFarmProblem(
    design_vars=dict(zip("xy", initial.T), xs=x_ss_init, ys=y_ss_init),
    cost_comp=cost_comp,
    constraints=[XYBoundaryConstraint(boundary), SpacingConstraint(500)],
    driver=driver,
    plot_comp=plot_comp,
    expected_cost=1e1,
)

Smart Start

x = np.linspace(500, 5500, 41)
y = np.linspace(-9500, -500, 41)
YY, XX = np.meshgrid(y, x)
tf.smart_start(
    XX, YY, tf.cost_comp.comp_0.get_aep4smart_start(ws=8, wd=270), random_pct=20
)
tf.evaluate()

Smartstart: 100%|██████████| 30/30 [00:04<00:00,  6.13it/s]

1681 possible points, 30 wt, 56.0 points pr wt, 1057(63%) unused points

(np.float64(18400091.896197386),
 {'x': array([4500., 3000.,  750., 4750., 1250., 3375., 2500.,  875., 5500.,
         1125.,  625., 1250., 1750.,  750.,  875., 3500., 2875., 1500.,
         3750.,  625., 1625., 1000., 1500., 2500., 1125., 3875., 1125.,
         2250., 1500., 2500.]),
  'y': array([-5675., -2300., -8150., -4775., -6350., -1850., -1175., -1625.,
         -3650., -5000., -5450., -9050., -7025., -4550., -5900., -3425.,
         -3425., -7475.,  -500., -2300.,  -950., -7250., -9500., -2975.,
         -3200., -3875.,  -950., -3650., -3875., -6800.]),
  'xs': array([3081.11677185]),
  'ys': array([-5262.48585166])})

Optimize

cost, state, recorder = tf.optimize()

Iteration limit reached    (Exit mode 9)
            Current function value: -2325385.2204741337
            Iterations: 100
            Function evaluations: 140
            Gradient evaluations: 101
Optimization FAILED.
Iteration limit reached
-----------------------------------

Plot

plot_list_recorder(recorder, dont_plot=["v0", "v1"])
tf.model.cost_comp.comp_2.wfn.plot()

<Axes: >