如何使用workbuddy实现能源站仿真的运行优化，实现运行费用最优（纯代码小白开始）

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import linprog

# ========== 1. 基础参数 ==========
H = 24
dt = 1.0                       # 小时
COP = 4.0                      # 制冷机能效比
# 上海商业分时电价(元/kWh)：谷(0.30) / 平(0.80) / 峰(1.30)
price = np.array([
    0.30,0.30,0.30,0.30,0.30,0.30,   # 0-6  谷
    0.80,0.80,                        # 6-8  平
    1.30,1.30,1.30,                   # 8-11 峰
    0.80,0.80,0.80,0.80,0.80,0.80,0.80, # 11-18 平
    1.30,1.30,1.30,                   # 18-21 峰
    0.80,                             # 21-22 平
    0.30,0.30                         # 22-24 谷
])
# 设备约束(统一用冷量 kW 表示)
chiller_cap       = 1200.0   # 制冷机最大冷量出力
ice_charge_max    = 800.0    # 最大蓄冰速率
ice_discharge_max = 1000.0   # 最大融冰速率
S_max             = 4000.0   # 蓄冰槽容量

# ========== 2. 模拟 24h 冷负荷(随机波动) ==========
rng = np.random.default_rng(42)
base = 300 + 900 * np.exp(-((np.arange(H)-14)**2)/(2*5**2))  # 午后峰值
noise = rng.normal(0, 60, H)
load = np.clip(base + noise, 100, None)
load = np.round(load, 1)

# ========== 3. 建模：最小化全天购电成本 ==========
# 变量顺序: q(24)制冷机冷量, c(24)蓄冰, d(24)融冰, s(24)蓄冰槽存量
n = 4 * H
def idx(k, h): return k * H + h

c_obj = np.zeros(n)
for h in range(H):
    c_obj[idx(0, h)] = price[h] / COP          # 电费 = 电价 × 电功率(=冷量/COP)

A_eq, b_eq = [], []
# 能量平衡: q + d - c = load
for h in range(H):
    row = np.zeros(n)
    row[idx(0, h)] = 1; row[idx(2, h)] = 1; row[idx(1, h)] = -1
    A_eq.append(row); b_eq.append(load[h])
# 蓄冰动态: s[h] = s[h-1] + c[h] - d[h]
for h in range(H):
    row = np.zeros(n)
    row[idx(3, h)] = 1
    if h > 0:
        row[idx(3, h-1)] = -1
    row[idx(1, h)] = -1
    row[idx(2, h)] = 1
    A_eq.append(row); b_eq.append(0.0)
# 初始蓄冰槽为空
row0 = np.zeros(n); row0[idx(3, 0)] = 1
A_eq.append(row0); b_eq.append(0.0)

bounds = []
for k in range(4):
    for _ in range(H):
        if k == 0: lo, hi = 0, chiller_cap
        elif k == 1: lo, hi = 0, ice_charge_max
        elif k == 2: lo, hi = 0, ice_discharge_max
        else:        lo, hi = 0, S_max
        bounds.append((lo, hi))

res = linprog(c_obj, A_eq=A_eq, b_eq=b_eq, bounds=bounds, method='highs')
assert res.success, res.message

q = res.x[idx(0, 0):idx(0, 24)]   # 制冷机冷量出力
c_ = res.x[idx(1, 0):idx(1, 24)]  # 蓄冰
d = res.x[idx(2, 0):idx(2, 24)]   # 融冰
s = res.x[idx(3, 0):idx(3, 24)]   # 蓄冰槽存量
elec = q / COP                     # 对应电功率
total_cost = float(np.sum(elec * price * dt))
peak_cost  = float(np.sum(elec[price == price.max()] * price.max()))

print(f"全天最低购电成本: {total_cost:,.0f} 元")
print(f"峰段电费占比: {peak_cost/total_cost*100:.1f}%")
print(f"蓄冰槽峰值存量: {s.max():.0f} / {S_max:.0f}")

# ========== 4. 可视化 ==========
hours = np.arange(H)
fig, ax = plt.subplots(3, 1, figsize=(10, 9), sharex=True)
ax[0].plot(hours, load, 'k--', label='冷负荷')
ax[0].plot(hours, q, 'b', label='制冷机出力')
ax[0].plot(hours, d, 'g', label='融冰供冷')
ax[0].plot(hours, c_, 'orange', label='蓄冰')
ax[0].legend(); ax[0].set_ylabel('冷量 kW'); ax[0].set_title('制冷机出力与蓄冰动作')

ax[1].plot(hours, s, 'purple')
ax[1].fill_between(hours, s, alpha=0.3, color='purple')
ax[1].set_ylabel('蓄冰量 kWh冷'); ax[1].set_title('蓄冰槽存量变化')

cost_h = elec * price
colors = ['green' if p == price.min() else ('red' if p == price.max() else 'gray') for p in price]
ax[2].bar(hours, cost_h, color=colors)
ax[2].set_ylabel('电费 元'); ax[2].set_xlabel('小时')
ax[2].set_title('逐时电费(绿=谷 红=峰 灰=平)')
plt.tight_layout(); plt.show()