简介面向综合能源优化调度与需求响应研究场景提供一套基于Matlab的MPC蓄热控制仿真程序代码结构清晰适合作为算法复现与教学示例。程序以纽约市一栋含被动与主动蓄热系统的办公楼为对象结合三级需求费与日前调度电价构建锥形需求响应程序并纳入人员热舒适约束以随机最优控制语言建模、模型预测控制近似求解可帮助电力系统、建筑节能方向的初学者与研究者在真实场景中理解蓄热型需求响应的建模流程。压缩包共11个文件以9个.m主程序与函数文件为核心涵盖参数配置、优化求解、结果对比等环节另有1个.mat数据文件存储仿真场景数据以及1篇配套PDF参考文档说明控制原理整体仅1.59MB轻量易部署。已有1727人学习。资源提供完整可运行源码、示例数据与结果绘图脚本可复现日前调度、三级需求费下的蓄热需求响应与MPC控制效果帮助系统理解热舒适约束建模、锥形需求响应程序构建及CVX工具箱调用流程。1. 热储能需求响应的MPC仿真这套Matlab程序值得先跑通再看论文做综合能源优化调度的人第一次拿到这套热储能需求响应的Matlab仿真程序通常会被一堆名字相近的TSDR脚本劝退。反直觉的结论是代码结构比论文更直白——把main.m跑通、把CVX工具箱路径配好所有结果图都会自动生成真正要读懂的是TSDRMPC.m里的那几十行约束。这套程序复现的是某国外城市一栋带有被动与主动蓄热的办公楼费率结构为三级需求费加日前调度逐时电价目标是让蓄热系统在满足人员热舒适的前提下把整栋楼的总电费压到最低。适合正在做MPC、需求响应或建筑能源管理的从业者也适合刚入门想找一个完整闭环仿真做参照的研究者。2. 先看懂模型锥规划、热储能状态方程和MPC的边界条件2.1 MPC在这里优化什么状态、控制量、扰动与约束这套代码里的MPC不是一个通用控制包而是针对“建筑热储能”这个具体对象写的优化控制器。先看它数学上在干什么否则改参数就是盲改。状态选取了两个室内温度Tz和蓄热罐温度Ts。室内温度代表被动蓄热因为墙体、家具本身有热容电价低谷时把室温拉高或拉低相当于把能量存在建筑结构里蓄热罐温度代表主动蓄热它才是真正可调度的储能设备。两个状态串在一起构成一个二阶线性状态方程X(k1) A * X(k) B * u(k) W * d(k)其中X是2维列向量[室温; 罐温]u是蓄热罐的充放电功率充为正、放为负d是外部扰动序列——包括室外温度、内部设备散热、人员产热这些不可控输入。MPC的任务就是在每个采样时刻解一个有限时域优化问题算出未来N步的最优充放电策略但只执行第一步然后滚动更新。目标函数分两块。第一块是能源费用即逐时电价乘以从电网购电的功率第二块是需求费用这是国内很多做电力市场的人不熟悉的地方。需求费不是按电量算的而是按“结算周期内最大功率”算的——哪怕你一个月只出现一次峰值这一个点就能决定整月的需求费。代码里把需求费的数学形式写成了锥形约束本质是用一个辅助变量D去夹住“月最大功率”这个不可导的max函数从而让整个问题变成标准的锥规划CVX才能直接求解。约束条件里有热舒适区间夏季制冷场景下室温一般限制在20~24度之间这给了MPC利用建筑热容做预冷的自由度。蓄热罐的约束包括容量上下限和充放电速率上限。注意这些约束在代码里是硬约束一旦参数设得太紧或者预测的扰动数据不合理CVX会直接报不可行。2.2 十个文件各管哪一段先看结构再动手这套代码的文件命名其实很清晰拿到手先把每个文件在做什么摸清楚避免一上来就乱改。下面的表格按“我从哪个文件入手”排序文件作用我一般怎么用main.m总入口串起参数、求解、绘图先跑一遍看有没有报错fixedTSDRparams.m物理固定参数建筑热容、蓄热罐容量、舒适区间、预测时域改模型时动这里dailyTSDRparams.m每日变化参数电价序列、三级需求费率、天气数据复现不同日期时改这里TSDRdata.mat数据文件历史电价、负荷、室外温度序列load进来供MPC滚动使用TSDRMPC.mMPC核心求解函数CVX建模在这里读懂它就等于读懂整个包TSDROLOC.m离线最优控制即事后全知视角的优化用于和MPC求出的结果做下界对比TSDRgD.m需求费处理的辅助函数一般不用改TSDRplots.m主结果绘图温度、功率、费用曲线跑完直接看输出TSDRcomparisonPlots.mMPC与离线最优的对比图汇报和写论文时用TSDRparamPlots.m参数敏感性扫描绘图做调参分析时用注意TSDRdata.mat和dailyTSDRparams.m里很可能包含同一套价格数据运行时如果两个都加载后加载的会覆盖先加载的变量。这个问题在第5章会展开说但你现在先有个印象这套代码的变量管理比较依赖工作区顺序。2.3 fixedTSDRparams.m里藏着物理模型的核心参数打开fixedTSDRparams.m你会看到一串赋值语句。不同版本的实现略有差异但核心参数逃不出这些% 建筑与蓄热系统物理参数 params.building_capacitance 1.2e7; % 建筑热容单位 J/K越大蓄热能力越强 params.building_r 0.02; % 建筑围护热阻单位 K/W越小保温越好 params.tank_capacity 3000; % 蓄热罐等效水容量单位 kg params.tank_init_temp 15; % 蓄热罐初始温度单位 摄氏度 params.tank_min_temp 7; % 蓄热罐最低允许温度 params.tank_max_temp 30; % 蓄热罐最高允许温度 params.Tz_min 20; % 室内温度下限舒适约束 params.Tz_max 24; % 室内温度上限舒适约束 params.predict_horizon 24; % 预测时域 N对应日前电价的完整周期 params.control_horizon 12; % 控制时域 M一般小于等于 N params.charge_rate_max 500; % 蓄热罐最大充热功率单位 kW params.discharge_rate_max 800; % 蓄热罐最大放热功率单位 kW建筑热容这个参数最容易被新手误解。它不是一个可以直接查表的物性值而是“整栋楼等效热容”——把所有墙、地板、家具的比热容乘以质量再累加单位是J/K。这个值设得越大被动蓄热的潜力就越大MPC就越倾向于在电价低时把室温拉偏、电价高时让室温自然漂移。反之如果你把这个值缩小十倍建筑的“储能”效果基本消失MPC只能依赖蓄热罐需求响应的降费空间大幅缩水。蓄热罐的温度上下限直接影响状态可行域。这里有一个常见的工程权衡罐温下限不能设太低否则热泵或锅炉的制热效率会恶化上限不能设太高否则热量损耗变大同时罐体保温成本上升。代码里默认的7到30度算是一个冷水蓄冷的典型区间如果你要改成蓄热场景记得同步修改状态方程的系数矩阵A和B而不是只改上下限——这一点README里一般不会强调。预测时域N24是这套代码里最值得琢磨的一个参数。它和日前调度电价的发布周期严格对应电网提前一天发布24小时逐时电价MPC在凌晨时能看到完整一天的价格曲线所以N24能保证优化问题覆盖“低谷充电、高峰放电”这个完整套利周期。如果你把N改成12MPC在早晨6点开始预测时只能看到白天的电价看不到凌晨的谷价蓄热罐就不会充电需求响应直接就失效了。2.4 TSDRdata.mat和TSDRplots.m配合出图的顺序TSDRdata.mat不是直接run的脚本而是一个数据文件。在main.m里通常这样加载% 加载历史数据电价序列、室外温度、内部热增益 load(TSDRdata.mat); % 变量名常见为 price、temperature、internal_gain具体以whos查看为准 whos加载之后建议立刻执行一次whos确认当前工作区里的变量名。我遇到过两次因为变量名不匹配导致的“图是空白的”问题——脚本里引用的是price_data但数据文件里变量名叫priceMatlab不会自动做模糊匹配结果是绘图时把空数组传进去了。TSDRplots.m的工作逻辑很简单把TSDRMPC求出的室温、罐温和充放电功率画在这套数据的时间轴上同时叠加电价曲线方便对照。如果你发现图出来了但曲线和论文顺序不一致先检查TSDRcomparisonPlots.m里画的是不是同一组变量。3. 把MPC跑起来CVX环境配置与main.m主流程拆解3.1 前置环境CVX不是装了就完事这套代码的运行前提是CVX工具箱而且默认依赖的是CVX自带的求解器。很多人在第一步就栽了——在Matlab命令窗口直接敲main.m结果报Undefined function or variable cvx_begin。这个报错八成是因为CVX的文件夹还没加入Matlab路径。CVX的安装不算复杂但有几个顺序细节需要注意。解压CVX压缩包后不要直接双击cvx_setup.m而是先在Matlab里切到CVX所在目录再运行cvx_setup。它会自动把子目录加入搜索路径并保存路径配置。安装完成之后在脚本开头加一段环境自检是值得养成的习惯%% 环境自检确认CVX是否可用 if ~exist(cvx_begin, file) error(CVX工具箱未加入MATLAB路径请先运行cvx_setup); end % 小规模优化推荐用SeDuMi大规模问题才考虑SDPT3 cvx_solver sedumiexist(cvx_begin, file)能检测到cvx_begin这个函数是否存在如果CVX路径没配好这里会直接报错而不是等到优化求解时才报一个难懂的异常。cvx_solver这一行不是必须的CVX默认会用SDPT3但对于这套代码的变量规模——2个状态、1~2个控制量、24步预测——SeDuMi通常快不少而且数值稳定性更好。如果机器内存紧张这行能显著降低求解时间。3.2 main.m主流程先跑参数再跑求解最后出图main.m的结构在标准版本里很清晰基本是三步加载参数、调用求解器、绘图。下面是一个简化但保留了关键调用顺序的版本clear; clc; close all; %% 第一步加载物理参数和每日参数 run(fixedTSDRparams.m); % 物理模型参数一般只改一次 run(dailyTSDRparams.m); % 电价、需求费、负荷序列按天更新 %% 第二步加载历史数据注意与dailyTSDRparams的变量覆盖关系 load(TSDRdata.mat); %% 第三步调用MPC求解 x0 [params.Tz_init; params.tank_init_temp]; % 初始状态室温罐温 [U_mpc, X_mpc, J_mpc] TSDRMPC(x0, price, params); %% 第四步调用离线最优控制做对比下界 [U_oloc, X_oloc, J_oloc] TSDROLOC(x0, price, params); %% 第五步出对比图 TSDRplots(X_mpc, U_mpc, X_oloc, U_oloc, price); TSDRcomparisonPlots(X_mpc, X_oloc, U_mpc, U_oloc, price); close all;这里有一个值得强调的顺序细节run(dailyTSDRparams.m)和load(TSDRdata.mat)之间极易发生变量覆盖。如果dailyTSDRparams.m里定义了price向量而TSDRdata.mat里也定义了同名变量谁在后面谁生效。代码的逻辑是“越靠后的数据越接近当前调度场景”但在实际复现时往往因为变量名不一致导致MPC拿到错误的价格序列而不自知。我会在加载完数据后立刻打印一句话确认价格长度和首尾值比如fprintf(price length: %d, first: %.2f, last: %.2f\n, ... length(price), price(1), price(end));3.3 TSDRMPC.m核心CVX建模的完整逻辑TSDRMPC.m是整套代码的心脏函数签名大概是这样的function [U, X, J] TSDRMPC(x0, price, params)它接收初始状态x0、未来N步电价price序列、结构体params返回控制序列U、状态轨迹X和最优目标值J。核心的CVX建模代码大致可以抽象为下面这个结构%% 提取参数 N params.predict_horizon; A params.A; % 2x2状态矩阵 B params.B; % 2x1控制矩阵 W params.W; % 2x1扰动矩阵 %% 定义CVX优化问题 cvx_begin quiet variable X(2, N1) % 状态轨迹第1行为室温第2行为罐温 variable Uch(1, N) % 蓄热罐充电功率 variable Udis(1, N) % 蓄热罐放电功率 variable D(1, N) % 需求辅助变量用于锥形需求约束 % 目标函数能源费用 需求费用 minimize( price(:) * (Udis(:) - Uch(:)) * params.dt ... params.demand_rate * max(D(:)) ) subject to % 状态方程约束 for t 1:N X(:, t1) A * X(:, t) B * (Udis(t) - Uch(t)) W * params.weather(t); % 热舒适约束 X(1, t) params.Tz_min; X(1, t) params.Tz_max; % 蓄热罐约束 X(2, t) params.tank_min_temp; X(2, t) params.tank_max_temp; % 充放功率约束 Uch(t) 0; Uch(t) params.charge_rate_max; Udis(t) 0; Udis(t) params.discharge_rate_max; % 锥形需求约束D要夹住净负荷 D(t) params.base_load(t) Uch(t) - Udis(t); D(t) 0; end % 初始状态约束 X(:, 1) x0; cvx_end这段代码里的逻辑说明很重要。max(D(:))在CVX里是不允许直接写在目标函数里的——严格来说CVX支持minimize中的max表达式但会给求解器增加不必要的负担。更规范的做法是把需求费用也用锥约束表达即在约束里定义D(t)大于等于净负荷同时让目标函数对D做线性求和。这样做的数学效果是一样的在最小化目标的推动下D会恰好等于净负荷的峰值而不会比它大。这就是“锥形需求响应程序”这个名字的来源——它不是一个玄学概念而是把max函数用一组线性不等式等价替换的建模技巧。求解结束之后检查cvx_status。如果输出是Solved说明问题收敛如果是Infeasible说明约束自相矛盾最常见的是热舒适区间太紧或者充放电功率上限和建筑热负荷不匹配。我一般在cvx_end之后加一行断言if ~strcmp(cvx_status, Solved) error(MPC求解失败%s请检查约束参数, cvx_status); end这个习惯值十次debug的时间。cvx_status除了Solved之外还可能出现Failed、Inaccurate/Solved等状态后者意味着“求出了解但精度不足”对蓄热调度这种不需要高精度的问题通常可以接受但如果你在做论文对比实验建议调高CVX的精度设置或者改用更稳定的求解器。4. 三级需求费、热舒适边界与数据文件三个关键参数怎么调4.1 三级需求费的费率结构和滚动折算很多第一次接触美国电力市场计费的人会困惑为什么电费要分“能源费”和“需求费”两块而且需求费的计算方式如此反直觉。这套代码针对的是某国外城市的商业电价结构能源费就是普通的分时电价按kWh结算需求费按结算周期内出现的最大30分钟平均功率乘以费率按kW结算。商业建筑的月账单里需求费经常占到总电费的30%以上所以MPC把需求费写进目标函数是非常有实际意义的。代码里的三级需求费通常对应三个时段高峰时段比如工作日下午、平时段、低谷时段。每个时段有不同的费率而且各自独立计算该时段内的最大功率。这样做会导致一个问题MPC可能会把高峰时段的部分负荷“转移”到平时段因为平时段需求费低但平时段的峰值也会产生费用。更麻烦的是MPC的预测时域只有24小时而需求费实际是按月结算的——在一个24小时的窗口里你根本不知道这个月已经出现的峰值是多少所以无法精确计算“再增加1kW要多交多少需求费”。代码里的常见做法是把“月度结算”折算成“窗口内滑动峰值”在每次滚动优化时把历史已发生的最大功率作为一个额外约束目标函数里只优化当前窗口内的增量需求。这个近似在工程上完全够用但你要明白它的边界预测时域外的用电行为会影响真实的需求费但MPC看不到。如果你想把代码改成更贴近真实账单的版本需要在状态里增加一个“本月迄今最大功率”的记忆变量并把它作为扰动传入优化。4.2 热舒适区间越大越容易算但需求响应就没意义热舒适约束是这套代码里影响求解难度的最大变量。Tz_min和Tz_max之间的带宽决定了两件事一是MPC有多少自由度去“预冷”或“预热”建筑二是热舒适区间越宽状态方程的可行域越大优化越不容易无解。默认的20到24度是一个典型的制冷季舒适区间。在这个带宽下MPC可以在凌晨电价低时把室温从24度预冷到20度把“冷量”储存在建筑结构里白天气温上升时让室温在允许范围内自然漂移减少电制冷机在高电价时段的出力。这相当于把电价波动转化为室温的小幅波动。如果你把区间改成22正负0.5度舒适度是好了但被动蓄热的调度空间几乎消失——MPC只能在21.5到22.5度之间折腾这点温差能储存的能量有限需求响应的削峰效果会大打折扣。反过来说如果你为了追求降费把区间放宽到18到26度MPC会倾向于在电价低谷时过度预冷导致实际体感温度过低人员投诉率上升。我的调参建议是先保持默认区间跑一遍记录总电费和需求费然后每次把Tz_min降低0.5度重新跑观察需求费的边际变化。当你发现继续放宽舒适区间不再显著降低需求费时就说明被动蓄热的潜力已经用尽再放宽只会牺牲舒适度换不来钱。4.3 从TSDRdata.mat看数据准备价格、温度和负荷的时间对齐TSDRdata.mat里存放的是历史时间序列数据结构一般是结构体或者多个列向量。写一个简单的检查脚本把所有变量的长度看一眼%% 检查数据对齐 load(TSDRdata.mat); len_price length(price); len_temp length(temperature); len_load length(base_load); fprintf(price%d, temp%d, load%d\n, len_price, len_temp, len_load); if len_price ~ len_temp || len_price ~ len_load warning(数据长度不一致请检查时间轴对齐); end这里最容易踩的坑是时间步长。如果TSDRdata.mat里的数据是15分钟一个点而dailyTSDRparams.m里设置的采样时间是1小时MPC拿到的状态方程和时间序列就错位了。A矩阵和B矩阵是基于1小时采样离散化得到的喂15分钟的数据进去物理含义完全变掉。解决办法是统一所有时间序列的采样间隔确保与fixedTSDRparams.m里的离散化时间步长一致。这套代码默认是1小时一个点对应日前调度96点中的24点你无需改采样但要确认TSDRdata.mat里的序列确实是24点一天而不是48点。5. 避坑与排查CVX报错、空图、算不动的五个实战问题5.1 报Undefined function or variable cvx_begin现象运行main.mMatlab直接报错指向cvx_begin这一行。原因CVX没有安装或者安装了但没把路径加入Matlab。另一个常见原因是CVX版本太旧不支持当前Matlab版本——比如CVX 2.1在R2022b之后的某些版本上会提示证书错误或直接无法启动。解决先确认安装目录在Matlab命令窗口执行cvx_setup如果报证书问题检查CVX的许可文件。更快的验证方式是直接敲cvx_version如果返回版本号说明安装成功。我习惯在main.m顶部加一段自检代码第3.1节里的exist判断让新人拿到代码能直接定位问题而不是去翻一堆readme。5.2 求解器返回Infeasible最优值为NaN现象cvx_status输出InfeasibleU和X全是NaN绘图报错。原因约束之间矛盾了。最常见的组合是Tz_max设成了22Tz_min也设成了22而室外温度在预测时域内持续高于26蓄热罐容量又不足以把室温压到22度以下于是无解。说白了是“设备能力不够还硬要舒适”。解决把热舒适区间放宽或者增大蓄热罐的容量和充放功率上限。更稳妥的做法是给舒适约束加一个松弛变量——让室温可以轻微越界代价是目标函数里加一个大的惩罚项。这样优化问题永远有解只是会在特殊天气下牺牲一点舒适度。这也是工程上普遍采用的做法和代码里硬约束的默认写法相比更抗风险。5.3 能跑通但图是空白的现象程序不报错cvx_status也是Solved但TSDRplots.m画出来的图里没有任何曲线或者只有坐标轴。原因绘图脚本里引用的变量名和实际工作区不一致。比如TSDRplots.m里画的是X_mpc但你调用时传入的是X_opt变量名对不上Matlab不会报错只会画一个空图。另一个原因是load(TSDRdata.mat)之后工作区里的price变量被覆盖成了空数组。解决在运行出图前用whos和size检查关键变量的维度。我一般会写一个debug_plot.m专门打印所有要传给绘图函数的变量尺寸。在脚本里传给TSDRplots前加一行断言assert(size(X_mpc,2) length(price), 状态轨迹与价格序列长度不一致);这个断言能拦截90%以上的空白图问题。5.4 算得很慢一步要几分钟现象单次MPC求解耗时超过30秒滚动时域整个仿真跑下来要半小时以上。原因三个因素叠加一是用了默认的SDPT3求解器它对LP不如SeDuMi快二是预测时域N被改到48或72变量规模翻倍三是约束里写了非线性的写法比如直接在constraint里用max函数导致CVX把它当成锥规划来处理求解难度上升一个级别。解决先用cvx_solver sedumi换求解器这是最简单有效的一步。然后把预测时域从72降回24和电价发布周期匹配。最后审查TSDRMPC.m确保所有约束都是线性的——CVX虽然能处理一部分非线性凸表达式但每多一个非标准写法求解时间都会成倍增加。记住这套代码的规模和体量本质上就是一个线性规划或二次规划不应该让SDPT3去硬算。5.5 结果与论文复现差很多费用降幅远小于预期现象跑出来的总电费和论文图表里的结果差20%以上MPC的降费效果看起来微乎其微。原因大概率不是算法不对而是初始状态或预测时域没对齐。TSDRMPC的寻优结果高度依赖初始罐温和初始室温如果蓄热罐初始温度已经接近上限MPC就没有充电空间如果室温初始值已经接近Tz_maxMPC就无法预冷。另外预测时域内的电价序列必须和代码里的“日前调度”假设一致——如果你加载的是历史实测价格而不是日前预报价格MPC相当于开了上帝视角求出的结果反而不如离线最优控制因为两者看到的信息不对称。解决把初始状态设置成论文里的基准场景比如罐温15度、室温22度这样的“中间态”给MPC留足双向调节的余地。然后确认price里的前N个值确实对应“仿真开始时刻能获取到的预报值”而不是事后才整理出的真实值。这两个细节对齐之后MPC和离线最优的差距应该小到可以解释为滚动时域的预测误差。6. 一个进阶技巧把单次MPC改成滚动时域用参数扫描找最优预测时域单次求解MPC能出图但如果你要做真实的调度仿真必须把它放进滚动时域循环里每个时刻用当前实测状态做初始条件求解未来N步只执行第一步然后推进到下一个时刻。这样才符合MPC在实际系统中的工作方式。%% 滚动时域仿真模拟48小时 sim_len 48; x [params.Tz_init; params.tank_init_temp]; x_log zeros(2, sim_len); u_log zeros(1, sim_len); cost_total 0; for t 1:sim_len % 取当前位置开始的N步价格作为预测 if t params.predict_horizon length(price) p_forecast price(t : t params.predict_horizon - 1); else p_forecast [price(t:end), zeros(1, t params.predict_horizon - 1 - length(price))]; end % 调用MPC求解 [U, ~, ~] TSDRMPC(x, p_forecast, params); % 只执行第一步 u_apply U(:, 1); % 状态更新用实际扰动而不是预测扰动 x params.A * x params.B * u_apply params.W * weather(:, t); % 记录 x_log(:, t) x; u_log(:, t) u_apply; cost_total cost_total price(t) * u_apply * params.dt; end这段代码里有几个必须注意的细节。第一个是价格序列不足时的补零策略滚动到仿真末尾预测窗口会超出已知价格范围补零在数学上没问题但会让MPC在最后几步把蓄热罐里的能量“贱卖”掉因为它以为后续电价是零。更合理的做法是让预测窗口超出部分沿用最后一个已知电价或者缩短有效预测时域。第二个是状态更新用实际weather序列而不是预测值——这是滚动时域仿真与单次求解的本质区别预测和实际扰动之间存在偏差MPC通过每个时刻的实测状态来修正这个偏差。跑完这48小时之后你可以用TSDRparamPlots.m做一个简单的参数敏感性分析把params.predict_horizon从12扫到48每个值跑一遍滚动仿真统计总电费和最大需求。你会发现N24附近有一个明显的拐点小于24时蓄热罐不会在低谷时段充分充电大于36时对预测误差的敏感度上升结果反而变差。这是因为日前调度电价的信用结构只承诺24小时你预测得更远那些远期电价在结算时根本不会兑现相当于MPC在根据一个假的未来做决策。我最早一次跑这套程序的时候图出来了电费却没降多少。查到最后发现预测时域N被设成了6连一个完整的峰谷周期都没覆盖蓄热罐在谷价时段根本没有充电机会所有的优化自由度全被浪费了。从那以后我每次拿到类似MPC蓄热的工程包第一件事就是确认预测时域长度至少覆盖一个完整电价周期再谈算力优化和求解器选择。希望这套代码能帮你在需求响应这条路上少走几个弯路。本文还有配套的精品资源点击获取