前段时间我在做微电网日前经济调度相关课题的时候把风电、光伏、储能和需求响应全部塞进了同一个24小时优化模型里用Matlab完成建模和求解。说实话刚接触这类题目时很多人会觉得无非就是列约束、调求解器但真正把所有环节的时间耦合处理好尤其是需求响应的两种典型机制和储能充放电的二进制互斥变量一起放进优化里坑远比想象中多。这篇文章把我这版Matlab实现的完整思路、建模细节、代码骨架和踩坑记录整理出来给正在做微电网优化、综合能源调度或者想给毕业设计加一点需求响应亮点的同学做一个可以直接参考的脚手架。1. 调度问题的规模感先想清楚这个模型里有哪些单元1.1 微电网拓扑与潮流简化假设我搭的参考模型是一个典型的单母线交流微电网包含一台微型燃气轮机、一个风电场、一个光伏电站、一组电池储能系统、与主网的联络线以及一块带需求响应能力的本地负荷。拓扑结构虽然看着有位置远近之分但在日前经济调度这个时间尺度下我做了一个很重要的简化假设忽略线路潮流和节点电压约束把整个微电网视为一个功率平衡节点。这个假设不是偷懒而是日前调度问题的核心矛盾在于24小时时间轴上“什么时候充、什么时候放、什么时候买电、什么时候削减负荷”空间维度的潮流约束通常放在日内滚动优化或者潮流校验阶段再处理。单母线模型足以回答最核心的问题储能参与峰谷套利的充放电策略是什么、需求响应在哪些时段被调用、微型燃机和购电之间的经济分配怎么走。1.2 时间尺度为什么是24小时、1小时一个点日前调度的时间尺度取24小时步长取1小时这个几乎是行业惯例。更短的时间尺度比如15分钟会大幅增加变量数量让混合整数规划求解速度变差更长的时间尺度又无法刻画负荷早晚高峰和分时电价的切换。24个时段在这种规模下决策变量大概在一百到两百个之间加上二进制变量Gurobi或CPLEX这类求解器基本秒内就能出结果非常适合做算例研究和毕设。在这个时间尺度下所有功率变量都默认是每小时的恒定平均功率单位用MW或kW能量用MWh或kWh储能SOC递推式里的电量变化就等于功率乘以1小时。理解和把握好这个离散化粒度后面写约束时就不会把量纲搞乱。2. 风光出力和储能约束最容易写错的地方都在这2.1 风电、光伏作为“可削减电源”的建模风电和光伏在日前调度里最典型的处理方式是把预测出力曲线作为调度时段内的最大可用出力实际调度出力可以在0到预测值之间任意选择也就是说允许弃风弃光。我这里说的“任意选择”在代码里体现为一条上下限约束P_w sdpvar(1, T); % 风电调度出力 P_pv sdpvar(1, T); % 光伏调度出力 C [C, 0 P_w P_w_forecast]; C [C, 0 P_pv P_pv_forecast];这里有个很关键的细节如果目标函数里没有给弃风弃光加惩罚优化器一定会尽量削减新能源出力因为这样可以减少燃气轮机出力和购电甚至让储能少充电。这显然违背新能源消纳优先的初衷。所以我习惯在目标函数里加一个很小的弃风弃光惩罚项比如每弃1MWh惩罚0.05元这里的价格系数是相对值具体看你的成本量纲让优化器在可行的情况下尽量消纳新能源。这个惩罚系数很小不会主导经济性但作用类似于“软约束”防止出现无意义的弃风弃光。2.2 储能SOC递推与充放电互斥的三种写法储能是日前调度里最容易写错的部分。核心约束包括SOC递推、充放电功率上下限、充放电状态互斥、初始SOC和末值SOC。SOC递推式是SOC(t1) SOC(t) (eta_ch * P_ch(t) - P_dis(t) / eta_dis) / E_rated注意充电时要乘充电效率放电时要把放电功率除以放电效率也就是说储能从电网取1MWh真正存进电池的只有0.95MWh左右放电时想送出1MWh电池内部要消耗约1.05MWh。效率放对位置能量守恒才严格成立。在Matlab代码里我通常用一个长度为25的SOC向量SOC(1)是初始值SOC(241)是调度结束时刻的SOCSOC sdpvar(1, T 1); C [C, SOC(1) SOC0]; C [C, SOC(T 1) SOC0]; % 日循环约束 for t 1:T C [C, SOC(t 1) SOC(t) (eta_ch * P_ch(t) - P_dis(t) / eta_dis) / E_rated]; end C [C, SOC_min SOC SOC_max];初始SOC等于末值SOC这条约束我建议一定要加。日前调度是一个滚动运行的周期过程今天调度结束时电池里的电量会直接影响第二天的初始状态。如果不加这条优化器会把电池在低谷时段充满、在高峰时段彻底放空看起来成本很低但第二天一开局就无电可用结果根本没法实际落地。充放电互斥可以通过二进制变量实现z_ch binvar(1, T); z_dis binvar(1, T); C [C, 0 P_ch P_ch_max * z_ch]; C [C, 0 P_dis P_dis_max * z_dis]; C [C, z_ch z_dis 1];这个“z_ch z_dis 1”保证了同一个时段不会同时充电和放电。有人会问能不能不加二进制变量只在功率平衡里让储能的净功率自由变化如果你的模型里没有储能运行损耗项而且充放电效率又被简化成1那确实可以不区分充放电用一个净功率变量就行。但只要涉及效率不对称、电池寿命损耗成本或者充放电电价套利的经济性分析就必须区分状态二进制互斥几乎是标准做法。2.3 与主网联络线的功率边界联络线购电功率我设定了一个上限常见做法是P_buy sdpvar(1, T); C [C, 0 P_buy P_buy_max];如果考虑余电上网卖电就要引入两个非负变量P_buy和P_sell并加二进制互斥防止同一时段既买又卖。我这版模型暂不考虑卖电因为新能源装机相对负荷来说不大多余电量通过弃风弃光处理更贴合模型设定的运行场景。3. 需求响应进优化模型价格型与激励型我分别怎么处理3.1 价格型需求响应的弹性矩阵线性近似价格型需求响应PDR的核心思想是用户根据电价高低调整用电行为。最常用的数理表达是价格弹性矩阵ΔL_i / L_i E_ij * (ρ_j - ρ_ref_j) / ρ_ref_j这个式子看着简单放进优化模型却很麻烦。因为电价如果是决策变量ΔL和ρ相乘就是非线性项如果电价是固定分时电价那么(ρ_j - ρ_ref_j)/ρ_ref_j就是个常数ΔL可以做线性化处理。在初版模型里我建议用简化的“高峰削减低谷回补”模型。具体做法是引入一个可削减量变量P_pdr_cut(t)让它在高峰时段有削减上限同时在全天尺度上加一个“削减的总电量不超过某一比例”的约束模拟用户把部分用电从高峰挪到其他时段的效果P_pdr_cut sdpvar(1, T); alpha 0.15; % 高峰时段可削减比例 for t 1:T C [C, 0 P_pdr_cut(t) alpha * P_load_raw(t)]; end这个模型的好处是直观、好调参而且不会出现非线性项。弹性矩阵方法更适合做敏感性分析或者把电价内生化的进阶模型但初版代码建议先用这种线性近似跑通全流程。3.2 激励型可中断负荷的二进制建模激励型需求响应IL指的是调度中心提前和用户签合同约定在特定时段可以中断一部分负荷并支付补偿费用。在优化模型里这天然适合用二进制变量表达“是否调用”这个动作P_il sdpvar(1, T); % 实际中断负荷量 z_il binvar(1, T); % 是否调用可中断负荷 IL_max 2.0; % 单时段允许中断上限单位MW c_il 0.95; % 中断补偿价格单位元/kWh C [C, 0 P_il IL_max * z_il];目标函数里加上补偿成本cost_il c_il * sum(P_il);这里有一个经济性上的细节可中断负荷的补偿价格如果设定得比高峰购电价还高优化器就不会调用它因为直接购电更便宜。所以c_il的取值一般要略低于高峰时段的购电价让模型在“购电还是切负荷”之间做一个真实的经济权衡。3.3 两类DR一起用的耦合边界我现在这个版本同时考虑了价格型和激励型需求响应把它们都体现在负荷侧功率平衡里P_g(t) P_w(t) P_pv(t) P_dis(t) P_buy(t) P_load_raw(t) - P_pdr_cut(t) - P_il(t) P_ch(t)如果把PDR的弹性矩阵做得更复杂PDR自身还要考虑削减和回补的守恒关系此时可以加一个总电量守恒约束sum(P_pdr_cut) beta * sum(P_load_raw)其中beta是需求响应比例上限比如全天的负荷转移比例不超过总负荷的5%。这个约束虽然让模型保守了一点但更接近真实用户行为也更不容易被优化器拿去钻空子。4. 目标函数、求解器选型和代码主体结构4.1 目标函数各项成本怎么拼起来日前经济调度的目标函数其实是“运行成本最小化”我这一版的目标函数包含四项微型燃气轮机运行成本a * P_g^2 b * P_g c其中用二次项近似燃料成本随出力增加而边际上升的特性向上级电网购电成本分时电价乘以各时段购电功率可中断负荷补偿成本补偿单价乘以中断量弃风弃光惩罚成本一个很小的惩罚系数乘以被削减的新能源电量写成Matlab片段cost_g a * sum(P_g.^2) b * sum(P_g) c * T; cost_buy sum(price_buy .* P_buy); cost_il c_il * sum(P_il); cost_curl curtail_penalty * (sum(P_w_forecast - P_w) sum(P_pv_forecast - P_pv)); objective cost_g cost_buy cost_il cost_curl;需要注意的是PDR削减负荷本身没有直接出现在目标函数里——它的价值体现在“负荷减少→购电和机组出力下降”的间接经济性上。如果某些文献里给PDR也加了补偿价格那就在目标函数里多一项即可不影响整体结构。4.2 为什么用YALMIPGurobi而不用手写linprog这个模型里有二进制变量、二次目标函数和耦合约束属于混合整数二次规划MIQP或混合整数线性规划MILP。Matlab自带的linprog只能解线性规划fmincon解非线性规划时碰到整数变量又很痛苦所以最可靠的方案是用YALMIP建模后端接Gurobi或CPLEX这类商业求解器。我做过一个简单对比在24时段这个规模下纯linprog思路得靠分段线性化手动处理代码量翻倍而且容易引入误差YALMIP建模的优势是约束可以用几乎和数学表达式一模一样的方式写出来排错非常直观。方案支持整数变量求解MIQP建模直观度适合场景matlab linprog否否低纯LP教学matlab intlinprog是否低小规模MILPYALMIP Gurobi是是高本文这种调度模型手写拉格朗日法很难否极低学术研究不推荐如果你没有Gurobi license可以申请学术版学生和科研用途是免费的。实在不方便也可以用YALMIP配开源的SCIP求解器处理这个规模也够用只是速度慢一些。安装好之后求解代码是ops sdpsettings(solver, gurobi, verbose, 1); diagn optimize(constraints, objective, ops);用YALMIP还有一个好处目标函数里的二次项P_g.^2可以直接写Gurobi会用MIQP求解器处理如果后续想换CPLEX只需要改sdpsettings里的solver名字模型代码一行都不用动。4.3 代码骨架与关键约束生成片段我习惯把整个代码分成四段参数定义、变量定义、约束组装、求解与结果提取。下面给出一个浓缩版骨架参数可以根据自己的算例调整%% 参数 T 24; eta_ch 0.95; eta_dis 0.95; E_rated 4; % MWh SOC0 0.5; SOC_min 0.2; SOC_max 0.9; P_ch_max 2; P_dis_max 2; % MW %% 决策变量 P_g sdpvar(1, T); P_buy sdpvar(1, T); P_w sdpvar(1, T); P_pv sdpvar(1, T); P_ch sdpvar(1, T); P_dis sdpvar(1, T); SOC sdpvar(1, T 1); P_il sdpvar(1, T); z_il binvar(1, T); z_ch binvar(1, T); z_dis binvar(1, T); %% 约束 C []; C [C, P_g_min P_g P_g_max]; C [C, 0 P_w P_w_forecast]; C [C, 0 P_pv P_pv_forecast]; C [C, 0 P_buy P_buy_max]; C [C, SOC(1) SOC0, SOC(T 1) SOC0]; C [C, SOC_min SOC SOC_max]; C [C, 0 P_ch P_ch_max * z_ch]; C [C, 0 P_dis P_dis_max * z_dis]; C [C, z_ch z_dis 1]; C [C, 0 P_il IL_max * z_il]; for t 1:T C [C, SOC(t 1) SOC(t) (eta_ch * P_ch(t) - P_dis(t) / eta_dis) / E_rated]; C [C, P_g(t) P_w(t) P_pv(t) P_dis(t) P_buy(t) ... P_load_raw(t) - P_pdr_cut(t) - P_il(t) P_ch(t)]; end在完整版代码里我还会加微型燃气轮机的爬坡约束形如for t 2:T C [C, -ramp_down P_g(t) - P_g(t - 1) ramp_up]; end爬坡约束在日前调度里是很容易漏掉的一条我见过很多模型跑出来燃气轮机出力在相邻时段剧烈跳动其实工程上根本做不到。加上爬坡约束后机组的运行曲线会平滑很多。5. 一组典型算例的搭建与调度结果解读5.1 参数配置全景为了验证这个模型我搭建了一组典型参数组合。负荷曲线取一个夏季工作日的典型形状凌晨负荷低约8MW中午一段小高峰晚上18点到21点达到全天峰值约18MW。风电预测曲线是夜间大、白天小光伏预测曲线是中午大、早晚几乎为零。具体参数如下参数取值说明燃气轮机容量15 MW下限3MW联络线购电上限8 MW与主网交换功率约束储能容量4 MWh最大充放电功率2MW风电预测峰值6 MW夜间接近满发光伏预测峰值4 MW中午接近满发分时电价谷0.4、平0.8、峰1.3元/kWh简化三时段可中断负荷上限2 MW单时段最多调用2MW可中断补偿价0.95元/kWh低于峰值购电价PDR削减比例15%高峰负荷可削减比例弃风弃光惩罚0.05元/kWh软约束尽量消纳这些数值不一定对应某个真实工程但量级是合理的跑出来的调度策略非常有解释性。5.2 24小时调度曲线怎么解读这段是我觉得整个课题最有意思的部分。优化结果出来后把24个时段的调度曲线画在同一张图里你会看到几个非常清晰的规律凌晨0点到6点电价是低谷负荷也低储能会在这段时间充电SOC从0.5一路充到接近上限。风电此时出力较大燃气轮机压低到最小技术出力附近优先用风电满足负荷多余电量和电网谷电存进电池。白天8点到11点光伏出力爬升电价进入平段储能开始停止充电甚至小幅放电燃气轮机出力小幅抬升购电功率维持在一个合理水平。傍晚18点到21点这是全天的重头戏。负荷冲到最高电价也进入峰值段储能在18点前已经完成充电此时开始放电顶峰同时可中断负荷在19点和20点被调用PDR削减了一部分高峰负荷。燃气轮机满发如果还不够继续从主网购电但购电功率会被联络线上限卡住。这个调度策略的本质是储能做“日内能量搬移”DR做“峰值负荷软削”燃气轮机做“基荷和腰荷平衡”联络线做“硬边界”。四种灵活性资源各司其职最终让整个微电网以最低成本度过24小时。5.3 成本与峰谷对比在我的这组参数下加了需求响应和储能之后微电网的日运行成本大约降低12%左右负荷峰值也从18MW降到约15.5MW。需要说明的是这个降幅高度依赖参数设置储能容量越大、峰谷电价差越大、可中断负荷补偿价越低降幅就越明显。如果把储能和DR的收益拆开看储能主要通过低谷充电高峰放电赚取峰谷价差DR主要通过削减高峰购电来降低电费。两者叠加在一起还能产生协同效应储能放电降低了高峰期的购电需求DR削减又进一步减少了联络线压力联合起来的效果往往大于单一措施的叠加这就是多类型柔性资源协同调度的价值。6. 调试中容易被坑到的五个细节6.1 无解的常见原因与排查顺序我调试这个模型时遇到的第一类坑就是约束过紧导致无解。最典型的情况是SOC终值约束设定为必须回到初始值但储能容量太小或者峰谷时段分布不合理导致无论怎么充放都无法完成循环求解器直接报infeasible。排查无解问题时我的习惯是从简到繁逐步加约束先不加DR、不加储能、不加爬坡约束跑一个纯出力的线性平衡模型确认基础可行然后依次加储能、加DR、加爬坡每加一类约束就重跑一次这样一旦出现无解能很快定位到是哪类约束出了问题。YALMIP有个很实用的函数checkdiagn optimize(C, objective, ops); if diagn.problem ~ 0 check(C); endcheck会把每条约束的残差和违反量打印出来违反最严重的那条往往就是问题所在。6.2 二进制变量和连续变量相乘的坑这个模型里有个隐藏的非线性陷阱如果我想模拟“储能以某个效率参与调频”或者“可中断负荷的补偿价格和中断时长有关系”很容易写出连续变量乘以二进制变量的表达式比如P_il * z_il。这在YALMIP里会变成双线性项Gurobi没法直接处理求解会非常慢甚至不收敛。处理办法是不要把这种乘法写进约束而是像前面代码那样把二进制变量放进不等式的上限里P_il IL_max * z_il。这样既表达了“调用时才能有中断量”又保持了线性Gurobi秒解。6.3 二次目标导致求解时间陡增的处理我一开始的目标函数直接保留了燃气轮机的二次成本项优化器跑MIQP。在24时段尺度下还好但如果把时间步长细化到96个时段或者微电网数增加到十几个MIQP的求解时间会明显上升。这时候可以考虑把二次成本曲线分段线性化做成MILP把每一段写成一组线性约束加二进制变量求解速度会快很多。YALMIP里可以用sdpvar和binary的组合自己实现分段线性化也可以用内置的blackbox函数处理但自己写更可控。我的建议是先跑通MIQP理解调度逻辑再考虑优化求解速度的问题不要一开始就上分段线性化否则排错时你根本分不清是建模错误还是求解数值问题。6.4 结果提取时最容易忘的value()YALMIP的变量是sdpvar对象它不直接保存数值。很多新手跑完optimize直接plot(P_load)结果发现是空图或报错就是因为忘了加value()P_g_opt value(P_g); P_ch_opt value(P_ch); P_dis_opt value(P_dis); P_il_opt value(P_il); soc_opt value(SOC);取了value之后才能正常绘图和统计。这个看起来特别基础但我至少在这里被绊过两三次尤其是用bar(P_g_opt)和stairs绘24小时曲线的时候一定要保证传入的是数值数组而不是符号对象。6.5 电价序列和负荷序列的维度对齐分时电价我通常这样生成price_buy [repmat(0.4, 1, 8), repmat(0.8, 1, 4), repmat(1.3, 1, 6), repmat(0.8, 1, 6)];这个序列一共24个数顺序是0点到23点。看起来简单但特别容易在拼接时段边界时搞错。我的建议是把电价、负荷、风电预测、光伏预测都放在同一个脚本顶部集中定义每个向量都写注释标明时段起止跑结果前先用plot看一眼24小时曲线是否连续合理。如果光看数字边界错一位LOAD和PRICE对不上调度结果就会很怪比如低谷时段却在买峰价电。最后的一点个人体会把风光储和需求响应全部耦合进日前调度最深的感受是模型本身并不复杂复杂的是对各种资源时间耦合关系的理解。储能把能量从低谷搬到高峰需求响应把负荷从高峰挪到低谷燃气轮机在整个24小时轴线上做兜底本质上都是在成本信号的引导下实现能量在时间和空间上的重新分配。Matlab在这个场景下最大的优势是建模灵活、图表方便配合YALMIP和商业求解器几乎可以把数学优化模型一比一翻译成代码。如果你正准备上手类似的课题建议先把储能SOC循环、DR调用逻辑这两块理解透再动手改代码进度会快很多。