做综合能源系统优化调度这些年，我见过太多代码跑得飞起、图表画得漂亮，但实际落地时用户根本不买账的项目。原因很简单，绝大多数调度模型只盯着电费和气费的账本，把室内温度当成一个可以随便推向约束边界的“软柿子”——温度贴着下限跑，运行成本是好看了，人在屋里冻得直哆嗦。这个MATLAB项目有意思的地方在于，它把PMV（预测平均投票）这个人体热舒适度指标塞进了冷热电多能互补优化调度的模型里，让系统在省钱的同时还要保证人待着舒服，这正好踩中了综合能源系统从“设备侧最优”走向“用户侧满意”的大趋势。
这篇文章适合正在做综合能源系统优化调度、楼宇能量管理、或者想用MATLAB把热舒适度指标落地成代码的同学参考。我会从PMV模型的工程简化讲起，把冷热电联供系统的物理架构和调度模型约束一条条拆开，再给出一套基于YALMIP的可行代码方案和参数配置经验。全文不会贴大段完整源码，但核心建模思路、关键约束怎么写、求解器怎么配，都会讲到位。
1. 为什么非要把“体感舒适度”写进调度模型
综合能源系统调度本质上是一个数学优化问题：给定设备特性、能源价格、负荷曲线，求一组设备出力方案，让运行成本最低。传统模型里，“用户需求”通常被简化成一条固定的电负荷、热负荷、冷负荷曲线，系统只需要“按需生产”就行。但实际场景里，热负荷和冷负荷恰恰是最具备调节弹性的那部分需求——你设定26℃空调和设定22℃空调，耗电量能差出30%。如果调度模型完全不碰温度，只被动地满足“设计负荷”，等于放弃了最大的省钱空间。
反过来，如果只为了省钱而把温度往极限推，又会出现前面说的“冻得哆嗦”问题。这个矛盾在冷热电联供系统里尤其突出：燃气轮机发电会产生大量余热，余热既可以供热也可以驱动吸收式制冷机供冷。当电价特别低的时候，系统倾向于多买电、用高效电制冷机供冷；当气价相对便宜时，则更愿意让燃气轮机多发电、余热多供能。这种“换能源品种”的决策会直接改变送到房间里的冷量或热量，从而影响室内温度。如果模型里没有舒适度约束，优化算法就一定会找出一个把室内温度压到约束边界的解，因为那才是成本最低点。
PMV的价值就在于，它比单纯的温度区间约束更贴近人体的真实感知。同样25℃的室温，梅雨季湿度80%、空气静止不动，人会感觉闷热难受；而湿度40%、有轻微吹风感时，反而会觉得凉爽舒适。PMV把这几个因素全部折算成一个从-3到+3的标量，0代表“中性”（最舒适），-3代表“很冷”，+3代表“很热”。工程上通常认为PMV在-0.5到+0.5之间，90%以上的人体感是满意的。
这个项目把PMV放进调度模型，等于给系统装了一个“以人为本”的质量监测器。你不再问“室温是不是在18到26℃之间”，而是问“这24小时里，人是不是始终处于热中性状态”。这显然是一个更高级、也更难实现的优化目标。
2. PMV模型的原理与工程化改造
2.1 PMV的物理含义和原始公式
PMV是丹麦学者Fanger在上世纪八十年代提出的热舒适评价指标，后来被ISO 7730标准采纳。它不是一个可以直接测量出来的物理量，而是综合了六个参数的函数：空气温度ta、平均辐射温度tr、相对湿度RH、空气流速v、人体代谢率M和服装热阻Icl。为了让读者有个直观概念，我给出它的基本形式：
PMV = [0.303·exp(-0.036·M) + 0.028] · { (M - W) - 3.05e-3·[5733 - 6.99·(M - W) - Pa] - 0.42·[(M - W) - 58.15] - 1.7e-5·M·(5867 - Pa) - 0.0014·M·(34 - ta) - 3.96e-8·fcl·[(tcl + 273)⁴ - (tr + 273)⁴] - fcl·hc·(tcl - ta) }
这式子看一眼就头疼，但真正让程序猿崩溃的不是公式长度，而是里面藏着两处迭代计算：服装表面温度tcl和对流换热系数hc互相耦合，必须联立求解。也就是说，PMV本身是一个隐函数，没法直接以解析表达式的方式写进约束条件。
好在工程上有成熟的处理办法。在冷热电联供建筑场景中，代谢率M可以近似固定（办公室取1.0到1.2 met），风速v在空调房间里也基本稳定（0.1到0.15 m/s），平均辐射温度tr可以近似等于室内空气温度ta。这样PMV就变成了一个只跟室内温度ta、相对湿度RH、服装热阻Icl有关的函数。再进一步，如果你在模型中把湿度假设为固定值（很多建筑模型都这么做），那PMV就退化成室温和服装热阻的单变量函数。
2.2 固定参数下的线性化近似
我先说结论：在舒适区附近，PMV与室内温度的关系几乎是一条直线。以冬季工况为例：
- M = 1.0 met（静坐办公），W = 0
- Icl = 1.0 clo（厚外套加毛衣）
- v = 0.1 m/s，RH = 50%，tr = ta
用完整公式算几个关键点：ta=22℃时PMV约-0.5，ta=24℃时PMV接近0，ta=26℃时PMV约+0.5。于是可以拟合出：
PMV ≈ 0.25 × (ta - 24)
这个式子说明，在这个特定工况下，室温每升高1℃，PMV大约上升0.25。夏季工况换薄衣服（Icl=0.5 clo），中性温度会上移，类似地可以拟合出PMV ≈ 0.25 × (ta - 26)。不同文献里系数会有一点差异，但线性关系的形式是稳定的。
有人可能担心这种近似会不会误差太大。说实话，在[-0.5, +0.5]这个舒适区间内，线性近似的误差一般在0.1以内，作为优化约束足够用了。真要追求更高精度，可以用MATLAB的cftool工具箱对PMV做个分段线性拟合，把室内温度范围分成几段，每段用一条直线逼近。这种分段线性函数在YALMIP里可以用二元变量轻松建模，求解效率完全不受影响。
2.3 舒适区间的取值策略
舒适区间取多宽，直接决定优化问题的可行域和经济性。取值建议如下：
| 区间设置 | 对应PMV范围 | 适用场景 | 备注 |
|---|---|---|---|
| 严格中性 | [-0.2, +0.2] | 高性能办公、医院病房 | 模型可行域极小，成本显著上升 |
| 标准舒适 | [-0.5, +0.5] | 普通办公楼 | ISO 7730推荐的90%满意区间，推荐默认 |
| 宽松舒适 | [-0.8, +0.8] | 工厂车间、仓储 | 成本控制优先，用户容忍度较高 |
| 季节动态 | 冬季[-0.5, 0]，夏季[0, +0.5] | 节能优先的公共建筑 | 允许冬令略冷、夏令略热，节能效果明显 |
我个人的经验是，第一版模型务必先取[-0.5, +0.5]跑通，确认无解排除逻辑问题后，再根据需要放宽或收紧。千万别一上来就用[-0.2, +0.2]，否则大概率会得到一个infeasible problem，还得回头排查是哪个约束打架。
3. 冷热电多能互补系统结构与调度模型构建
3.1 系统架构与能量流
要建立优化调度模型，首先得把系统的能量流画清楚。这里以一个典型的楼宇级冷热电联供系统为例，主要设备包括：
- 燃气轮机：消耗天然气发电，同时产出高温余热
- 余热回收装置：把燃气轮机排烟余热转换成热水或蒸汽
- 吸收式制冷机：利用余热驱动制冷，热力系数COP_ac约1.2到1.4
- 电制冷机：消耗电力制冷，COP_ec约3.5到5.0
- 燃气锅炉：补充供热，效率约0.9
- 电储能/蓄热/蓄冷装置：平移能量供需
能量输入有两路：从电网购电、从管网购天然气。能量输出有三路：供电负荷、供热负荷、供冷负荷。
这个系统里最有意思的权衡发生在制冷环节。表面看电制冷机COP高达4以上，吸收式制冷机才1.3，似乎电制冷完胜。但别忽略电的来源：如果电来自燃气轮机，假设发电效率η_gt=0.35，那么“天然气→电→冷”的一次能源效率只有0.35×4=1.4；而余热回收效率约0.45，“天然气→余热→冷”的效率是0.45×1.3≈0.585。这么一看余热制冷似乎又很差。关键点在于：余热是燃气轮机的“副产品”，如果供热负荷不足时不用来制冷，就只能白白排掉，边际成本是零。所以实际决策依赖热负荷大小、电价气价比值、设备特性等多个因素——这正是优化调度要解决的权衡。
3.2 目标函数与决策变量
目标是日运行总成本最小，主要包括购电费用、购气费用和可能的售电收益：
min Σ [P_grid(k)·price_grid(k) + V_gas(k)·price_gas - P_sell(k)·price_sell(k)]
这里k是时段序号，取1到24。P_grid是购电功率，V_gas是天然气消耗量，P_sell是余电上网功率。
决策变量包括：
- 各时段燃气轮机发电功率P_gt(k)
- 燃气轮机启停状态（0/1二进制变量）
- 电网购电功率P_grid(k)
- 燃气锅炉供热量Q_gb(k)
- 吸收式制冷机供冷量Q_ac(k)和电制冷机耗电功率P_ec(k)
- 储能的充放功率和状态
- 室内温度T_in(k)
- 供热/供冷功率Q_hvac(k)
3.3 约束条件的几个关键块
第一块是能量平衡约束。电力平衡要求购电、燃气轮机发电、储能放电之和，等于电负荷加上电制冷机耗电、储能充电。热平衡要求余热回收热量、燃气锅炉供热、蓄热放热之和，等于热负荷加上蓄热充热。冷平衡类似。这几条约束是任何综合能源调度模型的骨架。
第二块是设备运行约束。每个设备都有出力上下限、爬坡速率限制，燃气轮机还有最小启停时间约束。储能设备需要额外加上SOC状态更新方程以及周期始末SOC一致的约束，否则模型会在一个调度周期里把储能“榨干”，算出来的结果第二天根本无法延续。
第三块是PMV舒适度约束和建筑热动态方程的耦合，这是本项目的核心难点。建筑热动态采用一阶等效热参数模型，把房间等效成一个热容C_b加一个热阻R。微分方程形式是：
C_b × dT_in/dt = Q_hvac + (T_out - T_in) / R
离散化成递推式后变成：
T_in(k+1) = a × T_in(k) + b × (R × Q_hvac + T_out)
其中a = exp(-Δt/τ)，b = 1 - a，τ = R×C_b是建筑时间常数。室内温度T_in再通过PMV线性近似进入舒适度约束：
-0.5 ≤ 0.25×(T_in - 24) ≤ 0.5（冬季工况举例）
这一串约束把“设备出力→室内温度→人体舒适度”完整串起来了，优化算法在决定功率分配时，必须同时考虑对室温的连锁影响。
3.4 建筑热惯性：被低估的“免费储能”
我特意把建筑热惯性单独拿出来讲，因为它在实际项目中价值非常大，但很多代码里根本没有体现。借用上面的递推式做一个实际估算：假设100m²的普通办公室，综合传热系数UA为80 W/K，等效热容C_b取20000 kJ/K（约200 kJ/(m²·K)）。
热阻R = 1/80 = 0.0125 K/W，时间常数τ = R×C_b = 0.0125×20000000 = 250000秒，约69小时。这个数字意味着，你把供热功率降下来之后，室温不会马上掉，而是以十小时为尺度缓慢变化。
我实际算过：白天电价高峰时段把供热功率从2kW降到0.5kW，持续4小时。热功率少了1.5kW，稳态温差会减少18.75℃，但因为时间常数是69小时，4小时只完成了约5.6%的过渡过程，室温实际只下降了大约1℃。PMV只掉0.25，用户几乎察觉不到，但系统却实打实地躲过了4小时的高价电。这就是建筑热惯性的“免费储能”价值。调度模型如果能把这个机制写进去，经济性提升空间非常可观。
4. MATLAB代码实现与求解配置
4.1 代码工程结构规划
拿到这个项目，别急着写代码。先把文件结构定了，后面调试会省很多事。我习惯的划分方式如下：
| 文件/文件夹 | 作用 |
|---|---|
| main.m | 主程序：参数初始化、调用建模、求解、后处理 |
| data_case.m | 输入数据：分时电价、气价、负荷曲线、设备参数 |
| pmv_calc.m | PMV计算函数（线性近似或完整公式） |
| build_model.m | 用YALMIP定义决策变量、目标函数和约束 |
| solve_model.m | 调用CPLEX/GUROBI求解并检查收敛状态 |
| plot_results.m | 绘制设备出力、温度曲线、PMV曲线、成本对比图 |
| results/ | 存放不同场景的求解结果，方便后续对比分析 |
4.2 核心建模代码思路
建模推荐用YALMIP工具箱加外部求解器。YALMIP只是一个建模语言，它能把高层的变量描述自动转换成求解器需要的标准格式。核心代码框架如下：
%% 定义决策变量 P_gt = binvar(1, 24); % 燃气轮机启停（若是连续功率则用sdpvar） P_gtd = sdpvar(1, 24); % 燃气轮机发电功率 P_ec = sdpvar(1, 24); % 电制冷机耗电功率 Q_gb = sdpvar(1, 24); % 燃气锅炉供热量 Q_ac = sdpvar(1, 24); % 吸收式制冷机制冷量 P_grid = sdpvar(1, 24); % 电网购电功率 T_in = sdpvar(1, 24); % 室内温度 Q_hvac = sdpvar(1, 24); % 供热/供冷功率 %% 目标函数 objective = sum(P_grid .* price_grid) + sum(gas_consume .* price_gas); %% 约束集合 Constraints = []; % 电力平衡约束 Constraints = [Constraints, P_grid + P_gtd*eta_gt ... == P_load + P_ec + P_bat_ch - P_bat_dis]; % 热平衡约束 Constraints = [Constraints, Q_gt_hr + Q_gb + Q_hs_dis ... == Q_load + Q_hs_ch]; % 冷平衡约束 Constraints = [Constraints, Q_ac + P_ec*COP_ec == Q_cool_load]; % 建筑热动态方程（核心耦合约束） for k = 1:24 if k == 1 Constraints = [Constraints, T_in(1) == T_in_0]; else Constraints = [Constraints, T_in(k+1) == a*T_in(k) ... + b*(R*Q_hvac(k) + T_out(k))]; end end % PMV舒适度约束（用线性近似） PMV = pmv_linear(T_in, season); % 返回1x24的PMV序列 Constraints = [Constraints, -0.5 <= PMV <= 0.5]; %% 求解 options = sdpsettings('solver', 'cplex', 'verbose', 2); optimize(Constraints, objective, options);这段代码有好几个细节值得展开。
电力平衡等式里，为什么电制冷机的耗电被放在等式右边和电负荷并列？因为电制冷机是耗电设备，不是“负的发电”。这样写有利于后续做灵敏度分析，也符合能量平衡的物理含义。
建筑热动态约束用了24个等式，把T_in和Q_hvac绑定起来。这里特别要注意，T_in(k)与Q_hvac(k)的时序关系不能搞反：Q_hvac(k)是本时段施加的供能功率，它影响的是T_in(k+1)。很多初写者在这里把下标错位，导致结果出现一个小时的相位偏移。
PMV约束那行看起来简单，但它实质上把T_in的上限和下限都隐含进去了。线性近似的斜率0.25，意味着PMV区间[-0.5, 0.5]对应的温度区间只有4℃，比普通温度约束要严格得多。这也是为什么必须配合建筑热模型一起用——没有热惯性缓冲，这4℃的窄窗口会把设备出力逼得非常生硬。
4.3 求解器配置与数值稳定
求解器方面，如果模型是混合整数线性规划（MILP），首选CPLEX或GUROBI。YALMIP本身不带求解器，必须在电脑上单独安装IBM ILOG CPLEX Optimization Studio或GUROBI，然后确保安装目录加入MATLAB路径。装好之后用yalmiptest命令检查求解器识别情况。
数值缩放是一个经常被忽略但极其重要的问题。假设功率单位用kW，购电成本一算就是几千上万，而温度约束里的数值只有20℃左右，两者数量级差了3到4个量级。很多求解器在这种条件下容易出现数值病态，表现为目标函数迭代不收敛或者约束被莫名其妙地违反。建议所有功率量统一用MW，价格单位统一用元/MWh，成本自然落到每小时几百元的量级，跟温度数值的差距缩小，求解稳定得多。
还有一个经验：不要一开始就把所有约束一次怼上去。先把能量平衡和设备约束跑通，确认目标函数有合理解；然后加储能约束；最后再加PMV和建筑热耦合，这样哪一步出了问题都能快速定位。
5. 结果分析：经济性和舒适度到底怎么权衡
5.1 三组对比实验设计
为了验证PMV舒适度约束的价值，建议做三组场景对比：
- Case 1：无舒适度约束，只有设备运行下限和能量平衡（纯经济调度）
- Case 2：固定温度约束，要求室温始终落在18℃到26℃
- Case 3：PMV舒适度约束，冬季工况PMV落在[-0.5, +0.5]
三组场景共用同一套负荷曲线、电价气价和设备参数，只改变约束条件。这样对比出来的成本差异完全由约束策略引起。
5.2 结果解读：成本只涨一点点，舒适度飞跃式提升
基于我自己的项目经验，结果通常会呈现这样几个特征：
Case 1的运行成本最低，但室内温度的模拟值会长时间贴着安全下限跑。在冬季，系统倾向于把室温压到15到16℃，PMV跌到-1.5以下，体感已经是“冷得难受”了。这个场景虽然省了钱，但完全不具备实际可操作性。
Case 2的成本比Case 1高大约8%到10%，因为18℃的下限和26℃的上限在冷热负荷高峰期会强制设备进入高成本运行区间。更麻烦的是，固定温度约束在过渡季节特别吃亏：明明湿度低、体感暖和，系统还是机械地维持18℃下限，白白浪费燃气。
Case 3的成本比Case 2低5%到12%，而PMV基本控制在±0.5以内。为什么PMV约束反而比固定温度省钱？关键在于冬季工况下PMV允许的温度下限其实低于18℃——比如湿度40%、风速0.15m/s、穿厚外套时，16℃室温对应的PMV可能也只有-0.4左右。PMV不会像温度约束那样一刀切地限制“温度必须高于18℃”，它给了调度系统在特定气象条件、特定着装条件下的合理“偷工减料”空间。
从工程角度这样表述：PMV约束的本质是把房间的“舒适度资源”动态化、场景化，而固定温度是一个静态的、僵化的判断。这就像限速60km/h的公路和根据天气、车流动态建议车速的智能交通系统相比，后者一定更高效。
5.3 图形输出要点
结果图建议至少包含四张：
- 各设备出力堆叠图（电、热、冷三个子图）
- 室内温度与室外温度对比曲线，叠加舒适温度带
- PMV逐时变化曲线，叠加±0.5的舒适带
- 三种场景的成本对比柱状图，或累计成本曲线
画出PMV曲线之后有个小技巧：把采暖季和供冷季分开画，或者用不同颜色区分过渡季节。因为冬季PMV和夏季PMV对应的温度区间差异很大，混在一张图里看容易误判系统是否存在越限。
6. 实操避坑与调试心得
6.1 高频问题速查表
| 问题现象 | 可能原因 | 解决方案 |
|---|---|---|
| 求解器报Infeasible problem | PMV约束过紧，或建筑热模型参数和负荷数据不匹配 | 先放宽PMV区间到[-1, 1]，跑通后再逐步收紧；检查T_in初值是否在合理范围 |
| YALMIP报No solver found | CPLEX/GUROBI没有正确安装或路径未加入 | 运行yalmiptest看求解器列表；重新添加路径 |
| 温度曲线剧烈振荡 | C_b取值过小，建筑热模型时间常数太短 | 增大C_b，保证时间常数在20小时以上 |
| 目标函数数值跳跃或NaN | 功率和成本量级不一致 | 统一单位：功率用MW，价格用元/MWh |
| SOC在周期末不为0 | 缺少SOC(24) == SOC(0)的周期约束 | 在约束集合里显式加上 |
| PMV线性近似导致温度超限 | 线性拟合的斜率或截距与实际工况不符 | 用pmv_calc函数校准关键点，更新拟合系数 |
6.2 我的几条独家调试经验
第一，PMV约束不要做成硬约束，最好加上松弛变量。具体做法是把目标函数里加一项“PMV偏离惩罚”：
f_penalty = λ × Σ max(0, |PMV| - 0.5)
YALMIP里可以用临时的非负辅助变量表示这个max项。带松弛的软约束求解稳定性远好于硬约束，而且当气象数据预测误差导致系统本无解时，软约束能自动退让到“尽量接近舒适”，而不是直接报无解。这个技巧在处理真实气象数据时价值极大。
第二，建筑热模型的参数标定千万别想当然。C_b和UA的取值对结果影响巨大，但它们又完全取决于具体建筑的围护结构。最稳妥的做法是先做一次开环仿真，只给固定的Q_hvac序列，看T_in的仿真曲线是否符合直觉。比如冬季夜间不供热时，室温应该缓慢下降而不是瞬间跌落。如果参数设置合理，这个仿真曲线本身就能讲出一个让人信服的“热惯性”故事。
第三，留意热负荷基数与PMV约束的匹配问题。冬季工况下PMV约束的温度下限可能低于传统的18℃下限，这会显著降低系统的热负荷需求，燃气轮机的余热可能因此供大于求。很多新手在这里会发现吸收式制冷机在寒冬也从“制冷”模式变成了“余热消纳”通道——这其实不是Bug，反而说明模型抓住了多能互补的精髓：余热不是非得用来供热，也可以驱动制冷机在冬季给数据中心或者特定区域供冷。优化算法会自动探索这种跨季节的能量分配策略，前提是你在约束里给了它足够的自由度。
做这个项目到最后，我最大的体会是：综合能源系统优化调度的难点从来不是设备建模或求解器调参，而是怎么定义“好东西”。把用户舒适度量化成PMV并放进约束，看似只是多写了一行不等式，实质上是把调度目标从“系统省钱”扩展到了“用户满意且系统高效”。这种思路在越来越多强调体验的智慧园区、零碳建筑项目中正在变成刚需。而随着空调负荷柔性调节、建筑虚拟储能这些概念落地，PMV这类人体热舒适指标与优化调度的结合点只会越来越多，这个方向值得深入研究。
