Cost Models

Try this yourself (requires google account)

TopFarm can plug into the DTU’s costmodels package, containing a collection of economical cost models of wind turbines and wind farms. The package provides a convenient interface to compute common economical metrics which can be included as a objective function in the optimization. Moreover, additional user-defined cost models can easily be integrated in the TOPFARM optimization problems as well. Please see the documentation of the package in the repository.

# 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 numpy as np
from costmodels.finance import Product, Technology
from costmodels.models import DTUOffshoreCostModel
from costmodels.models import BatteryCostModel
from costmodels.project import Project
from py_wake.deficit_models.gaussian import IEA37SimpleBastankhahGaussian
from py_wake.examples.data.iea37._iea37 import IEA37_WindTurbines, IEA37Site

from topfarm import TopFarmGroup, TopFarmProblem
from topfarm.constraint_components.boundary import CircleBoundaryConstraint
from topfarm.constraint_components.spacing import SpacingConstraint
from topfarm.cost_models.cost_model_wrappers import CostModelComponent
from topfarm.drivers.random_search_driver import RandomizeTurbinePosition_Circle
from topfarm.cost_models.economic_models.finance_wrapper import (
    OffshoreWindPlantFinanceWrapper,
)
from topfarm.easy_drivers import (
    EasyRandomSearchDriver,
    EasyScipyOptimizeDriver,
    EasySimpleGADriver,
    ScipyOptimizeDriver,
)
from topfarm.plotting import NoPlot, XYPlotComp

Set up plotting capability

try:
    import matplotlib.pyplot as plt
    plt.gcf()
    plot_comp = XYPlotComp()
    plot = True
except RuntimeError:
    plot_comp = NoPlot()
    plot = False

<Figure size 640x480 with 0 Axes>

Set up IEA Wind Task 37 case study site with 16 turbines.

n_wt = 16 # number of wind turbines
site = IEA37Site(n_wt) # site is the IEA Wind Task 37 site with a circle boundary
windTurbines = IEA37_WindTurbines() # wind turbines are the IEA Wind Task 37 3.4 MW reference turbine
wake_model = IEA37SimpleBastankhahGaussian(site, windTurbines) # select the Gaussian wake model

# vectors for turbine properties: diameter, rated power and hub height. these are inputs to the cost model
Drotor_vector = [windTurbines.diameter()] * n_wt
power_rated_vector = [float(windTurbines.power(20)/1000)] * n_wt
hub_height_vector = [windTurbines.hub_height()] * n_wt

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

Setup cost model with DTU Offshore cost model case. Starting with bathymetry of the site

# Bathymetry grid
from topfarm.examples.bathymetry_ficticio import (
    get_bathymetry_func_circle,
    plot_bathymetry_circle,
)

RADIUS = 1300.1
f = get_bathymetry_func_circle(
    sigma=RADIUS, mu=0.0, x_peak_1=0, y_peak_1=1000, x_peak_2=3000, y_peak_2=300
)
plot_bathymetry_circle(f, RADIUS)


# Water depth function
def water_depth_func(x, y, **kwargs):
    xnew, ynew = np.meshgrid(x, y)
    points = np.array([xnew.flatten(), ynew.flatten()]).T
    return -np.diag(f(points).reshape(n_wt, n_wt).T)

Setup economic model based on water depth

For single technology offshore wind plants you can use the pre-confiugred economic component. Further down you can explore using the cost model package directly with the possibility to change the techno-economical models and add or change the technologies present in the plant by adding e.g. solar or batteries and change or add energy vectors like e.g. hydrogen.

# Economic
LIFETIME = 25
el_price = np.random.uniform(0, 6, LIFETIME)  # time series
el_price = 6  # fixed ppa price
npv_comp = OffshoreWindPlantFinanceWrapper(windTurbines, n_wt, el_price, LIFETIME)
cost_model = npv_comp.cost_model
out = cost_model.run(aep=373206.64521435613, water_depth=20.0)
print(out)

CostOutput(capex=Array(103.73955334, dtype=float64), opex=Array(0.00137346, dtype=float64))

Set up functions for the AEP and cost calculations.

# function for calculating aep as a function of x,y positions of the wind turbiens
def aep_func(x, y, **kwargs):
    return wake_model(x, y).aep().sum().values * 1e3

Now set up a problem to run an optimization using NPV as the objective function.

Note that the turbines are fixed so the main driver changing the NPV will be the AEP as the turbine positions change. Here you can select different drivers to see how the optimization result changes.

# create an openmdao component for aep and NPV to add to the problem
aep_comp = CostModelComponent(
    input_keys=["x", "y"],
    n_wt=n_wt,
    cost_function=aep_func,
    output_keys="aep",
    output_unit="MWh",
    objective=False,
    output_vals=0.0,
)

# Water Depth
water_depth_component = CostModelComponent(
    input_keys=[("x", site.initial_position.T[0]), ("y", site.initial_position.T[1])],
    n_wt=n_wt,
    cost_function=water_depth_func,
    objective=False,
    output_keys=[("water_depth", np.zeros(n_wt))],
)

# create a group for the aep and npv components that links their common input/output (aep)
npv_group = TopFarmGroup([aep_comp, water_depth_component, npv_comp])

# add the group to an optimization problem and specify the design variables (turbine positions),
# cost function (npv_group), and constraints (circular boundary and spacing)
problem = TopFarmProblem(
    design_vars=dict(zip("xy", site.initial_position.T)),
    n_wt=n_wt,
    cost_comp=npv_group,
    # specify driver to use: SLSQP (gradient-based), random search (gradient-free), COBYLA (gradient-free), genetic algorithm, GA (gradient-free)
    # driver=ScipyOptimizeDriver(optimizer="SLSQP", maxiter=30, tol=1e-6, disp=True),
    driver=EasyRandomSearchDriver(
        randomize_func=RandomizeTurbinePosition_Circle(max_step=50), max_iter=50
    ),
    # driver=EasyScipyOptimizeDriver(optimizer='COBYLA', maxiter=200, tol=1e-6, disp=False),
    # driver=EasySimpleGADriver(max_gen=100, pop_size=5, Pm=None, Pc=.5, elitism=True, bits={}),
    constraints=[SpacingConstraint(200), CircleBoundaryConstraint([0, 0], RADIUS)],
    plot_comp=plot_comp,
)

# assign data from optimizationn to a set of accessible variables and run the optimization
cost, state, recorder = problem.optimize()

from topfarm.utils import plot_list_recorder
plot_list_recorder(recorder)

Detailed Economic model

If you have a hybrid plant comprised of several technologies and/or energy vectors you can make a bespoke financial component

LIFETIME = 25
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=0.4,
    nwt=n_wt,
)

wind_technology = Technology(
    name="wind",
    lifetime=LIFETIME,
    product=Product.SPOT_ELECTRICITY,
    cost_model=cost_model,
)

battery_cost_model = BatteryCostModel(
    battery_power=27.0,
    battery_energy=108.0,
    state_of_health=np.hstack(
        [
            -1.7e-6 * np.arange(1.8e5) + 1,
            -2.5e-6 * np.arange(25 * 365 * 24 - 1.8e5) + 1,
        ]
    ).ravel(),
)
battery_technology = Technology(
    name="battery",
    lifetime=LIFETIME,
    cost_model=battery_cost_model,
)

wind_farm_project = Project(
    technologies=[wind_technology, battery_technology],
    product_prices={Product.SPOT_ELECTRICITY: np.random.uniform(0, 6, LIFETIME)},
)


def npv_func(aep, water_depth, **kwargs):
    return np.asarray(
        wind_farm_project.npv(
            productions={
                wind_technology.name: aep,
            },
            cost_model_args={
                wind_technology.name: {
                    "water_depth": water_depth,
                    "aep": aep,
                },
            },
        )
    )


def npv_grad_func(aep, water_depth, **kwargs):
    grads = wind_farm_project.npv_grad(
        productions={
            wind_technology.name: aep,
        },
        cost_model_args={
            wind_technology.name: {
                "water_depth": water_depth,
                "aep": aep,
            },
        },
    )
    prod_grad = grads[0][wind_technology.name]
    water_depth_grad = grads[1][wind_technology.name]["water_depth"]
    return np.asarray(prod_grad), np.asarray(water_depth_grad)

# # Economy
npv_comp = CostModelComponent(
    input_keys=[
        ("aep", 0),
        ("water_depth", 30 * np.ones(n_wt)),
    ],
    n_wt=n_wt,
    cost_function=npv_func,
    cost_gradient_function=npv_grad_func,
    objective=True,
    maximize=True,
    output_keys=[("npv", 0)],
    output_unit="EUR",
)
npv_group = TopFarmGroup([aep_comp, water_depth_component, npv_comp])

# add the group to an optimization problem and specify the design variables (turbine positions),
# cost function (npv_group), and constraints (circular boundary and spacing)
problem = TopFarmProblem(
    design_vars=dict(zip("xy", site.initial_position.T)),
    n_wt=n_wt,
    cost_comp=npv_group,
    # specify driver to use: SLSQP (gradient-based), random search (gradient-free), COBYLA (gradient-free), genetic algorithm, GA (gradient-free)
    driver=ScipyOptimizeDriver(optimizer="SLSQP", maxiter=30, tol=1e-6, disp=True),
    # driver=EasyRandomSearchDriver(
    #     randomize_func=RandomizeTurbinePosition_Circle(max_step=50), max_iter=50
    # ),
    # driver=EasyScipyOptimizeDriver(optimizer='COBYLA', maxiter=200, tol=1e-6, disp=False),
    # driver=EasySimpleGADriver(max_gen=100, pop_size=5, Pm=None, Pc=.5, elitism=True, bits={}),
    constraints=[SpacingConstraint(200), CircleBoundaryConstraint([0, 0], RADIUS)],
    plot_comp=plot_comp,
)

# assign data from optimizationn to a set of accessible variables and run the optimization
cost, state, recorder = problem.optimize()

E0211 11:03:43.561468    5364 slow_operation_alarm.cc:73] Constant folding an instruction is taking > 1s:

  %reduce-window.18 = s64[13688,16]{0,1} reduce-window(%constant.1653, %constant.346), window={size=1x16 pad=0_0x15_0}, to_apply=%region_3.7

This isn't necessarily a bug; constant-folding is inherently a trade-off between compilation time and speed at runtime. XLA has some guards that attempt to keep constant folding from taking too long, but fundamentally you'll always be able to come up with an input program that takes a long time.

If you'd like to file a bug, run with envvar XLA_FLAGS=--xla_dump_to=/tmp/foo and attach the results.
E0211 11:03:43.884196    5279 slow_operation_alarm.cc:140] The operation took 1.327804216s
Constant folding an instruction is taking > 1s:

  %reduce-window.18 = s64[13688,16]{0,1} reduce-window(%constant.1653, %constant.346), window={size=1x16 pad=0_0x15_0}, to_apply=%region_3.7

This isn't necessarily a bug; constant-folding is inherently a trade-off between compilation time and speed at runtime. XLA has some guards that attempt to keep constant folding from taking too long, but fundamentally you'll always be able to come up with an input program that takes a long time.

If you'd like to file a bug, run with envvar XLA_FLAGS=--xla_dump_to=/tmp/foo and attach the results.

Iteration limit reached    (Exit mode 9)
            Current function value: 123810488.5671965
            Iterations: 30
            Function evaluations: 59
            Gradient evaluations: 31
Optimization FAILED.
Iteration limit reached
-----------------------------------

plot_list_recorder(recorder)

Exercise

Interpret results and see if you can improve the result by giving optimizer more iteration or change the driver in the topfarm problem above to see if an improved objective function can be obtained.

Exercise

Manipulate the DTU Offshore Cost Model instantiation inputs to see how this influences the optimal NPV found.