做综合能源系统优化调度的人应该都有过这种经历优化程序跑出结果运行成本压得很低可一进入供暖季或空调季用户侧的投诉就来了——办公室冷得坐不住下午房间热得发闷。问题出在调度模型里只有经济账和设备约束没有人这个变量。我完整调试过一套MATLAB代码主题是考虑用户舒适度的冷热电多能互补综合能源系统优化调度核心用PMV预测平均投票指标衡量用户热舒适度再把舒适度和经济性放进同一个优化框架里权衡。这篇文章把整套代码从头拆一遍PMV模型怎么落地、冷热电多能互补的调度模型怎么搭、MATLAB实现里哪些细节最容易卡住人。不管你是做IES/CCHP课题的研究生还是搞园区综合能源系统优化的工程师只要被舒适度怎么进优化模型这个问题绊住过这篇都能帮上忙。1. 从省电要紧到人舒服要紧PMV指标进入调度模型的底层逻辑1.1 PMV到底算的是什么六参数与Fanger方程PMV全称Predicted Mean Vote预测平均投票最早由Fanger教授在20世纪70年代提出后来被ISO 7730标准采纳成为室内热环境评价的主流指标。它的输出是一根-3到3的标尺0代表中性热感觉最舒服正值偏热负值偏冷。ISO 7730给出的推荐范围是-0.5到0.5在这个区间内可以认为绝大多数使用者对热环境满意。Fanger方程的形式相当绕但拆开看并不复杂PMV [0.303×e^(-0.036M) 0.028] × [ (M - W) - 3.05×10⁻³×(5733 - 6.99(M - W) - pa) - 0.42×((M - W) - 58.15) - 1.7×10⁻⁵×M×(5867 - pa) - 0.0014×M×(34 - ta) - 3.96×10⁻⁸×fcl×((tcl 273)⁴ - (tr 273)⁴) - fcl×hc×(tcl - ta) ]其中M是人体代谢率单位W/m²1 met等于58.2 W/m²W是人体做机械功坐姿办公基本取0pa是环境水蒸气分压单位kPata是空气温度tr是平均辐射温度tcl是服装外表面温度由另一个隐式方程迭代求出fcl是服装面积系数hc是表面对流换热系数。再加上服装热阻Icl、空气流速va一共六个主要变量空气温度、平均辐射温度、相对湿度、风速、代谢率、服装热阻。这里有个非常关键的点tcl并不是直接从物理参数算出来的它同时出现在方程左右两边必须用迭代法求解。我后面写MATLAB函数时会专门处理这个循环依赖很多初学者在这里翻车算出来的PMV跳变、非单调一查全是tcl迭代没收敛。1.2 调度模型里为什么不直接套原版PMV理论公式很漂亮但把它直接搬进优化调度模型麻烦就大了。首先平均辐射温度tr在工程现场几乎不可能做到逐时准确测量它受窗面太阳辐射、墙体内表面温度、室内设备散热影响很大调度模型里如果把它当成一个变量那模型又多了一层复杂耦合其次风速va在空调房间内本来就不是恒定值送风方式一变局部风速就不同再次代谢率M和服装热阻Icl随人员活动状态、季节变化属于典型的不确定参数。如果你把这些全变量原封不动塞进一个以15分钟为时间步长、全天96个时段的大规模优化问题里那么每一个时段都要算一次嵌套迭代、还要考虑六参数之间的乘除耦合结果就是模型非线性程度极高MATLAB里fmincon这类求解器跑起来非常痛苦初始点稍微偏一点就收敛到局部最优甚至是不可行解。所以在工程落地时大家普遍的做法是做一次面向调度的降维把M固定成1.2 met对应办公室轻度坐着工作的状态Icl按季节取典型值夏季0.5 clo、冬季1.0 clo风速va假设为0.1 m/s的室内静风环境平均辐射温度tr近似等于空气温度ta。这样六个参数就压缩到了一个——空气温度ta成了唯一的主变量PMV从六维函数退化成单变量函数可以拟合、可以分段线性化、可以稳定地嵌进优化模型。1.3 工程化降维把PMV变成室温的舒适度曲线降维之后的PMV(ta)曲线仍然是非线性的但形状非常规律在15℃到35℃这个常见室内温度区间内PMV基本随ta单调上升而且近似线性。我实测过一组夏季典型工况M1.2 metIcl0.5 cloRH60%va0.1 m/strta用标准Fanger方程逐点算出来后做最小二乘拟合得到了一个很好用的线性关系PMV ≈ 0.28 × (ta - 25.7)这个式子的含义很直观夏季工况下中性舒适温度大约是25.7℃每升高1℃PMV上升约0.28。用它算几个典型值室温24℃时PMV约-0.48贴着ISO舒适区下边界26℃时PMV约0.08非常舒服28℃时PMV约0.64已经偏热了29℃时约0.92这已经是让人明显感到闷热的程度。有人可能会问既然最后都做成温度约束了为什么不直接用室温上下限非要绕一圈用PMV区别在于PMV把湿度、风速、服装热阻、代谢率都归一进了一个指标。同样是26℃的室温南方梅雨季节RH 85%和北方干燥环境RH 30%下人的体感完全不同同样一件T恤在空调房里和在没有风的房间里感受也不同。纯温度区间约束是设备侧思维PMV约束才是人体侧思维。而且在写论文、出报告、对接暖通标准时PMV是有ISO 7730背书的指标说服力比一句我们把室温控制在26℃强得多。2. 冷热电多能系统调度模型的数学骨架从母线平衡到目标函数2.1 系统拓扑与能量母线谁在由电生冷、由热生冷开始写代码之前先把系统拓扑想清楚。我这里讨论的冷热电多能互补综合能源系统至少包含这些部分供能设备燃气轮机或燃气内燃机CHP机组同时发电和热燃气锅炉作为补热源光伏作为可再生能源电源。制冷设备电制冷机耗电产冷吸收式制冷机耗热产冷利用CHP余热或者锅炉热。储能设备蓄电池、蓄热罐有条件的还会加蓄冷槽。负荷侧电负荷、热负荷采暖/生活热水、冷负荷空调。整个系统内部有三个能量母线电母线、热母线、冷母线。电母线上有市电购电、CHP发电、光伏发电、电池充放电对应的消费方是电负荷、电制冷机、电池充电热母线上有CHP余热、锅炉产热、蓄热罐放热消费方是热负荷和吸收式制冷机冷母线上只有电制冷机和吸收式制冷机两个产冷源加上蓄冷设备共同满足冷负荷。这套拓扑最容易让人绕晕的是热和冷的耦合关系吸收式制冷机不是独立的制冷设备它是挂在热母线上的一个热负荷一边消耗热一边产冷。所以热平衡方程里必须把供吸收式制冷机的热单列出来不能和热负荷混在一起。我见过好几个初学者把热负荷写成了纯采暖负荷结果夏季工况下热母线全部不平衡模型跑出来一堆奇怪结果。2.2 目标函数运行成本加上舒适度软约束优化调度的目标是让系统在满足负荷的前提下总成本最小同时尽量让用户待在舒适区。我采用的日运行目标函数是min C Σt [ C_elec(t)×P_buy(t) C_gas×(F_chp(t) F_boiler(t)) Σi C_om_i×P_i(t) ] C_pmv前半部分是常规经济成本购电费、燃气费、设备运维费。P_buy(t)是t时段从电网购电功率C_elec(t)是分时电价F_chp(t)和F_boiler(t)分别是CHP和锅炉的天然气耗量单位m³/hC_om_i是第i台设备的单位运维成本按发电/供能量计费。后半部分C_pmv是舒适度惩罚项我用的形式是C_pmv Σt α × ( a(t) b(t) )其中a(t)、b(t)是辅助变量满足PMV(t) - PMV_max ≤ a(t) PMV_min - PMV(t) ≤ b(t) a(t) ≥ 0b(t) ≥ 0这是一个典型的线性软约束写法。PMV超过上限PMV_max时a(t)被迫大于0惩罚项增加低于下限时b(t)大于0。α是舒适度惩罚权重α越大系统越愿意牺牲经济性来维持舒适。这里有个值得说透的设计决策为什么不直接写硬约束PMV_min ≤ PMV(t) ≤ PMV_max因为舒适区其实是个可以偶尔突破的区域用户并不是每时每刻都要求PMV严格落在±0.5以内——中午最热时室温冲到27.8℃半小时多数人可以接受但持续一个下午就不行了。如果写硬约束相当于要求全天96个时段全部严格达标可行域会被显著压缩有时甚至无解同时硬约束会把一些稍越界但成本极低的调度方案直接砍掉经济性损失很大。软约束让优化器自己权衡多花多少钱换多少舒适度这才是工程上真实的选择逻辑。2.3 约束条件里容易写错的几个关系式调度约束比较多我列几个最容易写错、也最影响求解质量的第一CHP机组的热电产出与燃气消耗关系。天然气热值H_gas取9.7 kWh/m³那么P_chp(t) η_e × H_gas × F_chp(t) H_chp(t) η_h × H_gas × F_chp(t)注意η_e和η_h不是独立的两者之和约等于CHP总效率。我代码里取η_e0.35、η_h0.45总效率80%已经是中小型燃气轮机比较乐观的水平。这里常见错误是直接把F_chp当成本单位写进购气费忘了乘H_gas结果成本低估好几倍。第二热母线的平衡必须带弃热变量。很多文献写H_chp(t) H_boiler(t) H_load(t)这在冬季勉强成立但夏季热负荷很低时CHP为了在峰电时段多发电必然产生大量副产品热而热负荷加吸收式制冷机的热需求都消耗不完多出来的热怎么办必须允许弃热。所以我写成H_chp(t) H_boiler(t) H_store_out(t) H_heat(t) H_abs(t) H_dump(t)H_dump(t)就是弃热功率是连续变量不设上限也行但为了数值稳定通常给一个很大上限。第三室温动态方程这是舒适度模型和系统模型之间的桥梁。建筑等效热容模型T_room(t1) T_room(t) (Δt / C_bld) × [ Q_hvac(t) - UA_wall × (T_room(t) - T_out(t)) - Q_solar_in(t) ]Δt是调度步长单位小时C_bld是建筑等效热容单位kWh/K可以理解成房间的热惯性——墙体越厚、面积越大C_bld越大温度变化越慢UA_wall是围护结构等效传热系数kW/KQ_solar_in是太阳得热夏季白天会让室温上升。Q_hvac(t)是暖通系统实际注入室内的热功率夏季制冷时它是负值从室内抽热。这个方程单独看没问题但和系统模型一耦合就容易出错供冷功率Q_ec(t)Q_abs(t)是冷负荷侧值而Q_hvac(t)是室内侧值两者之间还要考虑风机盘管送风效率、水管冷损代码里我会加一个0.9左右的折减系数。很多人直接画等号导致室温方程和能量平衡方程互相矛盾求解器报不可行。3. MATLAB代码实现从PMV函数到YALMIP求解框架3.1 代码目录与输入数据结构一套能反复改参数的调度代码目录结构我建议这样组织IES_PMV/ ├── main_dispatch.m # 主程序 ├── load_data.m # 读取负荷/电价/天气数据 ├── pmv_calc.m # 原版Fanger方程PMV计算 ├── pmv_fit_linear.m # PMV单变量拟合与分段线性化 ├── build_model.m # 建立优化变量、目标、约束 ├── solve_milp.m # 求解与结果诊断 ├── plot_results.m # 绘图 └── data/ ├── load_24h.xlsx # 电、热、冷负荷 ├── price_24h.xlsx # 分时电价 └── weather_24h.xlsx # 室外温度、太阳得热load_data.m里有一个关键习惯所有物理量在进入模型前统一换算成调度模型标准单位。我用的标准是功率kW、能量kWh、天然气流量m³/h、价格元/kWh和元/m³、温度℃。很多坑都是单位混乱导致的比如Excel里负荷是MW代码里默认powers是kW中间差1000倍电池容量算出来上百兆瓦时结果还恰好可行这种问题最难查。3.2 PMV计算函数精度与单位陷阱pmv_calc.m的原版Fanger方程实现核心是tcl的迭代。我贴一段能用起来的版本function pmv pmv_calc(ta, tr, rh, vel, met, clo) % ta: 空气温度 ℃ % tr: 平均辐射温度 ℃ % rh: 相对湿度 %不是0-1 % vel: 空气流速 m/s % met: 代谢率 met % clo: 服装热阻 clo M met * 58.2; % W/m2 W 0; Icl clo * 0.155; % m2*K/W pa (rh / 100) * exp(16.6536 - 4030.183 ./ (ta 235)); % kPa tcl ta 2; % 服装表面温度初值 for k 1:30 hc max(2.38 * abs(tcl - ta)^0.25, 12.1 * sqrt(vel)); if Icl 0.078 fcl 1 1.29 * Icl; else fcl 1.05 0.645 * Icl; end tcl_new 35.7 - 0.028 * (M - W) ... - Icl * (3.96e-8 * fcl * ((tcl 273)^4 - (tr 273)^4) ... fcl * hc * (tcl - ta)); if abs(tcl_new - tcl) 1e-4 tcl tcl_new; break; end tcl tcl_new; end 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)^4 - (tr 273)^4) ... - fcl * hc * (tcl - ta) ); end这个函数本身不难但有三处必须注意。第一相对湿度rh传入的是百分数0-100不是0-1pa公式里除以100第二tcl迭代初值取ta2℃一般20步内收敛若初值离真实值太远比如冬季给ta18℃却初值设成35℃中间可能出现负浮点导致NaN需要做保护第三所有温度在辐射项里要加273换算成开尔文但对流换热项里保持摄氏度混用是常见错误。写好这个函数后第一件事不是接优化模型而是做单元测试。Fanger原始文献里有一个经典算例ta25℃tr25℃rh50%vel0.1m/smet1.2clo0.5PMV应该接近0.1以内。肉眼看着差不多就放进模型后面报错了你都不知道是舒适度模块还是求解器的问题。3.3 分段线性化让PMV跑进MILP求解器原版PMV函数即使降维成单变量也还是带指数和四次方的强非线性不能直接进MILP。我的做法是用pmv_calc.m在温度网格上预计算一系列节点值然后用SOS2特殊有序集第二类约束做分段线性插值。温度区间我取16℃到34℃每1℃一个节点共19个节点。思路是引入连续权重λ_k和二进制变量z_k约束保证ta只能在相邻两个节点之间插值不允许跳过区间。YALMIP里可以这样写x_node 16:1:34; % 温度节点 y_node arrayfun((xx) pmv_calc(xx, xx, rh, vel, met, clo), x_node); lambda sdpvar(length(x_node), 1); z binvar(length(x_node) - 1, 1); F [F, sum(lambda) 1, sum(z) 1]; F [F, lambda(1) z(1)]; for k 2:length(x_node) - 1 F [F, lambda(k) z(k - 1) z(k)]; end F [F, lambda(end) z(end)]; F [F, T_room x_node * lambda]; % 用插值表示室温 F [F, PMV_lin y_node * lambda]; % 用插值表示对应的PMV这里面lambda的第k个分量被约束成只能出现在相邻两个节点上所以PMV_lin是ta在相邻节点之间的线性插值结果。注意因为T_room被x_node*lambda锁定了它不可能超出16℃到34℃的范围这其实是个天然的温度上下限约束。如果你希望允许极端室温就得扩大节点范围我的建议是从一开始就把它限制在站得住的物理范围内反而有利于求解稳定性。有同学问我为什么不直接用YALMIP的pwf函数构建分段函数。pwf确实更简洁但不同YALMIP版本对pwf的兼容性和约束处理有差异跑大模型时偶尔会生成冗余的二进制变量我自己的经验是手写SOS2更可控也方便排查问题。3.4 主程序骨架YALMIP建模与求解器选择主程序里把设备变量、能量平衡、室温方程、舒适度软约束全部组合起来最后调用求解器。核心骨架如下T 24; ops sdpsettings(solver, gurobi, verbose, 2); P_buy sdpvar(T, 1); F_chp sdpvar(T, 1); F_boiler sdpvar(T, 1); P_chp sdpvar(T, 1); H_chp sdpvar(T, 1); H_boiler sdpvar(T, 1); P_ec sdpvar(T, 1); Q_ec sdpvar(T, 1); H_abs sdpvar(T, 1); Q_abs sdpvar(T, 1); H_dump sdpvar(T, 1); P_bat_ch sdpvar(T, 1); P_bat_dis sdpvar(T, 1); SOC sdpvar(T 1, 1); T_room sdpvar(T 1, 1); PMV_lin sdpvar(T, 1); a_pmv sdpvar(T, 1); b_pmv sdpvar(T, 1); F []; % 电平衡 F [F, P_buy P_chp P_pv P_bat_dis P_elec P_ec P_bat_ch]; % 热平衡 F [F, H_chp H_boiler H_store_out H_heat H_abs H_dump]; % 冷平衡 F [F, Q_ec Q_abs Q_store_out Q_cool]; % CHP模型 F [F, P_chp 0.35 * 9.7 * F_chp]; F [F, H_chp 0.45 * 9.7 * F_chp]; % 制冷机模型 F [F, Q_ec 4.0 * P_ec]; F [F, Q_abs 1.2 * H_abs]; % 室温动态 for t 1:T F [F, T_room(t 1) T_room(t) (dt / C_bld) * ... (0.9 * (Q_ec(t) Q_abs(t)) - UA * (T_room(t) - T_out(t)) - Q_solar(t))]; end % SOS2分段线性化 F [F, sos2_pmv(T_room(1:T), PMV_lin)]; % 舒适度软约束 F [F, PMV_lin - 0.5 a_pmv, -0.5 - PMV_lin b_pmv]; F [F, a_pmv 0, b_pmv 0]; % 电池SOC SOC(1) 50; F [F, SOC(t 1) SOC(t) - P_bat_dis / 0.95 0.95 * P_bat_ch, ... SOC(end) 50, 0 SOC 100]; Objective sum(price_elec .* P_buy) ... sum(gas_price .* (F_chp F_boiler)) ... sum(om_cost ... ) ... alpha * sum(a_pmv b_pmv); sol optimize(F, Objective, ops);求解器选择上我强烈建议用Gurobi或CPLEX尤其当T从24扩展到96甚至8760时两者差距会拉得很大。如果实验室没有商业求解器licenseMATLAB自带的intlinprog也能跑但YALMIP默认求解器在MILP问题上表现一般有时候会卡在预求解阶段半小时不动别怪代码是求解器不行。4. 舒适性与经济性的权衡仿真结果与权重系数标定4.1 典型日场景与基准工况为了说清楚α怎么起作用我构造了一个夏季典型日场景。园区最大电负荷350kW最大冷负荷380kW热负荷只有生活热水60到120kW。分时电价谷时0.38元/kWh平时0.70元/kWh峰时1.10元/kWh天然气价3.5元/m³。CHP容量300kW发电效率0.35、热效率0.45燃气锅炉400kW电制冷机400kWCOP4.0吸收式制冷机300kWCOP1.2蓄电池200kWh/100kW。先跑一个基准场景α0也就是完全不管用户舒适度纯经济调度。结果非常有代表性室温曲线在中午到傍晚爬到28℃以上最高接近28.7℃对应PMV约0.84已经明显偏热全天统计下来PMV超限时段占比62%也就是说一天里有一大半时间用户处在不舒适状态。日运行成本4380元。4.2 不同权重α下的调度结果对比保持其他参数不变把α从0逐渐增大重复求解得到一组结果α日运行成本元成本增幅PMV超限时段占比室温范围℃043800%62%25.8 - 28.70.0544682.0%35%25.1 - 27.90.2045844.7%11%24.6 - 26.90.5047458.3%3%24.3 - 26.41.00489711.8%0%24.0 - 26.1这个表很值得细看。α从0增加到0.2成本只涨了4.7%PMV超限占比却从62%降到了11%这就是典型的低成本换高舒适区间性价比极高。α再往上走边际效果快速衰减——从0.2加到1.0成本额外上涨7.1个百分点换来的只是把超限占比从11%降到0%。对多数工程场景来说这部分投入不划算。从调度策略上拆解原因更直观。α0时优化器会在谷时段给电池充满电、尽可能让制冷设备在谷时蓄冷峰时段则卡着CHP最大出力来赚高电价收益至于室温是否超限它根本不在乎。α0.2时下午最热时段电制冷机和吸收式制冷机必须同时满发蓄冷罐提前放冷夜间还会多用谷电给蓄冷罐补冷量冷侧出力节奏明显变得更保守但整体效率依然不错。α1.0时系统几乎全程把室温锁在26℃以内为了维持高舒适度CHP在晚间低负荷时段也保持最小技术出力来保证余热驱动吸收式制冷发电效率下降、燃料成本上升这就是最后那11%超限时段被买走的代价。4.3 帕累托前沿解读与α建议值把α继续扫到更大值把每组结果画成散点横轴是日运行成本纵轴是PMV超限占比就是一条单调递减的帕累托前沿。工程上挑点不需要追求极端我看这条曲线的经验是找到成本增幅5%左右、超限占比降到10%上下的位置那个α就是合理值。按这套场景α取0.15到0.3都说得过去。更讲究一点的做法是把舒适度惩罚和实际运营挂钩。比如园区物业和租户约定的温度投诉赔偿是每超限时段100元那你的α就不该拍脑袋定而是让α≈单位PMV超限对应的等效经济损失这样舒适度惩罚项在经济侧有了真实含义再做多目标分析时说服力也更强。5. 调试实录非线性PMV带来的收敛问题与工程化处理5.1 初始解与热启动别让求解器从空荡荡的可行域出发纯MILP模型对初始解不敏感但我们的模型里有一段室温动态方程和SOS2约束初始点不合适很容易触发预求解阶段的数值问题。我的第一个坑就是T_room(1)初值随手填了28℃但前一天的模拟数据来自25℃的工况首个时段瞬间产生巨大的温度偏差迭代几十轮后目标函数出现振荡。处理方法有两个一是把T_room(1)也作为优化变量在约束里限定它偏离合理初值不超过0.5℃二是更推荐的做法——先用α0的纯经济调度跑一遍把得到的室温曲线作为α0场景的初始点。YALMIP里可以用assign给变量赋初值然后调用optimize时开启热启动MILP的迭代次数能明显减少。我实测过α0.2场景热启动比冷启动平均快30%左右别小看这个差距等扩展到8760小时全年优化时这差距就是几十秒和几分钟的区别。5.2 室温动态方程的数值稳定性室温方程用的是显式欧拉格式数值稳定性条件要求Δt/τ不能太大其中τC_bld/UA_wall是建筑围护结构的时间常数。我们算一下C_bld20 kWh/K、UA_wall1.5 kW/K时τ约13.3小时Δt取1小时Δt/τ约0.075很稳定但如果C_bld取5 kWh/K、UA_wall取1.5 kW/Kτ只有3.3小时Δt/τ约0.3已经开始出现温度锯齿再极端点C_bld取2 kWh/K、Δt取1小时Δt/τ约0.75室温曲线会明显振荡优化器甚至能利用这种伪振荡免费压低峰值温度结果毫无物理意义。我的工程经验是保持Δt/τ≤0.1也就是时间常数至少要达到调度步长的10倍。如果建筑确实轻薄、时间常数小就别强行用大步长把调度步长缩到15分钟或者给室温变化率加一个显式约束比如|T_room(t1)-T_room(t)|≤0.5℃这样既保证了数值稳定也模拟了真实建筑不可能瞬变的基本物理常识。5.3 罚函数处不光滑与SOS2边界两个隐蔽的坑舒适度软约束里用了max函数虽然可以通过辅助变量线性化但很多人在写目标函数时还是习惯直接用max(0, PMV-PMV_max)结果引入一个不连续项fmincon这类非线性求解器根本算不出梯度卡在原地报NaN。解决办法就是我前面写的辅助变量法把max拆成两个不等式约束。这是标准技巧但在YALMIP里有一个连带坑如果PMV_lin本身来自SOS2插值同时又出现在max函数的约束里求解器预求解阶段可能会做大量不可行松弛导致变量数量虚高。我的处理是把PMV_lin明确声明为sdpvar并显式给区间边界比如-3≤PMV_lin≤3让它和温度插值解耦能省不少事。另一个隐蔽的问题是SOS2的温度节点范围。如果T_room被x_node*lambda锁定那它的值永远落在16℃到34℃之间这本身没问题但如果你后来又加了别的约束让T_room超过这个范围预求解器会直接报不可行而且报错信息晦涩得让人一头雾水。所以我在build_model.m里都会显式写一行16T_room34并在注释里标明这是PMV线性化的定义域约束不是室内温度的真实上下限避免三个月后自己回来看代码忘记当初为什么写它。5.4 求解器选型与不可行问题排查清单最后的求解器选型我用一个经历说明问题。同样一套模型T24时用intlinprog跑大概15秒用Gurobi大概3秒T96时intlinprog可能要十分钟Gurobi还是20秒内结束T8760时intlinprog基本无法当天出结果Gurobi也要几分钟到十几分钟不等。所以如果你的模型时间维度要做全年扩展投资一个商业求解器是值得的。如果暂时没有license至少可以用cbcCoin-OR的开源求解器作为YALMIP的备选表达效果比intlinprog稳定。真遇到不可行按这个顺序排查第一步只跑α0的纯经济调度如果这都不可行问题在设备平衡或SOC初值第二步加上舒适度软约束但把α设成0此时理论上结果应该和第一步相同若有差异说明约束写重了第三步加α0若不可行99%是温度定义域、T_room初值或SOS2节点范围问题第四步把所有变量和约束打印出来人工检查重点看SOC的最后一个时段约束是否漏写蓄热罐是否有初末能量平衡。这套排查流程我用了很多次基本能在五分钟内定位问题。很多时候所谓模型不可行根本不是模型错而是某个辅助变量忘了初始化或者数据单位差了1000倍又或者室温方程里Q_hvac的正负号搞反了——夏天制冷时它应该是负值有人写成正值等于一边给房间加热一边还要求室温不超过26℃不冲突才怪。最后分享一个我个人的调试习惯PMV模块永远单独做单元测试先拿25℃、RH50%、vel0.1m/s、met1.2、clo0.5的标准算例校一遍确认函数输出落在0.1附近再把模块接进优化模型。另外α、温度节点范围、舒适度阈值、设备容量全部写成配置文件或load_data.m里的结构体参数不要散落在主程序里。我做权重灵敏度分析时只需要改一个数字然后批量循环求解不用每次手动改代码省下的时间足够把帕累托前沿整个扫一遍。