做微网优化的同行应该都有体会：光搭一个纯电系统的优化模型练手，半天就能跑通；但一旦把燃气锅炉、储热罐、余热回收、CHP机组全部拉进来，"多能互补"四个字马上就把问题复杂度抬高一个量级。热电联供型微网的核心矛盾在于，电和热两套系统在物理上强耦合，但在需求侧又各自独立变化，冬天热负荷飙升时电负荷不一定大，白天电负荷高峰时热负荷可能反而很低。这种电-热"剪刀差"，正是优化调度要解决的主要矛盾。

这篇文章我会完整拆解一套基于多能互补思路的热电联供型微网优化运行方案，给出可复用的Matlab建模方法、约束方程写法和求解思路，并附上我在实际调模型中踩过的坑。内容上不追求面面俱到，但保证每一步都符合工程习惯，适合正在做微网调度、综合能源系统优化、或者准备相关方向毕业设计的同学参考。

1. 项目思路与优化问题怎么定

1.1 热电联供微网里的"多能互补"到底指什么

先捋清楚概念。我们说的热电联供型微网（Combined Heat and Power Microgrid），本质是一个同时承担电负荷与热负荷的小型能源系统，常见组成包括：

• 燃气轮机或内燃机组成的CHP机组，发电同时回收余热
• 燃气锅炉作为备用热源
• 电锅炉或热泵作为电转热环节，实现"以电补热"
• 储电、储热装置，起到削峰填谷作用
• 外部电网购电通道，以及可能接入的光伏、风电等可再生能源

多能互补体现在两个层面：一是能源输入侧的多品种组合，气、电、热、可再生能源协同供能；二是转换环节的灵活性，当热负荷紧缺时可以从电网购电驱动电锅炉补热，当电价高、气价低时燃气机组多发电的同时多供热，这就是"互补"的本质。

很多新手在建模时容易把电系统和热系统拆开来独立做，然后结果看似收敛，实际不具备物理可行性。原因很简单：CHP机组的电出力和热出力共用同一个原动机，你不能让燃气轮机只发电不产热，也不能只产热不发电。这种强耦合关系必须写进约束里，否则优化结果就是空中楼阁。

1.2 目标函数选最小运行成本还是最小碳排放

优化运行的核心是目标函数。我见到的绝大多数研究都采用最小化系统日运行成本，因为这是工程上最直接、最可解释的指标。典型成本项包括：

1. 购电成本：从上级电网购电的电量电费
2. 燃料成本：燃气轮机和燃气锅炉的天然气消耗费用
3. 设备启停成本：机组启停带来的寿命损耗费用
4. 运维成本：按出力比例折算的设备维护费用

也可以把碳排放量纳入，做多目标优化。但我要提醒：多目标并不是简单地把成本和排放加权，碳定价法是工程上更合理的做法——给每吨CO2设定一个价格，放进成本函数里统一优化。这样既避免了加权系数的主观性，又能直接输出单一最优解。

我常用的目标函数形式如下：

% 目标函数：购电成本 + 燃料成本 + 碳成本 + 运维成本 objective = sum(C_grid.*P_grid) ... % 购电成本 + sum(C_gas.*F_chp + C_gas.*F_boiler) ... % 燃料成本 + sum(C_carbon.*(E_chp + E_boiler)) ... % 碳排放成本 + sum(K_om.*P_chp + K_om.*P_grid); % 运维成本

1.3 为什么选混合整数线性规划(MILP)框架

优化运行问题的建模框架有很多种选择：动态规划、遗传算法、粒子群、混合整数线性规划。我的建议是，除非你有特殊需求，否则优先选MILP。

原因有三。第一，MILP的全局最优性有数学保证，智能算法只能给近似解，在学术对比里容易被质疑；第二，MILP模型的约束可以直接对应物理逻辑，调试直观；第三，现在的商业求解器对MILP的处理能力已经很强，一台普通电脑上求解几百个变量的调度问题通常只需几十秒到几分钟。

使用MILP的前提是，目标函数和约束是线性的。这个条件在实际中基本可以满足：成本曲线可以分段线性化，设备效率可以取常数，机组启停用0/1变量表示。光伏和负荷预测的波动可以通过场景法处理，而场景法本身也只是把问题规模变大，不改变MILP框架。

2. 核心约束条件建模：从物理规律到数学公式

2.1 电功率平衡和热功率平衡，两条线不能乱

优化模型的最基本约束是功率平衡。简单说，任何时刻系统发出的电功率必须等于电负荷加上损耗和储能充放电功率。写成方程就是：

% 电功率平衡约束（T为调度时段数） Constraints = [Constraints, P_grid + P_chp + P_pv + P_dis == P_load + P_ch + P_eb]; % 注意：P_eb是电锅炉消耗的电功率

热功率平衡类似，所有热源产生的热功率必须覆盖热负荷：

% 热功率平衡约束 Constraints = [Constraints, H_chp + H_boiler + H_dis == H_load + H_ch];

这里有个细节容易出错：储热装置也有充放热状态，所以热平衡里同样要有储热充放变量。而且储电、储热不能同时充放，这需要加约束：

% 储能不能同时充放（M为足够大的数，用于Big-M法） Constraints = [Constraints, P_ch <= M.*u_ch, P_dis <= M.*u_dis, u_ch + u_dis <= 1];

2.2 CHP机组电热耦合可行域：关键中的关键

CHP机组是整个模型里最有技术含量的约束部分。燃气轮机的电出力和热回收之间不是简单的线性比例关系，而是存在一个凸多边形可行域。

对于典型的抽汽式CHP机组，在电-热功率平面上，可行运行区域可以描述为几个不等式围成的区域：

% CHP机组可行域约束（以最小出力-最大出力-最大热出力围成的凸多边形为例） % 变量：P_chp为电出力，H_chp为热出力 Constraints = [Constraints, P_chp >= max(P_min, 0.4*H_chp + 15)]; % 最小电出力随热出力抬升 Constraints = [Constraints, P_chp <= P_max - 0.2*H_chp]; % 最大电出力随热出力下降 Constraints = [Constraints, H_chp <= H_max]; % 最大热出力限制

这里的系数0.4和0.2来自机组的热电特性曲线，实际工程数据需要通过机组厂家提供的工况图取值。可行域内的任意点代表一个可运行状态，优化算法会在其中寻找最优工作点。

要点是：约束边界必须能反映"发热时发电能力会受限"这个物理规律。很多初稿模型只写了P_chp在一定区间内、H_chp在一定区间内，但把两者写成相互独立的变量，这会让优化结果出现"又发电又大量发热"的不合理工况。

2.3 储能约束和机组爬坡约束：别忘了时间耦合

储能设备的约束要体现跨时段的状态转移。储电设备的电量（SOC）按如下方式递推：

for t = 2:T Constraints = [Constraints, SOC(t) == SOC(t-1) + P_ch(t)*eta_ch - P_dis(t)/eta_dis]; end % SOC要限制在允许范围 Constraints = [Constraints, SOC_min <= SOC <= SOC_max]; % 初始和末尾SOC保持一致（循环调度常用） Constraints = [Constraints, SOC(1) == SOC_init, SOC(T) == SOC_init];

爬坡约束刻画了机组出力调整的物理限制。燃气轮机每分钟能升多少负荷是有限的，在15分钟粒度下有如下约束：

% 爬坡约束（上调/下调速率限制） Constraints = [Constraints, P_chp(t+1) - P_chp(t) <= ramp_up]; Constraints = [Constraints, P_chp(t) - P_chp(t+1) <= ramp_down];

我补充一个容易被忽略的细节：如果做了机组启停优化，还要加入启停相关的逻辑约束，比如"停机状态下出力为0"和"启动后最小运行时间限制"。这一组约束才是MILP比连续LP复杂的地方。

2.4 可再生能源出力的不确定性：场景法处理

很多论文会加入光伏和风电，那么如何处理它们的间歇性就成了绕不开的问题。最常见做法是采用多场景随机优化：生成若干个典型日的出力曲线，给每个场景赋一个概率，目标函数变成所有场景下的期望成本最小化。这样模型从确定性优化变成两阶段随机优化，规模会成倍增长。

如果不想一开始就碰随机优化，我建议先跑确定性模型，把光伏出力当成已知曲线来处理，先跑通整个框架。等确定性版本稳定后，再扩展为场景法。这个循序渐进的做法在硕士论文和实际工程项目里都适用，可以避免一上来就被不确定性和求解规模困住。

3. Matlab实现细节：代码框架与核心模块

3.1 建模工具怎么选：YALMIP还是直接用求解器API

Matlab环境下做优化建模，主流方案是YALMIP工具箱+第三方求解器。YALMIP提供了高层的建模语法，让你能像写数学表达式一样写约束和目标函数，然后自动翻译成求解器可识别的格式。这样做的好处是：代码结构清晰，想换求解器只需改一行配置。

当然也可以用linprog和intlinprog等原生函数，但需要手动把约束写成矩阵形式。变量一多，矩阵维度很容易错乱，调试起来非常痛苦。除非你的模型极小，否则我不推荐手写矩阵形式。

求解器方面，我实测过的几款如下表：

我日常用Gurobi最多，求解速度快，数值稳定性好。如果是学生做课程设计，先装SCIP免费版也能满足大部分场景。

% 在Matlab中设置YALMIP调用Gurobi求解 options = sdpsettings('solver','gurobi','verbose',1,'showprogress',1);

3.2 代码结构设计：模块化拆分开，别写成一坨

我建议按下面的模块划分代码目录，这是我从几个落地项目里沉淀下来的习惯：

1. data.m：输入数据，包括负荷曲线、能源价格、设备参数、预测曲线
2. model_vars.m：定义所有优化变量，连续变量用sdpvar，0/1变量用binvar
3. constraints.m：按类别添加约束（功率平衡、机组、储能、爬坡、备用）
4. objective.m：构建目标函数
5. run_optimization.m：主程序，依次调用上述模块并输出结果
6. plot_results.m：结果可视化

变量命名保持"类型_设备_时段"的规则，例如P_chp(t)、H_boiler(t)、SOC_ees(t)。这种命名方式在约束写多了之后极大降低排查成本。

调度周期用1小时为步长，一天24个时段起步；如果做高精度调度，可以细化到15分钟96个时段。每个时段的决策变量包括CHP电出力、CHP热出力、锅炉热出力、电锅炉电功率、储电充放、储热充放、购电功率等，加起来大约10个变量乘以96个时段，即约1000个变量，其中一小部分是0/1变量。这个规模对现代求解器很轻松。

3.3 一个最小的可运行模型，直接套用

下面给一个简化版模型核心代码。它以一天24小时为周期，包含CHP机组、燃气锅炉、电锅炉、储能、购电通道，目标是成本最小化。复制即用，体验一下完整流程。

%% 最小可运行的热电联供微网优化模型 T = 24; P_load = [数据导入]; % 电负荷曲线，1x24 H_load = [数据导入]; % 热负荷曲线，1x24 Price = [数据导入]; % 分时电价，1x24 C_gas = 3.4; % 天然气价格，元/m3 % 1. 定义变量 P_chp = sdpvar(1, T); % CHP电出力 H_chp = sdpvar(1, T); % CHP热出力 P_boiler = sdpvar(1, T); % 锅炉热出力（kW） P_eb = sdpvar(1, T); % 电锅炉电功率 P_grid = sdpvar(1, T); % 购电功率 SOC_e = sdpvar(1, T); % 储电SOC SOC_h = sdpvar(1, T); % 储热SOC u_chp = binvar(1, T); % CHP启停变量 % 2. 约束集合 Constraints = []; % 电功率平衡 for t = 1:T Constraints = [Constraints, P_grid(t) + P_chp(t) == P_load(t) + P_eb(t)]; end % 热功率平衡：CHP热出力 + 锅炉补热 == 热负荷 + 电锅炉耗热? 实际电锅炉耗电生热 % 注意：电锅炉产热功率记入热平衡，本简化模型中H_eb = COP * P_eb for t = 1:T Constraints = [Constraints, H_chp(t) + P_boiler(t) + COP_eb*P_eb(t) == H_load(t)]; end % CHP可行域 for t = 1:T Constraints = [Constraints, P_chp(t) >= 20.*u_chp(t) + 0.4.*H_chp(t)]; Constraints = [Constraints, P_chp(t) <= 100.*u_chp(t) - 0.2.*H_chp(t)]; Constraints = [Constraints, H_chp(t) <= 50.*u_chp(t)]; end % 储能SOC递推和容量限制 for t = 2:T Constraints = [Constraints, SOC_e(t) == SOC_e(t-1) + 0.9*P_ch_e(t) - P_dis_e(t)/0.9]; Constraints = [Constraints, SOC_h(t) == SOC_h(t-1) + 0.95*P_ch_h(t) - P_dis_h(t)/0.95]; end Constraints = [Constraints, 0.2 <= SOC_e <= 0.9, 0.2 <= SOC_h <= 0.9]; % 电网购电上限 Constraints = [Constraints, 0 <= P_grid <= 200]; % 3. 目标：购电成本 + 燃料成本（CHP和锅炉） F_chp = (P_chp + 0.5*H_chp) ./ 9.7; % 天然气消耗量，m3/h F_boiler = P_boiler ./ (0.9 * 9.7); % 锅炉气耗 objective = sum(Price.*P_grid) + C_gas*(sum(F_chp) + sum(F_boiler)); % 4. 求解 options = sdpsettings('solver','gurobi','verbose',1); optimize(Constraints, objective, options);

这个模型省略了储能充放电变量的定义，但约束框架是完整可扩展的。想完整跑通的话，把储能的充放变量补上，再补齐SOC递推逻辑即可。这是最小可行的模板，扩展性好。

3.4 结果可视化：学会读优化结果

运行完模型后，第一件事是画负荷平衡图。优化结果必须满足每个时段的功率平衡，而画图能直观看出电负荷由谁供给、热负荷由谁承担。我常用两张图：一张堆叠面积图展示电功率构成，另一张展示热功率构成。

figure; area(1:T, [P_grid', P_chp', P_eb'], 'LineWidth', 1); legend('购电','CHP发电','电锅炉耗电','Location','best'); xlabel('时段/h'); ylabel('功率/kW'); title('电功率平衡构成'); saveas(gcf, 'electric_balance.png');

热功率图同理。第三张图建议画SOC曲线，看储电储热是否处在合理区间。好的调度结果SOC曲线应该是平滑的锯齿形，峰谷时段有充放动作，而不是一条直线或者剧烈振荡。

4. 常见问题与排查技巧实录

4.1 求解器加载失败或许可证报错

YALMIP调用外部求解器时最常见的报错是“No suitable solver installed”，或者是License错误。我的建议，装求解器的时候务必把许可证路径配置好。Gurobi在Windows下通常安装到C盘，YALMIP会自动识别，但Mac和Linux下可能需要在MATLAB里手动执行：

% 手动添加Gurobi路径 addpath('C:\gurobi1100\matlab'); % 以实际安装路径为准 gurobi_setup() % 初始化尽调

如果是intlinprog自带的求解器，则是许可证随Matlab，没这个烦恼。但如果模型变量较多，自带求解器速度会明显下滑。

4.2 模型不可行，怎么定位是哪个约束引发的

不可行是MILP建模最容易踩的坑。解决办法是逐个约束加"松弛变量"：

% 把松弛量加到可能出问题的约束上 s = sdpvar(1, T); % 非负松弛变量 Constraints = [Constraints, P_grid(t) + P_chp(t) + s(t) == P_load(t) + P_eb(t)]; % 然后检查最优解里的s值，哪个时段的s不为0，问题就出在哪个时段

把目标函数暂时改为minimize sum(s)，求解后观察哪一时段松弛量最大，就可以反过来检查该时段是负荷高得离谱还是机组出力上限设太小。我用这个方法排查过很多次模型bug，比盲猜高效得多。

4.3 几个很值得留意的实操细节

第一，所有变量的单位务必统一。我做过一个项目，CHP发热量的单位用了GJ，而电出力的单位用了kWh，结果热平衡约束怎么检查都不收敛。统一成kW（功率）作为基准单位，每条约束都代入单位复核一遍，能省去大量排查时间。

第二，Big-M法里的M不能取太大。虽然理论上M取无穷大就行，但数值上M过大会造成病态矩阵，让求解器精度崩塌。M只取该变量物理上限的10倍左右，数值上已经足够。

第三，供热网络本身的"热惯性"在标准模型里一般忽略。如果你想做更精细的控制，则需要考虑热网管道蓄热效应，这时不能再用静态热平衡，要引入动态热网方程。

第四，目标函数里的平方项能用线性逼近就尽量线性。如果不小心写了二次项，模型会变成MIQP，求解难度立刻上一个台阶。除非确实需要，否则保持线性。

4.4 关于数据和场景的两点经验补充

运行调度模型，数据质量比模型本身更影响结果。我建议从公开数据平台抓取典型日数据，比如美国能源信息署(EIA)的部分公开负荷数据，或者国内高校开源数据仓库。自行构造假数据做演示可以，但如果后续要对标真实场景，务必找有据可查的数据源。

对于场景法处理不确定性，我建议先用“聚类法”把365天的数据聚成3-5个典型场景，每个场景对应一个概率。这样既保留了主要天气类型，又不至于让问题规模失控。随机优化在论文里的说服力比确定性优化明显更高，而且聚类实现起来不难。

5. 完整扩展方向的几条参考思路

前面把主干流程都讲完了，如果你还想把这套模型继续延展，我给出三个比较现实的方向。

第一是加入需求响应机制。通过引入可转移负荷和可削减负荷，让电负荷曲线从"刚性的"变成"柔性的"。做法是在约束里增加负荷转移变量和补偿成本，模型规模多出几十个变量，但能明显提升微网运行的经济性。

第二是引入热网管道动态模型。传统做法只做能源站内的调度，如果把管网的传输延迟、热损失、储热特性也建模，系统灵活性会再上一个台阶。这对大区域供热系统尤其有效，代价是模型的微分方程更多，求解时间变长。

第三是探索日前-日内两阶段调度。日前用场景法做随机优化确定机组启停和储能计划，日内用滚动优化跟踪实际负荷偏差。这种“先计划后调整”的结构是工程落地的主流做法，也比单阶段模型更有说服力。

这几个方向上，我踩过不少坑，比如日前计划用太激进会导致日内滚动优化频繁反推、热负荷预测不准导致储热罐放空等。实际工程里，计划一定要留裕度，运行边界宁可窄一点也不要在临界点跳舞。

6. 写给自己，也算写给正在做这个方向的你

我从第一次接触微网优化到现在，最深的感受是：优化模型的建成只是第一步，真正的功夫在“约束写得好不好”“结果能不能落回物理现场”。你写出的每一条不等式，都应该能对应到现场设备的某个调节阀或者某个功率上限，否则这个模型就只是个数学游戏。

热电联供微网优化运行这个课题，Matlab代码实现其实只是一个载体，核心还是对能源转换物理过程的理解深度。建议你做任何改动时都反问自己：如果把这个约束删掉，现场会不会出问题？如果把这个参数放大一倍，设备会不会扛不住？带着这种工程直觉去写模型，才能做出真正有价值的东西。

最后分享一个小技巧：每次跑完优化，把结果里的调度曲线导出来和实际运行曲线叠在一起看，不仅是看趋势像不像，还要看每个设备在峰值时有没有超出预期出力。有一次我发现CHP机组在优化结果中频繁启停，而实际根本不允许这么操作，就是加了最小开关机时间约束才解决问题。这类问题建模时想不起来，等对图横竖不对的时候自然而然就想到了。