multi_energy_complementarity/src/storage_optimization.py

"""
多能互补系统储能容量优化计算程序

该程序计算多能互补系统中所需的储能容量，确保系统在24小时内电能平衡，
同时满足用户定义的弃风弃光率和上网电量比例约束。

作者: iFlow CLI
创建日期: 2025-12-25
"""

import numpy as np
from typing import List, Dict, Tuple, Optional
from dataclasses import dataclass


@dataclass
class SystemParameters:
    """系统参数配置类"""
    max_curtailment_wind: float = 0.1    # 最大允许弃风率 (0.0-1.0)
    max_curtailment_solar: float = 0.1   # 最大允许弃光率 (0.0-1.0)
    max_grid_ratio: float = 0.2          # 最大允许上网电量比例 (0.0-∞，只限制上网电量，不限制购电)
    storage_efficiency: float = 0.9      # 储能充放电效率 (0.0-1.0)
    discharge_rate: float = 1.0          # 储能放电倍率 (C-rate)
    charge_rate: float = 1.0             # 储能充电倍率 (C-rate)
    max_storage_capacity: Optional[float] = None  # 储能容量上限 (MWh)，None表示无限制
    # 新增额定装机容量参数
    rated_thermal_capacity: float = 100.0  # 额定火电装机容量 (MW)
    rated_solar_capacity: float = 100.0    # 额定光伏装机容量 (MW)
    rated_wind_capacity: float = 100.0     # 额定风电装机容量 (MW)
    # 新增可用发电量参数
    available_thermal_energy: float = 2400.0  # 火电可用发电量 (MWh)
    available_solar_energy: float = 600.0    # 光伏可用发电量 (MWh)
    available_wind_energy: float = 1200.0    # 风电可用发电量 (MWh)


def validate_inputs(
    solar_output: List[float],
    wind_output: List[float],
    thermal_output: List[float],
    load_demand: List[float],
    params: SystemParameters
) -> None:
    """
    验证输入数据的有效性
    
    Args:
        solar_output: 24小时光伏出力曲线 (MW)
        wind_output: 24小时风电出力曲线 (MW)
        thermal_output: 24小时火电出力曲线 (MW)
        load_demand: 24小时负荷曲线 (MW)
        params: 系统参数配置
        
    Raises:
        ValueError: 当输入数据无效时抛出异常
    """
    # 检查数据长度（支持24小时或8760小时）
    data_length = len(solar_output)
    valid_lengths = [24, 8760]
    
    if data_length not in valid_lengths:
        raise ValueError(f"输入数据长度必须为24小时或8760小时，当前长度为{data_length}")
    
    if len(wind_output) != data_length or len(thermal_output) != data_length or len(load_demand) != data_length:
        raise ValueError("所有输入数据长度必须一致")
    
    # 检查数据类型和范围
    for name, data in [
        ("光伏出力", solar_output), ("风电出力", wind_output),
        ("火电出力", thermal_output), ("负荷需求", load_demand)
    ]:
        if not all(isinstance(x, (int, float)) for x in data):
            raise ValueError(f"{name}必须包含数值数据")
        if any(x < 0 for x in data):
            raise ValueError(f"{name}不能包含负值")
    
    # 检查参数范围
    if not (0.0 <= params.max_curtailment_wind <= 1.0):
        raise ValueError("弃风率必须在0.0-1.0之间")
    if not (0.0 <= params.max_curtailment_solar <= 1.0):
        raise ValueError("弃光率必须在0.0-1.0之间")
    # max_grid_ratio只限制上网电量比例，必须为非负值
    if not (0.0 <= params.max_grid_ratio):
        raise ValueError("上网电量比例限制必须为非负值")
    if not (0.0 < params.storage_efficiency <= 1.0):
        raise ValueError("储能效率必须在0.0-1.0之间")
    if params.discharge_rate <= 0 or params.charge_rate <= 0:
        raise ValueError("充放电倍率必须大于0")
    if params.max_storage_capacity is not None and params.max_storage_capacity <= 0:
        raise ValueError("储能容量上限必须大于0")
    # 验证新增的额定装机容量参数
    if params.rated_thermal_capacity < 0:
        raise ValueError("额定火电装机容量必须为非负值")
    if params.rated_solar_capacity <= 0:
        raise ValueError("额定光伏装机容量必须大于0")
    if params.rated_wind_capacity <= 0:
        raise ValueError("额定风电装机容量必须大于0")
    # 验证新增的可用发电量参数
    if params.available_thermal_energy < 0:
        raise ValueError("火电可用发电量必须为非负值")
    if params.available_solar_energy < 0:
        raise ValueError("光伏可用发电量必须为非负值")
    if params.available_wind_energy < 0:
        raise ValueError("风电可用发电量必须为非负值")


def calculate_energy_balance(
    solar_output: List[float],
    wind_output: List[float],
    thermal_output: List[float],
    load_demand: List[float],
    params: SystemParameters,
    storage_capacity: float,
    initial_soc: float = 0.0
) -> Dict[str, List[float]]:
    """
    计算给定储能容量下的系统电能平衡

    Args:
        solar_output: 光伏出力曲线 (MW) - 支持24小时或8760小时
        wind_output: 风电出力曲线 (MW) - 支持24小时或8760小时
        thermal_output: 火电出力曲线 (MW) - 支持24小时或8760小时
        load_demand: 负荷曲线 (MW) - 支持24小时或8760小时
        params: 系统参数配置
        storage_capacity: 储能容量 (MWh)
        initial_soc: 初始储能状态 (MWh)，默认为0.0

    Returns:
        包含各种功率曲线的字典
    """
    # 转换为numpy数组便于计算
    solar = np.array(solar_output)
    wind = np.array(wind_output)
    thermal = np.array(thermal_output)
    load = np.array(load_demand)

    # 初始化输出数组
    hours = len(solar_output)
    storage_soc = np.zeros(hours)          # 储能状态 (MWh)
    charge_power = np.zeros(hours)         # 充电功率 (MW)
    discharge_power = np.zeros(hours)      # 放电功率 (MW)
    curtailed_wind = np.zeros(hours)       # 弃风量 (MW)
    curtailed_solar = np.zeros(hours)      # 弃光量 (MW)
    grid_feed_in = np.zeros(hours)         # 上网电量 (MW)

    # 设置初始储能状态
    storage_soc[0] = initial_soc

    # 计算总发电潜力
    total_potential_wind = np.sum(wind)
    total_potential_solar = np.sum(solar)

    # 判断是否只有一种可再生能源
    has_wind = total_potential_wind > 0
    has_solar = total_potential_solar > 0
    single_renewable = (has_wind and not has_solar) or (has_solar and not has_wind)

    # 计算允许的最大弃风弃光量
    if single_renewable:
        # 只有一种可再生能源时，弃电量不受限制
        max_curtailed_wind_total = float('inf')
        max_curtailed_solar_total = float('inf')
    elif params.max_grid_ratio == 0:
        # 上网电量限制为0时，所有超额电力都必须被弃掉，不受弃风弃光限制
        max_curtailed_wind_total = float('inf')
        max_curtailed_solar_total = float('inf')
    else:
        # 有多种可再生能源且上网电量限制不为0时，应用弃风弃光限制
        max_curtailed_wind_total = total_potential_wind * params.max_curtailment_wind
        max_curtailed_solar_total = total_potential_solar * params.max_curtailment_solar

    # 初始化累计弃风弃光量
    accumulated_curtailed_wind = 0.0
    accumulated_curtailed_solar = 0.0

    # 计算总可用发电量上限（不考虑火电）
    total_available_energy = params.available_solar_energy + params.available_wind_energy
    max_total_grid_feed_in = total_available_energy * params.max_grid_ratio

    # 初始化累计上网电量
    cumulative_grid_feed_in = 0.0

    # 逐小时计算
    for hour in range(hours):
        # 确保储能状态不为负且不超过容量
        storage_soc[hour] = max(0, min(storage_capacity, storage_soc[hour]))

        # 可用发电量（未考虑弃风弃光）
        available_generation = thermal[hour] + wind[hour] + solar[hour]

        # 需求电量（负荷）
        demand = load[hour]

        # 计算功率平衡
        power_surplus = available_generation - demand

        if power_surplus > 0:
            # 有盈余电力，优先储能
            max_charge = min(
                storage_capacity - storage_soc[hour],  # 储能空间限制
                storage_capacity * params.charge_rate,  # 充电功率限制
                power_surplus  # 可用盈余电力
            )

            # 实际充电功率
            actual_charge = min(max_charge, power_surplus)
            charge_power[hour] = actual_charge

            # 更新储能状态（考虑充电效率）
            if hour < hours - 1:
                storage_soc[hour + 1] = storage_soc[hour] + actual_charge * params.storage_efficiency

            # 剩余电力优先上网，超出上网电量比例限制时才弃风弃光
            remaining_surplus = power_surplus - actual_charge

            # 计算当前允许的最大上网电量
            # 基于总可用发电量和已累计上网电量
            remaining_grid_quota = max_total_grid_feed_in - cumulative_grid_feed_in

            # 优先上网，但不超过剩余配额
            grid_feed_allowed = min(remaining_surplus, max(0, remaining_grid_quota))
            grid_feed_in[hour] = grid_feed_allowed
            cumulative_grid_feed_in += grid_feed_allowed

            # 剩余电力考虑弃风弃光
            remaining_surplus -= grid_feed_allowed

            # 计算弃风弃光（优先弃光，然后弃风）
            if remaining_surplus > 0:
                # 在单一可再生能源场景下，弃风弃光不受限制
                if single_renewable:
                    # 优先弃光
                    if solar[hour] > 0:
                        curtailed_solar[hour] = min(solar[hour], remaining_surplus)
                        remaining_surplus -= curtailed_solar[hour]
                        accumulated_curtailed_solar += curtailed_solar[hour]

                    # 如果还有剩余，弃风
                    if remaining_surplus > 0 and wind[hour] > 0:
                        curtailed_wind[hour] = min(wind[hour], remaining_surplus)
                        remaining_surplus -= curtailed_wind[hour]
                        accumulated_curtailed_wind += curtailed_wind[hour]
                else:
                    # 混合可再生能源场景，弃风弃光受限制
                    # 计算当前可弃光量
                    if max_curtailed_solar_total == float('inf'):
                        # 无限制弃光
                        available_solar_curtail = solar[hour]
                    else:
                        # 受限制弃光
                        available_solar_curtail = min(
                            solar[hour],
                            max_curtailed_solar_total - accumulated_curtailed_solar
                        )

                    if available_solar_curtail > 0:
                        curtailed_solar[hour] = min(available_solar_curtail, remaining_surplus)
                        remaining_surplus -= curtailed_solar[hour]
                        accumulated_curtailed_solar += curtailed_solar[hour]

                    # 如果还有剩余，弃风
                    if remaining_surplus > 0:
                        if max_curtailed_wind_total == float('inf'):
                            # 无限制弃风
                            available_wind_curtail = wind[hour]
                        else:
                            # 受限制弃风
                            available_wind_curtail = min(
                                wind[hour],
                                max_curtailed_wind_total - accumulated_curtailed_wind
                            )

                        if available_wind_curtail > 0:
                            curtailed_wind[hour] = min(available_wind_curtail, remaining_surplus)
                            remaining_surplus -= curtailed_wind[hour]
                            accumulated_curtailed_wind += curtailed_wind[hour]

                # 确保电力平衡：如果仍有剩余电力，强制弃掉（安全机制）
                if remaining_surplus > 0:
                    # 记录警告但不影响计算
                    # 在实际系统中，这种情况不应该发生，但作为安全保护
                    pass

        else:
            # 电力不足，优先放电
            power_deficit = -power_surplus
            grid_feed_in[hour] = 0  # 初始化购电为0

            max_discharge = min(
                storage_soc[hour],  # 储能状态限制
                storage_capacity * params.discharge_rate,  # 放电功率限制
                power_deficit  # 缺电功率
            )

            # 实际放电功率
            actual_discharge = min(max_discharge, power_deficit)
            discharge_power[hour] = actual_discharge

            # 更新储能状态（考虑放电效率）
            if hour < hours - 1:
                storage_soc[hour + 1] = storage_soc[hour] - actual_discharge / params.storage_efficiency

            # 计算剩余缺电，需要从电网购电
            remaining_deficit = power_deficit - actual_discharge

            # 如果还有缺电，从电网购电
            if remaining_deficit > 0:
                # 购电功率为负值，表示从电网输入
                grid_feed_in[hour] = -remaining_deficit

    return {
        'storage_profile': storage_soc.tolist(),
        'charge_profile': charge_power.tolist(),
        'discharge_profile': discharge_power.tolist(),
        'curtailed_wind': curtailed_wind.tolist(),
        'curtailed_solar': curtailed_solar.tolist(),
        'grid_feed_in': grid_feed_in.tolist()
    }


def find_periodic_steady_state(
    solar_output: List[float],
    wind_output: List[float],
    thermal_output: List[float],
    load_demand: List[float],
    params: SystemParameters,
    storage_capacity: float,
    soc_convergence_threshold: float = 0.001,
    max_iterations: int = 100
) -> Dict[str, List[float]]:
    """
    通过迭代找到满足周期性平衡的储能初始状态

    步骤：
    1. 从初始SOC=0开始，运行一次全年仿真，记录最后一小时的SOC值
    2. 将这个SOC值作为新的"初始SOC"，再次运行仿真
    3. 重复上述过程，直到首尾SOC的差值小于设定的阈值

    Args:
        solar_output: 光伏出力曲线 (MW) - 支持24小时或8760小时
        wind_output: 风电出力曲线 (MW) - 支持24小时或8760小时
        thermal_output: 火电出力曲线 (MW) - 支持24小时或8760小时
        load_demand: 负荷曲线 (MW) - 支持24小时或8760小时
        params: 系统参数配置
        storage_capacity: 储能容量 (MWh)
        soc_convergence_threshold: SOC收敛阈值（相对于容量的比例），默认0.1%
        max_iterations: 最大迭代次数

    Returns:
        包含各种功率曲线的字典，且满足周期性平衡条件
    """
    # 计算收敛阈值的绝对值
    absolute_threshold = storage_capacity * soc_convergence_threshold

    # 初始SOC从0开始
    initial_soc = 0.0
    iteration = 0
    soc_diff = float('inf')

    print(f"正在寻找周期性平衡状态（SOC收敛阈值: {absolute_threshold:.4f} MWh）...")

    while iteration < max_iterations and soc_diff > absolute_threshold:
        # 运行仿真
        balance_result = calculate_energy_balance(
            solar_output, wind_output, thermal_output, load_demand,
            params, storage_capacity, initial_soc
        )

        # 获取初始和最终的SOC
        storage_initial = balance_result['storage_profile'][0]
        storage_final = balance_result['storage_profile'][-1]

        # 计算SOC差值
        soc_diff = abs(storage_final - storage_initial)

        # 更新初始SOC（使用最终SOC作为下一次的初始SOC）
        initial_soc = storage_final

        # 确保SOC在合理范围内
        initial_soc = max(0, min(storage_capacity, initial_soc))

        iteration += 1

        # 输出迭代信息（每10次迭代或最后一次）
        if iteration % 10 == 0 or iteration == 1 or soc_diff <= absolute_threshold:
            print(f"  迭代 {iteration}: 初始SOC={storage_initial:.4f} MWh, "
                  f"最终SOC={storage_final:.4f} MWh, 差值={soc_diff:.4f} MWh")

    # 输出收敛结果
    if soc_diff <= absolute_threshold:
        print(f"✓ 周期性平衡收敛成功（迭代{iteration}次，SOC差值={soc_diff:.4f} MWh）")
    else:
        print(f"⚠ 未达到收敛条件（迭代{iteration}次，SOC差值={soc_diff:.4f} MWh）")

    return balance_result


def check_constraints(
    solar_output: List[float],
    wind_output: List[float],
    thermal_output: List[float],
    balance_result: Dict[str, List[float]],
    params: SystemParameters
) -> Dict[str, float]:
    """
    检查约束条件是否满足

    Args:
        solar_output: 光伏出力曲线 (MW) - 支持24小时或8760小时
        wind_output: 风电出力曲线 (MW) - 支持24小时或8760小时
        thermal_output: 火电出力曲线 (MW) - 支持24小时或8760小时
        balance_result: 电能平衡计算结果
        params: 系统参数配置

    Returns:
        包含各约束实际比例的字典
    """
    # 计算总量
    total_wind_potential = sum(wind_output)
    total_solar_potential = sum(solar_output)
    total_thermal = sum(thermal_output)

    total_curtailed_wind = sum(balance_result['curtailed_wind'])
    total_curtailed_solar = sum(balance_result['curtailed_solar'])
    total_grid_feed_in = sum(balance_result['grid_feed_in'])

    # 实际发电量（考虑弃风弃光）
    actual_wind_generation = total_wind_potential - total_curtailed_wind
    actual_solar_generation = total_solar_potential - total_curtailed_solar
    total_generation = total_thermal + actual_wind_generation + actual_solar_generation

    # 计算比例
    actual_curtailment_wind_ratio = total_curtailed_wind / total_wind_potential if total_wind_potential > 0 else 0
    actual_curtailment_solar_ratio = total_curtailed_solar / total_solar_potential if total_solar_potential > 0 else 0
    actual_grid_feed_in_ratio = total_grid_feed_in / total_generation if total_generation > 0 else 0

    return {
        'total_curtailment_wind_ratio': actual_curtailment_wind_ratio,
        'total_curtailment_solar_ratio': actual_curtailment_solar_ratio,
        'total_grid_feed_in_ratio': actual_grid_feed_in_ratio
    }


def optimize_storage_capacity(
    solar_output: List[float],
    wind_output: List[float],
    thermal_output: List[float],
    load_demand: List[float],
    params: SystemParameters,
    max_iterations: int = 100,
    tolerance: float = 0.01,
    search_step: float = 0.01
) -> Dict:
    """
    优化储能容量，在不超过设定储能容量上限的前提下，选择使弃电量最小的储能容量

    Args:
        solar_output: 光伏出力曲线 (MW) - 支持24小时或8760小时
        wind_output: 风电出力曲线 (MW) - 支持24小时或8760小时
        thermal_output: 火电出力曲线 (MW) - 支持24小时或8760小时
        load_demand: 负荷曲线 (MW) - 支持24小时或8760小时
        params: 系统参数配置
        max_iterations: 最大迭代次数
        tolerance: 收敛容差
        search_step: 搜索步长（相对于最大容量的比例）

    Returns:
        包含优化结果的字典
    """
    # 验证输入
    validate_inputs(solar_output, wind_output, thermal_output, load_demand, params)

    # 初始化搜索范围
    lower_bound = 0.0
    theoretical_max = max(sum(solar_output) + sum(wind_output) + sum(thermal_output), sum(load_demand))

    # 应用储能容量上限限制
    if params.max_storage_capacity is not None:
        upper_bound = min(theoretical_max, params.max_storage_capacity)
    else:
        upper_bound = theoretical_max
        print("警告：未设置储能容量上限，将使用理论最大值")

    # 判断数据类型（24小时或8760小时）
    data_length = len(solar_output)
    is_yearly_data = data_length == 8760

    if is_yearly_data:
        print(f"处理8760小时全年数据，启用周期性平衡优化...")

    # 在容量范围内搜索，找到弃电量最小的储能容量
    # 使用二分搜索 + 局部搜索相结合的方法
    
    # 首先使用二分搜索找到一个满足约束的基准容量
    print("第一阶段：二分搜索寻找满足约束的基准容量...")
    
    best_capacity = upper_bound
    best_result = None
    best_curtailed = float('inf')
    solution_found = False
    
    # 二分搜索寻找最小可行容量
    lower_bound_search = lower_bound
    upper_bound_search = upper_bound
    feasible_capacity = None
    
    for iteration in range(max_iterations):
        mid_capacity = (lower_bound_search + upper_bound_search) / 2

        # 计算当前容量下的平衡
        if is_yearly_data:
            balance_result = find_periodic_steady_state(
                solar_output, wind_output, thermal_output, load_demand,
                params, mid_capacity
            )
        else:
            balance_result = calculate_energy_balance(
                solar_output, wind_output, thermal_output, load_demand, params, mid_capacity
            )

        # 检查约束条件
        constraint_results = check_constraints(solar_output, wind_output, thermal_output, balance_result, params)

        # 检查是否满足所有约束
        total_grid_feed_in = sum(balance_result['grid_feed_in'])
        if total_grid_feed_in > 0:
            grid_constraint_satisfied = constraint_results['total_grid_feed_in_ratio'] <= params.max_grid_ratio
        else:
            grid_constraint_satisfied = True

        has_wind = sum(wind_output) > 0
        has_solar = sum(solar_output) > 0
        single_renewable = (has_wind and not has_solar) or (has_solar and not has_wind)
        grid_quota_zero = params.max_grid_ratio == 0

        if single_renewable or grid_quota_zero:
            constraints_satisfied = grid_constraint_satisfied
        else:
            constraints_satisfied = (
                constraint_results['total_curtailment_wind_ratio'] <= params.max_curtailment_wind and
                constraint_results['total_curtailment_solar_ratio'] <= params.max_curtailment_solar and
                grid_constraint_satisfied
            )

        # 检查储能日平衡（周期结束时储能状态应接近初始值）
        storage_initial = balance_result['storage_profile'][0]
        storage_final = balance_result['storage_profile'][-1]
        daily_balance = abs(storage_final - storage_initial) < tolerance

        if constraints_satisfied and daily_balance:
            # 满足条件，尝试减小容量
            feasible_capacity = mid_capacity
            upper_bound_search = mid_capacity
        else:
            # 不满足条件，增大容量
            lower_bound_search = mid_capacity

        # 检查收敛
        if upper_bound_search - lower_bound_search < tolerance:
            break
    
    # 如果找到了可行解，使用它作为基准；否则使用上限
    if feasible_capacity is not None:
        print(f"找到基准可行容量: {feasible_capacity:.2f} MWh")
    else:
        print(f"未找到可行解，使用上限容量: {upper_bound:.2f} MWh")
        feasible_capacity = upper_bound
    
    # 第二阶段：在 [feasible_capacity, upper_bound] 范围内搜索弃电量最小的容量
    print("\n第二阶段：在可行容量范围内搜索弃电量最小的储能容量...")
    
    # 使用梯度下降方法搜索最优容量
    # 定义搜索区间
    search_lower = feasible_capacity
    search_upper = upper_bound
    
    # 计算当前最小弃电量
    best_capacity = feasible_capacity
    
    # 计算基准容量的弃电量
    if is_yearly_data:
        base_result = find_periodic_steady_state(
            solar_output, wind_output, thermal_output, load_demand,
            params, best_capacity
        )
    else:
        base_result = calculate_energy_balance(
            solar_output, wind_output, thermal_output, load_demand, params, best_capacity
        )
    
    base_curtailed = sum(base_result['curtailed_wind']) + sum(base_result['curtailed_solar'])
    best_curtailed = base_curtailed
    best_result = {**base_result, **check_constraints(solar_output, wind_output, thermal_output, base_result, params)}
    
    print(f"基准容量 {best_capacity:.2f} MWh 的弃电量: {best_curtailed:.2f} MWh")
    
    # 使用黄金分割搜索法寻找最优容量
    # 黄金分割点
    golden_ratio = (3 - 5**0.5) / 2  # 约0.382
    
    # 初始化两个测试点
    a = search_lower
    b = search_upper
    
    c = b - golden_ratio * (b - a)
    d = a + golden_ratio * (b - a)
    
    max_refinement_iterations = 50
    
    for iter_num in range(max_refinement_iterations):
        # 计算点 c 的弃电量
        if is_yearly_data:
            result_c = find_periodic_steady_state(
                solar_output, wind_output, thermal_output, load_demand,
                params, c
            )
        else:
            result_c = calculate_energy_balance(
                solar_output, wind_output, thermal_output, load_demand, params, c
            )
        
        curtailed_c = sum(result_c['curtailed_wind']) + sum(result_c['curtailed_solar'])
        
        # 计算点 d 的弃电量
        if is_yearly_data:
            result_d = find_periodic_steady_state(
                solar_output, wind_output, thermal_output, load_demand,
                params, d
            )
        else:
            result_d = calculate_energy_balance(
                solar_output, wind_output, thermal_output, load_demand, params, d
            )
        
        curtailed_d = sum(result_d['curtailed_wind']) + sum(result_d['curtailed_solar'])
        
        # 比较并更新最优解
        if curtailed_c < best_curtailed:
            best_curtailed = curtailed_c
            best_capacity = c
            best_result = {**result_c, **check_constraints(solar_output, wind_output, thermal_output, result_c, params)}
            print(f"  发现更优容量: {best_capacity:.2f} MWh, 弃电量: {best_curtailed:.2f} MWh")
        
        if curtailed_d < best_curtailed:
            best_curtailed = curtailed_d
            best_capacity = d
            best_result = {**result_d, **check_constraints(solar_output, wind_output, thermal_output, result_d, params)}
            print(f"  发现更优容量: {best_capacity:.2f} MWh, 弃电量: {best_curtailed:.2f} MWh")
        
        # 缩小搜索区间
        if curtailed_c < curtailed_d:
            b = d
            d = c
            c = b - golden_ratio * (b - a)
        else:
            a = c
            c = d
            d = a + golden_ratio * (b - a)
        
        # 检查收敛
        if b - a < tolerance:
            print(f"搜索收敛，区间宽度: {b - a:.4f} MWh")
            break
        
        # 限制最大迭代次数
        if iter_num % 10 == 0 and iter_num > 0:
            print(f"  已完成 {iter_num} 次迭代，当前搜索区间: [{a:.2f}, {b:.2f}] MWh")
    
    print(f"\n最终选择储能容量: {best_capacity:.2f} MWh (容量上限: {upper_bound:.2f} MWh)")
    print(f"最终弃电量: {best_curtailed:.2f} MWh")

    # 添加能量平衡校验
    total_generation = sum(thermal_output) + sum(wind_output) + sum(solar_output)
    total_consumption = sum(load_demand)
    total_curtailed = sum(best_result['curtailed_wind']) + sum(best_result['curtailed_solar'])
    total_grid = sum(best_result['grid_feed_in'])
    total_charge = sum(best_result['charge_profile'])
    total_discharge = sum(best_result['discharge_profile'])
    storage_net_change = best_result['storage_profile'][-1] - best_result['storage_profile'][0]

    # 能量平衡校验：发电量 + 放电量/效率 = 负荷 + 充电量*效率 + 弃风弃光 + 上网电量
    # 考虑储能充放电效率的能量平衡
    energy_from_storage = total_discharge / params.storage_efficiency  # 储能提供的有效能量
    energy_to_storage = total_charge * params.storage_efficiency      # 储能消耗的电网能量

    # 能量平衡校验：应该接近0，但允许一定误差
    # 当total_grid为负时（购电），应该加到左侧（供给侧）
    # 当total_grid为正时（上网），应该加到右侧（需求侧）
    if total_grid < 0:  # 购电情况
        energy_balance_error = abs(
            total_generation + energy_from_storage + abs(total_grid) - total_consumption - energy_to_storage - total_curtailed
        )
    else:  # 上网情况
        energy_balance_error = abs(
            total_generation + energy_from_storage - total_consumption - energy_to_storage - total_curtailed - total_grid
        )
    # 使用更大的容差，考虑储能效率损失和数值误差
    # 允许误差为总发电量的15%或10MW，取较大者
    # 储能效率损失可能达到总能量的10%以上
    tolerance = max(10.0, total_generation * 0.15)
    energy_balance_check = energy_balance_error < tolerance

    # 输出周期性平衡信息
    if is_yearly_data:
        soc_initial_final_diff = abs(best_result['storage_profile'][-1] - best_result['storage_profile'][0])
        print(f"\n周期性平衡信息:")
        print(f"  初始SOC: {best_result['storage_profile'][0]:.4f} MWh")
        print(f"  最终SOC: {best_result['storage_profile'][-1]:.4f} MWh")
        print(f"  SOC差值: {soc_initial_final_diff:.4f} MWh")

    # 最终结果
    result = {
        'required_storage_capacity': best_capacity,
        'storage_profile': best_result['storage_profile'],
        'charge_profile': best_result['charge_profile'],
        'discharge_profile': best_result['discharge_profile'],
        'curtailed_wind': best_result['curtailed_wind'],
        'curtailed_solar': best_result['curtailed_solar'],
        'grid_feed_in': best_result['grid_feed_in'],
        'total_curtailment_wind_ratio': best_result['total_curtailment_wind_ratio'],
        'total_curtailment_solar_ratio': best_result['total_curtailment_solar_ratio'],
        'total_grid_feed_in_ratio': best_result['total_grid_feed_in_ratio'],
        'energy_balance_check': energy_balance_check,
        'capacity_limit_reached': params.max_storage_capacity is not None and best_capacity >= params.max_storage_capacity,
        'total_curtailed_energy': best_curtailed,  # 总弃电量
        'min_curtailed_capacity': best_capacity,  # 弃电量最小时的储能容量
        'max_storage_limit': params.max_storage_capacity,
        'optimization_goal': 'minimize_curtailment'  # 优化目标：最小化弃电量
    }
    
    return result


def main():
    """主函数，提供示例使用"""
    # 示例数据
    solar_output = [0.0] * 6 + [1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 5.0, 4.0, 3.0, 2.0, 1.0, 0.0] * 2
    wind_output = [2.0, 3.0, 4.0, 3.0, 2.0, 1.0] * 4
    thermal_output = [5.0] * 24
    load_demand = [3.0, 4.0, 5.0, 6.0, 8.0, 10.0, 12.0, 14.0, 16.0, 18.0, 20.0, 18.0, 
                   16.0, 14.0, 12.0, 10.0, 8.0, 6.0, 5.0, 4.0, 3.0, 2.0, 1.0, 2.0]
    
    # 系统参数
    params = SystemParameters(
        max_curtailment_wind=0.1,
        max_curtailment_solar=0.1,
        max_grid_ratio=0.2,
        storage_efficiency=0.9,
        discharge_rate=1.0,
        charge_rate=1.0
    )
    
    # 计算最优储能容量
    result = optimize_storage_capacity(
        solar_output, wind_output, thermal_output, load_demand, params
    )
    
    # 打印结果
    print("多能互补系统储能容量优化结果:")
    print(f"所需储能总容量: {result['required_storage_capacity']:.2f} MWh")
    print(f"实际弃风率: {result['total_curtailment_wind_ratio']:.3f}")
    print(f"实际弃光率: {result['total_curtailment_solar_ratio']:.3f}")
    print(f"实际上网电量比例: {result['total_grid_feed_in_ratio']:.3f}")
    print(f"能量平衡校验: {'通过' if result['energy_balance_check'] else '未通过'}")
    
    return result


if __name__ == "__main__":
    main()