做微电网运行优化，绕不开一个朴素的问题：电负荷和热负荷明明共存于同一个园区，为什么大多数优化模型只盯着“电”在算？我最近在某园区微电网项目的经济调度分析里，把基于电热联合调度的微电网运行优化模型在 MATLAB 里完整跑通了一遍。今天把这个模型的变量设计、约束写法、求解代码和调试过程一起整理出来，适合正在做微电网调度、综合能源系统优化的研究生和工程师参考，尤其是手里有燃气轮机、热电联产机组、电锅炉、蓄电和蓄热设备，想算清楚“一套设备怎样组合运行一天最省钱”的这类场景。

先说明白这个模型能回答什么问题：光伏和风电出力有波动，电价有峰谷，天然气也有价格，园区既有电负荷又有热负荷，电锅炉可以把电转成热，热电联产机组（CHP）可以同时产电产热，蓄电罐和蓄热罐可以把能量跨越时间搬运。电热联合调度就是在满足电、热两类负荷的前提下，决定每一台设备每个时段开多大功率、充多少放多少，让全天总运行成本最低。这本质是一个混合整数线性规划问题，MATLAB 加 YALMIP 工具箱处理起来非常顺手。

1. 电热联合调度的优化空间：为什么单算“电”会吃亏

1.1 微电网里的“热”角色：从热电联产到电锅炉

大多数微电网项目在设计初期，最容易犯的错是把热负荷单独扔给燃气锅炉去扛，电负荷单独靠光伏、储能和电网。这样做逻辑清晰，但经济性很差。因为园区里大量热负荷存在的地方，往往同时存在电负荷，而 CHP 机组本身就是“一份燃气投入，同时产出电和热”的设备。如果你只让它发电、不让它供热，相当于花了同样的燃料钱，白白扔掉了一半以上的可用能量。

反过来，电锅炉也是一种耦合设备，它消耗电能，产出热能。电锅炉本身不一定省钱，但在低谷电价时段，或者在光伏大发导致弃光时，用电制热可能比用天然气制热便宜得多。于是“电”和“热”会通过 CHP 和电锅炉形成两条耦合通道。联合调度模型的价值，就是把这两条通道放进同一个优化问题里，让系统自己决定何时走“气到电到热”的路径，何时走“电网电到热”的路径。

CAT在一个典型的园区场景里加入蓄热罐后，热负荷甚至可以像电储能一样实现时移：白天光伏多的时候多产热存起来，晚上热负荷上来再放热。这个过程单看电力系统模型是完全看不见的，但它在实际运行中能显著降低顶峰时段购电压力。

1.2 一个直观案例：分产调度和联合调度的成本差在哪

我在调试模型时用了两个对照方案。第一个方案是分产调度：CHP 按“以热定电”的方式运行，也就是热负荷多少，CHP 就发对应比例的电，缺电直接从电网买，电锅炉和蓄热罐不参与联动。第二个方案是电热联合调度：所有设备统一优化。

同样一组典型日负荷数据、同样一套设备参数，分产调度一天运行成本大约是 14230 元，联合调度是 12350 元，节省了约 13.2%。这 13.2% 主要来自三个地方。

一是低谷电价时段，电锅炉启动储热，替代了高峰时段部分燃气消耗。二是 CHP 不再被“以热定电”绑死，而是更多在峰电时段满发，在低谷时段压低出力或者停机，热缺口由蓄热罐释放的热量补上。三是光伏大发时段，系统主动减少购电甚至向电网卖电，把多余电能转化为蓄热罐里的热能，避免弃光。这三个收益来源，离开电热联合调度模型都很难被量化。

2. 模型的数学骨架：目标函数、约束与变量设计

2.1 决策变量的分层设计：连续功率变量与启停0-1变量

优化模型的第一步是确定决策变量，我的做法是分成连续变量和二进制变量两类。连续变量包括各设备每个时段的功率和能量状态，二进制变量主要用于 CHP 的启停状态。

CHP 的启停变量不是可选项，而是必须项。因为燃气轮机存在最小技术出力，比如 50kW，如果只建模成连续变量，优化器为了省钱会把 CHP 压到 1kW、2kW 这种根本不存在的运行点。引入 on_chp 后，功率下限写成P_chp_min * on_chp <= P_chp <= P_chp_max * on_chp，机组要么不开，要么至少开到 50kW，这才符合实际。

一个需要提前想清楚的点：蓄电和蓄热设备的“充放互斥”要不要加二进制变量。从理论严谨性角度，应该避免同时充电和放电。但在实际模型中，由于充放都有损耗，且损耗已经进了目标函数，优化器在经济上没有动机让同一台储能同时充放。所以我通常不加这组二进制变量，后续求解速度会快很多。

2.2 电热平衡约束与CHP可行域约束的建模细节

电功率平衡是微电网模型的地基，任何一个时段都必须满足：

pv(t) + wt(t) + P_chp(t) + P_dis(t) + P_buy(t) = P_load(t) + P_eb(t) + P_ch(t) + P_sell(t)

左边是电源侧，右边是负荷侧。注意电锅炉的耗电功率 P_eb 必须放在负荷侧，这是初学者最容易漏掉的一项。光伏、风电数据如果来自预测曲线，可以直接作为已知量放进等式。如果考虑弃光，就需要给光伏增加一个“可削减”变量。

热功率平衡相对简单：

H_chp(t) + H_eb(t) + H_tes_dis(t) = H_load(t) + H_tes_ch(t)

蓄热罐充热时算是热负荷，放热时算是热源。

CHP 的电气特性是电热联合调度里最关键的约束。我用固定热电比模型，也就是：

H_chp(t) = rho_chp * P_chp(t)

这种模型用起来最省事，适合做可行性分析。它意味着 CHP 发多少电就附带产生多少热。当热负荷很小、电负荷很大时，多余的热可以进蓄热罐；当热负荷很大、电负荷很小时，则只能通过降低 CHP 出力来减少热产出，缺电部分由电网补。也可以把热电比关系改成可行域约束，允许热电比在一定范围内变化，但那样会引入更复杂的耦合不等式，对线性求解器不友好。

CHP 的天然气消耗我采用的是线性近似：

gas_consume(t) = a_chp * P_chp(t) + b_chp * on_chp(t)

其中 a_chp 是边际气耗系数，单位是 m³/kWh，b_chp 是空载气耗系数，单位是 m³/h。这两个系数由 CHP 性能曲线拟合得到。这种方法比直接用效率除法更贴近实际，也保留了 MILP 的线性结构。

2.3 储能SOC、爬坡、购售电约束怎么列才不踩坑

储能约束最容易踩坑的是 SOC 的索引长度。如果优化周期是 T=24，充放电功率变量长度是 24，那么 SOC 状态变量长度应该是 25，因为要从初始时刻一直递推到 T 时刻结束。我习惯用 E_bess(1) 表示初始电量，E_bess(T+1) 表示最终电量，中间递推关系写成：

E_bess(t+1) = E_bess(t) + dt * (eta_ch * P_ch(t) - P_dis(t) / eta_dis)

充电时乘以效率，放电时除以效率。这个细节特别重要，如果方向写反，优化器会发现“充进去 1 度电，放出来 1.2 度电”，相当于凭空创造能量，结果会出现一堆反直觉的充放行为。

容量约束要加在电量状态变量上，而不是加在充放功率上：

E_bess_min <= E_bess(t) <= E_bess_max

同时还要给充放电功率分别设置上限。初始和终值电量可以设定成同一个值，保证日内调度的可持续性。如果允许终值高于初始值，调度结果其实就是“一天之内把储能充满”，成本会被人为压低，没有实际参考意义。

CHP 的爬坡约束是另一个容易出问题的地方。爬坡定义是相邻时段功率变化量：

-ramp_chp * dt <= P_chp(t+1) - P_chp(t) <= ramp_chp * dt

注意在 MATLAB 中的索引是从 1 开始，所以代码里会有P_chp(2:T) - P_chp(1:T-1)的写法。

购售电约束方面，购电和售电功率各设上限，要不要加“不同时购售”的互斥二进制变量，看实际需求。我一般先不加，因为分时电价下基本不会出现同时购售的荒谬解。如果最终结果异常，再补一组互斥变量也不迟。

3. MATLAB代码实现：数据组织、建模与求解器对接

3.1 参数与负荷曲线如何整理成矩阵

MATLAB 做优化建模，第一步就是统一数据格式。我强烈建议所有时序数据都用1 x T的行向量组织，电价、负荷、光伏、风电都是如此。这样后续与 YALMIP 的变量做点乘运算时可以直接.*，不需要来回转置。

单位也要统一。功率单位用 kW，能量单位用 kWh，电价单位用 元/kWh，天然气价格单位用 元/m³。如果功率数据是 15 分钟一个点，那么 T=96，dt=0.25，所有 SOC 递推式都要乘以 dt。很多模型跑出来的结果时间对不上，最后发现是 dt 漏乘了。

一段典型的数据初始化代码：

T = 24; % 优化时段数 dt = 1; % 时间步长，单位小时 % 电负荷和热负荷，单位 kW P_load = [105 98 92 88 90 96 110 135 150 148 132 128 ... 115 108 120 138 155 160 148 136 120 112 108 100]; H_load = [80 75 72 68 66 70 76 82 86 88 80 75 ... 70 68 72 80 90 95 92 88 82 78 74 70]; % 光伏和风电预测出力，单位 kW pv = [0 0 0 0 0 5 25 60 90 120 140 145 130 110 90 70 40 10 0 0 0 0 0 0]; wt = [18 20 22 24 22 20 18 15 12 10 9 8 10 14 16 18 20 22 24 22 18 16 14 12]; % 分时电价，单位 元/kWh price_buy = [0.40 0.40 0.40 0.40 0.40 0.40 1.20 1.20 1.20 1.20 ... 0.75 0.75 0.75 0.75 1.20 1.20 1.20 1.20 1.20 1.20 ... 0.75 0.75 0.75 0.75]; price_sell = repmat(0.40, 1, T); % 上网电价 % CHP 参数 Pchp_max = 200; Pchp_min = 50; % kW rho_chp = 1.2; % 热电比 a_chp = 0.25; % 边际气耗，m3/kWh b_chp = 2.5; % 空载气耗，m3/h gas_price = 2.5; % 元/m3 % 电锅炉 P_eb_max = 150; eta_eb = 0.95; % 储电 P_ch_max = 150; P_dis_max = 150; E_bess_max = 600; E_bess_min = 60; E_bess_init = 120; eta_bess_ch = 0.95; eta_bess_dis = 0.95; % 储热 H_tes_ch_max = 200; H_tes_dis_max = 200; E_tes_max = 1000; E_tes_min = 100; E_tes_init = 200; eta_tes_ch = 0.90; eta_tes_dis = 0.90; % 联络线 P_buy_max = 800; P_sell_max = 500;

这些参数不需要一次到位，第一次跑通后再根据设备手册细化。重点是让数据格式统一、单位一致。

3.2 基于YALMIP的核心建模代码

变量声明和约束整合可以这样写。YALMIP 的sdpvar定义连续变量，binvar定义 0-1 变量，optimize负责调用底层求解器。

% 决策变量 P_chp = sdpvar(1, T, 'full'); H_chp = sdpvar(1, T, 'full'); on_chp = binvar(1, T, 'full'); P_eb = sdpvar(1, T, 'full'); H_eb = sdpvar(1, T, 'full'); P_ch = sdpvar(1, T, 'full'); P_dis = sdpvar(1, T, 'full'); E_bess = sdpvar(1, T+1, 'full'); H_tes_ch = sdpvar(1, T, 'full'); H_tes_dis = sdpvar(1, T, 'full'); E_tes = sdpvar(1, T+1, 'full'); P_buy = sdpvar(1, T, 'full'); P_sell = sdpvar(1, T, 'full'); Cons = []; % CHP 出力区间与热电耦合 Cons = [Cons, Pchp_min * on_chp <= P_chp <= Pchp_max * on_chp]; Cons = [Cons, H_chp == rho_chp * P_chp]; % 电锅炉 Cons = [Cons, 0 <= P_eb <= P_eb_max]; Cons = [Cons, H_eb == eta_eb * P_eb]; % 储电约束 Cons = [Cons, 0 <= P_ch <= P_ch_max]; Cons = [Cons, 0 <= P_dis <= P_dis_max]; Cons = [Cons, E_bess_min <= E_bess <= E_bess_max]; Cons = [Cons, E_bess(1) == E_bess_init]; Cons = [Cons, E_bess(T+1) == E_bess_init]; Cons = [Cons, E_bess(2:T+1) == E_bess(1:T) + dt * (eta_bess_ch * P_ch - P_dis / eta_bess_dis)]; % 储热约束 Cons = [Cons, 0 <= H_tes_ch <= H_tes_ch_max]; Cons = [Cons, 0 <= H_tes_dis <= H_tes_dis_max]; Cons = [Cons, E_tes_min <= E_tes <= E_tes_max]; Cons = [Cons, E_tes(1) == E_tes_init]; Cons = [Cons, E_tes(T+1) == E_tes_init]; Cons = [Cons, E_tes(2:T+1) == E_tes(1:T) + dt * (eta_tes_ch * H_tes_ch - H_tes_dis / eta_tes_dis)]; % 购售电 Cons = [Cons, 0 <= P_buy <= P_buy_max]; Cons = [Cons, 0 <= P_sell <= P_sell_max]; % 电功率平衡 Cons = [Cons, pv + wt + P_chp + P_dis + P_buy == P_load + P_eb + P_ch + P_sell]; % 热功率平衡 Cons = [Cons, H_chp + H_eb + H_tes_dis == H_load + H_tes_ch]; % CHP 爬坡约束 Cons = [Cons, -60 * dt <= P_chp(2:T) - P_chp(1:T-1) <= 60 * dt];

目标函数里包含购电成本、售电收益、天然气成本、运维成本。这里的光伏运维成本我按实际出力pv算，模型里可以当成常数，不会影响优化变量最优解，但会改变总成本数值。

% 目标函数：购电成本 - 售电收益 + 天然气成本 + 运维成本 gas_cost = sum((a_chp * P_chp + b_chp * on_chp) .* gas_price) * dt; om_cost = sum(0.02 * P_chp + 0.01 * P_eb + 0.008 * (P_ch + P_dis) ... + 0.005 * (H_tes_ch + H_tes_dis) + 0.01 * pv) * dt; Objective = sum(P_buy .* price_buy) * dt - sum(P_sell .* price_sell) * dt ... + gas_cost + om_cost; % 求解 options = sdpsettings('solver', 'gurobi', 'mipgap', 0.001, 'verbose', 1); tic; sol = optimize(Cons, Objective, options); toc;

如果没有 Gurobi 或者 Cplex，把solver改成intlinprog也可以运行。只是大型 MILP 问题用内置的 intlinprog 求解速度会慢不少，我建议至少配置 Gurobi 试一下。

3.3 求解器配置与收敛经验

YALMIP 是建模层，真正算 MILP 的是底层求解器。Gurobi 和 Cplex 对线性整数规划的求解性能很强，适合 24 时段这种规模。对纯线性规划问题，用linprog就够，但一旦出现binvar，必须走整数规划路线。

sdpsettings里我惯用的几个参数：

options = sdpsettings('solver', 'gurobi', 'mipgap', 0.0005, 'verbose', 1); options.gurobi.TimeLimit = 300; options.gurobi.MIPFocus = 1;

mipgap设成 0.0005 意思是允许最优解的相对误差在 0.05% 以内。对于运行成本优化，这个精度足够。TimeLimit是防止模型规模过大时无限算下去。如果发现求解时间和最优性之间不平衡，可以适当放宽 MIPGap，比如 1%，结果往往只差几十块钱，但求解时间可能缩短一个数量级。

调试时建议先用optimize(Cons, [], options)只做可行性求解，不优化目标函数。如果这个可行解都找不到，说明约束之间有矛盾。接下来用check(Cons)查看约束残差，YALMIP 会返回每个约束的 residual，一眼能看出哪类约束冲突。

4. 仿真结果怎么看：典型日场景下的方案对比

4.1 结果输出与绘图思路

求解完成后，需要把 YALMIP 的变量取值提取出来：

P_chp_opt = value(P_chp); H_chp_opt = value(H_chp); on_chp_opt = value(on_chp); P_eb_opt = value(P_eb); P_ch_opt = value(P_ch); P_dis_opt = value(P_dis); E_bess_opt = value(E_bess); H_tes_ch_opt = value(H_tes_ch); H_tes_dis_opt = value(H_tes_dis); E_tes_opt = value(E_tes); P_buy_opt = value(P_buy); P_sell_opt = value(P_sell);

绘图时我喜欢用两张堆叠图看全天电平衡和热平衡。电平衡图横轴是小时，纵轴是功率，堆叠项包括光伏、风电、CHP、储能放电、购电；热平衡图堆叠项包括 CHP 供热、电锅炉供热、蓄热罐放热。一张图能直观看出各时段谁在出力、谁在充电、谁在买电。

figure; bar(1:T, [pv; wt; P_chp_opt; P_dis_opt; P_buy_opt]', 'stacked'); hold on; plot(1:T, P_load, 'k-', 'LineWidth', 2); legend('光伏', '风电', 'CHP', '储电放电', '购电', '电负荷'); xlabel('时段/h'); ylabel('功率/kW');

4.2 联合调度结果解读：何时充电、何时储热、何时购电

跑通后的典型日优化结果大致遵循这样一个规律：夜间低谷电价时段，系统倾向于多购电，一部分直接供电负荷，一部分给蓄电充电，同时电锅炉在低谷电价时段启动制热并存入蓄热罐。白天光伏大发、电负荷不是特别高时，CHP 可能降出力，避免产生多余的热，蓄电在光伏顶峰时充电、在晚高峰放电。晚高峰电价最高，CHP 满发，储电放电，蓄热罐放热，购电最少。

一个比较典型的运行安排：

这里需要特别理解“储热”的价值。白天如果光伏较多，而电负荷不够消化，CHP 再满发就会导致多余电力。但如果你把 CHP 在白天多发，一部分电可以卖给电网，附带的热量存入蓄热罐留到晚上用，整体账可能更划算。这种跨时段套利关系，只有在电热联合模型里才能被算出来。

5. 调试过程中最常踩的坑：规模失控与热负荷预测误差

5.1 0-1变量爆炸：从573秒到2.8秒的重构经验

第一次完整建模时，我抱着“约束越多越严谨”的想法，给很多关系都加了二进制变量。除了 CHP 启停，还加了储电充放互斥、储热充放互斥、购电售电互斥，结果二进制变量数量变成 4*T=96 个。看起来不多，但求解时间长达 573 秒，MIPGap 还卡在 1% 下不去。

后来我逐个删掉了“经济上不会发生”的互斥变量。储存设备同时充放会同时产生两笔运维成本，还要承担效率损耗，优化器没有任何理由这么做；购售电同时发生在分时电价下也不会省成本。保留 CHP 启停变量后，二进制变量降为 24 个，加上约束结构从大 M 型变成更紧凑的边界约束，最终求解时间降到 2.8 秒。

这个经验可以总结成一条原则：每一个二进制变量都要问自己，它在真实物理场景中是否必须存在、删掉它是否真的会出荒谬解。如果单纯为了理论严谨而加一堆 0-1 变量，最终拖累的是求解性能和模型实用性。先用线性规划跑通，再加少量整数变量，是更省时间的路径。

5.2 热负荷预测误差导致的结果敏感性

另一件让我印象很深的事是热负荷预测误差。某天早上我看当地气象台温度和历史热负荷数据，估了一条偏低 20% 的热负荷曲线。优化结果出来，蓄热罐夜间蓄的热量明显偏少。到了下午热负荷实测比预测高，系统只能临时让 CHP 多发电、多供热，还把储电提前放出，晚高峰购电骤增，当天实际成本比用准确热负荷预测跑出来的结果高了 12%。

电负荷预测误差对结果的影响被讲得很多，但热负荷预测误差其实更隐蔽。因为热系统惯性大，储热时间常数长，一旦蓄热不足，后面的补救手段非常有限。

我在工程上用了两个办法。第一个是滚动时域优化，每隔一小时用最新实测数据重新优化未来 24 小时，把预测误差的影响限制在单步之内。第二个是在热负荷预测值上人为加一个 5% 到 10% 的安全系数，或者强制蓄热罐末端 SOC 不低于 20%。这两种做法都会让日运行成本略微上升，但换回来的运行可靠性提升远大于那点成本。

6. 从这个模型还能延伸出什么

6.1 加入不确定性：鲁棒优化和场景随机优化

确定性的电热联合调度模型可以直接用，但园区微电网里光伏、风电、热负荷都有不确定性。如果想更进一步，一个方向是做鲁棒优化。对光伏出力和热负荷分别设置盒式不确定集，模型在“最坏情况”下仍然满足电热平衡。代价是最优成本会比确定性模型高，因为系统要为最坏情况预留更多备用容量。

另一个更贴近实际的方向是场景法随机优化。从历史数据抽 N 个典型日场景，每个场景分配一个概率，目标函数改成所有场景总成本的期望值。决策变量分成第一阶段和第二阶段：第一阶段的 CHP 启停状态、储电充放决策在场景之间保持一致，第二阶段功率变量随场景调整。这个模型比鲁棒优化更精细，但计算量也翻倍。

用 YALMIP 做随机优化时，可以把场景展开成一个超大的时间序列，用向量索引对应不同场景的位置，这样仍然可以保持单次optimize求解。只是约束矩阵规模会随着场景数量线性膨胀。

6.2 热网时空特性建模的扩展方向

目前的模型把热力系统看成单节点，所有热负荷、产热设备、蓄热罐都汇聚到一个热母线。如果园区里存在几百米甚至更长的供热管网，水温传输延迟和管网热损失就不能再被忽略。

一种实用的扩展是把热网简化为带延迟的蓄热模型：热源到热用户之间的管道相当于一个天然的大蓄热罐，具有管存效应。另一种更精确的做法是建模供水温度和回水温度，用二阶锥约束逼近管网水力特性，模型从 MILP 变成 MISOCP，求解器则需要支持二阶锥整数规划。

从工程落地角度看，单节点电热联合调度模型最适合项目前期的容量配置分析和运行策略研究。等真正到了系统设计阶段，再把热网网络约束加进去，两者在思路上一脉相承。

最后分享一个我自己的土办法：跑完优化之后，不要直接信任结果，先把各变量取值手动带进电平衡和热平衡等式，逐时段算一遍残差，再用check(Cons)核对。我遇到过的大部分异常结果，最后追查下来都是效率参数写反了、SOC 索引错位了，或者某台设备的单位没换算。先固定 CHP 启停状态只做线性规划，检查结果合理后再放回二进制变量，是排查模型错误最省力的顺序。这个习惯帮我省下过大量时间。