import numpy as np import matplotlib.pyplot as plt # 设置中文字体 plt.rcParams['font.sans-serif'] = ['SimHei'] plt.rcParams['axes.unicode_minus'] = False class GreenAmmoniaAnalyzer: """绿电制氨系统分析器""" def __init__(self): """初始化系统参数""" # 系统参数 self.P_fixed = 20.75 # 固定负荷 (MW) self.P_base_con = 6 # 常规负荷基准 (MW) self.P_base_wind = 40 # 风电基准 (MW) self.P_base_pv = 64 # 光伏基准 (MW) self.M_daily = 36 # 日产量 (吨) # 标幺值曲线 (24小时) self.p_con_pu = np.array([ 0.2283, 0.1833, 0.1333, 0.1650, 0.1967, 0.2717, 0.4167, 0.5617, 0.6683, 0.7000, 0.6867, 0.6500, 0.5433, 0.4483, 0.6050, 0.6367, 0.6250, 0.5233, 0.4350, 0.3533, 0.3100, 0.2717, 0.2600, 0.2467 ]) self.p_wind_pu = np.array([ 0.2726, 0.3788, 0.3328, 0.3768, 0.4681, 0.2407, 0.5701, 0.4463, 0.3582, 0.2976, 0.1934, 0.2139, 0.3630, 0.1677, 0.2794, 0.2527, 0.1596, 0.1607, 0.1434, 0.0938, 0.0923, 0.1220, 0.0625, 0.0798 ]) self.p_pv_pu = np.array([ 0.00, 0.00, 0.00, 0.00, 0.00, 0.02, 0.12, 0.28, 0.45, 0.58, 0.67, 0.70, 0.72, 0.67, 0.58, 0.45, 0.28, 0.08, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00 ]) # 分时电价 (元/kWh) self.price_buy = np.array([0.3424]*7 + [0.6074]*3 + [0.8024]*5 + [0.6074]*3 + [0.8024]*3 + [0.6074]*2 + [0.3424]) # 计算结果存储 self.results = {} def calculate_power_balance(self): """计算功率平衡""" # 计算实际功率 P_con = self.P_base_con * self.p_con_pu P_wind = self.P_base_wind * self.p_wind_pu P_pv = self.P_base_pv * self.p_pv_pu P_ren = P_wind + P_pv P_load = P_con + self.P_fixed delta = P_load - P_ren P_buy = np.maximum(delta, 0) P_sell = np.maximum(-delta, 0) # 存储结果 self.results.update({ 'P_con': P_con, 'P_wind': P_wind, 'P_pv': P_pv, 'P_ren': P_ren, 'P_load': P_load, 'P_buy': P_buy, 'P_sell': P_sell }) return P_con, P_wind, P_pv, P_load, P_buy, P_sell def calculate_energy_metrics(self): """计算电量指标""" P_buy = self.results['P_buy'] P_sell = self.results['P_sell'] P_load = self.results['P_load'] P_ren = self.results['P_ren'] P_wind = self.results['P_wind'] P_pv = self.results['P_pv'] # 累计电量 E_load = np.sum(P_load) E_ren = np.sum(P_ren) E_wind_total = np.sum(P_wind) E_pv_total = np.sum(P_pv) E_buy = np.sum(P_buy) E_sell = np.sum(P_sell) # 计算指标 lambda_self = (E_load - E_buy - E_sell) / E_ren * 100 lambda_green = (E_ren - E_sell) / E_load * 100 lambda_grid = E_sell / E_ren * 100 lambda_self_intuitive = (E_ren - E_sell) / E_ren * 100 self.results.update({ 'E_load': E_load, 'E_ren': E_ren, 'E_wind_total': E_wind_total, 'E_pv_total': E_pv_total, 'E_buy': E_buy, 'E_sell': E_sell, 'lambda_self': lambda_self, 'lambda_green': lambda_green, 'lambda_grid': lambda_grid, 'lambda_self_intuitive': lambda_self_intuitive }) return E_load, E_ren, E_buy, E_sell def calculate_cost(self): """计算运行成本""" P_buy = self.results['P_buy'] E_wind_total = self.results['E_wind_total'] E_pv_total = self.results['E_pv_total'] E_sell = self.results['E_sell'] # 成本计算 cost_wind = E_wind_total * 1000 * 0.15 # 风电成本 0.15元/kWh cost_pv = E_pv_total * 1000 * 0.12 # 光伏成本 0.12元/kWh cost_buy = np.sum(P_buy * self.price_buy * 1000) # 购电成本 income_sell = E_sell * 1000 * 0.3779 # 售电收入 om_total = 24000 + 36000 + 36 # 运维总成本 total_cost = cost_wind + cost_pv + cost_buy - income_sell + om_total c_nh3 = total_cost / self.M_daily self.results.update({ 'cost_wind': cost_wind, 'cost_pv': cost_pv, 'cost_buy': cost_buy, 'income_sell': income_sell, 'om_total': om_total, 'total_cost': total_cost, 'c_nh3': c_nh3 }) return total_cost, c_nh3 def print_power_balance(self): """打印功率平衡表""" P_con = self.results['P_con'] P_wind = self.results['P_wind'] P_pv = self.results['P_pv'] P_load = self.results['P_load'] P_buy = self.results['P_buy'] P_sell = self.results['P_sell'] print("\n" + "="*70) print("表7.1 典型日逐时功率平衡（MW）") print("-"*70) print(f"{'时刻':<6} {'常规负荷':>10} {'风电':>10} {'光伏':>10} " f"{'总负荷':>10} {'购电':>10} {'售电':>10}") print("-"*70) for t in range(24): print(f"{t:02d}:00 {P_con[t]:10.2f} {P_wind[t]:10.2f} " f"{P_pv[t]:10.2f} {P_load[t]:10.2f} " f"{P_buy[t]:10.2f} {P_sell[t]:10.2f}") print("="*70) def print_energy_metrics(self): """打印电量指标""" print("\n表7.2 典型日电量") print("-"*50) print(f"总用电量 E_load: {self.results['E_load']:.2f} MWh") print(f"新能源发电量 E_ren: {self.results['E_ren']:.2f} MWh " f"(风 {self.results['E_wind_total']:.2f}, 光 {self.results['E_pv_total']:.2f})") print(f"网购电量 E_buy: {self.results['E_buy']:.2f} MWh") print(f"上网电量 E_sell: {self.results['E_sell']:.2f} MWh") print("-"*50) def print_compliance(self): """打印达标判定""" print("\n表7.3 绿电直连指标及达标判定") print("-"*50) metrics = [ ('λ_self', self.results['lambda_self'], 60, '>'), ('λ_green', self.results['lambda_green'], 30, '>'), ('λ_grid', self.results['lambda_grid'], 20, '<') ] for name, value, threshold, op in metrics: is_compliant = (value > threshold) if op == '>' else (value < threshold) status = "✓ 符合" if is_compliant else "✗ 不符合" print(f"{name}: {value:6.2f}% (阈值{op}{threshold}%) {status}") print(f"\n提示：λ_self直观值 = {self.results['lambda_self_intuitive']:.2f}%") print("-"*50) def print_cost_breakdown(self): """打印成本构成""" print("\n表7.4 吨氨运行成本构成") print("-"*50) total = self.results['total_cost'] wind_pv_cost = self.results['cost_wind'] + self.results['cost_pv'] items = [ ('风光发电成本', wind_pv_cost), (' 风电', self.results['cost_wind']), (' 光伏', self.results['cost_pv']), ('网购电成本', self.results['cost_buy']), ('售电收入', -self.results['income_sell']), ('电解槽运维', 24000 + 36000), ('合成氨运维', 36) ] for name, value in items: if name.startswith(' '): # 子项缩进显示 print(f"{name}: {value:10.0f} 元") else: pct = value / total * 100 if value > 0 else -value / total * 100 print(f"{name}: {value:10.0f} 元 ({pct:.1f}%)") print("-"*50) print(f"总运行成本: {self.results['total_cost']:10.0f} 元") print(f"吨氨成本: {self.results['c_nh3']:10.0f} 元/吨") print("="*50) def plot_power_balance(self): """绘制功率平衡曲线""" fig, ax1 = plt.subplots(figsize=(12, 5)) # 负荷侧 ax1.stackplot(range(24), self.results['P_con'], [self.P_fixed] * 24, labels=['常规负荷', '制氢制氨固定负荷'], alpha=0.8) ax1.plot(range(24), self.results['P_load'], 'k-', lw=2, label='总负荷') ax1.set_ylabel('功率 (MW)', fontsize=12) ax1.set_xlabel('时刻 (h)', fontsize=12) # 电源侧 ax2 = ax1.twinx() ax2.stackplot(range(24), self.results['P_wind'], self.results['P_pv'], labels=['风电', '光伏'], alpha=0.8) ax2.plot(range(24), self.results['P_ren'], 'b--', lw=2, label='可再生总出力') ax2.set_ylabel('功率 (MW)', fontsize=12) # 合并图例 lines1, labels1 = ax1.get_legend_handles_labels() lines2, labels2 = ax2.get_legend_handles_labels() ax1.legend(lines1 + lines2, labels1 + labels2, loc='upper left', fontsize=10) plt.title('典型日功率平衡曲线', fontsize=14, fontweight='bold') plt.tight_layout() plt.show() def plot_purchase_sale(self): """绘制购售电曲线""" fig, ax = plt.subplots(figsize=(12, 4)) ax.bar(range(24), self.results['P_sell'], color='#2E8B57', label='售电（上网）', alpha=0.7) ax.bar(range(24), -self.results['P_buy'], color='#CD5C5C', label='购电（网购）', alpha=0.7) ax.axhline(0, color='black', linewidth=0.8) ax.set_xlabel('时刻 (h)', fontsize=12) ax.set_ylabel('功率 (MW)', fontsize=12) ax.legend(loc='upper right', fontsize=10) ax.grid(True, alpha=0.3) plt.title('典型日购售电功率曲线', fontsize=14, fontweight='bold') plt.tight_layout() plt.show() def plot_compliance(self): """绘制达标情况图""" fig, ax = plt.subplots(figsize=(8, 5)) metrics = ['λ_self', 'λ_green', 'λ_grid'] values = [self.results['lambda_self'], self.results['lambda_green'], self.results['lambda_grid']] thresholds = [60, 30, 20] bars = ax.bar(metrics, values, color='steelblue', alpha=0.7) for i, (bar, v, th) in enumerate(zip(bars, values, thresholds)): # 绘制阈值线 ax.plot([i - 0.3, i + 0.3], [th, th], 'r--', lw=2, label='政策阈值' if i == 0 else '') # 达标判定 is_compliant = (v > th) if i < 2 else (v < th) color = 'green' if is_compliant else 'red' ax.text(i, v + 1, '✓' if is_compliant else '✗', ha='center', fontsize=18, color=color, fontweight='bold') # 添加直觉值标注 ax.text(0, self.results['lambda_self_intuitive'] + 2, f'直觉值: {self.results["lambda_self_intuitive"]:.1f}%', ha='center', fontsize=10, color='purple', style='italic') ax.set_ylabel('比例 (%)', fontsize=12) ax.set_ylim(0, max(values + thresholds) * 1.15) ax.legend(['政策阈值'], loc='upper right', fontsize=10) ax.grid(True, alpha=0.3, axis='y') plt.title('绿电直连指标达标情况', fontsize=14, fontweight='bold') plt.tight_layout() plt.show() def plot_cost_pie(self): """绘制成本构成饼图""" fig, ax = plt.subplots(figsize=(7, 6)) wind_pv_cost = self.results['cost_wind'] + self.results['cost_pv'] labels = [f'风光发电\n{wind_pv_cost:.0f}元', f'网购电\n{self.results["cost_buy"]:.0f}元', f'运维\n{self.results["om_total"]:.0f}元'] sizes = [wind_pv_cost, self.results['cost_buy'], self.results['om_total']] colors = ['#4682B4', '#CD5C5C', '#3CB371'] explode = (0, 0.05, 0) wedges, texts, autotexts = ax.pie(sizes, explode=explode, labels=labels, colors=colors, autopct='%1.1f%%', startangle=90, textprops={'fontsize': 11}) # 美化百分比文字 for autotext in autotexts: autotext.set_color('white') autotext.set_fontweight('bold') plt.title('吨氨运行成本构成', fontsize=14, fontweight='bold', pad=20) plt.tight_layout() plt.show() def run_analysis(self): """运行完整分析""" print("\n" + "="*70) print(" "*20 + "绿电制氨系统分析报告") print("="*70) # 执行计算 self.calculate_power_balance() self.calculate_energy_metrics() self.calculate_cost() # 输出结果 self.print_power_balance() self.print_energy_metrics() self.print_compliance() self.print_cost_breakdown() # 绘制图表 print("\n正在生成图表...") self.plot_power_balance() self.plot_purchase_sale() self.plot_compliance() self.plot_cost_pie() # 主程序入口 if __name__ == "__main__": try: analyzer = GreenAmmoniaAnalyzer() analyzer.run_analysis() except Exception as e: print(f"程序运行出错: {e}") import traceback traceback.print_exc()