1.3.0

下载 Notebook 下载完整案例

风电场布局优化问题 (WFLOP)¶

类别： 网络设计

使用 OptAgent 的 Python 接口描述变量、约束与目标。

问题与原始示例来源：Hexaly Code Templates。

问题描述¶

风电场布局优化问题 (WFLOP) 涉及在指定区域内确定固定数量的风力涡轮机的最佳布置位置。每台涡轮机必须位于该区域内部,并且与其他涡轮机之间保持最小安全距离。目标是最大化年发电量 (AEP),这是衡量风电场将风能转化为电能效率的关键指标。该优化问题的核心挑战在于最小化尾流损失,即涡轮机之间相互干扰其风流所产生的损失。

虽然许多方法处理该问题的连续版本,但我们的工作聚焦于一种离散近似,从而能够在各种可行域类型上进优化。在本示例中,我们关注圆形区域,但该方法同样适用于任何几何形状。

学到的建模技巧¶

• 使用 Python JSON 模块读取输入文件
• 使用 OptAgent 非线性运算符计算尾流损失和功率
• 最大化非线性年发电量目标 (AEP)

数据¶

我们使用 IEA Task 37 on System Engineering in Wind Energy 提供的实例。其中包含三种场景,分别涉及 16、36 和 64 台风力涡轮机,部署在半径分别为 1300 m、2000 m 和 3000 m 的圆形区域内。我们已将原始的 .yaml 数据文件转换为 .json 格式以便处理。

每个 .json 文件包含以下参数:

• "radius": 部署区域半径

• "D": 风轮直径

• "nT": 涡轮机数量

• "prated": 单台涡轮机的额定功率

• "urated": 额定风速

• "ucut-in": 切入风速

• "ucut-out": 切出风速

• "wind": 风数据,包括:

• "degrees": 风向角度的离散化

• "probs": 风向的概率分布

• "J": 考虑的风向总数

• "speed": 来流风速

• "cT": 推力系数

• "kY": 基于湍流强度 0.075 的动态系数

可行区域以间距 1.7 × D 进行均匀离散化,涡轮机之间的最小距离为 2 × D。涡轮机既可放置在圆形区域内的离散网格点上,也可放置在边界上,其中针对每个角度评估一个可行位置(详见下文)。

建模思路¶

用于风电场布局优化问题 (WFLOP) 的 OptAgent 模型使用一个布尔决策向量,其中每个元素表示对应离散位置上是否放置风力涡轮机(1 表示放置,0 表示不放置)。

第一个约束确保可行区域内放置的涡轮机总数等于 nT。

为强制最小距离要求,我们对任意一对距离小于允许阈值的位置施加约束,使得这两个位置中至多只能放置一台涡轮机。

每个潜在位置的风条件采用 Bastankhah 高斯尾流模型的简化版本 进行计算,涡轮机的功率输出采用一条平滑、凸且单调递增的简化功率曲线得到(参见 IEA Task 37 Anouncements p2)。

目标函数通过聚合所有风向下的期望功率输出(以各自概率加权)来最大化年发电量 (AEP)。

Python 实现¶

import json
import math
from pathlib import Path

from optagent import OptModel, solve


def main(instance_file, output_file=None, time_limit=60):
    data = read_data(instance_file)
    points_x, points_y = build_discretization(data["radius"], data["disc_delta"])
    nb_locations = len(points_x)
    incompatibility_sets = compute_incompatibilities(points_x, points_y, data["min_distance"])

    model = OptModel()
    chosen_turbines = [model.bool(name=f"location_{location}") for location in range(nb_locations)]
    model.constraint(model.sum(chosen_turbines) == data["nb_turbines"], name="turbine_count")

    for location in range(nb_locations):
        for incompatible_location in incompatibility_sets[location]:
            if incompatible_location <= location:
                continue
            model.constraint(
                chosen_turbines[location] + chosen_turbines[incompatible_location] <= 1,
                name=f"spacing_{location}_{incompatible_location}",
            )

    # Store sparse wake rows as offsets into shared index/coefficient arrays.
    # Reuse one power formula across all location/direction pairs.
    nb_degrees = len(data["degrees"])
    wake_offsets = [0]
    wake_indices = []
    wake_coefficients = []
    for location_1 in range(nb_locations):
        turbine_1 = (points_x[location_1], points_y[location_1])
        for angle in range(nb_degrees):
            for location_2 in range(nb_locations):
                loss = wake_loss(
                    turbine_1,
                    (points_x[location_2], points_y[location_2]),
                    data["degrees"][angle],
                    data["dyn_coeff"],
                    data["turbine_diameter"],
                    data["thrust_coeff"],
                )
                if loss != 0:
                    wake_indices.append(location_2)
                    wake_coefficients.append(loss**2)
            wake_offsets.append(len(wake_indices))

    turbine_choices = model.array(chosen_turbines)
    offsets = model.array(wake_offsets)
    # A sentinel keeps the element type explicit even when every row is empty.
    indices = model.array(wake_indices + [0])
    coefficients = model.array(wake_coefficients + [0.0])
    probabilities = model.array(data["probabilities"])
    wake_term = model.lambda_function(lambda i: turbine_choices[indices[i]] * coefficients[i])

    def weighted_power(row):
        wake_loss_sum = model.sum(model.range(offsets[row], offsets[row + 1]), wake_term)
        turbine_wind = data["inflow_speed"] * (turbine_choices[row // nb_degrees] - model.sqrt(wake_loss_sum))
        polynomial_case = (data["min_speed"] <= turbine_wind) * (turbine_wind < data["optimal_speed"])
        polynomial_value = ((turbine_wind - data["min_speed"]) / (data["optimal_speed"] - data["min_speed"])) ** 3
        constant_case = (data["optimal_speed"] <= turbine_wind) * (turbine_wind < data["max_speed"])
        power = (
            polynomial_case * data["turbine_nominal_power"] * polynomial_value
            + constant_case * data["turbine_nominal_power"]
        )
        # With nonnegative cut-in speed an absent turbine has zero power.
        # Keep the original formula for data outside that physical assumption.
        if data["optimal_speed"] > data["min_speed"] >= 0 and data["inflow_speed"] >= 0:
            power = model.iif(turbine_choices[row // nb_degrees], power, 0.0)
        return probabilities[row % nb_degrees] * power

    annual_energy = 8760 * model.sum(
        model.range(0, nb_locations * nb_degrees), model.lambda_function(weighted_power)
    )
    model.maximize(annual_energy, name="annual_energy_production")

    solution = solve(model, time_limit_s=float(time_limit))
    if not solution.feasible:
        print(f"No feasible solution found; Status = {solution.status}")
        return solution

    points = [
        (points_x[location], points_y[location])
        for location in range(nb_locations)
        if chosen_turbines[location].value == 1
    ]
    print(f"Annual Energy Production = {annual_energy.value} Wh; Status = {solution.status}")

    if output_file is not None:
        lines = [
            "*************************************",
            "*** WIND FARM LAYOUT OPTIMIZATION ***",
            "*************************************",
            "",
            (
                f"Annual Energy Production : {annual_energy.value} Wh "
                f"({round(annual_energy.value / 1e9 * 100) / 100} GWh)."
            ),
            "",
            "Locations of the Wind Turbines : ",
        ]
        lines.extend(f"   - ({point[0]}, {point[1]})" for point in points)
        Path(output_file).write_text("\n".join(lines) + "\n", encoding="utf-8")
        print("Solution written in file", output_file)
    return solution


def read_data(instance_file):
    """
    Reads the input files of the problem.
    """
    data = json.loads(Path(instance_file).read_text(encoding="utf-8"))

    # Wind information
    degrees = data["wind"]["degrees"]  # Discretization of the degrees rose
    probabilities = data["wind"]["probs"]  # Probabilities of wind in every direction
    inflow_speed = data["wind"]["speed"]  # Inflow wind speed
    thrust_coeff = data["wind"]["cT"]  # Thrust coefficient
    dyn_coeff = data["wind"]["kY"]  # Dynamic coefficient

    # Information about the problem
    radius = data["radius"]  # Radius of the perimeter to fullfil
    nb_turbines = data["nT"]  # Number of Wind Turbines
    turbine_diameter = data["D"]  # Diameter of each Turbine
    turbine_nominal_power = data["p_rated"] * 1e6  # Power of the Turbine
    min_distance = 2 * turbine_diameter  # Minimal distance between two turbines
    optimal_speed = data["u_rated"]  # Optimal wind speed (maximum efficiency)
    min_speed = data["u_cut_in"]  # Minimal wind speed producing energy
    max_speed = data["u_cut_out"]  # Maximal wind speed producing energy

    disc_delta = 1.7 * turbine_diameter  # Distance between two points in the discretization

    return {
        "degrees": degrees,
        "probabilities": probabilities,
        "inflow_speed": inflow_speed,
        "thrust_coeff": thrust_coeff,
        "dyn_coeff": dyn_coeff,
        "radius": radius,
        "nb_turbines": nb_turbines,
        "turbine_diameter": turbine_diameter,
        "turbine_nominal_power": turbine_nominal_power,
        "min_distance": min_distance,
        "optimal_speed": optimal_speed,
        "min_speed": min_speed,
        "max_speed": max_speed,
        "disc_delta": disc_delta,
    }


def build_discretization(radius, disc_delta):
    """
    Builds the discretization in a circular field of radius {radius} and with a
    distance between points of {disc_delta}.

    The discretization follows the one described in 'eawe' paper, i.e. building
    a regular uniform mesh inside the circle, and a point every degree on the
    frontier of the circle.
    """
    points_x = []
    points_y = []
    # Maximum number of points in any direction from the center
    max_onedir_points = int(radius / disc_delta + 1)
    # Interior points
    for c_x in range(max_onedir_points):
        for side_x in [-1, 1]:
            if not (c_x == 0 and side_x == 1):
                new_x = side_x * c_x * disc_delta
                for c_y in range(max_onedir_points):
                    for side_y in [-1, 1]:
                        new_y = side_y * c_y * disc_delta
                        if math.sqrt(new_x**2 + new_y**2) < radius and not (c_y == 0 and side_y == 1):
                            points_x.append(new_x)
                            points_y.append(new_y)

    # Points on the border of the circle
    for deg in range(360):
        points_x.append(radius * math.sin(deg * math.pi / 180))
        points_y.append(radius * math.cos(deg * math.pi / 180))

    return points_x, points_y


def compute_incompatibilities(points_x, points_y, min_distance):
    """
    Compute incompatibility sets for every location on the field.

    N_i[i] := {l in 0...nb_locations : l != i && ||l - i|| < min_distance}
    """
    # Number of available points
    nb_locations = len(points_x)
    incompatibilities = [[] for location in range(nb_locations)]
    # Fill the set for every location
    for location_1 in range(nb_locations):
        for location_2 in range(location_1 + 1, nb_locations):
            if (
                math.sqrt(
                    (points_x[location_1] - points_x[location_2]) ** 2
                    + (points_y[location_1] - points_y[location_2]) ** 2
                )
                < min_distance
            ):
                incompatibilities[location_1].append(location_2)
                incompatibilities[location_2].append(location_1)

    return incompatibilities


def distances(location_1, location_2, theta):
    """
    This function calculates the distance between two points, considering a
    frame of reference with angle theta.
    """
    # Rotation of the reference
    theta_deg = 270 - theta
    theta_rad = theta_deg * math.pi / 180
    cos_wind = math.cos(-theta_rad)
    sin_wind = math.sin(-theta_rad)

    # Change the reference of the coordinates
    location_2_x = (location_2[0] * cos_wind) - (location_2[1] * sin_wind)
    location_2_y = (location_2[0] * sin_wind) + (location_2[1] * cos_wind)

    location_1_x = (location_1[0] * cos_wind) - (location_1[1] * sin_wind)
    location_1_y = (location_1[0] * sin_wind) + (location_1[1] * cos_wind)

    return {"parallel": location_1_x - location_2_x, "perpendicular": location_1_y - location_2_y}


def wake_loss(turbine_1, turbine_2, theta, dyn_coeff, turbine_diameter, thrust_coeff):
    """
    Wake loss evaluated at i, caused by a Wind Turbine at l, considering a wind
    of angle {theta}.
    This loss is a simplified version of Bastankhah's Gaussian model.
    """
    # We calculate parallel and perpendicular distances
    d = distances(turbine_1, turbine_2, theta)
    # If x_i - x_l < 0, the wake loss is zero
    if d["parallel"] > 0:
        # Standard deviation of the wake deficit
        s_y = dyn_coeff * d["parallel"] + turbine_diameter / math.sqrt(8)
        # Separation of the product in two components
        coeff = thrust_coeff / (8 * (s_y / turbine_diameter) ** 2)
        pdt_1 = 1 - math.sqrt(1 - coeff)
        pdt_2 = math.exp(-1 / 2 * (d["perpendicular"] / s_y) ** 2)

        return pdt_1 * pdt_2

    return 0

INSTANCE_DIR = Path.cwd() / "instances"
print("Instances:", INSTANCE_DIR)

solution_nt02 = main(INSTANCE_DIR / "nT02.json", time_limit=10)