"""
Home Energy Optimization — Simplified LP Model (CVXPY)
========================================================
基于题设模型直接实现，无效率损耗，无自放电。

目标函数：
    min sum_t(P_buy,t * pi_buy,t - P_sell,t * pi_sell,t + C_cycle * P_discharge,t) * dt

约束：
    功率平衡 / SOC 动态 / SOC 边界 / 初始 SOC / 逆变器限制 / 并网限制
"""

# pip install numpy
# pip install matplotlib
# pip install cvxpy

import cvxpy as cp
import numpy as np
import matplotlib.pyplot as plt

# ═══════════════════════════════════════════════════════════
# 1. 参数配置（替换为你的实际数据）
# ═══════════════════════════════════════════════════════════

T = 24          # 优化时域总时段数
dt = 1.0        # 时间步长 [h]（1 = 1小时，0.25 = 15分钟）

# 预测数据（示例值，替换为你的预测模型输出）
P_pv     = np.array([0, 0, 0, 0, 0, 0,              0.2, 0.8, 2.1, 3.5, 4.2, 4.8,
                     4.5, 3.8, 2.5, 1.2, 0.3, 0,    0, 0, 0, 0, 0, 0])  # kW
P_rigid  = np.array([0.4, 0.3, 0.3, 0.3, 0.4, 0.6,  0.8, 1.0, 0.9, 0.7, 0.6, 0.7,
                     0.8, 0.9, 0.8, 0.9, 1.2, 1.8,  2.5, 2.8, 2.2, 1.5, 0.9, 0.6])  # kW (18-22点为高峰期)

# 分时电价 [元/kWh]（示例广东电网分时电价 TOU（Time of Use））
pi_buy  = np.array([0.3, 0.3, 0.3, 0.3, 0.3, 0.3,  0.3, 0.5, 0.8, 0.8, 0.8, 0.8,
                    0.8, 0.8, 0.5, 0.5, 0.5, 0.8,  0.8, 0.5, 0.3, 0.3, 0.3, 0.3])
pi_sell = pi_buy * 0.4   # 上网电价 ≈ 买价的 40%

# 电池参数
C_batt    = 13.5    # 电池总容量 [kWh]
SOC_min   = 0.10    # 最低 SOC（10% DOD）
SOC_max   = 0.90    # 最高 SOC（90% DOD）
SOC_init  = 0.50    # 初始 SOC [比例]，来自 BMS 实时读取

# 功率限制
P_inv_max  = 5.0    # 逆变器功率上限 [kW]
P_grid_max = 5.0    # 并网功率上限 [kW]

# 经济参数
C_cycle = 0.15      # 电池循环成本 [元/kWh]，仅惩罚放电

# ═══════════════════════════════════════════════════════════
# 2. 决策变量
# ═══════════════════════════════════════════════════════════

P_buy       = cp.Variable(T, nonneg=True)       # 从电网购电 [kW]
P_sell      = cp.Variable(T, nonneg=True)       # 向电网售电 [kW]
P_charge    = cp.Variable(T, nonneg=True)       # 电池充电功率 [kW]
P_discharge = cp.Variable(T, nonneg=True)       # 电池放电功率 [kW]
SOC         = cp.Variable(T + 1)                # SOC[0] = 初始, SOC[1..T] = 优化结果

# ═══════════════════════════════════════════════════════════
# 3. 目标函数
# ═══════════════════════════════════════════════════════════

cost = (
    cp.sum(P_buy @ pi_buy)               # 购电成本
    - cp.sum(P_sell @ pi_sell)           # 售电收益
    + C_cycle * cp.sum(P_discharge)      # 电池循环成本（仅放电）
) * dt

objective = cp.Minimize(cost)

# ═══════════════════════════════════════════════════════════
# 4. 约束条件
# ═══════════════════════════════════════════════════════════

constraints = []

# ── 4.1 功率平衡 ───────────────────────────────────────────
for t in range(T):
    # 功率平衡：光伏 + 购电 + 放电 = 刚性负荷 + 充电 + 售电
    constraints.append(
        P_pv[t] + P_buy[t] + P_discharge[t]
        == P_rigid[t] + P_charge[t] + P_sell[t]
    )

    # SOC 动态：SOC_t = SOC_{t-1} + (P_charge - P_discharge) * dt
    constraints.append(
        SOC[t + 1] == SOC[t] + (P_charge[t] - P_discharge[t]) * dt
    )

# ── 4.2 初始 SOC ──────────────────────────────────────────
constraints.append(SOC[0] == SOC_init * C_batt)

# ── 4.3 SOC 动态 ──────────────────────────────────────────
for t in range(T):
    # SOC 动态：SOC_t = SOC_{t-1} + (P_charge - P_discharge) * dt
    constraints.append(
        SOC[t + 1] == SOC[t] + (P_charge[t] - P_discharge[t]) * dt
    )

# ── 4.4 SOC 边界 ──────────────────────────────────────────
constraints.append(SOC >= SOC_min * C_batt)
constraints.append(SOC <= SOC_max * C_batt)

# ── 4.5 逆变器功率限制 ────────────────────────────────────
constraints.append(P_charge   <= P_inv_max)
constraints.append(P_discharge <= P_inv_max)

# ── 4.6 并网功率限制 ──────────────────────────────────────
constraints.append(P_buy  <= P_grid_max)
constraints.append(P_sell <= P_grid_max)

# ═══════════════════════════════════════════════════════════
# 5. 求解
# ═══════════════════════════════════════════════════════════

problem = cp.Problem(objective, constraints)
result = problem.solve(solver=cp.SCS, verbose=True)

# ═══════════════════════════════════════════════════════════
# 6. 结果输出
# ═══════════════════════════════════════════════════════════

if problem.status in ["optimal", "optimal_inaccurate"]:
    print(f"\n求解状态：{problem.status}")
    print(f"最优电费：¥{problem.value:.2f}")
    print("+" + "-" * 87 + "+")
    print(f"| {'时':>3} | {'PV':>6} | {'负荷':>4} | {'【购电】':>4} | {'【售电】':>4} | {'【充电】':>4} | {'【放电】':>4} | {'SOC':>6} | {'成本':>5} |")
    print("+" + "-" * 87 + "+")

    total_cost = 0.0
    for t in range(T):
        soc_pct = SOC[t + 1].value / C_batt * 100
        step_cost = (
            P_buy.value[t] * pi_buy[t]
            - P_sell.value[t] * pi_sell[t]
            + C_cycle * P_discharge.value[t]
        ) * dt
        total_cost += step_cost

        print(
            f"| "
            f"{t:>3}  | "
            f"{P_pv[t]:6.2f} | "
            f"{P_rigid[t]:6.2f} | "
            f"{P_buy.value[t]:6.2f}   | "
            f"{P_sell.value[t]:6.2f}   | "
            f"{P_charge.value[t]:6.2f}   | "
            f"{P_discharge.value[t]:6.2f}   | "
            f"{soc_pct:5.1f}% | "
            f"¥{step_cost:6.2f} |"
        )

    print("+" + "-" * 87 + "+")
    print(f"| {'合计':>2} | {'':>6} | {'':>6} | {'':>8} | {'':>8} | {'':>8} | {'':>8} | {'':>6} | {'¥':>1} {total_cost:.2f} |")
    print(f"\n初始 SOC：{SOC_init*100:.0f}% → 最终 SOC：{SOC[T].value / C_batt * 100:.1f}%")

    # ── plot ──────────────────────────────────────────────
    P_buy_val       = P_buy.value
    P_sell_val      = P_sell.value
    P_charge_val    = P_charge.value
    P_discharge_val = P_discharge.value
    SOC_val         = SOC.value

    hours = np.arange(T) * dt

    fig, axes = plt.subplots(3, 1, figsize=(14, 10), sharex=False)
    plt.rcParams['font.size'] = 10
    plt.rcParams['font.sans-serif'] = ['Noto Sans CJK JP', 'Noto Sans CJK SC', 'Noto Sans CJK TC', 'Noto Sans CJK KR', 'DejaVu Sans']
    plt.rcParams['axes.unicode_minus'] = False

    # ========== 图1: 功率平衡（面积图，正负分列）==========
    ax1 = axes[0]

    # 正向：发电/买电/放电
    ax1.fill_between(hours, 0, P_pv, label=r'$P_{pv}$ 光伏', alpha=0.75, color='#2ecc71', linewidth=0)
    ax1.fill_between(hours, P_pv, P_pv+P_buy_val, label=r'$P_{buy}$ 购电', alpha=0.75, color='#3498db', linewidth=0)
    ax1.fill_between(hours, P_pv+P_buy_val, P_pv+P_buy_val+P_discharge_val, label=r'$P_{discharge}$ 放电', alpha=0.75, color='#e67e22', linewidth=0)

    # 负向：售电/充电（向下）
    ax1.fill_between(hours, 0, -P_rigid, label=r'$P_{rigid}$ 用电', alpha=0.75, color='#e74c3c', linewidth=0)
    ax1.fill_between(hours, -P_rigid, -P_rigid-P_sell_val, label=r'$P_{sell}$ 售电', alpha=0.75, color='#bdc3c7', linewidth=0)
    ax1.fill_between(hours, -P_rigid-P_sell_val, -P_rigid-P_sell_val-P_charge_val, label=r'$P_{charge}$ 充电', alpha=0.75, color='#9b59b6', linewidth=0)

    ax1.set_ylabel('功率 (kW)')
    ax1.set_xlabel('时间 (h)')
    ax1.set_title('功率平衡')
    ax1.legend(loc='upper right', ncol=6, fontsize=10, framealpha=0.9)
    ax1.grid(True, alpha=0.3, linestyle='--')
    ax1.axhline(0, color='black', linewidth=0.8)
    ax1.set_ylim(-6, 10)
    ax1.set_xticks(range(T))

    # ========== 图2: 电池SOC ==========
    ax2 = axes[1]

    ax2.bar(hours, SOC_val[1:] / C_batt * 100, width=dt * 0.8, label=r'电池 $SOC$', color='#3498db', alpha=0.85, edgecolor='white')
    ax2.axhline(SOC_min * 100, color='red',    linestyle='--', linewidth=1, label=f'SOC_min={SOC_min*100:.0f}%')
    ax2.axhline(SOC_max * 100, color='green',  linestyle='--', linewidth=1, label=f'SOC_max={SOC_max*100:.0f}%')
    ax2.set_ylabel('SOC (%)')
    ax2.set_xlabel('时间 (h)')
    ax2.set_title('电池荷电状态')
    ax2.legend(loc='upper right', fontsize=10)
    ax2.grid(True, alpha=0.3, linestyle='--', axis='y')
    ax2.set_ylim(0, 100)
    ax2.set_xticks(range(T))

    # ========== 图3: 分时电价 ==========
    ax3 = axes[2]

    ax3.step(hours, pi_buy,  'r-', where='post', linewidth=1.5, label=r'$\pi_{buy}$ 买电价格')
    ax3.step(hours, pi_sell, 'g-', where='post', linewidth=1.5, label=r'$\pi_{sell}$ 卖电价格')

    ax3.set_ylabel('电价 (元/kWh)')
    ax3.set_xlabel('时间 (h)')
    ax3.set_title('分时电价')
    ax3.legend(loc='upper right', fontsize=10)
    ax3.grid(True, alpha=0.3, linestyle='--')
    ax3.axhline(0, color='black', linewidth=0.8)
    ax3.set_ylim(0, max(pi_buy) * 1.3)
    ax3.set_xticks(range(T))

    plt.tight_layout()
    plt.savefig('home_energy_optimization.png', dpi=150, bbox_inches='tight')
    plt.show()
    print("图片已保存为 home_energy_optimization.png")

else:
    print(f"求解失败：{problem.status}")
    if problem.status == "infeasible":
        print("→ 约束无解")