考虑风电经济调度的储能运行优化

在“双碳”目标下，低碳排放的风力发电广泛接入电力系统，装机比例逐年显著提升，但风电具有的反调峰性加剧了电力系统的调峰压力，也不得不产生弃风，储能的存在可有效解决风电的弃风问题，实现弃风消纳，降低电力系统的生产成本以及减少碳排放。因此，在推文 经济调度问题的建模及求解——考虑储能约束 的基础上，增加对风电场弃风惩罚，实现风电高比例消纳的风火储经济调度。

Part 1、风电场功率约束介绍

考虑风电的利用率问题，其并网功率和弃风功率需要满足：

P t w + P t w c ≤ Λ t w \begin{aligned} P_{t}^{w}+P_{t}^{\mathrm{wc}} \leq \Lambda_{t}^{w} \end{aligned} Ptw+Ptwc≤Λtw

其中， Λ t w \Lambda_{t}^{w} Λtw为风电场实际可发功率（MW）， P t w c P_{t}^{\mathrm{wc}} Ptwc为弃风功率（MW）， P t w P_{t}^{w} Ptw为风力发电并网功率（MW）。在此基础上，结合火电机组发电以及储能相关约束，含风火储的经济调度模型可以表示为：

min ⁡ D V T C = ∑ g , t a g P g , t 2 + b g P g , t + c g + ∑ t V W C × P t w c D V = { P g , t , S O C t , P t d , P t c , P t w , P t w c } P g min ⁡ ≤ P g , t ≤ P g max ⁡ P g , t − P g , t − 1 ≤ R U g P g , t − 1 − P g , t ≤ R D g S O C t = S O C t − 1 + ( P t c η c − P t d / η d ) Δ t P min ⁡ c ≤ P t c ≤ P max ⁡ c P min ⁡ d ≤ P t d ≤ P max ⁡ d S O C min ⁡ ≤ S O C t ≤ S O C max ⁡ P t w + ∑ g P g , t + P t d ≥ L t − P t c P t w + P t w c ≤ Λ t w \begin{aligned} \min _{\mathrm{DV}} \mathrm{TC}=& \sum_{g, t} a_{g} P_{g, t}^{2}+b_{g} P_{g, t}+c_{g}+\sum_{t} \mathrm{VWC} \times P_{t}^{\mathrm{wc}} \\ & \mathrm{DV}=\left\{P_{g, t}, \mathrm{SOC}_{t}, P_{t}^{d}, P_{t}^{c}, P_{t}^{w}, P_{t}^{\mathrm{wc}}\right\} \\ & P_{g}^{\min } \leq P_{g, t} \leq P_{g}^{\max } \\& P_{g, t}-P_{g, t-1} \leq \mathrm{RU}_{g} \\& P_{g, t-1}-P_{g, t} \leq \mathrm{RD}_{g} \\ & \mathrm{SOC}_{t}=\mathrm{SOC}_{t-1}+\left(P_{t}^{c} \eta_{c}-P_{t}^{d} / \eta_{d}\right) \Delta_{t} \\ & P_{\min }^{c} \leq P_{t}^{c} \leq P_{\max }^{c} \\& P_{\min }^{d} \leq P_{t}^{d} \leq P_{\max }^{d} \\& \mathrm{SOC}_{\min } \leq \mathrm{SOC}_{t} \leq \mathrm{SOC}_{\max } \\ & P_{t}^{w}+\sum_{g} P_{g, t}+P_{t}^{d} \geq L_{t}-P_{t}^{c} \\ & P_{t}^{w}+P_{t}^{w c} \leq \Lambda_{t}^{w} \end{aligned} DVminTC=g,t∑agPg,t2+bgPg,t+cg+t∑VWC×PtwcDV={Pg,t,SOCt,Ptd,Ptc,Ptw,Ptwc}Pgmin≤Pg,t≤PgmaxPg,t−Pg,t−1≤RUgPg,t−1−Pg,t≤RDgSOCt=SOCt−1+(Ptcηc−Ptd/ηd)ΔtPminc≤Ptc≤PmaxcPmind≤Ptd≤PmaxdSOCmin≤SOCt≤SOCmaxPtw+g∑Pg,t+Ptd≥Lt−PtcPtw+Ptwc≤Λtw

其中， a g a_{g} ag、 b g b_{g} bg、 c g c_{g} cg表示发电机组的成本系数， P g , t P_{g,t} Pg,t表示每个发电机组的输出功率，且 P g , t P_{g,t} Pg,t满足发电上、下限约束和爬坡约束； L t L_{t} Lt表示负载； P t c P_{t}^c Ptc表示储能系统的充电功率； P t d P_t^d Ptd表示储能系统的放电功率； S O C t SOC_t SOCt为储能系统的荷电状态； V W C VWC VWC为风电场弃风惩罚价格。

Part 2、Python+Gurobi 代码实现

# 导入包
import gurobipy as grb
import matplotlib.pyplot as plt
import numpy as np# 建立模型
M = grb.Model()a = np.array([0.12, 0.17, 0.15, 0.19])
b = np.array([14.8, 16.57, 15.55, 16.21])
c = np.array([89, 83, 100, 70])eta_c = 0.95
eta_d = 0.9
SOC_ini = 100
SOC_max = 300
SOC_min = 60N_g = 4
N_t = 24
VWC = 50
Pg_min = np.array([28, 20, 30, 20])
Pg_max = np.array([200, 290, 190, 260])
RU = np.array([40, 30, 30, 50])
RD = np.array([40, 30, 30, 50])Load = np.array([510, 530, 516, 510, 515, 544, 646, 686, 741, 734, 748, 760, 754, 700, 686, 720, 714, 761, 727, 714, 618, 584, 578, 544])
W = np.array([44.1, 48.5, 65.7, 144.9, 202.3, 317.3, 364.4, 317.3, 271, 306.9, 424.1, 398, 487.6, 521.9, 541.3, 560, 486.8, 372.6, 367.4, 314.3, 316.6, 311.4, 405.4, 470.4])# 声明变量
PG = M.addVars(N_g, N_t, vtype = grb.GRB.CONTINUOUS, lb=0, ub=300, name='PG')
SOC = M.addVars(N_t, vtype = grb.GRB.CONTINUOUS, lb=0.2*SOC_max, ub=SOC_max, name='SOC')
Pc = M.addVars(N_t, vtype = grb.GRB.CONTINUOUS, lb=0, ub=0.2*SOC_max, name='Pc')
Pd = M.addVars(N_t, vtype = grb.GRB.CONTINUOUS, lb=0, ub=0.2*SOC_max, name='Pd')
PW = M.addVars(N_t, vtype = grb.GRB.CONTINUOUS, lb=0, ub=W, name='PW')
PWc = M.addVars(N_t, vtype = grb.GRB.CONTINUOUS, lb=0, ub=W, name='PWc')# 给出目标函数
M.setObjective(sum(a[g] * PG[g, t] * PG[g, t] + b[g] * PG[g, t] + c[g]for t in range(N_t) for g in range(N_g)) \+ VWC * sum(PWc[t] for t in range(N_t)), grb.GRB.MINIMIZE)# 添加约束
M.addConstrs((PG[g, t] <= Pg_max[g] for t in range(N_t) for g in range(N_g)), 'Pg_max')
M.addConstrs((PG[g, t] >= Pg_min[g] for t in range(N_t) for g in range(N_g)), 'Pg_min')
M.addConstrs((PG[g, t+1]-PG[g, t] <= RU[g] for t in range(N_t-1) for g in range(N_g)), 'RU')
M.addConstrs((PG[g, t-1]-PG[g, t] <= RD[g] for t in range(1,N_t) for g in range(N_g)), 'RD')
# M.addConstrs((PG[g, t]-PG[g, t+1] <= RD[g] for t in range(N_t-1) for g in range(N_g)), 'RD')M.addConstr(SOC[0] == SOC_ini+(Pc[0]*eta_c-Pd[0]/eta_d), name='SOC_1')
M.addConstrs((SOC[t] == SOC[t-1]+(Pc[t]*eta_c-Pd[t]/eta_d) for t in range(1, N_t)), 'SOC')
M.addConstr(SOC[23] == SOC_ini, name='SOC_24')M.addConstrs(((PW[t] + sum(PG[g, t] for g in range(N_g)) + Pd[t] >= Load[t] + Pc[t]) for t in range(N_t)), 'Load')
M.addConstrs((PW[t] + PWc[t] == W[t] for t in range(N_t)), 'Wind')# 求解目标函数
M.optimize()# 输出结果
print('Obj=', M.objVal)
M.write('M.lp')demand = np.zeros((25))
Soc = np.zeros((25))
pw = np.zeros((25))
T_PG = np.zeros((25))
Pec = np.zeros((25))
Ped = np.zeros((25))
for t in range(1,25):demand[t] = Load[t-1]Soc[t] = SOC[t - 1].Xpw[t] = PW[t - 1].XT_PG[t] = sum(PG[g, t - 1].X for g in range(N_g))Pec[t] = Pc[t - 1].XPed[t] = Pd[t - 1].X# 绘图
plt.figure(1)
plt.plot(np.arange(1, 25), demand[1:], color='red', label='Demand')
plt.plot(np.arange(1, 25), T_PG[1:], color='black', label='Total power by thermal units')
plt.plot(np.arange(1, 25), Soc[1:], color='pink', label='SOC')
plt.plot(np.arange(1, 25), pw[1:], color='green', label='PW')
plt.legend()
plt.xlim(0, 24 + 1)
plt.ylim(0, 800)
plt.show()

运行结果展示：
Obj= 223360.06

注意事项：

SOC(t=1)时刻表达式的书写
爬坡约束中 t 的区间
绘图中，x轴为25个点，以方便展示

Part 3、创作人员名单

作者| 李雄，李同学
代码| 李雄，李同学
审核| 李同学，脱同学

欢迎大家关注公众号“菜鸟的能源优化之路”获取最新资讯

参考文献

A. Soroudi, Power System Optimization Modeling in GAMS. 2017. DOI: 10.1007/978-3-319-62350-4.

往期内容回顾

15. 列与约束生成(Column and Constraint Generation, C&CG/CCG)算法
14. 考虑负荷灵敏度分析的火电机组经济调度
13. 拉格朗日乘子法与KKT条件
12. 经济调度问题的建模及求解——考虑爬坡约束
11. 火电机组经济调度建模及求解——基础篇
10. 基于Distflow的最优潮流模型(OPF)–模型推导篇
9. 火电机组经济调度建模及求解——基础篇
8. Python|Gurobi——零基础学优化建模-终章
7. Python|Gurobi——零基础学优化建模-压轴篇：多目标优化
6. Python|Gurobi——零基础学优化建模-分段模型线性化(PWL)
5. Python|Gurobi——零基础学优化建模-QCP
4. Python|Gurobi——零基础学优化建模-NLP
3. Python|Gurobi——零基础学优化建模-MIP
2. Python|Gurobi——零基础学优化建模-LP
1. Python|Gurobi——零基础学优化建模