简介一份基于内点法求解电力市场实时最优电价的MATLAB实现资源面向电力系统优化、智能电网定价领域的工程师与研究人员。项目以30节点电力网络为算例完整呈现从系统建模、供需平衡约束到实时电价计算的算法流程适合用于研究分时电价TOU策略、理解内点法在非线性规划问题中的实际应用。压缩包共7个文件包含4个MATLAB源码文件主程序main.m、预测-校正模块mainaf.m、30节点输入数据inputcase30.m及初始化程序initial.m、2个fig格式迭代过程分析图以及1个说明文档整体仅40KB体量轻巧。已有119人学习下载。通过源码研读与图表分析读者可以掌握内点法求解实时电价的核心思路了解不同向心参数SIGMA及预测-校正环节对迭代收敛过程的影响为后续开展电力市场价格机制研究或算法改进提供可直接运行的基础代码。1. 现货电价不能直接用作分时电价实时最优电价为什么绕不开内点法现货市场spot market的电价每五分钟到每小时跳一次直接把现货价格透传给零售用户结算会乱用户也受不了这种波动。反过来按峰谷平拍脑袋三档又丢掉了现货价格里最有价值的激励信号。所以售电公司和电网调度侧常见的做法是把现货电价当作输入用内点法求解一组“实时最优电价”在购电成本、用户价格承受力、峰谷价差约束之间找最优折中。这条链路从现货曲线到分时费率核心求解器就是内点法。适合做售电套餐定价、微电网能量管理、分时电价政策测算的工程师。2. 搭建分时电价优化模型目标函数、决策变量与约束怎么落到算式分时电价优化的关键不是先选求解器而是先把“电价”从业务语言翻译成数学语言。这一步做不对后面用内点法再熟练也是白搭。我拿到一个现货电价背景的分时电价项目一般先问三个问题目标函数优化的是什么决策变量细分到什么时间粒度哪些约束必须硬性满足。下面按这个顺序展开。2.1 目标函数从购售电利润到社会福利选哪个更稳现货背景下的分时电价本质上是售电公司或电网在替用户承担现货价格波动风险。最常见的做法是把目标函数定为购售电利润最大化。设p_t是用户在t时段的零售电价s_t是该时段的现货购电成本q_t是该时段的用电量目标函数写为max Σ (p_t - s_t) * q_t这个式子很直观但实际项目里有个问题如果没有约束优化器会把价格推到用户接受上限所有时段价格全贴上界分时电价变成“天花板价”政策上没法看。所以我在实际项目中一般先用利润最大化形式但必须配合价格上下限、峰谷价差这类硬约束把电价拉回合理区间。另一种做法是用户福利最大化在目标函数里加一个用户用电效用的惩罚项把用户的“不满意度”也建模进去。目标变成更平滑解出来的价格也更温和但用户效用函数的系数很难标定一不小心就成了另一个玄学参数。我的做法是第一版先用利润最大化把用户响应写进需求函数跑通之后如果要考虑政策平缓性再切到福利最大化。这里的需求响应模型行业内最常用的是线性形式q_t q_base_t - ε_t * (p_t - p_ref_t)其中q_base_t是现货价格透明时该时段的自然负荷p_ref_t是用户预期的参考电价ε_t是需求价格弹性。把q_t代回利润目标函数后目标变成关于价格的二次函数这正好是后文内点法求解器最喜欢处理的凸二次规划形态。弹性系数ε_t的取值建议用历史日负荷与历史电价做回归拟合回归时把温度、星期几、节假日作为控制变量。没有历史数据时常见做法是取归一化弹性 0.10.3但要注意弹性在不同时段差异很大夏季午高峰工商业用户产线不能停弹性接近 0夜间谷段用户可转移负荷多弹性明显更大。我一般按峰、平、谷三段设不同的ε_t而不是整条曲线用一个值。2.2 决策变量与约束时段价格、峰谷价差和爬坡约束怎么设置决策变量是每个时段的零售电价。常见的时间粒度是 24 点小时级或 96 点15 分钟级。15 分钟粒度能更好跟随现货波动但约束数量和求解耗时也随之增加。对 96 点规模内点法完全能处理但模型维护难度上升我通常建议第一版先做 24 点跑通后再加密。约束方面分时电价模型至少要包含四类它们也是内点法不等式约束的主要来源第一类是价格上下限p_min ≤ p_t ≤ p_max。下限用来保护售电公司收益上限用来保护用户利益也是政策合规的硬要求。第二类是峰谷价差约束既要有上限防止峰谷价差过大也要有下限防止分时退化成单一电价。价差约束看起来简单做起来有个小坑max(p)-min(p)这个形式本身不可导内点法处理不了常见做法是把约束展开成任意两个时段之间的线性不等式p_i - p_j ≤ Δ_max对所有峰谷时段对生效。96 点规模下这种约束大约 4500 条对内点法来说仍然是线性代数的常规操作。第三类是爬坡约束|p_t - p_{t-1}| ≤ ramp。这是我个人很看重的一条。现货价格可以在半小时内跳涨三倍但零售电价不能跟着一步跳否则用户侧设备会集体响应形成新的负荷冲击。爬坡约束等价于两条线性不等式可以很好地压住相邻时段的价格突变。第四类是可选的总用电量守恒约束Σ q_t(p_t) Σ q_base_t它隐含地模拟了用户把负荷从高峰转移到低谷的行为比单纯加交叉弹性项更容易实现后面避坑章节会详细讲。这一组约束整理下来整个问题就是一个带线性不等式和线性等式约束的凸二次规划。求解这个问题的经典算法路线就是内点法。表 2.1 是常见约束的参数取值范围参考具体数值取决于当地现货价格水平和政策要求。约束名常见取值备注价格下限p_min0.30.5 倍现货均价保护售电公司收入价格上限p_max1.52.0 倍现货均价防止电价过高最大峰谷价差Δ_max0.81.2 倍现货均价时间粒度越细取值越小爬坡速率ramp0.20.5 倍现货均价每时段小时级取大值15 分钟级取小值2.3 内点法为什么比单纯形和罚函数法更适合电价问题分时电价优化用内点法不是因为它听起来高级而是问题结构决定的。早期做电价优化有人用单纯形法但单纯形法只能处理线性目标需求响应写进目标函数后问题变成二次型单纯形就不适用了。也有人用罚函数法处理不等式约束把约束乘一个大系数加到目标函数里做法简单但数值上容易振荡惩罚系数小了约束不满足大了目标函数变形解出来的电价在边界上来回跳。内点法处理不等式约束的思路完全不同。它把每个不等式约束变成对数障碍项加进目标函数例如- (1/μ) * Σ log(-g_i(x))其中g_i(x) ≤ 0是原约束μ是障碍参数。迭代点被障碍函数“推”在可行域内部μ从大到小逐渐逼近 0解沿中心路径收敛到原问题的最优点。这个机制的直观效果是迭代过程稳定不会像罚函数那样在约束边界来回弹跳而且对约束数量和维度的扩张非常友好。实际使用中不需要自己写障碍函数CVXOPT 的二次规划求解器内部就是大尺度内点法只需要把模型整理成标准的二次规划格式。但理解内点法的这条路径有一个实际价值当看到求解器输出“对偶间隙不收敛”或“迭代数超过上限”时能判断是μ衰减策略的问题还是模型本身病态而不是盲目调大迭代次数。这是调优化问题的重要能力。3. 用内点法求解实时最优电价CVXOPT 标准型与 Python 实现模型搭好之后下一步是把公式落成可以跑的代码。我习惯用 CVXOPT因为它输入格式明确内部实现的内点法比较成熟不需要去手写 KKT 迭代。但 CVXOPT 要求先把模型写成它规定的标准型才能调用这一步是新手最容易卡住的地方。3.1 把利润目标展开成二次规划标准型CVXOPT 二次规划求解器solvers.qp解决的是如下形式的问题min 0.5 * x^T P x q^T x约束: Gx ≤ h, Ax b所以要把目标函数max Σ(p_t - s_t) * q_t变成这个格式。把q_t q_base_t - ε_t*(p_t-p_ref_t)代入并展开先写成最小化形式对每个时段 t目标项展开为ε_t * p_t^2 - p_t * (q_base_t ε_t*p_ref_t ε_t*s_t)加上一个与p_t无关的常数项。常数项不影响最优解直接丢弃。因此二次项矩阵P的主对角元素是2*ε_t一次项向量q_lin的第 t 个元素是-(q_base_t ε_t*(p_ref_t s_t))。这里有个容易错的系数0.5 * x^T P x里已经带了 0.5所以P对角元要设为2*ε_t而不是ε_t。我见过不止一次因为这个系数少乘 2导致求出来的最优价格和手工推导结果对不上。线性约束部分相对直接价格上下限、峰谷价差、爬坡约束都可以写成Gx ≤ h的行形式具体见下面的代码。3.2 基于 CVXOPT 内点法的 Python 实现与参数说明为了便于阅读我用 8 个时段的简化数据演示完整求解流程扩展到 96 点只需要把矩阵规模改大代码结构不变。先准备输入数据。import numpy as np from cvxopt import matrix, solvers T 8 # 时段数演示用 8 点实际可扩到 24 或 96 # 现货价格 s_t单位元/MWh spot np.array([120., 110., 105., 150., 280., 350., 300., 180.]) # 基础负荷 q_base_t单位MWh q_base np.array([55., 50., 48., 52., 78., 89., 82., 65.]) # 需求弹性峰时段小午高峰最大负荷对应时段弹性最低 eps np.array([0.25, 0.25, 0.25, 0.15, 0.08, 0.05, 0.08, 0.15]) # 参考电价 p_ref_t通常取现货当日均价或历史零售均价 p_ref np.full(T, spot.mean()) # 价格上下限、峰谷价差上限、爬坡速率 p_min, p_max 0.4 * spot.mean(), 1.6 * spot.mean() delta_max 0.9 * spot.mean() ramp 0.3 * spot.mean()这段代码定义了问题的输入。spot和q_base是当日现货电价和预测负荷eps是需求弹性。注意弹性较小的时段对应午高峰这符合工商业负荷对电价迟钝的真实特征。p_ref用现货均价是做第一版时的简化实际项目中应该用用户当前的零售电价因为那是用户已经形成的价格心理锚点。接下来构造二次规划矩阵并求解。# 目标函数二次项0.5 * x^T P xP 为对角矩阵 P matrix(np.diag(2.0 * eps).astype(float)) # 目标函数一次项q_lin[t] -(q_base_t eps_t * (p_ref_t spot_t)) q_lin matrix(-(q_base eps * (p_ref spot)).astype(float)) # 约束矩阵 G 和 h构造 Gx h G_rows, h_vals [], [] # 1) 价格上下限p_min x_i p_max for i in range(T): row np.zeros(T); row[i] -1.0 G_rows.append(row); h_vals.append(-p_min) row np.zeros(T); row[i] 1.0 G_rows.append(row); h_vals.append(p_max) # 2) 峰谷价差上限p_i - p_j delta_max # 峰时段集合与谷时段集合由现货价格排序或业务规则决定 peak_idx [4, 5, 6] # 高价时段 valley_idx [0, 1, 2] # 低价时段 for i in peak_idx: for j in valley_idx: row np.zeros(T); row[i] 1.0; row[j] -1.0 G_rows.append(row); h_vals.append(delta_max) # 3) 爬坡约束|p_t - p_{t-1}| ramp for t in range(1, T): row np.zeros(T); row[t] 1.0; row[t-1] -1.0 G_rows.append(row); h_vals.append(ramp) row np.zeros(T); row[t] -1.0; row[t-1] 1.0 G_rows.append(row); h_vals.append(ramp) G matrix(np.array(G_rows).astype(float)) h matrix(np.array(h_vals).astype(float)) # 内点法求解器参数 solvers.options[show_progress] False solvers.options[maxiters] 100 solvers.options[feastol] 1e-6 # 求解凸二次规划 sol solvers.qp(P, q_lin, G, h) p_opt np.array(sol[x]).reshape(-1) print(最优分时电价(元/MWh):, np.round(p_opt, 2)) print(迭代次数:, sol[iterations]) print(求解状态:, sol[status])代码里的核心逻辑是三步。第一步构建P和q_lin对应 3.1 节推出的目标函数标准型第二步把四类约束按行组装进G和h第三步调用求解器。这里解释几个关键参数solvers.options[maxiters]100是内点法最大迭代次数8 点规模通常 10 次以内收敛96 点规模需要 50 次左右100 是安全上限feastol是 KKT 条件里可行性残差的容忍度默认 1e-7如果模型约束多、量纲差异大放宽到 1e-5 往往比调大迭代次数更有效。sol返回字典里sol[x]是最优解向量sol[status]是求解状态。状态为optimal表示正常收敛状态为unknown意味着数值问题或迭代超限这时不要急着改代码先检查量纲和约束是否有冲突。这一点后面避坑章会展开。3.3 内点法关键参数障碍系数、收敛阈值与迭代上限怎么定用 CVXOPT 时并不直接设置障碍系数μ但理解μ的衰减机制对判断求解行为很有帮助。内点法在每轮迭代中会同时缩小障碍参数和推进牛顿方向μ减太快对数障碍项把迭代路径“推”向边界的力过猛步长会变小甚至振荡μ减太慢迭代次数明显增加。成熟求解器会自适应调节μ所以手动造轮子时最容易翻车的地方其实是三个参数的选择。第一个是可行性容忍度feastol。它控制的是约束违反程度的可接受水平分时电价这类工程问题并不需要机器精度1e-4到1e-6足够。设成1e-9会让内点法在最后几轮迭代里反复修整浪费计算时间。第二个是对偶间隙容忍度abstol/relatol它控制最优解和理论最优值的距离。电价连续量保留两位小数就已经很精确1e-4是合理选择。第三个是maxiters。很多新手一看到unknown就把它调到 500但 96 点模型迭代 150 步还不收敛多半是模型病态或约束冲突不是迭代次数不够。手写内点法时才需要直接选择障碍系数初值和衰减策略常见做法是初值μ1每步乘以 0.1 到 0.5 的衰减因子步长则用线搜索保证目标函数单调下降。如果后续你想在实时系统里绕过 CVXOPT自己实现一个轻量内点法可以把μ的衰减过程打印出来观察对偶间隙是否稳定下降。只要对偶间隙在连续下降说明路径没问题一旦它停滞或反弹问题大概率出在矩阵条件数上也就是下一章要讲的量纲问题。4. 分时电价优化避坑指南四类翻车现场与参数修正内点法不是调通一次就一劳永逸的。现货输入一换、负荷曲线一换后面几个坑几乎必踩。我把最常遇到的四类问题按“现象 → 原因 → 解决”的格式写出来每一条都来自实际调试经验。4.1 解出来的价格全在下限目标函数符号与矩阵正定性现象内点法正常返回statusoptimal但p_opt全部等于p_min分时电价变成“白菜价”售电公司收益为负。原因通常不是内点法的问题而是目标函数符号写反了。利润最大化max Σ(p_t-s_t)q_t转化为最小化形式时原目标是凹函数取负后应变成凸函数二次项系数应该为正。如果在这里漏了变号或者P矩阵对角元直接用了-ε_t那么目标函数仍然凹内点法沿下降方向找点时会把解推向约束边界价格全部抵到下界。解决先用小例子做手推验证。检查P的最小特征值是否大于 0用np.linalg.eigvalsh(np.diag(2*eps))看第一列结果如果出现负数或接近 0 的数值回查 3.1 节的推导。另外打印q_lin的符号正常情况下它应该是负数因为q_base_t ε_t*(p_ref_ts_t)为正。另一个隐蔽的同类问题是eps值本身为负这会让整个模型的物理意义反转表现为价格越高负荷越大怎么看怎么不对。4.2 分时退化为单一电价缺少最小峰谷价差约束现象内点法收敛但p_opt在 24 个时段里几乎都是同一个值分时电价名存实亡。原因目标函数里的需求响应项天然会把价格向现货价格的中枢拉。如果ε_t设得偏大优化器发现“抬高峰时段电价的结果是用户不用电利润反而下降”于是一方面压低峰价一方面抬高谷价最后殊途同归。这个现象在现货价格本身波动较小的日子里尤其明显。解决加最小峰谷价差约束。3.2 节代码里只有p_i-p_j delta_max这是上限约束要避免退化成单一电价还需要一组下限约束p_i-p_j delta_min写法是反过来在G里加一行j位置为 1、i位置为 -1h值设为-delta_min。同时注意epsilon的取值范围要和价差约束匹配。如果弹性系数标定得过大比如归一化后超过 0.8优化器在数学上就有动机拉平价格这不是求解器毛病而是模型参数不合理。建议弹性系数的第一版取值不要超过 0.3调试稳定后再逐步放宽。4.3 削峰填谷变成峰谷倒挂忽略跨时段替代弹性现象优化后的电价确实是“峰时段贵、谷时段便宜”但用户负荷曲线反而出现新的晚高峰整体峰谷差没有变小甚至在原本的低谷时段出现了二次尖峰。原因线性需求响应模型q_t q_base_t - ε_t*(p_t-p_ref_t)只表示同一个时段内价格对负荷的抑制效应这在经济学上叫自价格弹性。但真实用户会在不同时段之间做用电替代比如电动汽车充电从晚高峰挪到深夜热水器错峰启动。这部分行为需要交叉价格弹性来描述如果模型里完全没有跨时段的负荷迁移优化的价格会诱导用户把一部分负荷从峰时段挪到紧挨着的平段甚至谷段结果谷段被填出一个新的峰来。解决最省事的办法是加总用电量守恒约束Σ q_t(p_t) Σ q_base_t。把这个等式展开后得到Σ ε_t * p_t Σ ε_t * p_ref_t这是一条线性等式约束在 CVXOPT 里用A和b参数传入即可。加了这条约束后系统不再可能把所有时段价格都压低来规避需求响应损失跨时段的负荷转移被隐含地考虑进来。条件允许的情况下更精细的做法是显式加入交叉弹性矩阵把q_t写成对所有时段价格的线性函数代价是约束规模从 T 涨到 T*T但 96 点规模下仍然可行。这一条踩坑记录提醒我们内点法只能保证数学上最优不能保证物理上合理模型结构的完整度要靠工程师把关。4.4 内点法迭代震荡不收敛量纲与初值问题现象扩展到 96 点、加上峰谷价差和爬坡约束后求解器返回statusunknown迭代在 80 到 100 步之间反复震荡sol[iterations]显示已经跑到maxiters上限。原因这是数值问题。现货价格s_t在 100 到 1000 元/MWh 量级而需求弹性ε_t可能是 0.0001 级别如果单位没归一化的话P矩阵元素和G矩阵元素的量纲差距能到 1e7。内点法每一步都要解一个牛顿线性方程组矩阵条件数这么差数值误差会直接放大KKT 残差很难降到容忍度以内。解决全部变量先归一化。常见做法是让价格除以当日现货均价得到标幺值负荷除以日最大负荷弹性系数也按标幺值体系重新标定。求解得到归一化最优电价后再乘回基准价格。归一化后P、G、h的元素量纲都落在 0.1 到 10 之间矩阵条件数下降几个数量级迭代次数通常能减少一半。我在 96 点模型里遇到过一次迭代 200 步不收敛的问题归一化后只迭代了 34 步就拿到optimal状态。另外如果发现feastol从 1e-6 放宽到 1e-4 后收敛但价格结果与业务预期偏差大不要强行用更严格的容差而是先检查是不是模型约束漏了这是另一个维度的坑。5. 从离线测算到实时下发滚动优化的初值策略与验证方法实时最优电价不是解一次就完事的。现货价格每小时在变负荷预测也在滚动更新实际系统里普遍采用滚动优化每 15 分钟或每 1 小时重新求解一次只把当前窗口里最近一个时段的电价真正下发。这样能让分时电价始终跟随最新现货信息同时避免一次优化结果被后续扰动彻底击穿。滚动优化里有个非常实用的技巧是热启动。把上一次求解得到的p_opt作为下一次迭代的初始点内点法在可行域内部沿中心路径找解热启动能省掉从任意初值逼近中心路径的步骤。CVXOPT 的qp接口没有直接传x0的参数我一般用两个替代方案一是把滚动窗口缩小每步只求解 4 个时段每步构建矩阵的开销很小二是换用支持初值参数的求解器。实测里冷启动的内点法迭代 20 到 30 次热启动能降到 5 到 10 次实时性差距明显。接入真实系统之前一定要做回放验证。我的习惯是取历史两周的现货价格和负荷数据前 14 天做弹性系数回归第 15 天做优化回放对比基线方案把现货均价直接切成峰谷平三段的购电成本和峰谷差。回放不通过就不要接实时系统。特别要提醒的是回放数据要覆盖价格突变的样例比如现货价格跳涨三倍的日子否则模型里少一个价格突变抑制逻辑上线第一天就会翻车。我在第一次做实时电价时跳过回放直接上了 96 点模型结果在现货跳涨那天优化器给出了一个负电价调度侧直接拒绝执行。后来排查发现是现货输入里有一个极端值没做预处理而这天的数据根本没出现在回放集里。从那以后我在数据预处理里加了价格突变抑制超过上一时刻三倍的现货价格先做截断再进入优化器这个逻辑也成了模型输入检测的一部分。最后提醒一句上线第一周最好保留手工覆盖开关把内点法算出的电价与现货曲线并排画出来用肉眼校一遍趋势。这套从模型到代码再到回放的流程能帮你把分时电价优化做成一个可维护的方案希望帮到你少走点弯路。本文还有配套的精品资源点击获取