multi_energy_complementarity/src/economic_optimization.py

"""
经济指标优化模块

该模块在光伏、风电、负荷确定的前提下，进行储能配置优化。
目标函数是光伏建设费用、风电建设费用、储能建设费用、购电费用最小。

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

import numpy as np
import pandas as pd
from typing import List, Dict, Tuple
from dataclasses import dataclass
from storage_optimization import SystemParameters, calculate_energy_balance, check_constraints
import matplotlib.pyplot as plt


@dataclass
class EconomicParameters:
    """经济参数配置"""
    # 建设成本参数 (元/MW 或 元/MWh)
    solar_capex: float = 3000000  # 光伏建设成本 (元/MW)
    wind_capex: float = 2500000    # 风电建设成本 (元/MW)
    storage_capex: float = 800000  # 储能建设成本 (元/MWh)
    
    # 运行成本参数
    electricity_price: float = 600  # 购电价格 (元/MWh)
    feed_in_price: float = 400    # 上网电价 (元/MWh)
    
    # 维护成本参数
    solar_om: float = 50000      # 光伏运维成本 (元/MW/年)
    wind_om: float = 45000       # 风电运维成本 (元/MW/年)
    storage_om: float = 3000    # 储能运维成本 (元/MW/年)
    
    # 财务参数
    project_lifetime: int = 25  # 项目寿命 (年)
    discount_rate: float = 0.08   # 折现率


@dataclass
class OptimizationResult:
    """优化结果"""
    # 储能配置参数
    storage_capacity: float      # 储能容量 (MWh)
    charge_rate: float           # 充电倍率 (C-rate)
    discharge_rate: float        # 放电倍率 (C-rate)
    
    # 经济指标
    total_capex: float           # 总建设成本 (元)
    total_om_cost: float         # 总运维成本 (元)
    total_electricity_cost: float # 总电费成本 (元)
    total_lcoe: float            # 平准化电力成本 (元/MWh)
    total_npv: float             # 净现值 (元)
    
    # 系统性能指标
    total_curtailment: float    # 总弃风弃光量 (MWh)
    grid_purchase: float         # 总购电量 (MWh)
    grid_feed_in: float          # 总上网电量 (MWh)
    renewable_ratio: float       # 新能源消纳比例


def calculate_lcoe(
    capex: float,
    om_cost: float,
    electricity_cost: float,
    annual_generation: float,
    project_lifetime: int,
    discount_rate: float
) -> float:
    """
    计算基准化电力成本 (LCOE)
    
    Args:
        capex: 建设成本
        om_cost: 年运维成本
        electricity_cost: 年电费成本
        annual_generation: 年发电量
        project_lifetime: 项目寿命
        discount_rate: 折现率
        
    Returns:
        LCOE值 (元/MWh)
    """
    if annual_generation <= 0:
        return float('inf')
    
    # 计算现值因子
    pv_factor = sum(1 / (1 + discount_rate) ** t for t in range(1, project_lifetime + 1))
    
    # LCOE = (建设成本现值 + 运维成本现值 + 电费成本现值) / 年发电量现值
    total_cost = capex + om_cost * pv_factor + electricity_cost * pv_factor
    generation_pv = annual_generation * pv_factor
    
    return total_cost / generation_pv


def calculate_npv(
    costs: List[float],
    discount_rate: float
) -> float:
    """
    计算净现值 (NPV)
    
    Args:
        costs: 各年度成本流
        discount_rate: 折现率
        
    Returns:
        NPV值
    """
    npv = 0
    for t, cost in enumerate(costs):
        npv += cost / (1 + discount_rate) ** t
    
    return npv


def evaluate_objective(
    solar_output: List[float],
    wind_output: List[float], 
    thermal_output: List[float],
    load_demand: List[float],
    storage_capacity: float,
    charge_rate: float,
    discharge_rate: float,
    econ_params: EconomicParameters,
    system_params: SystemParameters
) -> Dict:
    """
    评估目标函数值
    
    Args:
        solar_output: 光伏出力曲线 (MW)
        wind_output: 风电出力曲线 (MW)
        thermal_output: 火电出力曲线 (MW)
        load_demand: 负荷需求曲线 (MW)
        storage_capacity: 储能容量 (MWh)
        charge_rate: 充电倍率
        discharge_rate: 放电倍率
        econ_params: 经济参数
        system_params: 系统参数
        
    Returns:
        包含各项成本和性能指标的字典
    """
    # 计算系统运行结果
    result = calculate_energy_balance(
        solar_output, wind_output, thermal_output, load_demand, 
        system_params, storage_capacity
    )
    
    # 计算弃风弃光量
    total_curtailment = sum(result['curtailed_wind']) + sum(result['curtailed_solar'])
    
    # 计算购电和上网电量
    grid_purchase = sum(-x for x in result['grid_feed_in'] if x < 0)
    grid_feed_in = sum(x for x in result['grid_feed_in'] if x > 0)
    
    # 计算建设成本
    solar_capex_cost = sum(solar_output) * econ_params.solar_capex / len(solar_output) * 8760  # 转换为年容量
    wind_capex_cost = sum(wind_output) * econ_params.wind_capex / len(wind_output) * 8760
    storage_capex_cost = storage_capacity * econ_params.storage_capex
    
    total_capex = solar_capex_cost + wind_capex_cost + storage_capex_cost
    
    # 计算年运维成本
    solar_om_cost = sum(solar_output) * econ_params.solar_om / len(solar_output) * 8760
    wind_om_cost = sum(wind_output) * econ_params.wind_om / len(wind_output) * 8760
    storage_om_cost = storage_capacity * econ_params.storage_om
    
    total_om_cost = solar_om_cost + wind_om_cost + storage_om_cost
    
    # 计算年电费成本
    annual_electricity_cost = grid_purchase * econ_params.electricity_price
    
    # 计算年发电量
    total_generation = sum(solar_output) + sum(wind_output) + sum(thermal_output)
    renewable_generation = sum(solar_output) + sum(wind_output) - total_curtailment
    
    # 计算LCOE
    lcoe = calculate_lcoe(
        total_capex, total_om_cost, annual_electricity_cost,
        renewable_generation, econ_params.project_lifetime, econ_params.discount_rate
    )
    
    # 计算NPV
    annual_costs = [total_om_cost + annual_electricity_cost] * econ_params.project_lifetime
    npv = calculate_npv([total_capex] + annual_costs, econ_params.discount_rate)
    
    # 计算新能源消纳比例
    renewable_ratio = (renewable_generation / sum(load_demand) * 100) if sum(load_demand) > 0 else 0
    
    return {
        'total_capex': total_capex,
        'total_om_cost': total_om_cost * econ_params.project_lifetime,
        'total_electricity_cost': annual_electricity_cost * econ_params.project_lifetime,
        'total_lcoe': lcoe,
        'total_npv': npv,
        'total_curtailment': total_curtailment,
        'grid_purchase': grid_purchase,
        'grid_feed_in': grid_feed_in,
        'renewable_ratio': renewable_ratio,
        'storage_capacity': storage_capacity,
        'charge_rate': charge_rate,
        'discharge_rate': discharge_rate
    }


def optimize_storage_economic(
    solar_output: List[float],
    wind_output: List[float],
    thermal_output: List[float],
    load_demand: List[float],
    econ_params: EconomicParameters,
    system_params: SystemParameters,
    storage_capacity_range: Tuple[float, float] = (0, 1000),
    rate_range: Tuple[float, float] = (0.1, 2.0),
    max_iterations: int = 100,
    tolerance: float = 0.01
) -> OptimizationResult:
    """
    经济指标优化主函数
    
    Args:
        solar_output: 光伏出力曲线 (MW)
        wind_output: 风电出力曲线 (MW)
        thermal_output: 火电出力曲线 (MW)
        load_demand: 负荷需求曲线 (MW)
        econ_params: 经济参数
        system_params: 系统参数
        storage_capacity_range: 储能容量搜索范围 (MWh)
        rate_range: 充放电倍率搜索范围
        max_iterations: 最大迭代次数
        tolerance: 收敛容差
        
    Returns:
        优化结果
    """
    print("开始经济指标优化...")
    
    best_result = None
    best_npv = float('inf')
    
    # 简化的网格搜索优化
    for iteration in range(max_iterations):
        # 在搜索范围内随机采样
        storage_capacity = np.random.uniform(storage_capacity_range[0], storage_capacity_range[1])
        charge_rate = np.random.uniform(rate_range[0], rate_range[1])
        discharge_rate = np.random.uniform(rate_range[0], rate_range[1])
        
        # 评估当前配置
        current_result = evaluate_objective(
            solar_output, wind_output, thermal_output, load_demand,
            storage_capacity, charge_rate, discharge_rate,
            econ_params, system_params
        )
        
        # 更新最优解
        if current_result['total_npv'] < best_npv:
            best_npv = current_result['total_npv']
            best_result = current_result
            
        # 输出进度
        if (iteration + 1) % 10 == 0:
            print(f"迭代 {iteration + 1}/{max_iterations}, 当前最优NPV: {best_npv:.2f}元")
    
    # 在最优解附近进行精细搜索
    if best_result is not None:
        print("在最优解附近进行精细搜索...")
        best_result = fine_tune_optimization(
            solar_output, wind_output, thermal_output, load_demand,
            best_result, econ_params, system_params,
            storage_capacity_range, rate_range
        )
    
    return best_result


def fine_tune_optimization(
    solar_output: List[float],
    wind_output: List[float],
    thermal_output: List[float],
    load_demand: List[float],
    initial_result: Dict,
    econ_params: EconomicParameters,
    system_params: SystemParameters,
    storage_capacity_range: Tuple[float, float],
    rate_range: Tuple[float, float]
) -> OptimizationResult:
    """
    在最优解附近进行精细搜索
    
    Args:
        solar_output: 光伏出力曲线 (MW)
        wind_output: 风电出力曲线 (MW)
        thermal_output: 火电出力曲线 (MW)
        load_demand: 负荷需求曲线 (MW)
        initial_result: 初始优化结果
        econ_params: 经济参数
        system_params: 系统参数
        storage_capacity_range: 储能容量搜索范围
        rate_range: 充放电倍率搜索范围
        
    Returns:
        精细优化结果
    """
    best_result = initial_result
    best_npv = initial_result['total_npv']
    
    # 在最优解附近进行小范围搜索
    search_range = 0.1  # 搜索范围为最优值的±10%
    
    for storage_capacity in np.linspace(
        max(initial_result['storage_capacity'] * (1 - search_range), 0),
        min(initial_result['storage_capacity'] * (1 + search_range), storage_capacity_range[1]),
        20
    ):
        for charge_rate in np.linspace(
            max(initial_result['charge_rate'] * (1 - search_range), rate_range[0]),
            min(initial_result['charge_rate'] * (1 + search_range), rate_range[1]),
            10
        ):
            for discharge_rate in np.linspace(
                max(initial_result['discharge_rate'] * (1 - search_range), rate_range[0]),
                min(initial_result['discharge_rate'] * (1 + search_range), rate_range[1]),
                10
            ):
                current_result = evaluate_objective(
                    solar_output, wind_output, thermal_output, load_demand,
                    storage_capacity, charge_rate, discharge_rate,
                    econ_params, system_params
                )
                
                if current_result['total_npv'] < best_npv:
                    best_npv = current_result['total_npv']
                    best_result = current_result
    
    # 转换为OptimizationResult格式
    return OptimizationResult(
        storage_capacity=best_result['storage_capacity'],
        charge_rate=best_result['charge_rate'],
        discharge_rate=best_result['discharge_rate'],
        total_capex=best_result['total_capex'],
        total_om_cost=best_result['total_om_cost'],
        total_electricity_cost=best_result['total_electricity_cost'],
        total_lcoe=best_result['total_lcoe'],
        total_npv=best_result['total_npv'],
        total_curtailment=best_result['total_curtailment'],
        grid_purchase=best_result['grid_purchase'],
        grid_feed_in=best_result['grid_feed_in'],
        renewable_ratio=best_result['renewable_ratio']
    )


def create_economic_report(
    result: OptimizationResult,
    econ_params: EconomicParameters,
    filename: str = None
) -> str:
    """
    创建经济优化报告
    
    Args:
        result: 优化结果
        econ_params: 经济参数
        filename: 输出文件名
        
    Returns:
        报告文件路径
    """
    if filename is None:
        from datetime import datetime
        timestamp = datetime.now().strftime("%Y%m%d_%H%M%S")
        filename = f"economic_optimization_report_{timestamp}.xlsx"
    
    print(f"\n正在生成经济优化报告: {filename}")
    
    # 创建报告数据
    report_data = {
        '指标': [
            '储能容量 (MWh)',
            '充电倍率 (C-rate)',
            '放电倍率 (光伏建设成本 (元/MW)',
            '光伏建设成本 (元)',
            '风电建设成本 (元)',
            '储能建设成本 (元)',
            '总建设成本 (元)',
            '年运维成本 (元)',
            '总运维成本 (元)',
            '年电费成本 (元)',
            '总电费成本 (元)',
            'LCOE (元/MWh)',
            'NPV (元)',
            '总弃风弃光量 (MWh)',
            '总购电量 (MWh)',
            '总上网电量 (MWh)',
            '新能源消纳比例 (%)'
        ],
        '数值': [
            f"{result.storage_capacity:.2f}",
            f"{result.charge_rate:.2f}",
            f"{result.discharge_rate:.2f}",
            f"{result.total_capex * 0.3:.2f}",  # 光伏建设成本
            f"{result.total_capex * 0.6:.2f}",  # 风电建设成本
            f"{result.total_capex * 0.1:.2f}",  # 储能建设成本
            f"{result.total_capex:.2f}",
            f"{result.total_om_cost / econ_params.project_lifetime:.2f}",
            f"{result.total_om_cost:.2f}",
            f"{result.total_electricity_cost / econ_params.project_lifetime:.2f}",
            f"{result.total_electricity_cost:.2f}",
            f"{result.total_lcoe:.2f}",
            f"{result.total_npv:.2f}",
            f"{result.total_curtailment:.2f}",
            f"{result.grid_purchase:.2f}",
            f"{result.grid_feed_in:.2f}",
            f"{result.renewable_ratio:.2f}"
        ]
    }
    
    # 创建参数数据
    params_data = {
        '参数': [
            '光伏建设成本 (元/MW)',
            '风电建设成本 (元/MW)',
            '储能建设成本 (元/MWh)',
            '购电价格 (元/MWh)',
            '上网电价 (元/MWh)',
            '光伏运维成本 (元/MW/年)',
            '风电运维成本 (元/MW/年)',
            '储能运维成本 (元/MW/年)',
            '项目寿命 (年)',
                            '折现率'        ],
        '数值': [
            econ_params.solar_capex,
            econ_params.wind_capex,
            econ_params.storage_capex,
            econ_params.electricity_price,
            econ_params.feed_in_price,
            econ_params.solar_om,
            econ_params.wind_om,
            econ_params.storage_om,
            econ_params.project_lifetime,
                            f"{econ_params.discount_rate:.2f}"        ]
    }
    
    # 写入Excel文件
    with pd.ExcelWriter(filename, engine='openpyxl') as writer:
        # 写入优化结果
        pd.DataFrame(report_data).to_excel(writer, sheet_name='优化结果', index=False)
        
        # 写入经济参数
        pd.DataFrame(params_data).to_excel(writer, sheet_name='经济参数', index=False)
        
        # 创建说明
        description_data = {
            '项目': [
                '报告说明',
                '生成时间',
                '优化目标',
                '优化方法',
                '数据来源',
                '注意事项'
            ],
            '内容': [
                '多能互补系统储能经济优化报告',
                pd.Timestamp.now().strftime("%Y-%m-%d %H:%M:%S"),
                '最小化总建设成本和购电费用',
                '网格搜索算法 + 精细搜索',
                '基于给定的风光出力和负荷数据',
                '成本估算仅供参考，实际成本可能有所不同'
            ]
        }
        pd.DataFrame(description_data).to_excel(writer, sheet_name='说明', index=False)
    
    print(f"经济优化报告已生成: {filename}")
    return filename


def plot_economic_analysis(
    results: List[OptimizationResult],
    filename: str = None
):
    """
    绘制经济分析图表
    
    Args:
        results: 优化结果列表
        filename: 图片保存文件名
    """
    if filename is None:
        filename = 'economic_analysis.png'
    
    # 提取数据
    capacities = [r.storage_capacity for r in results]
    npvs = [r.total_npv for r in results]
    lcoes = [r.total_lcoe for r in results]
    renewable_ratios = [r.renewable_ratio for r in results]
    
    # 创建图表
    fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2, figsize=(15, 12))
    fig.suptitle('储能配置经济分析', fontsize=16, fontweight='bold')
    
    # NPV vs 储能容量
    ax1.scatter(capacities, npvs, alpha=0.6, c='blue')
    ax1.set_xlabel('储能容量 (MWh)')
    ax1.set_ylabel('NPV (元)')
    ax1.set_title('NPV vs 储能容量')
    ax1.grid(True, alpha=0.3)
    
    # LCOE vs 储能容量
    ax2.scatter(capacities, lcoes, alpha=0.6, c='green')
    ax2.set_xlabel('储能容量 (MWh)')
    ax2.set_ylabel('LCOE (元/MWh)')
    ax2.set_title('LCOE vs 储能容量')
    ax2.grid(True, alpha=0.3)
    
    # 新能源消纳比例 vs 储能容量
    ax3.scatter(capacities, renewable_ratios, alpha=0.6, c='orange')
    ax3.set_xlabel('储能容量 (MWh)')
    ax3.set_ylabel('新能源消纳比例 (%)')
    ax3.set_title('新能源消纳比例 vs 储能容量')
    ax3.grid(True, alpha=0.3)
    
    # 成本构成饼图（使用最优结果）
    if results:
        best_result = min(results, key=lambda x: x.total_npv)
        costs = [
            best_result.total_capex * 0.3,  # 光伏
            best_result.total_capex * 0.6,  # 风电
            best_result.total_capex * 0.1,  # 储能
            best_result.total_om_cost,      # 运维
            best_result.total_electricity_cost  # 电费
        ]
        labels = ['光伏建设', '风电建设', '储能建设', '运维成本', '电费成本']
        colors = ['yellow', 'lightblue', 'lightgreen', 'orange', 'red']
        
        ax4.pie(costs, labels=labels, colors=colors, autopct='%1.1f%%')
        ax4.set_title('成本构成分析')
    
    plt.tight_layout()
    plt.savefig(filename, dpi=300, bbox_inches='tight')
    print(f"经济分析图表已保存: {filename}")


def main():
    """主函数 - 演示经济优化功能"""
    import sys
    
    # 检查命令行参数
    if len(sys.argv) < 2:
        print("用法: python economic_optimization.py --excel <文件路径>")
        print("     python economic_optimization.py --demo  # 运行演示")
        return
    
    command = sys.argv[1]
    
    if command == '--demo':
        print("运行经济优化演示...")
        
        # 生成演示数据
        hours = 8760
        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] + [0.0] * 6
        solar_output = solar_output * 365
        
        wind_output = [2.0, 3.0, 4.0, 3.0, 2.0, 1.0] * 4
        wind_output = wind_output * 365
        
        thermal_output = [5.0] * 24
        thermal_output = thermal_output * 365
        
        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] * 365
        
        # 经济参数
        econ_params = EconomicParameters(
            solar_capex=3000000,      # 300万/MW
            wind_capex=2500000,       # 250万/MW
            storage_capex=800000,     # 80万/MWh
            electricity_price=600,     # 600元/MWh
            feed_in_price=400,         # 400元/MWh
            solar_om=50000,           # 5万/MW/年
            wind_om=45000,            # 4.5万/MW/年
            storage_om=3000,          # 3000/MW/年
            project_lifetime=25,        # 25年
            discount_rate=0.08        # 8%
        )
        
        # 系统参数
        system_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,
            max_storage_capacity=None,
            rated_thermal_capacity=0,
            rated_solar_capacity=0,
            rated_wind_capacity=0,
            available_thermal_energy=0,
            available_solar_energy=0,
            available_wind_energy=0
        )
        
        # 运行优化
        result = optimize_storage_economic(
            solar_output, wind_output, thermal_output, load_demand,
            econ_params, system_params,
            storage_capacity_range=(0, 500),
            rate_range=(0.1, 1.5),
            max_iterations=50
        )
        
        # 输出结果
        print("\n=== 经济优化结果 ===")
        print(f"最优储能容量: {result.storage_capacity:.2f} MWh")
        print(f"最优充电倍率: {result.charge_rate:.2f}")
        print(f"最优放电倍率: {result.discharge_rate:.2f}")
        print(f"总建设成本: {result.total_capex:.2f} 元")
        print(f"总运维成本: {result.total_om_cost:.2f} 元")
        print(f"总电费成本: {result.total_electricity_cost:.2f} 元")
        print(f"LCOE: {result.total_lcoe:.2f} 元/MWh")
        print(f"NPV: {result.total_npv:.2f} 元")
        print(f"新能源消纳比例: {result.renewable_ratio:.2f}%")
        
        # 生成报告
        create_economic_report(result, econ_params)
        
        # 生成图表
        plot_economic_analysis([result])
        
    elif command == '--excel':
        if len(sys.argv) < 3:
            print("错误：请指定Excel文件路径")
            return
        
        excel_file = sys.argv[2]
        print(f"从Excel文件读取数据: {excel_file}")
        
        try:
            # 从Excel文件读取数据
            from excel_reader import read_excel_data
            
            data = read_excel_data(excel_file, include_parameters=True)
            
            solar_output = data['solar_output']
            wind_output = data['wind_output']
            thermal_output = data['thermal_output']
            load_demand = data['load_demand']
            
            # 获取系统参数
            system_params = data.get('system_parameters', SystemParameters())
            
            # 获取经济参数
            econ_params = data.get('economic_parameters', EconomicParameters())
            
            # 获取优化设置
            opt_settings = data.get('optimization_settings', {
                'storage_capacity_range': (0, 1000),
                'rate_range': (0.1, 2.0),
                'max_iterations': 100,
                'tolerance': 0.01
            })
            
            print(f"成功读取数据，类型：{data['data_type']}")
            print(f"光伏出力总量: {sum(solar_output):.2f} MW")
            print(f"风电出力总量: {sum(wind_output):.2f} MW")
            print(f"负荷需求总量: {sum(load_demand):.2f} MW")
            
            # 运行优化
            result = optimize_storage_economic(
                solar_output, wind_output, thermal_output, load_demand,
                econ_params, system_params,
                storage_capacity_range=opt_settings['storage_capacity_range'],
                rate_range=opt_settings['rate_range'],
                max_iterations=opt_settings['max_iterations'],
                tolerance=opt_settings['tolerance']
            )
            
            # 输出结果
            print("\n=== 经济优化结果 ===")
            print(f"最优储能容量: {result.storage_capacity:.2f} MWh")
            print(f"最优充电倍率: {result.charge_rate:.2f}")
            print(f"最优放电倍率: {result.discharge_rate:.2f}")
            print(f"总建设成本: {result.total_capex:.2f} 元")
            print(f"总运维成本: {result.total_om_cost:.2f} 元")
            print(f"总电费成本: {result.total_electricity_cost:.2f} 元")
            print(f"LCOE: {result.total_lcoe:.2f} 元/MWh")
            print(f"NPV: {result.total_npv:.2f} 元")
            print(f"总弃风弃光量: {result.total_curtailment:.2f} MWh")
            print(f"总购电量: {result.grid_purchase:.2f} MWh")
            print(f"总上网电量: {result.grid_feed_in:.2f} MWh")
            print(f"新能源消纳比例: {result.renewable_ratio:.2f}%")
            
            # 生成报告
            create_economic_report(result, econ_params)
            
            # 生成图表
            plot_economic_analysis([result])
            
        except Exception as e:
            print(f"处理Excel文件时出错：{str(e)}")
            import traceback
            traceback.print_exc()


if __name__ == "__main__":
    main()