当前位置： 首页 > news >正文

Electricity Market Optimization 探索系列（VII）- 最优潮流 与 节点边际电价

news 来源：原创 2025/5/9 3:38:24

本文参考链接link
\hspace{1.6em} 在上一节中，我们讨论了电力网络中潮流方程的由来，引出了 PTDF（Power Transmission Distribution Factor），本节我们将讨论最优潮流与节点边际电价

直流最优潮流问题的概念与基本假设

\hspace{1.6em} 最优潮流问题是指在一个特定的时段内，考虑了电力网络中潮流的分布情况的安全性约束的经济调度，其影响着整个电力系统运营。加入的这个约束是为了保证输电线路上的电流不会过大。我们考虑一下假设：

  • 假设1: 假设一个特定的时段内，每个节点的负荷是不变的
  • 假设2: 假设直流潮流方程和未经简化前的潮流方程的物理特性足够接近

\hspace{1.6em} 本文将延续之前所说的线路直流潮流方程的两种形式，稠密矩阵形式和稀疏矩阵形式，通过 Gurobi 求解三节点系统的直流最优潮流问题。

直流最优潮流问题的模型建立

\hspace{1.6em} 我们在之前的文章中 link 提到了经济调度问题，最优调度问题就是在经济调度问题的基础之上，添加一个直流潮流方程作为约束

带稠密矩阵形式的直流潮流方程约束的最优潮流问题

模型建立

min ⁡ ∑ k ∈ [ G ] c k ( p k ) s.t ( μ ) : ∑ k ∈ [ G ] p k = ∑ i ∈ [ N ] D i ( − ) : p k ≤ P ‾ k ∀ k ∈ [ G ] ( − ) : f l = ∑ i ∈ [ N ] T l , i ( ∑ j ∈ [ G ] i p j − D i ) ∀ l ∈ [ L ] ( λ l + ) : f l ≤ F ‾ l ∀ l ∈ [ L ] ( λ l − ) : − f l ≤ F ‾ l ∀ l ∈ [ L ] \begin{aligned} \min \quad & \sum_{k\in[G]} c_k(p_k) \\\\ \text{s.t} \quad (\mu):&\qquad \sum_{k\in[G]} p_k = \sum_{i\in[N]} D_i \\ (-):&\qquad p_k \le \overline{P}_k && \forall k\in[G] \\\\ (-):&\qquad f_l = \sum_{i\in[N]} T_{l,i}\Big(\sum_{j\in[G]_i} p_j - D_i\Big)&& \forall l\in[L] \\\\ (\lambda_l^+):&\qquad f_l \le \overline{F}_l && \forall l\in[L] \\ (\lambda_l^-):&\qquad -f_l \le \overline{F}_l && \forall l\in[L] \end{aligned} mins.t(μ):(−):(−):(λl+​):(λl−​):​k∈[G]∑​ck​(pk​)k∈[G]∑​pk​=i∈[N]∑​Di​pk​≤Pk​fl​=i∈[N]∑​Tl,i​(j∈[G]i​∑​pj​−Di​)fl​≤Fl​−fl​≤Fl​​​∀k∈[G]∀l∈[L]∀l∈[L]∀l∈[L]​
\hspace{1.6em} 其中 ( ∑ j ∈ [ G ] i p j ) − D i (\sum_{j\in[G]_i} p_j) - D_i (∑j∈[G]i​​pj​)−Di​ 是第 i i i 个节点的节点注入功率， [ N ] [N] [N] 是网络中的节点的编号集合， [ G ] i [G]_i [G]i​ 是网络中 i i i 节点的发电机的编号集合， [ L ] [L] [L] 是网络中的线路的编号的集合，对应 F = TP \mathbf F = \textbf T \textbf P F=TP 这个方程中 P \textbf P P 的第 i i i 个元素。为方便叙述，本文将在之后称这个模型为 模型(1)
\hspace{1.6em} 参数设定和求解代码均在下面的这段代码中，这段代码通过 Gurobi 进行求解，代码如下：

import pandas as pd
import numpy as np
import math
import gurobipy as gp
from gurobipy import *
import matplotlib.pyplot as plt# 设置潮流的上下界
c = [5, 20, 100] 
pmax = [102, 100, 100] 
fmax = [100, 100, 30] # Case 1
D = [0, 0, 130] # nodal demand
# Case 2
# D = [0, 30, 130] # nodal demandn_bus = 3 # number of buses/nodes
n_gen = 3 # number of generators
n_line = 3 # number of linesT = np.array([[-1/3, 1/3, 0], [-2/3, -1/3, 0],[1/3,  2/3, 0]])m = gp.Model()
m.setParam("OutputFlag", 0)p = m.addVars(n_gen, lb=0, ub=GRB.INFINITY, name="p") 
f = m.addVars(n_line, lb=-GRB.INFINITY, ub=GRB.INFINITY, name="f") # 供需平衡约束
m.addConstr(sum(p[i] for i in range(n_gen)) == sum(D[i] for i in range(n_bus)), name="enerbal")
m.addConstrs(p[i] <= pmax[i] for i in range(n_gen))# 稠密矩阵形式的直流潮流方程以及线路潮流约束
m.addConstrs((f[l] == sum(T[l,i]*(p[i] - D[i]) for i in range(n_bus)) for l in range(n_line)), name="flowdef")
m.addConstrs((f[l] <= fmax[l] for l in range(n_line)), name="flowlim_pos")
m.addConstrs((-f[l] <= fmax[l] for l in range(n_line)), name="flowlim_neg")m.setObjective(sum(p[i]*c[i] for i in range(n_gen)), GRB.MINIMIZE)
m.optimize()#  打印结果
print("System cost")
print(f"{m.ObjVal:.0f}")
print("Production:", end=" ")
p_ret = [m.getVarByName(f"p[{i}]").X for i in range(n_gen)]
print(p_ret)
print("Flow:", end=" ")
f_ret = [m.getVarByName(f"f[{i}]").X for i in range(n_line)]
print(f_ret)

代码的运行结果：

System cost
4450
Production: [90.0, 0.0, 40.0]
Flow: [-30.0, -60.0, 30.0]

\hspace{1.6em} 同样的，和本系列的其他诸多价格相关的问题一样，我们要考虑对偶问题以及其其对偶变量的最优值和原问题的关系，并引出我们的节点边际电价

带稀疏矩阵形式的直流潮流方程约束的最优潮流问题

模型建立

min ⁡ ∑ k ∈ [ G ] c k ( p k ) s.t ( λ i ) : ∑ k ′ ∈ [ G ] i p k ′ − D i = ∑ j : i j ∈ [ L ] f i j ∀ i ∈ [ N ] ( σ ‾ i ) : p k ≤ P ‾ k ∀ k ∈ [ G ] ( β i j ) : f i j = B i j ( θ i − θ j ) ∀ i j ∈ [ L ] θ i = s l a c k = 0 ( μ i j + ) : f i j ≤ F ‾ i j ∀ i j ∈ [ L ] ( μ i j − ) : − f i j ≤ F ‾ i j ∀ i j ∈ [ L ] \begin{aligned} \min \quad & \sum_{k\in[G]} c_k(p_k) \\\\ \text{s.t} \quad (\lambda_i):&\qquad \sum_{k^{'}\in[G]_i}p_{k^{'}} - D_i = \sum_{j:ij\in[L]}f_{ij} && \forall i \in [N] \\\\ (\overline{\sigma}_i):&\qquad p_k \le \overline{P}_k && \forall k\in[G] \\\\ (\beta_{ij}):&\qquad f_{ij} = B_{ij}(\theta_i - \theta_j )&& \forall ij\in[L] \\\\&\hspace{2em}\theta_{i=slack} = 0 \\\\ (\mu_{ij}^+):&\qquad f_{ij} \le \overline{F}_{ij} && \forall ij\in[L] \\\\ (\mu_{ij}^-):&\qquad -f_{ij} \le \overline{F}_{ij} && \forall ij\in[L] \end{aligned} mins.t(λi​):(σi​):(βij​):(μij+​):(μij−​):​k∈[G]∑​ck​(pk​)k′∈[G]i​∑​pk′​−Di​=j:ij∈[L]∑​fij​pk​≤Pk​fij​=Bij​(θi​−θj​)θi=slack​=0fij​≤Fij​−fij​≤Fij​​​∀i∈[N]∀k∈[G]∀ij∈[L]∀ij∈[L]∀ij∈[L]​

\hspace{1.6em} 还是上面的那三台电机，还是一样的三条线路，每个节点的需求还是那个时段的一个常数，只是现在给定这三条线路的阻抗值均为 1，并给定参考节点的编号，求解这三台电机的最优发电量。为方便叙述，本文将在之后称这个模型为 模型(2)
\hspace{1.6em} 参数设定和求解代码均在下面的这段代码中，这段代码通过 Gurobi 进行求解，代码如下：

import numpy as np
import math
import gurobipy as gp
from gurobipy import *
import matplotlib.pyplot as plt# 设置潮流的上下界
c = [5, 20, 100] 
pmax = [102, 100, 100] 
fmax = [100, 100, 30] # Case 1
D = [0, 0, 130] # nodal demand
# Case 2
# D = [0, 30, 130] # nodal demandn_bus = 3 # number of buses/nodes
n_gen = 3 # number of generators
n_line = 3 # number of lines#假设线路中的电纳恒定为1
B = 1
slack = 2
# 定义线路中的潮流限制矩阵如下, 其中fmax[i, j]表示的是从 i 节点到 j 节点的潮流最大值
fmax = np.array([[0, 100, 100], [100, 0, 30],[100, 30, 0]])m = gp.Model()
m.setParam("OutputFlag", 0)
p = m.addVars(n_gen, lb=0, ub=GRB.INFINITY, name="p") 
f = m.addVars(n_bus, n_bus, lb=-GRB.INFINITY, ub=GRB.INFINITY, name="f") 
theta = m.addVars(n_line, lb=-GRB.INFINITY, ub=GRB.INFINITY, name="f") 
for i in range(n_bus):m.addConstr(p[i] - D[i] == sum(f[i,j] for j in range(n_bus)), name="enerbal")
m.addConstrs(p[i] <= pmax[i] for i in range(n_gen))for i in range(n_bus):for j in range(n_bus):m.addConstr((f[i,j] == B*(theta[i] - theta[j])), name="flowdef")for i in range(n_bus): for j in range(n_bus):m.addConstr((f[i,j] <= fmax[i,j]), name="flowlim_pos")m.addConstr((-f[i,j] <= fmax[i,j]), name="flowlim_neg")m.addConstr(theta[slack] == 0,name="slack_bus")
m.setObjective(sum(p[k]*c[k] for k in range(n_gen)), GRB.MINIMIZE)
m.optimize()print("System cost")
print(f"{m.ObjVal:.0f}")print("Production:", end=" ")
p_ret = [m.getVarByName(f"p[{i}]").X for i in range(n_gen)]
print(p_ret)print("Flow:", end=" ")
f_ret = [f[i, j].x for i in range(n_bus) for j in range(n_bus)]
for i in range(0, len(f_ret), 3):print(f"[{f_ret[i]:.2f}, {f_ret[i + 1]:.2f}, {f_ret[i + 2]:.2f}]")

代码的运行结果：

System cost
4450
Production: [90.0, 0.0, 40.0]
Flow: [0.00, 30.00, 60.00]
[-30.00, 0.00, 30.00]
[-60.00, -30.00, 0.00]

\hspace{1.6em} 可以看出，模型(1) 和 模型(2) 是等价的，

节点边际电价

\hspace{1.6em} 根据经济学的定义，在有限的资源下，如果需求每增加一个单位，都要消耗更多的某一种资源，那么这种资源增加量所需要的成本即为边际价格。根据这一定义，在电力系统中，节点边际电价取决于系统向这一节点多提供的单位功率所增加的成本，定量的分析可以通过对 模型(1) 的对偶问题得到。下面是简要的推导。

节点边际电价的数学推导

模型(1) 的拉格朗日函数为:
L ( p → , λ L → + , λ L → − , μ ) = ∑ k ∈ [ G ] c k ( p k ) + ∑ l ∈ [ L ] λ l + ( ( f l − ∑ i ∈ [ N ] T l , i ( ∑ j ∈ [ G ] i p j − D i ) ) − F ‾ l ) + ∑ l ∈ [ L ] λ l − ( − ( f l − ∑ i ∈ [ N ] T l , i ( ∑ j ∈ [ G ] i p j − D i ) ) − F ‾ l ) + μ ( ∑ k ∈ [ G ] p k − ∑ i ∈ [ N ] D i ) \begin{aligned} L(\mathbf{\overrightarrow{p}}, \mathbf{\overrightarrow{\lambda_L}}^+,\mathbf{\overrightarrow{\lambda_L}}^-,\mu)=\\ &\sum_{k\in[G]} c_k(p_k) \\\\ &+\sum_{l\in[L]}\lambda_l^+( \Big(f_l -\sum_{i\in[N]} T_{l,i}\Big(\sum_{j\in[G]_i} p_j - D_i\Big)\Big)-\overline{F}_l)\\\\ &+\sum_{l\in[L]}\lambda_l^-(- \Big(f_l -\sum_{i\in[N]} T_{l,i}\Big(\sum_{j\in[G]_i} p_j - D_i\Big)\Big)-\overline{F}_l)\\\\ &+\mu\Big(\sum_{k\in[G]} p_k - \sum_{i\in[N]} D_i\Big) \end{aligned} L(p ​,λL​ ​+,λL​ ​−,μ)=​k∈[G]∑​ck​(pk​)+l∈[L]∑​λl+​((fl​−i∈[N]∑​Tl,i​(j∈[G]i​∑​pj​−Di​))−Fl​)+l∈[L]∑​λl−​(−(fl​−i∈[N]∑​Tl,i​(j∈[G]i​∑​pj​−Di​))−Fl​)+μ(k∈[G]∑​pk​−i∈[N]∑​Di​)​
L ( ⋅ ) L(\cdot) L(⋅) 对 p → \mathbf{\overrightarrow{p}} p ​ 求偏导，可以得到下述一系列式子： ∂ L ∂ p i = c i + μ − ∑ l ∈ [ L ] i λ l + T l , i + ∑ l ∈ [ L ] i λ l − T l , i ∀ i ∈ [ G ] ∀ l ∈ [ L ] \begin{align} \frac{\partial L}{\partial p_i}=c_i+\mu- \sum_{l\in[L]i} \lambda_l^+T_{l,i}+\sum_{l\in[L]i} \lambda_l^-T_{l,i}&&&&&\forall i\in[G] &&\forall l\in[L] \end{align} ∂pi​∂L​=ci​+μ−l∈[L]i∑​λl+​Tl,i​+l∈[L]i∑​λl−​Tl,i​​​​​​∀i∈[G]​​∀l∈[L]​​
\hspace{1.6em} 其中 [ L ] i [L]i [L]i 是与 i i i 节点直接相连的所有线路的编号的集合，还是以上面的代码中假设的参数为例，可以根据这三个偏导数在 λ L → + \mathbf{\overrightarrow{\lambda_L}}^+ λL​ ​+， λ L → − \mathbf{\overrightarrow{\lambda_L}}^- λL​ ​−， μ \mu μ 的不同取值范围内，会得到什么不同的对偶函数 g ( λ L → + , λ L → − , μ ) g(\mathbf{\overrightarrow{\lambda_L}}^+,\mathbf{\overrightarrow{\lambda_L}}^-,\mu) g(λL​ ​+,λL​ ​−,μ)，下面来分析这个以最大化对偶函数的对偶问题

定性分析

\hspace{1.6em} 可以通过本系列的其他博客得知，对偶函数是一个线性分段函数，而且是一个凹函数，其最大值就是节点边际电价。对于不同的发电机组，都有 ( 1 ) (1) (1) 式成立，对于那些在对偶问题中最优解带回 ( 1 ) (1) (1) 式，分为三种情况：

  • 编号带入后偏导数为 0 的机组，我们称为边际机组，这些机组的报价将会等于节点电价。
  • 编号带入后偏导数小于 0 的机组，他们的机组出力达到上限 P k ‾ \overline{{P}_k} Pk​​，这些机组的报价小于节点电价
  • 编号带入后偏导数大于 0 的机组，他们的机组出力达到下限 P k ‾ \underline{{P}_k} Pk​​，这些机组的报价大于节点电价

定量分析

\hspace{1.6em} 本文讨论的三台发电机组的边际成本分别是 5, 20, 100，发电机最大容量分别为 102, 100, 100，线路最大容量分别为 100, 100, 30，考虑这些参数就可以得到具体的对偶函数，人工得到这个对偶函数的过程比较复杂，幸运的是，Gurobi 提供了一个约束的对偶变量的最优解。(不过话又说回来，python 搭载的求解器求解一个新定义的模型挺慢的，不知道其他语言的怎么样)
{ ∂ L ∂ p 1 = 5 + μ + 1 3 ( λ 1 − − λ 1 + ) − 1 3 ( λ 1 − − λ 1 + ) = 5 + μ ∂ L ∂ p 2 = 20 + μ + 2 3 ( λ 2 − − λ 2 + ) + 1 3 ( λ 2 − − λ 2 + ) = 20 + μ + ( λ 2 − − λ 2 + ) ∂ L ∂ p 3 = 100 + μ + 1 3 ( λ 3 − − λ 3 + ) + 2 3 ( λ 3 − − λ 3 + ) = 100 + μ + ( λ 3 − − λ 3 + ) \begin{aligned} \begin{cases} \frac{\partial L}{\partial p_1}=5+\mu+\frac{1}{3}(\lambda_1^--\lambda_1^+)-\frac{1}{3}(\lambda_1^--\lambda_1^+)=5+\mu\\\\ \frac{\partial L}{\partial p_2}=20+\mu+\frac{2}{3}(\lambda_2^--\lambda_2^+)+\frac{1}{3}(\lambda_2^--\lambda_2^+) = 20+\mu+(\lambda_2^--\lambda_2^+)\\\\ \frac{\partial L}{\partial p_3}=100+\mu+\frac{1}{3}(\lambda_3^--\lambda_3^+)+\frac{2}{3}(\lambda_3^--\lambda_3^+)=100+\mu+(\lambda_3^--\lambda_3^+) \end{cases} \end{aligned} ⎩ ⎨ ⎧​∂p1​∂L​=5+μ+31​(λ1−​−λ1+​)−31​(λ1−​−λ1+​)=5+μ∂p2​∂L​=20+μ+32​(λ2−​−λ2+​)+31​(λ2−​−λ2+​)=20+μ+(λ2−​−λ2+​)∂p3​∂L​=100+μ+31​(λ3−​−λ3+​)+32​(λ3−​−λ3+​)=100+μ+(λ3−​−λ3+​)​​

# 阻塞电价
mu_pos = [m.getConstrByName(f"flowlim_pos[{l}]").Pi for l in range(n_line)]
mu_neg = [m.getConstrByName(f"flowlim_neg[{l}]").Pi for l in range(n_line)]
print("mu+", end="")
print(mu_pos)
print("mu-", end="")
print(mu_neg)# 节点电价：对于一个节点来说，节点电价取决于供需关系和与其相连的所有线路的阻塞电价
lam = m.getConstrByName("enerbal").Pi
lmp = [lam + sum(T[l,i]*(mu_neg[l]- mu_pos[l]) for l in range(n_line)) for i in range(n_bus)]
print("lambda: ", end="")
print(lam)
print("LMPs:", end="")
print(lmp)

代码的运行结果：

mu+：[0.0, 0.0, -285.0]
mu-：[0.0, 0.0, 0.0]
lambda: 100.0
LMPs: [195.0, 290.0, 100.0]

相关文章：

  • 微信小程序三种裁剪动画有效果图
  • 智能云图库-12-DDD重构
  • pycharm无法识别到本地python的conda环境解决方法
  • Allegro、Etsy、Joom、敦煌测评自养号技术优势
  • AI写论文哪个工具靠谱（四个AI工具对比分析）
  • ACL(访问控制列表）
  • 突破反爬限制的智能数据采集实战 —— 面向中小企业的高效信息监控方案
  • Redis面试问题详解2
  • 开源一体化白板工具Drawnix本地部署打造毫秒级响应的远程协作空间
  • win10和win11系统修复工具各类故障解决
  • DevOps-文章目录
  • Linux 防火墙（ iptables ）
  • 全球6G大会 | 紫光展锐用“芯”推动空天地一体创新纪元
  • 2025.4.17学习日记 初识JavaScript 以及Java和JavaScript有什么区别
  • [漏洞篇]CSRF漏洞详解
  • visual Studio+Qt插件检查内存泄漏
  • opencv HSV的具体描述
  • div（HTML标准元素）和view（微信小程序专用组件）的主要区别体
  • FPGA入门学习Day1——设计一个DDS信号发生器
  • Python开发环境打包迁移指南：离线与在线环境的完美解决方案
  • 85后清华博士黄佐财任湖北咸宁市咸安区委副书记、代区长
  • 一季度全国消协组织为消费者挽回经济损失23723万元
  • “用鲜血和生命凝结的深厚情谊”——习近平主席署名文章中的中俄友好故事
  • 江苏淮安优化村级资源配置：淮安区多个空心村拟并入邻村
  • 体坛联播｜国米淘汰巴萨晋级欧冠决赛，申花击败梅州避免连败
  • 巴基斯坦军方：印度导弹袭击已造成至少3人死亡