1. 项目背景与核心问题拆解1.1 热电联产机组为什么会“堵住”风电消纳做风电、火电联合调度的朋友应该都遇到过这个场景冬季夜间风电场满发但全网用不完这些电于是调度只能下令弃风。明明是可再生能源白白扔掉看着就心疼。问题出在哪很大程度出在那批“以热定电”的热电联产机组CHP身上。北方城市供热季是靠热电厂带基础热负荷的。供热管网需要稳定的热源而CHP机组的运行模式是“以热定电”——只要供热量确定了发电出力就被机组的电热特性曲线绑死了。冬季夜间恰好是供热的峰值时段CHP机组必须高负荷运转发出来的电却正好撞上电网负荷低谷和风电大发期。电多、用不掉风电就得让路。这个矛盾我前几年在北方某省份的实际调度数据里看得非常直观三月供热季结束前后弃风率能差出五六个百分点。传统思路是上电锅炉、蓄热罐这类设备来解耦热电耦合但这涉及一台台机组怎么搭配运行的问题。单纯加设备不解决机组之间的协同调度投资打了水漂风电照样弃。这篇项目做的就是一个“联合优化控制”框架把多台CHP机组、电加热装置、电能负荷、热负荷和风电出力放在同一个优化模型里用Matlab求解出未来若干个时段内每台机组的出力计划目标就是在保证供热和供电的前提下最大程度把风电消纳掉。1.2 这个优化控制的本质是在做什么剥开“联合优化控制”这层外壳本质就是一个带约束的多时段经济调度问题。我更愿意把它理解为“给一整片电热联合系统做一个未来24小时的动作编排计划”风电出力多出来的时段怎么压CHP电出力、怎么用储热来顶上供热缺口风电出力小的时段怎么让CHP多发电、给储热罐蓄热给后续时段留出调节空间所有动作都要满足电网功率平衡、热网平衡、机组爬坡速率、出力上下限等硬约束。项目实施路径分成四步先梳理CHP机组的电热运行可行域再搭目标函数和约束的数学表达式然后通过Matlab调用规划求解器完成优化计算最后做多场景仿真验证模型有效性。下面我把每一个环节的实操细节都展开说清楚包括我踩过的坑。2. 数学模型搭建方法2.1 目标函数怎么选才符合工程直觉目标函数是优化模型的“指挥棒”定得对不对直接决定结果是否符合调度预期。这个项目里有一个容易被忽略的细节最大化风电消纳和最小化系统运行成本表面上方向一致实际建模时处理方式不同。“最大化风电消纳”最直接的写法是让目标函数里包含一个风电出力最大化项比如最大化每个时段风电场并网出力之和。但在实际工程调度里单纯追求风电消纳会导致求解器利用一切手段压CHP出力哪怕某些CHP机组的运行效率很高、成本很有优势也被一刀切地压到最低这不符合经济性要求。我见过不少初学的朋友在这里翻车目标函数只写风电消纳最大化求解出来的方案是“不管成本多高都要弃煤保风”。等对上实际考核指标比如煤耗、供热保障率时方案根本没法落地。这个项目采用的替代做法是把目标函数改为最小化系统运行总成本同时在成本的各项组成中引入“弃风惩罚成本”弃风成本系数设置得比常规机组发电成本高一个量级这样求解器在权衡“压CHP出力节省的煤耗成本”和“弃风惩罚成本”时会主动选择优先消纳风电但当风电已经无法全额消纳时系统才会选择弃风避免无限制牺牲经济性。这种间接建模方式叫“惩罚系数法”是工程里最实用、最不容易出问题的目标表达方式。具体写作弃风惩罚成本 K_penalty * sum( Wind_available(t) - Wind_used(t) )其中 K_penalty 取值要根据系统具体情况调整通常设定在500~1000元/MWh量级远高于火电边际煤耗成本。热负荷侧如果也追求经济性目标函数里再加一个恒温控制惩罚项防止为了省电费导致供热温度波动过大。2.2 弃风弃在哪里风电并网约束的表达风电消纳问题在模型里的核心约束是风电场并网功率约束0 P_wind_used(t) P_wind_available(t) P_wind_used(t) P_load(t) P_electric_boiler(t) - P_chp_sum(t)第二条不等式实际上是电网功率平衡的变体它表达了一个直观逻辑风电能用的最大量取决于“用电负荷电锅炉额外吃掉的电CHP必须发的电”。风电用不完就是因为右边的可用空间被CHP的“以热定电”出力给挤占了。电锅炉在这里扮演的角色很关键。它本质上是把风电转化成热量供给供热管网相当于给风电“找个出口”。在约束表达里电锅炉出力算作电负荷的一部分同时计入热功率平衡方程电功率平衡风电出力 CHP电出力 电负荷 电锅炉耗电热功率平衡CHP热出力 电锅炉产热 储热罐放热 热负荷 储热罐蓄热这个“电转热”power-to-heatP2H的解耦思路是北方供热地区提高风电消纳能力最成熟的手段之一。2.3 CHP机组电热可行域的建模细节CHP机组不是随便一个“电出力和热出力”组合都能运行。它有一个电热运行可行域主要由以下几段构成背压工况线热电比固定电出力随热出力线性变化通常是可行域的上边界纯凝工况线热出力为零时的最小和最大电出力区间;抽汽调节区间电出力可以在一定范围内独立于热出力调节。把这四条边界围出来的多边形区域用一组线性不等式约束表示就是CHP机组的可行域。项目里我用的是最常见的四边形近似法P_chp_min P_chp(t) P_chp_max 0 H_chp(t) H_chp_max P_chp(t) - c1 * H_chp(t) P_min_pure_condensing P_chp(t) - c2 * H_chp(t) P_max_back_pressure其中 c1、c2 是机组的电热特性系数由厂家热力试验数据拟合得到。这个四边形近似在工程精度足够但要注意如果实际机组是抽凝式且运行区间包含多个工作模式四边形会高估调节能力需要做更精细的可行域聚类。这一点我后面在常见问题部分还会再提。3. Matlab代码实现要点3.1 工具箱选型我为什么选YALMIPGurobiMatlab平台做这类优化调度可选方案有好几套手写线性规划直接调linprog、用fmincon做非线性规划、用YALMIP建模后调Gurobi/CPLEX。这个项目由于约束多、变量维度大每台机组、每个时段都有变量linprog手动拼矩阵容易拼错且难维护fmincon处理大规模线性约束又不够稳最终我选的是YALMIPGurobi组合。原因是这个模型本质上是混合整数线性规划MILP——机组启停状态的0-1变量和功率连续变量混在一起Gurobi在求解速度上有碾压级优势。举个例子24时段、3台CHP机组、1台电锅炉、1个储热罐的模型变量数大约800个、约束约1500条Gurobi默认参数下一般十几秒就能收敛到0.1%的MIP gap。用linprog手搓基本要在建模阶段耗一半时间fmincon则可能因为非凸的可行域陷入局部最优。当然没有Gurobi license的朋友可以改用免费的CBC求解器YALMIP底层可以直接调用。速度和稳定性略差但中小规模场景完全够用。我建议初学先用CBC跑通流程再换Gurobi提升性能。3.2 数据结构的组织从杂乱数据到规范矩阵Matlab代码实现这个模型最大难点不是求解公式而是数据的组织。我强烈建议一开始就把所有输入数据封装成结构体struct不要散落着写在脚本里。项目里我的数据组织方式如下% 输入数据结构 sys.H 24; % 调度时段数 sys.nt 3; % CHP机组数 % 风电预测数据 wind.P_avail [34 32 30 ...]; % 各时段可用风电出力MW % 负荷数据 load.P_demand [120 115 ...]; % 各时段电负荷MW load.H_demand [210 205 ...]; % 各时段热负荷MW % CHP机组参数3台机组分别定义 chp(1).P_max 120; chp(1).P_min 40; chp(1).H_max 120; chp(1).c1 0.25; chp(1).c2 0.75; chp(1).ramp_rate 20; % MW/h爬坡速率 % 电锅炉参数 eb.P_max 60; % 最大耗电功率MW eb.eta 0.98; % 电转热效率 % 储热罐参数 ts.V_max 300; % 最大储热量MWh ts.charge_rate_max 30; % 最大充/放热功率MW把这些数据集中定义之后模型构建阶段反复引用同一个结构体字段不容易发生“变量名打错但语法合法”这类隐错。我见过太多同行把数据散落在一堆变量里调试的时候根本理不清哪些属于同一台机组。3.3 核心代码框架逐段拆解模型构建的核心代码分五大块变量定义、目标函数、等式约束、不等式约束、求解与结果回读。第一步定义决策变量% 定义一个27x24的优化问题3CHP电锅炉储热风电功率平衡相关 P_chp sdpvar(3, 24); % CHP电出力 H_chp sdpvar(3, 24); % CHP热出力 P_wind sdpvar(1, 24); % 风电并网出力 P_eb sdpvar(1, 24); % 电锅炉耗电功率 E_ts sdpvar(1, 25); % 储热罐储热量25个节点覆盖首尾 H_ts_charge sdpvar(1, 24);% 储热罐蓄热功率 H_ts_release sdpvar(1,24);% 储热罐放热功率 U_chp binvar(3, 24); % 机组启停状态特别注意这里储热罐的储热量 E_ts 定义成了25个点而不是24个这是为了表达储能状态在相邻时段之间的递推关系E_ts(t1) E_ts(t) 充热 - 放热。很多新手直接定义24个点写递推约束时发现 t24 时刻的状态没法闭环后期还得回头改维度不如一开始就预留一个节点。第二步目标函数% 成本项 objective 0; for t 1:24 % 煤耗成本线性近似也可以取二次函数的线性化分段 objective objective sum(coal_cost_coef .* P_chp(:,t)); % 启停成本 objective objective sum(start_cost .* max(0, U_chp(:,t)-U_chp(:,t-1))); % 弃风惩罚成本核心 objective objective K_penalty * (wind.P_avail(t) - P_wind(t)); end第三步约束条件组装constraints []; % 电功率平衡 for t 1:24 constraints [constraints, sum(P_chp(:,t)) P_wind(t) ... load.P_demand(t) P_eb(t)]; end % 热功率平衡 for t 1:24 constraints [constraints, sum(H_chp(:,t)) eb.eta * P_eb(t) ... H_ts_release(t) - H_ts_charge(t) load.H_demand(t)]; end % CHP可行域约束 for k 1:3 for t 1:24 constraints [constraints, chp(k).P_min * U_chp(k,t) P_chp(k,t) chp(k).P_max * U_chp(k,t)]; constraints [constraints, 0 H_chp(k,t) chp(k).H_max * U_chp(k,t)]; constraints [constraints, P_chp(k,t) - chp(k).c1 * H_chp(k,t) chp(k).P_min]; constraints [constraints, P_chp(k,t) - chp(k).c2 * H_chp(k,t) chp(k).P_max]; end end第四步求解与回读ops sdpsettings(solver, gurobi, verbose, 2, mipgap, 0.001); optimize(constraints, objective, ops); % 回读结果 P_chp_opt value(P_chp); H_chp_opt value(H_chp); P_wind_opt value(P_wind); P_eb_opt value(P_eb); E_ts_opt value(E_ts);3.4 爬坡约束和储热罐约束两类容易漏掉的物理约束除了上述基础约束还必须给机组加爬坡约束。CHP机组有功出力不能瞬间跳变通常用for k 1:3 for t 2:24 constraints [constraints, ... -chp(k).ramp_rate P_chp(k,t) - P_chp(k,t-1) chp(k).ramp_rate]; end end这里我吃过亏不加爬坡约束的模型优化结果很“漂亮”风电消纳率高得惊人但拿到现场一执行机组响应不了那么快计划全部作废。加上爬坡约束后消纳率会下降两三个百分点这才是真实可信的结果。储热罐方面有两组约束一是储热量上下限和初始/末尾能量约束防止调度周期结束把热全部放空二是充放热功率互补约束同一时段不能同时充和放。% 储热罐能量递推 for t 1:24 constraints [constraints, E_ts(t1) E_ts(t) ... (H_ts_charge(t) - H_ts_release(t))]; end % 储热罐容量约束 constraints [constraints, E_ts(1) ts.E_initial]; constraints [constraints, E_ts(25) ts.E_final]; constraints [constraints, 0 E_ts(2:25) ts.V_max]; % 充放热不能同时进行用二进制变量互斥 C_ts binvar(1, 24); constraints [constraints, H_ts_charge(t) ts.charge_rate_max * C_ts(t)]; constraints [constraints, H_ts_release(t) ts.charge_rate_max * (1-C_ts(t))];末能量约束这一条特别重要。如果把 E_ts(25) 放开不管求解器会把所有热都留在罐子里等于白给的成本或者把储热量压到零相当于调度吃光了所有储能余量两种都不是运营商想要的。通常做法是让始末储热量相等保证日循环可持续性。4. 案例仿真与结果分析4.1 典型场景设置我用一个基础算例来验证模型。场景设定如下3台CHP机组容量分别为120MW、100MW、80MW电锅炉1台最大耗电60MW效率0.98储热罐有效容量300MWh最大充放热功率30MW初始储热量100MWh调度周期24小时风电预测数据选用北方某风电场冬季典型日夜间23:00-07:00风电大发白天出力回落热负荷全天高位运行尤其夜间更大电负荷白天高、夜间低。这个场景基本复现了“风电大发时段恰好供热需求高、电负荷低”的典型困境。4.2 没有联合优化时发生了什么为了说明效果先跑一版约束“CHP按以热定电方式运行、无电锅炉无储热罐”的基准场景。结果一目了然夜间23:00至次日07:00CHP为了带供热总电出力最低仍有160MW左右而电负荷只有120MW硬生生占掉了风电的上网空间。风电可用功率在夜间峰值达到100MW实际并网只有30MW左右弃风率高达60%以上。这数据看着触目惊心但在北方供热期是真实存在的普遍现象。4.3 引入联合优化后的效果加入电锅炉和储热罐、并运行联合优化模型后结果发生了明显变化风电消纳比例从约52%提升到了约90%以上弃风集中时段明显缩短从原来夜间连续7小时弃风缩减到凌晨两三个小时的轻度弃风CHP总电出力曲线不再是完全跟着热负荷走而是在风电大发的时段主动压低出力把供热任务部分转移给电锅炉和储热罐放热系统总煤耗成本比基准场景有所上升因为电锅炉耗电增加但加上弃风惩罚成本后的“综合成本”明显下降。我额外做了一个 CHP储热罐无电锅炉的对照算例发现储热罐单独使用只能把风电消纳率提高到70%左右。原因在于储热罐只能平移热负荷的时序不能真正扩大用电空间电锅炉则直接把风电转化为热负荷相当于“新增了用电需求”两者配合才能最大化消纳效果。这个结论和业内多篇文献的结论一致也验证了模型的合理性。提示实际工程中不要指望消纳率能到100%。受机组最小技术出力、网络约束等因素影响轻微弃风往往是经济上最优的。模型的价值是把弃风控制到合理区间而不是数字上做到最好看。4.4 不同惩罚系数下的敏感性分析消纳效果和成本之间有一个平衡点而这个平衡点主要由弃风惩罚系数 K_penalty 控制。我跑了一组敏感性分析K_penalty (元/MWh)风电消纳率系统总成本万元/日机组启停次数20072%156650086%162880092%16412120094%16618注意看惩罚系数从800提高到1200消纳率只提高2个百分点成本却多出2万元机组启停次数也大幅增加——频繁启停对设备寿命非常不利。也就是说盲目追求“零弃风”并不理性。实际工程里我一般把这个系数定在600~900元/MWh之间让模型在“尽量消纳”和“避免过度调节”之间取一个平衡。5. 常见问题与调试心得5.1 求解器报“数值鲁棒性”问题怎么办YALMIP/Gurobi求解时经常遇到提示“Numerical problems / Ill-conditioned data”。这通常不是模型写错了而是量纲差异太大。比如储热罐蓄热量单位是MWh数值在几百而弃风成本单位是元数值在几万甚至几十万目标函数里同时出现不同量级的数值数值求解器容易出状况。我的处理方式是统一量纲所有功率量用MW能量量用MWh费用项全部除以1e3变成“千元”或者干脆把目标函数整体缩小1000倍保持求解数值在1e-3到1e5区间内。这类细节在大规模模型中影响更大数值扰动会导致收敛速度骤降甚至报错。还有个细节储热罐的储热量维度MWh和充放热功率维度MW在递推约束里相乘要匹配E_ts(t1) E_ts(t) [充放热功率 × 1小时]这里每个时段长度是1小时所以功率数值可以直接累加。如果改成15分钟一个时段递推式里必须乘以0.25否则能量不守恒的bug非常隐蔽查半天查不出来。5.2 风电消纳率算出来明显偏高检查约束有没有漏我发现一个很常见的排查路径模型刚建完跑出来风电几乎全额消纳高到不真实。这时候优先排查三类约束是否遗漏了CHP机组出力下限和可行域下边界约束如果 P_chp 的边界写成了 0 P_chp P_max那么求解器可以把 CHP 出力压到零来腾空间给风电这在物理上完全不合理是否遗漏了爬坡约束夜间风电大发时求解器可能让CHP从晚高峰的高出力直接跳到夜间的低出力一个时段内跨了太大幅度是否给电锅炉设置了过大的容量如果电锅炉容量超过热负荷需求模型会无脑把电锅炉拉到满发来消耗风电热平衡被“强行”满足这也不是真实可行方案。这三条我全踩过。排查手段也不复杂把优化得到的各变量曲线打印出来逐一核对每个时段功率平衡是否严格满足、设备出力是否落在可行域内半小时就能定位问题。5.3 整数变量导致的求解时间剧增模型里有机组启停状态和储热罐互斥状态两类二进制变量如果再加电锅炉的开关状态总整数变量会达到“机组数×时段数×2~3”个。以10台机组、24时段算就有至少480个整数变量直接跑默认参数可能要好几分钟。优化手段有三招第一把非必要整数变量替换成连续变量边界约束。比如电锅炉如果运行模式是“深度调峰期间才开”可以预判时段提前固定开关状态减少整数变量第二给Gurobi设置合理的MIP gap比如 0.001 或 0.005。对调度计划而言1%以内的次优解完全可接受没必要追求数学上的绝对最优第三用“滚动优化”替代全局24小时优化例如每4小时滚动一次每次只优化未来4小时并冻结后20小时的预测信息模型规模大幅下降求解速度从分钟级降到秒级。5.4 CHP可行域多边形过简化的后果最后专门提一个模型精度问题。前面用四边形近似CHP可行域其实隐性地“缩减”了机组实际调节空间。在实际项目中我遇到过抽汽式CHP机组的实际可行域接近五边形甚至六边形考虑低压缸最小冷却流量约束四边形近似可能导致计算结果偏向保守系统调节能力被低估风电消纳率低于实际可达到水平极端情况下优化的运行点在四边形内但不在真实可行域内下发给现场后机组无法跟踪执行。解决方法是把机组厂家的热力特性试验数据找出来用凸包法convhull拟合真实可行域的多边形顶点然后转化成多面体约束% 用凸包顶点构建可行域约束YALMIP支持polygon约束 Poly [P1 H1; P2 H2; P3 H3; ...]; % 从试验数据得到的顶点 constraints [constraints, ismember([P_chp(k,t); H_chp(k,t)], Poly)];如果没有试验数据至少也要把四边形的关键顶点坐标标定准确不要随便拍脑袋定系数。C1、C2这类系数偏差10%最终调度计划就能差出一个档次。6. 一些值得坚持的工程习惯项目做到这个程度模板代码其实已经完成核心验证。如果后续要往真实现场走我建议再补齐三个动作把优化结果做24小时功率曲线可视化分析每一张图都要能解释“为什么在这个时刻CHP压低出力、储热罐放热”做一个“风电预测误差”影响分析在风电预测值±20%扰动下重复求解观察消纳率的波动范围预留与SCADA/AGC系统的数据接口把优化结果输出成标准格式的调度指令文件。我个人在实际操作中有一个习惯值得分享模型跑通之后第一件事不是调参数刷高消纳率而是先把“基准场景”的每一条约束都对一遍物理意义再放开对比。Matlab里用YALMIP建模的最大便利是约束可读性高逐条打印出来就能当报告用。这套代码我后来又扩充过多区域版本、加了日前实时两阶段调度逻辑核心框架基本没动过。前期把模型和数据结构写规范后面改起来会顺手很多。