Scipy风电场AEP优化场景下添加湍流强度约束的实现方法咨询

Scipy的minimize函数支持传入自定义不等式约束，我们只需额外定义约束校验函数，同时增加仿真结果缓存避免重复计算，不需要改动原有AEP优化的核心逻辑，修改步骤如下：

• 新增缓存逻辑，存储相同入参c对应的仿真结果，避免目标函数和约束函数重复调用仿真接口浪费性能
• 定义约束函数：要求所有被选中运行的风机的全部风速、风向下的有效湍流强度TI_eff ≤ 阈值0.2，转换为Scipy要求的约束返回值 ≥ 0格式，返回0.2 - 各TI_eff值即可
• 调用minimize时指定支持约束的求解器（如SLSQP），同时传入约束配置

"""
References:
    scipy.optimize.minimize官方文档
    PyWake项目
"""

import time
import numpy as np
from scipy.optimize import minimize
from py_wake.examples.data.hornsrev1 import V80 
from py_wake.examples.data.hornsrev1 import Hornsrev1Site
from py_wake import BastankhahGaussian

# 全局缓存，避免重复仿真
_last_c = None
_last_sim_res = None
_last_neg_aep = None
TI_THRESHOLD = 0.2  # 湍流强度阈值

def funC(x, y, c):
    """
    根据c值筛选要运行的风机，c≥0.5的风机投入运行
    """
    mask = c >= 0.5
    x_selected = x[mask]
    y_selected = y[mask]
    return x_selected, y_selected, mask

def run_simulation(c):
    """统一仿真入口，带缓存逻辑"""
    global _last_c, _last_sim_res, _last_neg_aep
    # 入参未变化直接返回缓存结果
    if _last_c is not None and np.array_equal(c, _last_c):
        return _last_neg_aep, _last_sim_res
    # 入参变化重新仿真
    site = Hornsrev1Site()
    x, y = site.initial_position.T
    windTurbines = V80()
    wf_model = BastankhahGaussian(site, windTurbines)
    x_new, y_new, _ = funC(x, y, c)
    # 运行仿真
    sim_res = wf_model(
        x_new, y_new,
        h=None,
        type=0,
        wd=None,
        ws=None
    )
    aep_output = sim_res.aep().sum()
    neg_aep = -float(aep_output)
    # 更新缓存
    _last_c = c.copy()
    _last_sim_res = sim_res
    _last_neg_aep = neg_aep
    return neg_aep, sim_res

def wt_simulation(c):
    """目标函数，仅返回负AEP供优化器最小化"""
    neg_aep, _ = run_simulation(c)
    return neg_aep

def ti_constraint(c):
    """湍流强度约束函数，返回值要求全部≥0"""
    _, sim_res = run_simulation(c)
    # 取出所有被选中风机的全部wd、ws下的TI_eff值，展平为一维数组
    all_ti = sim_res.TI_eff.values.flatten()
    # 约束要求TI ≤ TI_THRESHOLD，即 TI_THRESHOLD - TI ≥ 0
    return TI_THRESHOLD - all_ti

def solve():
    t0 = time.perf_counter()
    wt = 80  # V80风机总数
    x0 = np.ones(wt)  # 初始值默认全部风机投入运行
    bounds = [(0, 1) for _ in range(wt)]
    # 定义约束配置，类型为不等式约束
    constraints = [{'type': 'ineq', 'fun': ti_constraint}]
    # 调用minimize，指定支持约束的SLSQP求解器
    res = minimize(
        wt_simulation,
        x0=x0,
        bounds=bounds,
        constraints=constraints,
        method='SLSQP'
    )
    print(f'success status: {res.success}')
    print(f'aep: {-res.fun} GWh')  # 取反得到真实最大AEP
    print(f'c values: {res.x}\n')
    print(f'elapse: {round(time.perf_counter() - t0)}s')

# 启动优化
solve()

注意：如果需要排除非运行风机的TI校验，可在约束函数中通过funC返回的选中掩码过滤，上述代码默认仅校验实际投入运行的风机的TI值，符合业务逻辑。

内容的提问来源于stack exchange，提问作者sadra