Electricity Market Optimization 探索系列（V）

本文参考链接link

\hspace{1.6em} 众所周知, 社会福利是指消费者剩余和生产者剩余的和，也等价于产品的市值减去产品的成本，在电力市场中也非常关注社会福利这一概念，基于电力商品的同质性的特点，我们引入反价格需求函数来形象地刻画电力市场中的社会福利这一概念

电力市场中的需求函数和供给函数

反价格需求函数：
\quad\space\space    而反价格需求函数则是将价格表示为需求量的函数，即$\pi (q^D)=\beta q^D+{\pi }^{max}$。它反映了在不同的需求量下，消费者对于每单位商品或服务所愿意支付的最高价格

其中已知的参数列表如下

符号	说明	值
$\beta$	需求曲线的斜率	-25
$p_{i,max}$	第$i$个发电机的最大发电量	[200, 100, 500, 200]
$c_i$	第$i$个发电机的成本系数	[0, 10, 50, 80]
${\pi }^{max}$	最大价格	11000
$q_{max}$	最大需求量	440

import numpy as np
import plotly.graph_objects as go# 反价格需求函数
# pi(q) = beta * q + pi_max
beta = -25
pi_max = 11e3
pi = lambda q: beta*q + pi_max
c = [0, 10, 50, 80]  # vector of cost [$/MW]
p_max = [200, 100, 500, 200]  # vector of maximum capacity [MW]
gen_type = ["RES", "cheap", "base", "peak"]
n_gen = len(c)  # number of generators
fig = go.Figure()hide_traces = [2, 3, 5, 6, 8, 9, 11, 12]
for i in range(n_gen):x1 = sum(p_max[:i]) x2 = sum(p_max[:i + 1])x = [x1, x2]y = [c[i], c[i]]trace_num = 3 * i + 1# 横轴show_legend = trace_num not in hide_tracesfig.add_trace(go.Scatter(x=x, y=y, mode='lines', name=gen_type[i], line=dict(width=3), showlegend=show_legend))# 纵轴show_legend = (trace_num + 1) not in hide_tracesfig.add_trace(go.Scatter(x=[x1, x1], y=[0, c[i]], mode='lines', line=dict(color='black', dash='dash', width=0.8), showlegend=show_legend))show_legend = (trace_num + 2) not in hide_tracesfig.add_trace(go.Scatter(x=[x2, x2], y=[0, c[i]], mode='lines', line=dict(color='black', dash='dash', width=0.8), showlegend=show_legend))xs = np.arange(0, 1000, 50)  
ys = [pi(x) for x in xs]  
fig.add_trace(go.Scatter(x=xs, y=ys, mode='lines', name="demand", line=dict(color='black')))fig.update_layout(yaxis_title="Price [$]",xaxis_title="Production [MW]",yaxis_range=[-1, 100],legend=dict(orientation="h",yanchor="bottom",y=1.02,xanchor="right",x=1)
)
fig.show()

程序结果图：

\quad\space\space    如果从供需曲线上看，社会福利还可以通过 三角形 + 矩形 - 去多个小矩形面积之和求得，其中三角形 + 矩形的面积就是产品市值，可以用$(1),(2),(3)式$进行总结
$$\begin{array}{rl}V& =\frac{1}{2}{q}^{\ast }({\pi }^{max}-{\pi }^{\ast })+q^{\ast }{\pi }^{\ast }=\frac{1}{2}{q}^{\ast }({\pi }^{\ast }+{\pi }^{max})\\ & =\frac{1}{2}{q}^{\ast }(\beta q^{\ast }+{\pi }^{max}+{\pi }^{max})=\frac{1}{2}\beta (q^{\ast }{)}^2+{\pi }^{max}{q}^{\ast }\\ C& =\sum _{i=1}^G{c}_i{p}_i^{\ast }\end{array}$$

那么社会福利最大化模型，就是在满足平衡的基础之上实现社会福利最大化，其模型可以简要的概括为如下数学模型

原问题

$$\begin{array}{rlrl}\underset{p,q^D}{\max }& V-C& & \\ \text{s.t.}& \sum _{i=1}^G{p}_i=q^D& & \\ & p_i\le p_{i,max}& & \\ & p_i,q^D\ge 0& & \forall i\in [G]\end{array}$$
为了方便数学分析，先讨论成以下模型，注意这和上述原文题并不完全等价，区别会在后面的分析中解释

$$\begin{array}{rlr}& \underset{\overrightarrow{\mathbf{p}},q^D}{\min }-\frac{1}{2}\beta (q^D{)}^2-{\pi }^{max}{q}^D+{\overrightarrow{\mathbf{c}}}^T\overrightarrow{\mathbf{p}}& \\ & \text{s.t.}\sum _{i=1}^G{p}_i=q^D& \overrightarrow{\mathbf{p}},q^D\in domf\end{array}$$

下面将使用优化理论分析如上模型，求出对偶问题

拉格朗日函数：
$$\begin{array}{rl}& L(\overrightarrow{\mathbf{p}},q^D,\mu )=-\frac{1}{2}\beta (q^D{)}^2-{\pi }^{max}{q}^D+{\overrightarrow{\mathbf{c}}}^T\overrightarrow{\mathbf{p}}+\mu \left(\sum _{i=1}^G{p}_i-q^D\right)\\ & \\ & 令{\overrightarrow{\mathbf{x}}}^T=[q^D,p_1,p_2,p_3,p_4]，并将各种参数带入\\ & \\ & L(\cdot )=\frac{1}{2}{\overrightarrow{\mathbf{x}}}^T\left[\begin{array}{ccccc}25& 0& 0& 0& 0\\ 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0\end{array}\right]\overrightarrow{\mathbf{x}}+[-\mu -11000,\mu ,10+\mu ,50+\mu ,80+\mu ]\cdot \overrightarrow{\mathbf{x}}\\ & \\ & 其中\left[\begin{array}{c}0\\ 0\\ 0\\ 0\\ 0\end{array}\right]\le \overrightarrow{\mathbf{x}}\le \left[\begin{array}{c}440\\ 200\\ 100\\ 500\\ 200\end{array}\right]\end{array}$$

多元函数求极值，对 $\overrightarrow{\mathbf{x}}$ 求偏导
$$\begin{array}{rl}\frac{\partial L}{\partial \overrightarrow{\mathbf{x}}}& =0\\ & \\ \Rightarrow \left[\begin{array}{ccccc}25& 0& 0& 0& 0\\ 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0\end{array}\right]\overrightarrow{\mathbf{x}}& =\left[\begin{array}{c}11000+\mu \\ -\mu \\ -\mu -10\\ -\mu -50\\ -\mu -80\end{array}\right]\end{array}$$

显然该线性方程无解，所以 $L(\cdot )$ 在边界处取到极值，此时要反过来考虑 $\mu$ 对于拉格朗日对偶函数 $g(\mu )$ 的影响，
$$\begin{array}{r}g(\mu )=\{\begin{array}{llcc}-\frac{1}{50}(11000+\mu {)}^2& & & \mu >0\\ -\frac{1}{50}(11000+\mu {)}^2+200\mu & & & -10<\mu <0\\ -\frac{1}{50}(11000+\mu {)}^2+300\mu +1000& & & -50<\mu <-10\\ -\frac{1}{50}(11000+\mu {)}^2+800\mu +26000& & & -80<\mu <-50\\ -\frac{1}{50}(11000+\mu {)}^2+1000\mu +42000& & & \mu <-80\end{array}\end{array}$$
最后，对于 $(8),(9)$ 两式确定的优化问题，其对偶问题可以写成
$$\begin{array}{rl}& \underset{\mu }{\max }g(\mu )\end{array}$$
用python求解：

import numpy as np
from scipy.optimize import minimize_scalar# 定义负的分段函数，用于求最大值
def negative_g(mu):if mu > 0:return (1 / 50) * (11000 + mu) ** 2elif -10 < mu < 0:return (1 / 50) * (11000 + mu) ** 2 - 200 * muelif -50 < mu < -10:return (1 / 50) * (11000 + mu) ** 2 - 300 * mu - 1000elif -80 < mu < -50:return (1 / 50) * (11000 + mu) ** 2 - 800 * mu - 26000elif mu < -80:return (1 / 50) * (11000 + mu) ** 2 - 1000 * mu - 42000else:return np.nan# 定义区间，将无穷区间替换为足够大的有限区间
intervals = [(-1e6, -80), (-80, -50), (-50, -10), (-10, 0), (0, 1e6)]max_values = []
max_args = []for interval in intervals:result = minimize_scalar(negative_g, bounds=interval, method='bounded')max_values.append(-result.fun)max_args.append(result.x)# 找到全局最大值及其对应的自变量值
global_max_index = np.argmax(max_values)
global_max = max_values[global_max_index]
global_max_arg = max_args[global_max_index]print(f"分段函数的最大值为: {global_max}，此时自变量 μ 的取值为: {global_max_arg}")

最终 $g(\mu )$ 在 $\mu =-50$ 处取到最大值，最大值为-2412050对于这个结果，需要做以下解释：

1. -2412050是由于 $maxf(x)$ 转换成 $min-f(x)$ 之后，这两个问题求得的最优值是互为相反数的， $(10)$ 是 $(8),(9)$ 的对偶问题，而且经验证，满足强对偶性质，所以2412050是原文题的最优值

2. $g(\mu )$ 这个函数是一个连续可微的凹函数，其图像如下
   

3. 原问题中 $\overrightarrow{\mathbf{x}}$ =$\left[\begin{array}{c}438\\ 200\\ 100\\ 138\\ 0\end{array}\right]$ 是最优解，那么这和转化问题之后的模型的拉格朗日函数有什么关系？
   当 $\mu =-50$ 时，$L(\overrightarrow{\mathbf{p}},q^D,-50)=$ 这个函数要取到最小值要取到最小值，分析 $L(\overrightarrow{\mathbf{p}},q^D,-50)=-\frac{1}{2}\beta (q^D{)}^2-{\pi }^{max}{q}^D+{\overrightarrow{\mathbf{c}}}^T\overrightarrow{\mathbf{p}}-50\left({\sum }_{i=1}^G{p}_i-q^D\right)$ 可知，取到最小值要求 $q^D=438$，且 $x_1$ 和 $x_2$ 都取到各自定义域的上限，$x_4$ 取到定义域的下限，对 $x_3$ 没有要求，说明只有 $x_3$ 不是在边界处取到，也就是约束对 $x_3$ 而言是松弛的，

4. 从 3. 中的分析还可以更深刻的理解文章一开始强调的电力市场中的需求函数和供给函数中供给函数为什么一定要那么画，从图形来说，可以说供给函数先将发电量按照价格从低到高排序是为了将面积变得更大，从数学分析来说，一旦需求提高，首先需要达到上限的一定是价格较低的发电机。

用 Gurobi 来求解社会福利最大化问题

在实际的电力市场中，需要先考虑社会福利最大化，但是根本不需要经过上面这样手算数学模型，只需要借助 gurobi 求解器这样的工具帮我们求解原问题即可，具体代码如下：

import numpy as np
import gurobipy as gp
from gurobipy import GRB#反需求函数
# pi(q) = beta * q + pi_max
beta = -25
pi_max = 11e3
pi = lambda q: beta*q + pi_max# 发电机的价格向量
c = [0, 10, 50, 80] 
# 发电机的容量向量
p_max = [200, 100, 500, 200]
gen_type = ["RES", "cheap", "base", "peak"]
n_gen = len(c)m = gp.Model()
m.setParam('OutputFlag', 0)
p = m.addVars(n_gen, lb=0, ub=GRB.INFINITY, name="p")
d = m.addVar(lb=0, ub=GRB.INFINITY, name="d")m.addConstr(p.sum() == d, name="enerbal")
for i in range(n_gen):m.addConstr(p[i] <= p_max[i])# 目标函数
V = 0.5 * beta * (d**2) + pi_max * d
C = sum(p[i]*c[i] for i in range(n_gen))
welfare = V - C
m.setObjective(welfare, GRB.MAXIMIZE)m.optimize()welfare_res = m.ObjVal
d_res = d.X
p_res = [p[i].X for i in range(n_gen)]
C_res = C.getValue()
print(f"Resulting welfare: {welfare_res:.1f} $")
print(f"Resulting generation cost: {C_res:.1f} $")
print(f"Resulting demand:  {d_res:.1f} MW")
print(f"Generator outputs:")
for i in range(n_gen):print(f"   {gen_type[i]}: {p_res[i]:.1f} MW")# 计算市场出清价格
pi_opt = pi(d_res)
print(f"The market clearing price is {pi_opt:.1f} $/MW")# 对偶问题的最优解
lam = m.getConstrByName("enerbal").Pi
print(f"The dual of the energy balance is {lam:.1f}")

运行结果：

Resulting welfare: 2412050.0 $
Resulting generation cost: 7900.0 $
Resulting demand: 438.0 MW
Generator outputs:
RES: 200.0 MW
cheap: 100.0 MW
base: 138.0 MW
peak: 0.0 MW
The market clearing price is 50.0 $/MW
The dual of the energy balance is -50.0

社会福利最大化之后的经济调度问题

由结果可知，与我们上面的分析一致，不难发现，社会福利最大化模型和经济调度模型的区别就是社会福利最大化模型考虑了发电量和需求，而经济调度模型只考虑了在给定用电需求的条件下的各个发电机的发电量问题，关于经济调度的具体问题，可见我的另外一篇博客：link