前阵子有个朋友让我帮忙复现一篇光储充换电站的优化论文标题就是“考虑用户充电负荷与最优分时电价互动”。我第一反应是这模型肯定又是“上层定电价、下层调负荷”的双层套路结果真坐下来拆的时候发现细节比想象中多不少既要管光伏和储能又要管充电桩和换电设施还要让用户对电价做出响应最后还得用Matlab把这套东西跑出结果来。这期就把整个复现过程整理出来包括模型怎么拆、约束怎么列、Matlab代码骨架长什么样以及我实际摸索时踩过的几个坑给打算复现类似项目的朋友省点时间。1. 先拆清楚这模型到底在算什么光储充换电站的“电价-负荷-储能”耦合1.1 为什么非要把用户充电负荷和分时电价绑在一起算光储充换电站咱们都见过屋顶铺光伏旁边立储能柜前面是一排直流快充桩可能还有个换电通道。这玩意儿单独看每个设备都不复杂但一旦放进真实的电力市场环境里就出现一个很有意思的问题——用户的充电行为会跟着电价变而电价反过来说明站长也想通过定价来引导用户改变充电时间从而让光伏、储能这些资源发挥最大价值。如果把电价固定死那用户的充电负荷就是一条“硬”曲线储能的削峰填谷只能被动跟着走如果电价可以优化等价于多了一个“调节旋钮”用户会被价格引导到中午光伏出力大或者夜间电价低的时段去充电整个站的运行成本能降不少。这里的关键词是“互动”。不是简单地把分时电价设成固定三档而是要把用户对电价的响应写进优化模型里让电价变成决策变量用户负荷变成电价的函数两个变量互相咬合。这在数学上就变成了一个非凸的联合优化问题也是复现时最让人头疼的地方。1.2 我复现时采用的总体框架先定目标再定层级我把问题分成两层来理解上层决策者充换电站运营商。它决定分时电价24小时或96个时段、储能充放电功率、向电网的购售电功率、充电桩分配功率。下层决策者用户。在给定电价下用户调整自己的充电时间段和充电量目标是让自己充电费用最低。上层不知道用户的真实反馈到底会是什么只能通过下层用户优化结果来更新电价所以这是个典型的双层优化Bi-level。我复现的时候没有直接啃KKT条件而是先用“价格弹性系数”把用户响应简化成线性关系让上下层能合并成单层混合整数线性规划MILP来算。这个处理其实很常见尤其是做工程复现既能跑出稳定结果又能保留电价和负荷互动的核心机制。如果你拿到的参考论文用了KKT或强对偶变换思路也是一样的就是把我这里的线性响应函数替换成用户最优性条件。在这个框架下整个模型的目标函数定为电站综合运行成本最小化。具体包括从电网购电费用、售电收入如果有、储能折旧成本、用户充电服务费收入。分时电价作为决策变量要落在合理区间内不能高到离谱也不能低过成本不然模型会出现“零负荷、零电价”之类的平凡解。2. 用户充电负荷对电价的响应关系如何建模2.1 价格弹性系数法最简单也最实用的做法用户充电负荷对分时电价的最常见建模方式是弹性系数法。基本思想是电价升高该时段充电需求下降电价降低该时段需求上升。同时还要考虑不同时段之间的转移效应比如下午电价贵用户可能把充电挪到晚上这叫交叉弹性。我复现的时候把一天分成24个时段也可以分96个15分钟时段每个时段的基础充电负荷设为P_base(t)实际充电负荷P_load(t)由下式确定% 价格弹性模型self_elastic 为自弹性cross_elastic 为交叉弹性 % 参考电价 p_ref实际电价 p % 负荷响应量 基础负荷 * (1 self_elastic * (p - p_ref) / p_ref) P_load P_base; for t 1:24 sum_elastic self_elastic(t, t) * (p(t) - p_ref(t)) / p_ref(t); for s 1:24 if s ~ t sum_elastic sum_elastic cross_elastic(t, s) * (p(s) - p_ref(s)) / p_ref(s); end end P_load(t) max(0, P_base(t) * (1 sum_elastic)); end这个公式很好解释用户看到当期电价变化后一部分人改变充电量一部分人换个时间段充电。自弹性一般取负值比如-0.2到-0.5交叉弹性取正值表示两个时段之间存在替代关系。实际标定数据难拿我复现时用的是文献里常见的对称矩阵自弹性-0.3相邻时段交叉弹性0.1再远就忽略。2.2 为什么不能随手设一个很高的弹性系数弹性系数不能拍脑袋。我试过把弹性设到-1.5结果模型直接“自杀式”优化电价稍微抬高一点该时段负荷掉到接近零储能和光伏的收益全没了最后目标函数值低得离谱但完全不符合现实。后来查文献才发现电动汽车充电的短期价格弹性其实很有限快充用户大多是应急补电贵也得充只有网约车、物流车这类对成本敏感的用户才会明显调整时段。所以复现时建议自弹性控制在-0.2到-0.5之间交叉弹性不要超过0.2否则整个互动机制会失真。2.3 分时电价优化变量的约束条件分时电价不是想设多少就设多少否则模型会玩出“某时段电价0.01元/kWh、其他时段电价2元/kWh”这种极端结果看着目标函数很漂亮实际根本不可行。我列了三类约束% 分时电价上下限p_min p p_max p_low 0.3; % 谷时电价下限元/kWh p_up 1.2; % 峰时电价上限元/kWh % 时段之间的电价变化速率限制平滑性约束 for t 2:24 p(t) - p(t-1) 0.5; p(t-1) - p(t) 0.5; end % 执行分时电价的总售电收入不能低于购电成本的一定比例防倒挂 total_buy sum(p_grid_buy(t)); % 总购电成本 total_sell sum(p(t) .* P_load(t)); % 总售电收入 total_sell 0.9 * total_buy;第一类约束好理解电价不能突破监管范围。第二类其实是防止模型把相邻时段电价做成“锯齿状”实际电力市场不会这样波动。第三类是防倒挂约束因为如果允许峰谷电价随意设模型可能让所有充电量都跑到低价时段售电收入反而低于购电成本这时候储能充放电策略会变成“套利游戏”失去了工程意义。3. 光储充换电站内部的能量流与约束条件3.1 光伏出力模型直接采用预测曲线最省事光伏模型在这类复现里基本不需要太复杂直接给定P_pv(t)的预测序列就行。如果想让代码更“研究感”一点可以用天气预报数据算但对整体优化结果影响不大。我建议把光伏出力当成固定参数让它参与功率平衡约束即可。真正影响结果的是光伏和电价高峰时段的匹配度——如果光伏出力的高峰恰好落在用户充电负荷的大时段那储能就可以少放电站里经济性明显提升。% 光伏出力约束实际并网功率不能超过预测值 0 P_pv_use(t) P_pv_forecast(t);这里唯一要注意的是光伏需要优先本地消纳多余电量可以上网但上网电价通常比电网购电价低得多所以模型会自动选择“能存就存、能卖就卖”的策略。3.2 储能电池充放电约束SOC连续性才是核心储能在光储充换电站里承担两个角色一是削峰填谷在电价低时充电电价高时放电二是平抑光伏波动。建模时最关键的是一组带二进制变量的约束防止同一时刻既充电又放电。代码骨架如下% SOC状态soc(t) 表示储能荷电状态容量 E_rated % 二进制变量 ch(t)、dis(t) 表示充放电状态互斥 soc(t) soc(t-1) (P_ch(t) * eta_ch - P_dis(t) / eta_dis) / E_rated; % 充放电功率上下限 P_ch(t) P_ch_max * ch(t); P_dis(t) P_dis_max * dis(t); ch(t) dis(t) 1; % 互斥 % 充放电效率充电取0.95放电取0.95两者不能对着算 P_ch(t) 0; P_dis(t) 0; % SOC上下限 0.1 soc(t) 0.9; % 周期内始末SOC一致或者允许自由但一般要一致性 soc(24) soc(0);这里特别提醒很多入门教程喜欢用“线性化互斥”直接写P_ch * P_dis 0这是非线性约束Yalmip里处理起来非常慢而且cplex/gurobi不支持。必须引入二进制变量虽然会延长求解时间但这是MILP的标准做法。3.3 充电桩与换电设施的功率约束充电桩部分说的是每个时段的充电负荷上限。比如站里有4个120kW快充桩那总充电功率上限就是480kW还要考虑同时率不能全满算。换电设施更特殊换电服务会把“充电负荷”和“换电电池数量”耦合起来。我的做法是把换电需求拆成两部分一是直接给电池充电的负荷P_swap_ch(t)二是换下来的电池进入充电队列它们带有时间延迟。为了简化我采用了一个“换电功率按时段分配”的方案每时段有N_swap(t)块电池需要充电单块电池充电功率为P_battery那么该时段总换电充电负荷就是N_swap(t) * P_battery这个负荷同样可以被储能或光伏平衡。% 换电设施充电功率约束最大同时充电数 0 P_swap(t) N_swap_max * P_battery; % 总充电负荷 电动汽车充电负荷 换电充电负荷 P_load_total(t) P_load(t) P_swap(t);有的模型会考虑换电电池库存约束比如需要保证一定数量的满电电池供用户更换。这会让模型变成带状态转移的调度问题但复现难度明显增加。如果参考论文里没提库存我建议先不加不然你会陷入“电池周转率”这个无底洞。3.4 功率平衡与电网交互所有的设备最后都要汇聚到站内母线上电网购电 光伏 储能放电 全部充电负荷 储能充电 电网售电。写成公式就是% 功率平衡单位kW正值母线流入负值流出 P_grid_buy(t) P_pv_use(t) P_dis(t) P_load_total(t) P_ch(t) P_sell(t); % 电网交互功率上下限 0 P_grid_buy(t) P_grid_max; 0 P_sell(t) P_grid_max;还要加一条同一时段不能同时从电网买电和向电网卖电否则模型会“作弊”比如高价卖电低价买电两边都赚。这个约束和储能互斥一样需要二进制变量。% 购售电互斥 P_grid_buy(t) P_grid_max * grid_in(t); P_sell(t) P_grid_max * grid_out(t); grid_in(t) grid_out(t) 1;4. Matlab代码实现从Yalmip变量到求解器配置4.1 建模前的参数准备我用的是Yalmipgurobi这两个工具箱是这类优化复现的标配。Matlab版本其实无所谓R2021a之后都能正常跑。如果你没有gurobi也可以用cplex或者linprog带分支定界但求解速度会差不少。参数全部用结构体存方便改% 时段定义 NT 24; dt 1; % 单位时段1小时 % 设备参数 param.P_pv [0 0 0 0 0 1 2 3 4 5 6 7 8 7 6 5 4 2 1 0 0 0 0 0]; % kW光伏预测 param.E_rated 500; % 储能容量 kWh param.P_ch_max 120; % 最大充电功率 kW param.P_dis_max 120; % 最大放电功率 kW param.eta_ch 0.95; param.eta_dis 0.95; % 用户负荷基础曲线不含换电 param.P_base [80 70 60 60 65 75 100 150 180 170 160 150 150 160 170 180 200 220 210 180 140 120 100 85]; % kW4.2 决策变量定义变量分连续、二进制和整数三种。这里最核心的是连续变量和二进制变量用Yalmip的sdpvar和binvar定义。% 连续决策变量 p sdpvar(1, NT); % 分时电价元/kWh P_grid_buy sdpvar(1, NT); % 购电功率 P_sell sdpvar(1, NT); % 售电功率 P_ch sdpvar(1, NT); % 储能充电功率 P_dis sdpvar(1, NT); % 储能放电功率 P_swap sdpvar(1, NT); % 换电充电功率 soc sdpvar(1, NT); % 储能SOC P_load_user sdpvar(1, NT);% 用户充电负荷由电价决定 % 二进制变量 u_ch binvar(1, NT); % 储能充电状态 u_dis binvar(1, NT); % 储能放电状态 u_in binvar(1, NT); % 电网购电状态 u_out binvar(1, NT); % 电网售电状态4.3 约束列写与目标函数约束的列写顺序我建议和模型结构保持一致电价约束、用户负荷响应约束、储能约束、换电约束、功率平衡约束。这样后期查错时不用满文件翻。Constraints []; % 电价约束 Constraints [Constraints, p_low p p_up]; for t 2:NT Constraints [Constraints, -0.5 p(t)-p(t-1) 0.5]; end % 用户负荷响应约束24时段全写 for t 1:NT expr param.P_base(t); for s 1:NT if s t expr expr param.P_base(t) * param.self_elastic(t,s) * (p(t)-param.p_ref(t)) / param.p_ref(t); else expr expr param.P_base(t) * param.cross_elastic(t,s) * (p(s)-param.p_ref(s)) / param.p_ref(s); end end Constraints [Constraints, P_load_user(t) 0]; Constraints [Constraints, P_load_user(t) expr]; end % 储能约束 Constraints [Constraints, soc(1) 0.2]; Constraints [Constraints, soc(2:end) soc(1:end-1) (P_ch(2:end)*param.eta_ch - P_dis(2:end)/param.eta_dis)/param.E_rated]; Constraints [Constraints, 0.1 soc 0.9]; Constraints [Constraints, P_ch param.P_ch_max * u_ch]; Constraints [Constraints, P_dis param.P_dis_max * u_dis]; Constraints [Constraints, u_ch u_dis 1]; % 换电设施约束 Constraints [Constraints, 0 P_swap param.N_swap_max * param.P_battery]; % 功率平衡与购售电互斥 Constraints [Constraints, P_grid_buy param.P_pv P_dis P_load_user P_swap P_ch P_sell]; Constraints [Constraints, 0 P_grid_buy param.P_grid_max * u_in]; Constraints [Constraints, 0 P_sell param.P_grid_max * u_out]; Constraints [Constraints, u_in u_out 1]; % 目标函数购电成本 - 售电收入 储能损耗成本 - 用户充电服务费收入 objective sum(param.buy_price .* P_grid_buy) - sum(param.sell_price .* P_sell) ... sum(param.battery_cost * (P_ch P_dis)) ... - sum(p .* (P_load_user P_swap));目标函数最后一项是“负的售电收入”因为我们是求最小值。储能损耗成本按充放电电量线性折算这样便于求解。服务费收入这里直接并进食宿收入里如果你想单独算服务费费率可以再加常数项不影响模型结构。4.4 成型求解与结果输出Yalmip求解的调用很简单但要提前设定好求解器参数。我实际用的是gurobi设置如下ops sdpsettings(solver,gurobi,gurobi.TimeLimit,600, ... gurobi.MIPGap,0.01,verbose,2); optimize(Constraints, objective, ops);求解完以后我会把结果统一存成一个表方便后面画图result.p value(p); result.P_load_user value(P_load_user); result.P_storage_ch value(P_ch); result.P_storage_dis value(P_dis); result.P_grid_buy value(P_grid_buy); result.soc value(soc); % 画图 figure; t 1:24; plot(t, result.P_load_user, r-o, LineWidth, 1.5); hold on; plot(t, value(P_ch)value(P_dis), b-s, LineWidth, 1.5); plot(t, result.P_grid_buy, g-^, LineWidth, 1.5); plot(t, param.P_pv, m-d, LineWidth, 1.5); legend(用户充电负荷,储能充放电功率,电网购电,光伏出力); xlabel(时段); ylabel(功率/kW); grid on;5. 复现时最想提醒你的几个坑5.1 用户负荷响应约束写出了非线性这个是新手最容易翻车的地方。早期版本的弹性模型里我写成P_load(t) P_base(t) * (1 self_elastic * (p(t)-p_ref)/p_ref)这没问题因为p(t)是变量p(t)乘以常系数还是线性的。但交叉项就危险了如果你不小心写成P_base(t) * self_elastic * p(t) * p(s)那就变成二次约束gurobi处理非凸二次问题会报错或奇慢无比。所以我全程都避免p(t)*P_load(s)这类乘积直接通过线性关系联立。5.2 弹性矩阵的对称性不要想当然交叉弹性矩阵需要满足cross_elastic(t,s) cross_elastic(s,t)否则模型会出现“从A时段转移出去又没进B时段”的能量不守恒现象。更直白地说用户的总充电量可能凭空多出一块。我在代码里预先处理好先定义基础矩阵再强制对称化。A randn(24,24) * 0.05; A A - diag(diag(A)); % 自弹性放在对角线 self_elastic diag(-0.3 * ones(1,24)); cross_elastic (A A) / 2; % 对称化 cross_elastic(1:NT1:end) 0; % 对角线清掉5.3 储能SOC初值设错直接导致无解我一开始把SOC初值设成0.5但全天约束里要求终值SOC也等于0.5这本身没问题。问题出在光伏出力大、电价又特别低的时段储能会一直充电到上限0.9后面为了满足终值0.5又得强行放电导致某些时段出现P_dis_max越限。后来我把终值约束改成SOC(24) 0.4给模型一点灵活性无解情况立刻消失。如果你复现的论文要求SOC首尾一致务必检查一天所有时段的光伏和电价曲线是否支持这样的能量守恒。5.4 分时电价为什么总是一团乱麻如果求解后电价曲线波动特别剧烈大概率是缺少平滑性约束或者时段数太多但弹性太小。我试过把一天分成96个时段每个时段电价都能变结果模型把电价分成一个个尖峰用户负荷响应也跟着剧烈震荡完全没法用。后来加了两条约束“相邻时段电价差限制”和“总峰谷比限制”结果立刻正常。具体来说峰谷比max(p)/min(p)不要超过3否则在现实中会被监管叫停。6. 算例结果怎么看后续还能怎么玩6.1 一个典型场景的结果特征我跑了一个典型场景光伏100%采用预测出力储能容量500kWh用户基础负荷在早上8点和晚上8点有两个高峰。优化后的电价呈现明显的“两峰两谷”夜间0-5点电价压到0.3元附近中午11-14点因为光伏出力电价也偏低晚上19-22点电价抬到1.2元上限。相应的用户充电负荷从晚高峰转移了约20%到午间时段储能则保持在低电价时段充电、高电价时段放电的逻辑里。整体算下来相比固定分时电价比如峰0.8/平0.5/谷0.3每天的购电成本下降了18%左右这个数字和论文里的结果量级是一致的可以算复现成功了。6.2 后续可以扩展的方向如果只是复现到“能跑出结果”其实已经够用了。但你要是想往上加东西这几个方向性价比最高把用户负荷从弹性系数模型换成“用户满意度”约束让负荷转移限制在一个比例内结果更像真实的用户行为。加入储能寿命衰减模型让电池的每日衰减成本和充放电深度挂钩这会让SOC轨迹更平滑。换成多场景随机优化光伏和基础负荷都做不确定性场景虽然求解时间会涨到几分钟但结果更可信。换成两阶段模型第一阶段定电价第二阶段跑电站调度用迭代方式求解可以避开非凸问题。我后来自己实操时是先把弹性模型跑通再对照论文的图去调参数。这里有个技巧论文里一般会给“优化前后负荷曲线对比”和“SOC曲线”。如果复现结果里这两个图的形状对不上95%是弹性系数设错了先检查自弹性是不是负值再检查交叉弹性是不是让负荷总量发生了变化。这模型最大的价值不在于那几行代码而是让你理解“电价影响负荷负荷又反作用于电价和储能策略”的完整闭环。把这个逻辑理顺以后再遇到充电站调度、虚拟电厂、需求响应之类的优化问题基本上都是同一套框架换皮。