Electricity Market Optimization(VI) - 机组组合模型以及 Gurobi 求解

本文参考链接：link

\hspace{1.6em} 机组组合问题在电力系统中非常重要，这个问题也是一个优化问题，研究的就是如何调度现有的机组，调度的对象是以煤炭、石油、天然气为燃料的火力发电机以及水力发电机等可预测处理的发电机组，调度的过程是指派某个发电机组在未来的某个特定的时间开始启动，在预定的时段开始发电，调度的目标是在满足约束的前提下尽量的是发电成本最小化。截止24年底，中国以煤炭、石油、天然气为燃料的火力发电机的装机容量占全部发电装机容量的 47.62%，所以我们这主要讨论火力发电机的启停的情况，下面是一个记录大型火电机组启动过程的甘特图

\hspace{1.6em} 可以看到一个大型火电机组的启动过程比较复杂，这也带来了机组组合问题的核心特征：

• 发电成本函数是一个在原点不连续的函数
• 需要提前确定哪些机组需要开，以保证第二天的电能实际交割

\hspace{1.6em} 为了描述机组的开关状态，引入如下方程：
v i , t − w i , t = u i , t − u i , t − 1 ∑ s = t − ( U T i − 1 ) t v i , s ≤ u i , t ∑ s = t − ( D T i − 1 ) t w i , s ≤ 1 − u i , t u i , t ∈ { 0 , 1 } 0 ≤ v i , t , w i , t ≤ 1 ∀ i ∈ [ G ] , t ∈ [ T ] \begin{aligned} & v_{i,t} - w_{i,t} = u_{i,t} - u_{i,t-1} \\ & \sum_{s=t-(UT_i-1)}^t v_{i,s} \leq u_{i,t} \\ & \sum_{s=t-(DT_i-1)}^t w_{i,s} \leq 1-u_{i,t} \\\\ & u_{i,t}\in\{0,1\} \\\\ & 0 \le v_{i,t}, w_{i,t} \le 1 &&&&&&&&&&&&& \forall i \in [G], t \in [T] \end{aligned} ​vi,t​−wi,t​=ui,t​−ui,t−1​s=t−(UTi​−1)∑t​vi,s​≤ui,t​s=t−(DTi​−1)∑t​wi,s​≤1−ui,t​ui,t​∈{0,1}0≤vi,t​,wi,t​≤1​​​​​​​​​​​​​∀i∈[G],t∈[T]​
\hspace{1.6em} 方程中 $u_{i,t}=1$ 表示在 $t$ 时刻，发电机处于运行状态，这种运行状态表示发电机正在准备或者已经处于发电状态，
$v_{i,t}>0$ 在 表示 $t$ 时刻，发电机正在准备发电
$w_{i,t}>0$ 在 表示 $t$ 时刻，发电机正在准备停机

机组组合问题的约束

\hspace{1.6em} 机组组合问题的约束进行分类，分为三类：

• 功率平衡约束
• 火电机组启停约束和爬坡约束
• 备用容量约束

\hspace{1.6em} 我在之前的博客中对功率平衡约束和爬坡约束进行了讨论，本文就机组启停约束和备用容量约束进行讨论

火电机组启停约束和爬坡约束

启动时功率变化和停机功率变化

\hspace{1.6em} 当机组在时刻 $t$ 启动时，即 $v_{i,t}=1$ 时，发电机的功率可以从 0 直接增加到最小发电功率 $\underline{P_i}$, 同样当机组在时刻 $t$ 停止时，即 $w_{i,t}=1$ 时，发电机的功率可以从最小发电功率直接减少到 0，所以可以用下面的式子来描述这一点：
$$\begin{array}{rlrlrlrlrlrlrlrlr}& p_{i,t}^{\mathrm{c}}-p_{i,t-1}^{\mathrm{c}}\le R_i^{60}+v_{i,t}{\underline{P}}_i& & & & & & & & & & & & & & & \\ & & & & & & & & & & & & & & & & \\ & p_{i,t-1}^{\mathrm{c}}-p_{i,t}^{\mathrm{c}}\le R_i^{60}+w_{i,t}{\underline{P}}_i& & & & & & & & & & & & & & & \forall i\in [G],t\in [T]\end{array}$$
其中 $R_i^{60}$ 是机组功率变化的上限

\begin{aligned}&&&&&&&&&&&&&&&&\end{aligned} ​​​​​​​​​​​​​​​​​至此，成功往模型中加入离散变量！

备用容量约束

\hspace{1.6em} 为了避免出现发电容量短缺的情况，在考虑发电机组组合的时候需要在日前就考虑到有一些火电机组本身就要留有备用容量，这里的备用容量是只针对处于开机状态的机组而言的，即热备用。这是为了防止在电能交割日出现的发电容量无法满足需求的情况。无法满足需求的原因包括：负荷突然升高，多台发电机组同时故障。站在发电厂的角度来看，解决短缺的措施之一就是保证有足够的机组备用容量，下面是 link 中介绍的几个关于备用容量的基本约束：

• 发电机只能提供其输出功率之下的备用容量，启动状态下无法提供非旋转备用
• 只有一部分的发电机才能提供旋转备用，只有具备启动可靠性才能提供
• 总备用容量应该大于或等于最大单发电机的输出功率
• 总备用容量应该大于或等于要求的最低备用容量 $r_t^{\mathrm{req}}$.
• 一半的最低备用容量 $r_t^{\mathrm{req}}$ 应该通过旋转备用来满足
  $$\begin{array}{rlrlrlrlrlrl}& p_{i,t}+r_{i,t}^{\mathrm{s}}\le {\overline{P}}_i{u}_{i,t}& & & & & & & & & & \\ & p_{i,t}+r_{i,t}^{\mathrm{s}}+r_{i,t}^{\mathrm{ns}}\le {\overline{P}}_i& & & & & & & & & & \\ & r_{i,t}^{\mathrm{s}}=0& & \forall i\in [G]\setminus {\mathcal{G}}^{\mathrm{S}},& & & & & & & & \\ & r_t^{\mathrm{req}}\ge p_{i,t}& & & & & & & & & & \\ & \sum _{i\in [G]}(r_{i,t}^{\mathrm{s}}+r_{i,t}^{\mathrm{ns}})\ge r_t^{\mathrm{req}}& & & & & & & & & & \\ & \sum _{i\in [G]}{r}_{i,t}^{\mathrm{s}}\ge 0.5r_t^{\mathrm{req}}& & & & & & & & & & \forall i\in [G],t\in [T]\end{array}$$

模型目标函数

\hspace{1.6em} 机组组合问题的目标函数中的成本分为两部分：启停成本和发电成本

启动成本和停机成本

\hspace{1.6em} 那么机组启停成本可以用下式表示：
$$\begin{array}{rlrlrlrlrlrlrlrlrlrlrl}{C}_t^{\mathrm{uc}}=\sum _{i\in [G]}(c_i^{\mathrm{su}}{v}_{i,t}+c_i^{\mathrm{sd}}{w}_{i,t})& & & & & & & & & & & & & & & & & & & & & \forall t\in [T]\end{array}$$

发电成本

\hspace{1.6em} 我们这里考虑了新能源机组的运行成本，新能源机组的运行成本构成已经在之前的博客中有说明link

对于数据的说明

启动可靠性

\hspace{1.6em} 这是评判一个机组是否能够正常进行启动的标准，一个机组是否具有启动可靠性需要看：

• 硬件性能是否正常，比如机械结构，电气系统是否正常运行
• 软件性能是否正常，比如检测和控制系统是否能精准测量并调节设备的各种性能参数，这些参数包括：机组的转速、负荷、温度、压力等运行参数
• 机组启动时的电网运行条件是否正常，发电机启动初期需要从电网获取电能建立足够的磁场和感应电动势，此时如果电网电压或者频率不稳定，会影响发电机的正常启动

利用Gurobi + python 建立模型并求解

import pandas as pd
import numpy as np
import math
import gurobipy as gp
from gurobipy import *
import matplotlib.pyplot as pltgen_colors = {"Hydro": "#1f78b4",  # Deep Blue"Nuclear": "#e31a1c",  # Red"Coal": "#8b4513",  # Dark Brown"Gas CC": "#8e44ad",  # Medium Purple"Gas CT": "#a569bd",  # Light Purple"Oil CT": "#4d4d4d",  # Dark Gray"Wind": "#6baed6",  # Light Sky Blue"PV": "#ff7f00",  # Bright Orange"Storage": "#33a02c" # Green
}def print_gp_status(m):status = m.Statusif status == GRB.OPTIMAL:print("The model is optimal.")elif status == GRB.INFEASIBLE:print("The model is infeasible.")elif status == GRB.UNBOUNDED:print("The model is unbounded.")else:print(f"Optimization ended with status {status}.")print()return status# 读入数据
data_file = "ts_and_gen_scuc.xlsx"
load_and_res_data = pd.read_excel(data_file, sheet_name=0)
gen_data = pd.read_excel(data_file, sheet_name=1)# 不同发电机的类型、编号、可变成本、容量、爬坡速率
gen_id = gen_data['Gen ID'].to_numpy()
gen_type = gen_data['Unit Type'].to_numpy()
mc = gen_data['Variable Cost USD/MWh'].to_numpy()
Pmax = gen_data['Installed in MW'].to_numpy()
R60 = gen_data['Ramp Rate in MW /hr'].to_numpy()# 不同发电机的固定成本、最小发电功率、最小启停时间、启动可靠性状态变量
su_cost = gen_data['SU cost $'].to_numpy()
Pmin = gen_data['PMin in MW'].to_numpy()
min_uptime = gen_data['Min Up Time Hr'].to_numpy()
min_downtime = gen_data['Min Down Time Hr'].to_numpy()
sr_elig = gen_data['SR eligible'].to_numpy()
n_gen = len(gen_id)# 负荷数据
load = load_and_res_data['Load MW'].to_numpy()
wind = load_and_res_data['Wind MW'].to_numpy()
pv = load_and_res_data['PV MW'].to_numpy()
n_t = len(load)# 建立模型# 弃风弃光的成本系数
wind_curt_pen = 0 # $/MWh
pv_curt_pen = 0 # $/MWhm = gp.Model()
m.setParam("OutputFlag", 0)# 向模型中添加变量，这些变量包括：
# 发电机的发电功率、运行状态、启动状态、停机状态、风力发电机的发电功率
# 光伏发电机的发电功率、对于电机的备用容量要求、旋转备用容量、非旋转备用容量
p = m.addVars(n_gen, n_t, lb=0, ub=GRB.INFINITY, name="p")
u = m.addVars(n_gen, n_t, vtype=GRB.BINARY, name="u") 
v = m.addVars(n_gen, n_t, lb=0, ub=1, name="v") 
w = m.addVars(n_gen, n_t, lb=0, ub=1, name="w")
p_wind = m.addVars(n_t, lb=0, ub=GRB.INFINITY, name="p_wind")
p_pv = m.addVars(n_t, lb=0, ub=GRB.INFINITY, name="p_pv")
r_req = m.addVars(n_t, lb=0, ub=GRB.INFINITY, name="r_req")
p_rs = m.addVars(n_gen, n_t, lb=0, ub=GRB.INFINITY, name="p_rs") 
p_rns = m.addVars(n_gen, n_t, lb=0, ub=GRB.INFINITY, name="p_rns") # 向模型中添加约束，这些约束包括：
pkey, pval = multidict(p)
# 1. 供需平衡约束
m.addConstrs((gp.quicksum(p[key] for key in pkey.select('*', t)) + p_wind[t] + p_pv[t] == load[t] for t in range(n_t)), name="energy_balance")for t in range(n_t):for i in range(n_gen):# 2. (传统能源)发电机容量约束m.addConstr(p[i,t] + p_rs[i,t] <= Pmax[i]*u[i,t]) # spinning reserve can be provided only from "on" generatorsm.addConstr(p[i,t] + p_rs[i,t] + p_rns[i,t] <= Pmax[i]) # non-spiining reserve can be provided also from "off" generatorsm.addConstr(p[i,t] >= Pmin[i]*u[i,t])# 3. 只有满足启动可靠性的那些(传统能源)发电机可以提供旋转容量if sr_elig[i]==0:m.addConstr(p_rs[i,t] == 0)# 4. 爬坡约束if t > 0:m.addConstr(p[i,t] - p[i,t-1] <= R60[i] + v[i,t]*Pmin[i])m.addConstr(p[i,t-1] - p[i,t] <= R60[i] + w[i,t]*Pmin[i])# 5. 5.1 机组启停状态受最小启停时间影响的约束ut = max(t-(math.ceil(min_uptime[i])-1), 0) dt = max(t-(math.ceil(min_downtime[i])-1), 0)m.addConstr(sum(v[i, s] for s in range(ut,t)) <= u[i,t])m.addConstr(sum(w[i, s] for s in range(dt,t)) <= 1 - u[i,t])# 5.2 描述机组启停状态的关系约束if t>0:m.addConstr(v[i,t] - w[i,t] == u[i,t] - u[i,t-1])else:m.addConstr(v[t,i] - w[t,i] == u[t,i])# 6. (新能源)发电机的容量约束
for t in range(n_t):m.addConstr(p_wind[t] <= wind[t])m.addConstr(p_pv[t] <= pv[t])# 7. 备用容量要求约束
for t in range(n_t):# 能够提供的总备用容量是决策变量，最少也要大于最大单机容量m.addConstrs(r_req[t] >= p[i,t] for i in range(n_gen))m.addConstr(gp.quicksum(p_rs[i,t] + p_rns[i,t] for i in range(n_gen)) >= r_req[t])# 要求总的备用容量中至少有一半是通过旋转备用容量实现的m.addConstr(gp.quicksum(p_rs[i,t] for i in range(n_gen)) >= 0.5 * r_req[t])# 此优化模型的目标函数包括几部分：
# 1. (传统能源)发电机在发电过程中产生的成本
gen_cost = gp.quicksum(LinExpr(mc, [p[key] for key in pkey.select('*', t)]) for t in range(n_t))
# 2. (新能源)发电机的机会成本
curt_cost = gp.quicksum((pv[t]-p_pv[t])*pv_curt_pen + (wind[t]-p_wind[t])*wind_curt_pen for t in range(n_t))
# 3. (传统能源)发电机在启动时产生的成本
startup_cost = gp.quicksum(LinExpr(su_cost, [v[key] for key in pkey.select('*', t)]) for t in range(n_t))# 利用 Gurobi 求解，并打印最优值
m.setObjective(gen_cost + curt_cost + startup_cost, GRB.MINIMIZE)
m.optimize()
print_gp_status(m)
objective = m.ObjVal
print(f"Objective value {objective/1e3:.2f} k$.\n")# 结果可视化，展示风、光、传统能源发电机在不同时段的发电量
# 其中传统能源发电机分为几种，也将打印其随时间变化的曲线
p_res = {type: np.array([m.getVarByName(f"p[{i},{t}]").X for t in range(n_t)]) for (i,type) in enumerate(gen_id)}
p_res['Wind'] = np.array([m.getVarByName(f"p_wind[{t}]").X for t in range(n_t)])
p_res['PV'] = np.array([m.getVarByName(f"p_pv[{t}]").X for t in range(n_t)])# 分别打印成本风、光电机的总机会成本
wind_curt_res = wind - np.array(p_res['Wind'])
pv_curt_res = pv - np.array(p_res['PV'])
print(f"A total of {sum(wind_curt_res):.2f} MWh wind were curtailed.")
print(f"A total of {sum(pv_curt_res):.2f} MWh PV were curtailed.")# 将不同种类的传统能源发电机发电量进行整合后可视化
gen_type_unique = np.array(list(dict.fromkeys(gen_type)))
gens_per_type = {}
for i, gt in enumerate(gen_type):if gt not in gens_per_type:gens_per_type[gt] = []gens_per_type[gt].append(i)p_res_type = {gt : sum([p_res[gen_id[sub]] for sub in sublist]) for gt, sublist in gens_per_type.items()}
p_res_type['Wind'] = np.array([m.getVarByName(f"p_wind[{t}]").X for t in range(n_t)])
p_res_type['PV'] = np.array([m.getVarByName(f"p_pv[{t}]").X for t in range(n_t)])

上述程序的运行结果如下：

The model is optimal.
Objective value 188.18 k$.
A total of 1475.71 MWh wind were curtailed.
A total of 2866.40 MWh PV were curtailed.

\hspace{1.6em} 现在我们想分别得到每个时段的电价，可以通过对偶问题的最优解这种供需平衡这个等式约束的对偶乘子的取值得到，这也被称为拉格朗日松弛方法(Lagrangian relaxation method)，关于如何通过 Gurobi 求优化问题的影子价格，我已经在我的其他博客link中进行了阐述，有兴趣的可以看看。
\hspace{1.6em} 但是机组组合问题的对偶函数是难以直接列写的，我们可以考虑与机组组合问题紧密联系的另一个问题：经济调度问题。这种转化只是为了得到电价，作为一种 “捷径”，必须承认在数学上并不严谨，没有直面机组组合这一 MIP 问题，只是为了和这个系列博客之前的内容进行联系。
\hspace{1.6em} 经济调度问题就是已经确定必开机组，设计一个优化问题来确定这些必开机组如何调度这些机组的发电来降低成本，下面就把上面求UC问题求出来的机组开关状态固定，以此去掉离散决策变量，在此之前，不妨来看看我们在之前的模型中，不同时段下都打开了哪些些机组，

# plot of aggregated production
fig, ax = plt.subplots(1,1)
fig.set_size_inches(10,6)
x = np.arange(24)
color_list = [gen_colors[g] for g in list(p_res_type.keys())]
ax.stackplot(x, list(p_res_type.values()), labels=list(p_res_type.keys()), colors=color_list);
ax.plot(load, linewidth=3, color='black', label='Load')
ax.set_xlabel("Hour")
ax.set_ylabel("Generation [MW]")
ax.legend(ncols=5, loc="lower center", bbox_to_anchor=(0.5, -0.23));

运行结果

\hspace{1.6em} Gurobi 中为优化面模型提供了一个 fixed() 方法，这个方法将在模型中把每个整数变量固定为它在MIP解或MIP起始值中取的值。这一方法将会去掉模型中的整数变量及其相关的约束，模型就退化成了 经济调度模型，关于经济调度模型，我在本系列博客link中已经说明。
\hspace{1.6em} 固定机组组合后，让我们重新求解模型并可视化结果

m_fixed = m.fixed() 
m_fixed.optimize()print_gp_status(m)
objective = m_fixed.ObjVal
print(f"Objective value {objective/1e3:.2f} k$.\n")
p_res = {type: np.array([m_fixed.getVarByName(f"p[{i},{t}]").X for t in range(n_t)]) for (i,type) in enumerate(gen_id)}
p_res['Wind'] = np.array([m_fixed.getVarByName(f"p_wind[{t}]").X for t in range(n_t)])
p_res['PV'] = np.array([m_fixed.getVarByName(f"p_pv[{t}]").X for t in range(n_t)])# 之后的代码就是在可视化p_res，上面已经写过，不再赘述

上述程序的运行结果如下：

The model is optimal.
Objective value 188.18 k$.
A total of 1870.84 MWh wind were curtailed.
A total of 2471.27 MWh PV were curtailed.

\hspace{1.6em} 那么现在就可以获得每个时段的对偶价格了，将其求出来之后放到原图中呈现

system_price = [m_fixed.getConstrByName(f"energy_balance[{t}]").Pi for t in range(n_t)]
print(system_price)
fig, ax = plt.subplots(1,1)
fig.set_size_inches(10,6)
x = np.arange(24)
color_list = [gen_colors[g] for g in list(p_res_type.keys())]
ax.stackplot(x, list(p_res_type.values()), labels=list(p_res_type.keys()), colors=color_list);
ax.plot(load, linewidth=3, color='black', label='Load')
ax2 = ax.twinx() 
ax2.plot(system_price, color='blue', linewidth=3, ls="--", label="Price")
ax2.set_ylim(0)
ax2.set_ylabel("System price [$]", color="blue")
ax.set_xlabel("Hour")
ax.set_ylabel("Generation [MW]")
ax.legend(ncols=5, loc="lower center", bbox_to_anchor=(0.5, -0.23));

上述程序的运行结果如下：

[25.0, 14.0, 14.0, 14.0, 16.0, 4.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 2.0, 10.0, 14.0, 14.0, 12.0, 10.0, 10.0, 10.0]

更多内容，请见我的其他博客