简介这是一套基于内点法求解实时最优电价的MATLAB实现方案面向电力市场研究者、电气工程专业学生以及从事电网调度优化的技术人员。资源以30节点电力网络为例完整展示了分时电价TOU策略下如何通过内点法迭代求解实时电价兼顾供需平衡与运营成本最小化。压缩包共7个文件包括4个MATLAB脚本主程序、数据输入、初始化及辅助计算、2个fig格式的迭代过程图分别展示不同向心参数及预测-校正环节的影响和1个说明文档整体仅40KB结构紧凑。该资源已有118人学习适合希望掌握最优潮流计算与电价建模的读者参考。通过运行模型可以直观对比不同参数下的收敛行为并结合实际电价数据验证算法效果为电力市场定价机制设计提供可复现的实验基础。1. 基于现货的实时最优电价为什么固定分时电价越来越不够用现货市场每 15 分钟出一个价格凌晨风光大发时可能只有几分钱晚高峰却能冲到一块五以上。而多数用户签的还是那种一年不变的峰谷分时电价Time of Use, TOU晚高峰照收高价深谷时段也照收固定低价。售电公司被现货价格和固定 TOU 价差两头挤用户也没拿到现货红利。更麻烦的是光伏渗透率越高的地区午间现货价被压得越低固定 TOU 却还在收平段价这等于把市场信号完全扭曲了。要做实时最优电价就得把现货价格、用户需求响应和电网约束一起扔进一个优化问题再用内点法在可接受的秒级时间内求出每个时点的最优电价。这套东西不算新真正让人头疼的是建模和调参——特别是 spot.rar 这类数据包里时间戳、缺失值和弹性系数任何一个没处理好解出来的电价曲线就是一张废纸。2. 从现货价格到 TOU 电价目标函数、约束与内点法选型2.1 先把问题写成数学形式需求响应函数与电价决策变量做实时最优电价第一步不是写代码而是把问题定义清楚。决策变量很简单未来一天 96 个时点15 分钟粒度的零售电价 (p_t)。输入是现货价格 (c_t)以及每个时点用户对电价的响应关系。用户需求不是固定的会随电价变化。常用的线性需求响应模型是[ D_t(p_t) A_t - b_t p_t ]其中 (A_t) 是用户在该时点的基线需求(b_t) 是价格弹性斜率。(b_t) 越大用户对电价越敏感。实际建模时(b_t) 往往由弹性系数 (\eta) 换算而来(b_t \eta \cdot A_t / p_{\text{ref}})(\eta) 一般取 0.10.3工业用户高一些居民低一些。目标函数我一般这么写让售电公司从现货市场购电的总成本最小同时尽量避免电价曲线剧烈跳动再叠加一个削峰惩罚项。[ \min \sum_t c_t D_t(p_t) \frac{\rho}{2} \sum_t D_t(p_t)^2 \frac{\alpha}{2} \sum_t (p_t - p_{t-1})^2 ]第一项是买电成本第二项是削峰软约束第三项是相邻时段电价平滑项。约束条件包括零售电价上下限 (p_{\min} \le p_t \le p_{\max})、平均电价不低到亏本、响应后的负荷不能超过基线负荷的 1.2 倍。这个模型里目标函数是非线性的因为 (p_t) 和 (D_t) 相乘约束又是一堆不等式用内点法处理最顺手。2.2 内点法 vs 线性规划、梯度下降实时电价场景为什么选它很多同行第一反应是用线性规划但把需求响应写进去目标不再是线性的单纯形法用不了。用普通梯度下降也能跑可碰上 (D_t \ge 0)、削峰限值这类不等式约束梯度下降没法保证解始终在可行域里罚函数法又得反复调惩罚系数时间都耗在调参上了。内点法的思路是把所有不等式约束 (g_i(x) \ge 0) 以对数障碍项的形式并入目标函数然后对一系列衰减的障碍参数 (\mu) 做牛顿迭代。它兼顾了非线性目标和大量不等式约束收敛速度和稳定性都远好于手动罚函数。方法非线性目标大量不等式约束秒级求解实现难度线性规划/单纯形不支持支持快低梯度下降罚函数支持勉强看调参低内点法支持支持快中实时电价的求解规模通常只有 96336 个变量一天 96 点或一周 336 点内点法在这种情况下几十次迭代就能收敛。cvxpy 自带的 CLARABEL 求解器本质就是内点法实现直接调用即可不需要自己造轮子除非你在做研究或者要嵌入 C 程序。2.3 拿到 spot.rar 先核四件事数据体检与时间戳对齐spot.rar 这类压缩包解压后通常是一堆 CSV 或 Excel 文件里面放的现货电价、负荷、温度数据。我拿到手不会急着跑模型先做四件事看目录结构、看字段名、看时间间隔、看缺失值。$ unar spot.rar -o ./spot_data/ $ ls -R ./spot_data解压后先扫一眼文件布局确认有没有说明文档。接着用 pandas 读取重点核对时间戳。import pandas as pd df pd.read_csv(./spot_data/spot_price.csv, parse_dates[time]) print(df.head()) print(df.info()) print(df.isna().sum()) # 缺失检测 print(df[time].diff().value_counts()) # 时间间隔是否均匀这一步能筛出很多暗坑时间戳是不是按 15 分钟对齐的、有没有重复行、现货电价单位是元/MWh 还是元/kWh。我见过把功率单位 MW 和电量单位 MWh 混在一起的数据直接在价格上差了两个数量级。单位不核对清楚后面算出来的最优电价谁都不敢用。3. 用内点法求解实时最优电价最小可复现代码3.1 用 cvxpy 三分钟跑通第一版建模代码与求解器选型把数学问题翻译成 cvxpy 代码逻辑非常直接。下面是一天 96 个时点的完整建模import cvxpy as cp import numpy as np T 96 p cp.Variable(T) # 待求的零售电价 c cp.Parameter(T, nonnegTrue) # 现货电价 A cp.Parameter(T, nonnegTrue) # 基线需求 b cp.Parameter(T, nonnegTrue) # 需求弹性斜率 rho 0.001 # 削峰惩罚系数 alpha 0.5 # 相邻时段平滑系数 p_min_s 0.10 # 零售电价下限单位元/kWh p_max_s 1.50 # 零售电价上限 D A - cp.multiply(b, p) # 需求响应后的负荷 cost_power cp.sum(cp.multiply(c, D)) # 现货购电成本 cost_peak 0.5 * rho * cp.sum_squares(D) cost_smooth 0.5 * alpha * cp.sum_squares(p[1:] - p[:-1]) objective cp.Minimize(cost_power cost_peak cost_smooth) constraints [ p p_min_s, p p_max_s, cp.mean(p) 0.35, # 平均电价下限防整体压价 D 0, # 负荷不能为负 D 1.2 * A # 削峰响应后负荷不超过基线1.2倍 ] prob cp.Problem(objective, constraints) prob.solve(solvercp.CLARABEL, verboseFalse) price_opt p.value print(price_opt[:10])这段代码的rhos和alpha值得仔细调。rho太大最优解会牺牲经济性强行压平负荷alpha太大电价曲线变成一条直线等于把实时信号抹掉了。我一般先跑一组固定参数rho0.001作为起点然后做灵敏度扫描再定。cp.multiply(b, p)是按位相乘D是 96 维的变量表达式。这里的关键是c、A、b用 Parameter 而不是直接的 numpy 数组这样后续做灵敏度分析和滚动优化时只需更新参数再重新solve()不用改模型结构。3.2 手写一个简化内点法KKT 条件、障碍参数与牛顿迭代cvxpy 适合落地但要调试内点法的收敛问题还是得理解底层在干嘛。最经典的是原始-对偶内点法这里给一个更直白的障碍法实现适合用来理解障碍参数 (\mu) 的作用。import numpy as np class BarrierIPM: 简化的障碍内点法适用于目标为二次、约束为线性的问题。 约束形式A x b即每个约束写成 g_i(x) b_i - (A_i)·x 0 def __init__(self, Q, q, A, b, mu00.1, mu_ratio0.5, tol1e-8): self.Q, self.q, self.A, self.b Q, q, A, b self.mu mu0 self.mu_ratio mu_ratio self.tol tol def g(self, x): return self.b - self.A x def phi(self, x): # 障碍目标原目标 - mu * sum(log(g_i)) return 0.5 * x self.Q x self.q x - self.mu * np.sum(np.log(self.g(x))) def grad(self, x): gv self.g(x) return self.Q x self.q - self.mu * (self.A.T (1.0 / gv)) def hess(self, x): gv self.g(x) return self.Q self.mu * (self.A.T np.diag(1.0 / gv**2) self.A) def solve(self, x0): x x0.copy() for outer in range(50): if self.mu 1e-8: break for inner in range(100): gv self.g(x) if np.min(gv) 0: raise RuntimeError(约束被冲破检查初始点或步长策略) dx -np.linalg.solve(self.hess(x), self.grad(x)) alpha 1.0 # 回溯线搜索先保证可行再保证目标下降 while np.min(self.g(x alpha * dx)) 0: alpha * 0.5 x_new x alpha * dx while self.phi(x_new) self.phi(x) 1e-4 * alpha * self.grad(x).dot(dx): alpha * 0.5 x_new x alpha * dx x x_new if np.linalg.norm(dx, np.inf) self.tol: break self.mu * self.mu_ratio # 障碍参数衰减逐步逼近原问题 return x核心逻辑就两件事外循环衰减 (\mu)内循环对当前障碍目标做牛顿迭代。grad和hess我用的都是解析式因为目标二次、约束线性求导不难。障碍项 (-\mu \ln g_i(x)) 在约束边界附近形成一个势垒无论怎么迭代解都不会越过边界。实际调试时mu00.1mu_ratio0.5是比较稳的组合。mu_ratio再小一点比如 0.2收敛更快但约束刚好在边界上时容易震荡甚至发散。这个手写版本因为没加原始-对偶残差校正对初始点要求较高必须严格可行也就是说所有约束初始值都大于 0。3.3 内点法必调的三个参数障碍因子、收敛阈值、回溯步长内点法的调参经验基本能直接套用到 cvxpy 那些黑匣子求解器上。三个参数最重要。参数经验值作用调参方向障碍初始值 (\mu_0)0.1决定最早迭代时障碍项的强度约束多时减小防止初始迭代太钝障碍衰减比 (\mu_{\text{ratio}})0.5外循环收敛速度0.20.5太大收敛慢太小易震荡收敛阈值1e-8牛顿步长小于该值时认为内循环收敛调度场景 1e-6 也够用回溯线搜索的系数我习惯用1e-4判断 Armijo 条件步长衰减用 0.5。这里有个翻车点如果初始点刚好落在边界附近(g_i(x)) 接近 0障碍项梯度会变得巨大第一轮牛顿步就可能直接冲出可行域。解决办法是初始化时尽量取可行域的中点或者先用一个很小的 (\mu_0) 起步。另外如果求解报告提示 Numerical problems 或者 Terminated because of numerical difficulties九成是数据量纲问题。现货电价用 0.2 和 800 这种不同量纲去喂同一个模型海森矩阵会病态。先把所有输入归一化到 01 之间求解完再映射回实际电价这是最省事的后悔药。4. 把最优解落成实时电价曲线时段划分与后处理4.1 把 96 个时点的电价聚成峰谷平KMeans 与连续时段合并实时最优电价求解出来是 96 个连续数值但用户侧套餐不能 15 分钟一个价否则根本没法对外解释。常见做法是聚成峰、平、谷三段。用 KMeans 把电价归成 3 类按类别均值排序再映射回时段。from sklearn.cluster import KMeans price_flat price_opt.reshape(-1, 1) km KMeans(n_clusters3, n_init10, random_state0) labels_raw km.fit_predict(price_flat) # 按类中心从小到大排序0谷1平2峰 center_rank np.argsort(km.cluster_centers_.flatten()) labels np.empty_like(labels_raw) for new_label, old_label in enumerate(center_rank): labels[labels_raw old_label] new_label for label in range(3): idx np.where(labels label)[0] print(f时段{label}: 时点 {idx.min()}-{idx.max()}, f平均电价 {price_opt[idx].mean():.3f})KMeans 的坑是标签号不是按价格高低排列的必须用cluster_centers_排序重贴标签否则你可能把峰段标成谷段。另外聚类结果常出现峰-谷-峰这种断裂时段一个小时高峰、两个小时平段、又一个小时高峰用户侧没法用。我的处理办法是少于 2 小时8 个时点的时段段落直接合并到相邻类别再重新计算边界。这种分段方法有个玄学点K 到底取几。有人用肘部法则有人直接按政府目录电价设峰谷平三段。我一般先跑 3 段看结果如果峰段和谷段价差小于 0.2 元/kWh就换 4 段说明现货波动已经细化到需要尖峰段了。4.2 电价修正三步平滑、削尖、约束回填聚类只是第一轮后处理真正能落地还得过三关。第一关是平滑KMeans 分类边界处两个时点的电价可能从 0.2 直接跳 1.1这种尖刺对用户不友好。第二关是削尖把超过现货价 2 倍或低于现货价 0.5 倍的时点电价往回收。第三关是校验约束回填后重新检查削峰约束是否被破坏。def smooth_and_clip(p, c, win5, k3): w np.ones(win) / win p_s np.convolve(p, w, modesame) # 滑动平均 p_s np.minimum(p_s, 2.0 * c) # 不高过现货2倍 p_s np.maximum(p_s, 0.5 * c) # 不低于现货0.5倍 p_s np.maximum(p_s, p_min_s) # 回填上下限 p_s np.minimum(p_s, p_max_s) return p_s滑动平均会削掉真正的现货尖峰信号所以c现货价参与截断很关键。边界处np.convolve(modesame)会引入边缘效应前两个点数值偏低我常手动回填这两个点为原始值或者干脆对首尾时点不做平滑。最后一步是把平滑后的电价拿回去跑一遍需求响应函数确认D 1.2 * A没有被破坏如果削峰约束还差一点就把rho调大重新求解这不是后处理能救回来的。4.3 输出格式给营销系统和账单系统的电价表输出电价表要同时满足两拨消费方一是营销系统要按峰谷平时段展示给用户二是账单系统按时间戳和价格计算电费。字段设计一句话说清一天 96 个时点、每个时点一个价格、标注时段类型、时间戳用 ISO8601。import pandas as pd df_out pd.DataFrame({ time: pd.date_range(2024-01-01, periods96, freq15min), price: price_opt, period: labels }) df_out[period_name] df_out[period].map({0: 谷, 1: 平, 2: 峰}) df_out.to_csv(tou_price_20240101.csv, indexFalse)账单系统最怕两件事时间戳没有时区、时段边界不连续。我建议文件里直接写2024-01-01 00:00:0008:00这种带时区的格式让下游系统自己解析。时段边界要检查没有空档比如 0:006:00 是谷那么 6:00 必须接上平中间不能漏时点。曾有同事输出文件里 23:45 和 00:00 之间断了 15 分钟账单系统当成无电价时段用户投诉一堆。5. 内点法做实时最优电价最常见的四个坑现象、原因、解决5.1 迭代发散或震荡目标函数出现 NaN现象求解器报Numerical problems或者手写内点法跑到第 10 轮目标函数变成 NaN。原因(\mu) 衰减太快牛顿方向在约束边界附近产生巨大的障碍梯度一步跨过可行域。另一个常见原因是初始点不可行比如用全零向量初始化而约束里有 (p \ge 0.1)零点在边界外。解决(\mu_0) 从 0.1 起步衰减比别小于 0.5初始点取可行域中点比如所有约束下界的 1.1 倍如果数据跨度大先归一化再求解。5.2 解出负电价或极端尖峰现货低价时模型让人白嫖现象现货某时段价格是 0.02 元/kWh模型解出的最优零售电价只有 0.001 元甚至直接触到 0。原因目标函数里购电成本是 (c_t D_t)现货越低模型越希望用低价刺激负荷去消纳电力这在数学上没错但零售侧价格倒挂用户会拼命加装充电桩套利售电公司每度亏一毛。解决给零售电价加一个和现货联动的下限约束至少覆盖现货价加合理附加。接线时改一行约束即可margin cp.Parameter(T, nonnegTrue) # 每度电的输配附加与利润留成 margin.value np.full(T, 0.02) constraints.append(p c margin)5.3 现货价时间戳错位峰谷时段整体算反现象模型跑出来的峰段是凌晨谷段是晚高峰和实际用电曲线完全相反。原因spot.rar 里的现货价格是按交易日期还是执行日期存的很多数据源两者差一天直接用行号对齐就整体平移了 24 小时。还有时区问题有些厂站存的是 UTC没转北京时间。解决拿现货价和实际负荷做互相关看滞后几小时相关性最高。正常情况下现货价和负荷同周期峰值同步才对。对齐代码不复杂# 假设df已有spot和load两列 corr_shift {} for lag in range(-6, 7): corr_shift[lag] df[spot].shift(lag).corr(df[load]) best_lag max(corr_shift, keycorr_shift.get) print(f最佳滞后: {best_lag}个时点)正数表示现货价领先负荷负数表示落后。发现错位后用shift(best_lag)重贴时间标签而不是重新下载数据。5.4 需求响应系数是拍脑袋设的结果对弹性系数极其敏感现象把弹性系数从 0.2 改成 0.25最优电价曲线完全变形峰段时长从 4 小时变成 9 小时。原因需求响应模型里 (A_t) 和 (b_t) 是强耦合的(b_t) 微小的变化会被削峰约束放大。尤其是把 (b_t) 当成常数填进去而不是按各时段基线负荷归一化结果就更不稳定。解决先归一化再进入模型(b_t \eta \cdot A_t / p_{\text{ref}})让 (\eta) 成为一个无量纲的平均弹性。上线前对 (\eta) 做一遍 0.050.3 的灵敏度扫描看峰谷价差和削峰率的变化曲线。如果某个 (\eta) 区间内结果剧烈跳变说明模型在那个区域病态要回头查约束是否合理而不是硬找一个最优弹性系数。这个教训百试不爽——先跑灵敏度再谈参数调优。6. 验证与进阶用历史现货数据回测你的实时最优电价6.1 回测看三个数账单、削峰率、求解耗时模型做完不能直接上生产拿过去一个月的现货价格跑一遍离线回测。我只看三个数用户总账单变化、最大负荷削峰率、单次求解耗时。账单降太多说明电价可能压到亏损线了削峰率超过 15% 基本不现实求解耗时超过 5 秒就不适合做滚动实时计算。6.2 灵敏度扫描把玄学参数变成验收依据rho、alpha、(\eta) 三个参数每个取 35 档做网格扫描把削峰率、账单、峰谷价差三张曲线画出来。运营团队要问为什么这个方案可行你直接甩灵敏度图比解释一小时模型都管用。这也是把内点法从黑匣子变成可验收交付物的关键一步。6.3 一个实用技巧滚动窗口的定时更新策略实时最优电价不等于每秒都算常见做法是每天 16:00 用最新出清的现货价格重算一次次日 96 点电价遇到极端天气或电网阻塞事件临时追加一次重算。滚动窗口的好处是不依赖预测精度太高的光伏出力现货价本身就是最新市场信息。我现在每次跑完都会习惯性看一眼灵敏度曲线里有没有突变点——某个参数从 0.1 变到 0.15 时削峰率跳了 5 个百分点这种位置一定要搞清楚原因再发版。这套流程走顺之后实时最优电价就不再是实验室里的玩具而是售电公司每天都能用的定价工具了。希望帮到你。本文还有配套的精品资源点击获取