做含光热电站的综合能源系统优化,说句实在话,最折磨人的不是数学公式难推,而是“模型逻辑理顺了,代码却跑不通”或者“跑通了,结果又不符合物理直觉”。这篇我把一个可以实际落地的“包含光热电站的综合能源系统优化运行规划”项目完整拆开,从物理建模、约束推导、Yalmip+Cplex求解实现,到结果分析和调参避坑,一次性讲透。适合正在做综合能源调度、光热电站建模仿真、或者刚接触MATLAB优化工具箱的研究生和从业者参考。
1. 项目整体思路:为什么光热电站是综合能源系统的“特殊杠杆”
1.1 这个项目到底解决什么问题
先给这个项目定个调。所谓“综合能源系统优化运行规划”,说白了就是在电、热、冷等多种能源耦合的园区或区域级系统里,回答两个问题:未来24小时(或更长时间尺度)每个机组出多少力、储能设备充多少放多少、从电网买多少电,以及这套运行方案对应的总成本是多少。
项目标题里有两个关键限定:一是“包含光热电站”,二是“优化运行规划”。前者决定了你要重点处理的是光热电站这种带储热环节的特殊电源,后者决定了整个项目的输出形式是一套多时段的优化调度方案。这套方案的决策变量通常是各机组的时序出力、储热罐的充放热功率、与外网购电功率,目标函数一般取系统总运行成本最小化或者碳排放最小化,约束条件包括功率平衡约束、设备容量约束、爬坡约束、储能状态约束等。
我觉得这个项目最值得讲透的一个点就在这:光热电站不是普通的光伏或者风电,它的特殊性在于“集热—储热—发电”三个环节是解开耦合的。光伏是“太阳来了就发,没太阳就没法发”,光热却可以通过储热罐把中午富余的太阳能存起来,到晚高峰再释放出来发电。这就让光热电站天然具备了时间搬移能力,是综合能源系统里一种非常优质的调节资源。
所以,这个项目不仅仅是套一个优化模型跑一遍结果,而是要通过数学模型把光热电站“源—储—机”解耦的物理特性表达精准。做这个项目的建议路径是:先搞清楚物理结构与能量流,再转成数学约束,最后写Yalmip代码用Cplex求解。这样理解出来的模型,改一改就能用在别的园区级IES项目上。
1.2 技术路线选型:MATLAB+Yalmip+Cplex凭什么能当主力
优化建模领域能用的组合其实不少:Python调Pyomo、Gurobi、Cplex,MATLAB直接调linprog/intlinprog,或者像项目标题这种用MATLAB配Yalmip工具箱再对接Cplex求解器。
MATLAB+Yalmip+Cplex这套组合能成为学术和工程项目的常见阵容,原因很实际。MATLAB的矩阵操作直观,做时序调度模型的变量声明和约束堆叠很方便;Yalmip提供了一套高层抽象语法,你不用手写CPLEX的LP/MILP矩阵格式,只需把约束和目标函数一行行“写出来”,它会在底层自动转换成求解器能识别的标准形式;Cplex是公认处理混合整数线性规划(MILP)效率极高的商用求解器之一,对大规模多时段优化问题的求解速度非常可靠。
有人会问,那Gurobi也很强,为什么这里选Cplex?我的经验是这两款求解器在中小规模MILP问题上差距并不悬殊,但Cplex在MATLAB环境下的集成路径成熟、学术许可申请方便,而且很多论文的对比基准都用的Cplex,复现对比起来更有参照性。你如果机器上同时装了Gurobi,用Yalmip也只需一行代码换求解器,不影响模型主体。
1.3 系统拓扑与能量流设计
为把这个项目讲清楚,我设定一个典型的园区级综合能源系统拓扑。系统内包含风电机组、光伏阵列、光热电站(CSP带储热)、燃气轮机(热电联产)、电储能电池、电锅炉,它们共同供应该区域24小时的电负荷和热负荷,允许从上级电网购电作为补充。
能量流的方向大致是这样的:电能侧由风电、光伏、光热发电模块、燃气轮机发电和电网购电共同供给,电负荷优先由可再生能源出力满足,缺额由燃气轮机和购电兜底;富余时由电储能吸收,必要时电锅炉消耗电能转成热能。热能侧由燃气轮机余热、光热电站的抽汽/换热输出、电锅炉产热共同满足热负荷需求。这套结构不算最复杂,但覆盖面足够:既有新能源随机性,又有储热储电双重灵活性,还有热电联产这种典型多能耦合单元。
明确系统拓扑后再去做优化建模,你会发现自己不会在写约束的时候“缺胳膊少腿”。这是整个项目最不该跳过的一步。
2. 光热电站建模:全系统最有技术含量的核心模块
2.1 光热电站的物理结构与能量流拆分
光热电站(Concentrating Solar Power, CSP)的物理结构一般包含三大块:集热场、储热系统(TES)和发电模块。集热场通过反射镜或塔式定日镜把太阳直射辐射(DNI)聚焦到吸热器,加热导热流体;储热系统通常采用双罐(冷罐和热罐)存储热能;发电模块从热罐抽取导热流体经换热产生蒸汽,推动汽轮机发电。
做数学建模的时候,不要把光热电站整个当成一个黑箱。最简洁有力的做法是把它的能量流拆成三段来建模:第一段是集热场吸收太阳热量,这部分取决于DNI、镜场面积和综合光热转换效率;第二段是热量的分配环节,集热场产出的热功率可以进入储热罐存储、也可以直接进入发电模块,当产热量有富余且储热罐已满时还会产生弃热;第三段是储热罐与发电模块之间的耦合,发电模块的出力由从罐中释放的热量决定。
这样拆开之后,光热电站变成了一个“输入太阳辐射、输出电功率和可能的热功率、中间夹着一个能量缓冲容器”的模块。这种表达方式的好处是,它天然地支持了时间搬移能力,也让Cplex求解时能清楚地看到每一个能量端口的状态变量。
2.2 关键约束的数学表达与推导
代码里最先要确定的是集热场热功率。给定时刻t的太阳直射辐照度 DNI_t、镜场面积 A_sf、综合效率 η_sf,集热场产出的热功率为:
H_sf(t) = η_sf × A_sf × DNI_t
这里的 η_sf 是镜场光学效率、吸热器效率和管路热损失的折合值,工程上通常在0.4到0.55之间,看具体镜场设计。
接下来是最核心的“热量分配”约束。任一时刻 t,集热场产热 H_sf(t) 等于进入储热罐的充热功率 H_ch(t)、直接进入发电模块的热功率 H_pb(t) 和弃热功率 H_spill(t) 之和:
H_sf(t) = H_ch(t) + H_pb(t) + H_spill(t)
这个约束是光热电站建模的第一关键。它保证能量守恒,同时引入了弃热(可以理解为光热版本的“弃光”)这种惩罚项,避免模型通过凭空消纳热量来骗优化结果。
储热罐的时序状态是整个模型的第二关键。设储热罐 t 时段末的储热量为 E_tes(t),充热和放热效率分别为 η_ch 和 η_dis,则:
E_tes(t+1) = E_tes(t) + η_ch × H_ch(t) - H_dis(t) / η_dis
这个递推式看起来简单,但有个容易踩的坑:储热罐的状态变量在时段上是“多一拍”的。如果模型跑24小时,储热容量变量一般要定义成25个点,否则边界时刻的递推会断开。代码调试时如果发现E_tes(t+1)索引越界或者结果里第一小时储热状态异常,先检查这个维度定义对不对。
发电模块约束则包括最大最小出力限制和爬坡约束:
P_pb_min ≤ P_pb(t) ≤ P_pb_max
-P_ramp ≤ P_pb(t) - P_pb(t-1) ≤ P_ramp
其中 P_pb(t) 是发电功率,它与消耗的热功率 H_pb(t) 之间用发电效率 η_pb 关联:P_pb(t) = η_pb × H_pb(t)。为了保持线性模型,通常把 η_pb 当作常数,虽然汽轮机的实际效率随负荷率变化,但做日前调度这种时间尺度,常效率假设在工程上是可接受的。
储热罐的物理限制也很关键:储热量有上下限,充放热功率有上限,一般还会加一个调度周期末储热量等于初始储热量的约束,保证储能不是“一次性用完”的不可持续方案:
E_min ≤ E_tes(t) ≤ E_max
0 ≤ H_ch(t) ≤ H_ch_max
0 ≤ H_dis(t) ≤ H_dis_max
E_tes(T) = E_tes(0)
这套约束写完后,光热电站的数学建模基本就完整了。
2.3 电热双输出模式的扩展建模
前面说的是光热电站只发电的情况。如果项目扩展到包含热负荷的综合能源系统,一个更实用的做法是把光热电站建模成热电联产模式(CSP-CHP),让它不仅能发电,还能抽出部分热能直接供给热负荷。
扩展思路也不复杂:在热量分配约束里增加一项供热功率 H_hx(t),把发电模块的冷凝余热或抽汽热量引到供热管网。对应的热量分配式变成:
H_sf(t) + H_dis(t) = H_pb(t) + H_hx(t) + H_spill(t)
这里的放热方向要重新梳理:实际运行时,进入发电模块的热量除了来自集热场直接产热,也可以来自储热罐放热。所以更严谨的表达是:t时段发电模块和供热模块的总输入热功率 = 集热场直接供热功率 + 储热罐放热功率。同时,储热罐充放热需要和集热场产热协调,避免出现“一边充一边放”的假性优化。
实操中,为了降低模型复杂度,很多人会把 H_ch(t) 和 H_dis(t) 设为互斥变量,但我的经验是:只要充放热功率上限和储热状态递推约束写完整,模型在目标函数成本驱动下自然不会有非物理的同时充放行为,不必强加整数变量来制造求解负担。
3. 优化模型构建与求解实现:Yalmip建模的完整链路
3.1 目标函数与成本项拆解
这个项目的目标函数我采用系统总运行成本最小化。成本项拆开来看包含五块:从上级电网购电费用、燃气轮机燃料成本、各设备运维成本、弃风弃光惩罚和弃热惩罚(如果2.2节的H_spill需要惩罚)。
目标函数数学表达为:
min Σ_t [ c_buy(t) × P_buy(t) + c_gas × F_gt(t) + Σ_i c_i × P_i(t) + c_curt × (P_wt_avail(t) - P_wt(t) + P_pv_avail(t) - P_pv(t)) ]
其中 P_buy(t) 是购电功率,F_gt(t) 是燃气轮机的燃料耗量,P_i(t) 代表各设备的输出功率,c_curt 是弃电惩罚系数。燃气轮机的燃料成本通常用二次函数逼近,但为了保持MILP模型线性,会做分段线性化处理,或者直接简化成线性效率模型。
有个经验要分享:目标函数里各项成本系数的量级不能差距太大,否则求解器会把优化重心全放在量级大的那项上,其他项失去调节意义。比如购电价格是0.8元/kWh时,运维成本如果是0.01元/kWh量级,建议在模型里统一单位的量纲,避免灵敏度失衡。
3.2 约束完整清单与Yalmip代码框架
完整约束除了2.2节光热电站的自身约束外,还需要系统层面的公共约束。电能平衡约束是每时段系统内各电源出力加总等于电负荷与电锅炉耗电之和;热能平衡约束是各热源产热加总等于热负荷。此外,风电和光伏的出力要小于等于预测可用出力,燃气轮机有出力和爬坡限制,电储能遵循与储热罐类似的时序递推约束,电锅炉出力受装机容量限制。
用Yalmip写这套模型的代码骨架,大致长这样:
%% 参数定义 T = 24; E0 = 0.5 * E_max; % 初始储热量 %% 决策变量 P_wt = sdpvar(1, T); % 风电出力 P_pv = sdpvar(1, T); % 光伏出力 P_pb = sdpvar(1, T); % 光热发电功率 H_ch = sdpvar(1, T); % 储热充热功率 H_dis = sdpvar(1, T); % 储热放热功率 E_tes = sdpvar(1, T+1);% 储热罐储热量 P_gt = sdpvar(1, T); % 燃气轮机发电 P_es = sdpvar(1, T); % 电储能净功率,正放负充 P_eb = sdpvar(1, T); % 电锅炉耗电 P_buy = sdpvar(1, T); % 购电功率 %% 约束集 C = []; % 光热约束 H_sf = eta_sf * A_sf * DNI; % 参数数组,非变量 C = [C, H_sf + H_dis == H_pb_to_th + H_spill]; % 热量分配 C = [C, P_pb == eta_pb * H_pb_to_th]; C = [C, E_tes(2:T+1) == E_tes(1:T) + eta_ch*H_ch - H_dis/eta_dis]; C = [C, E_tes >= E_min, E_tes <= E_max]; C = [C, E_tes(1) == E0, E_tes(T+1) == E0]; C = [C, 0 <= H_ch <= H_ch_max, 0 <= H_dis <= H_dis_max]; % 电平衡 C = [C, P_wt + P_pv + P_pb + P_gt + P_es + P_buy == P_load + P_eb]; % 热平衡 C = [C, eta_heat * (H_dis_from_csp) + eta_gt_h * P_gt + eta_eb * P_eb == H_load]; % 新能源上限 C = [C, 0 <= P_wt <= P_wt_forecast]; C = [C, 0 <= P_pv <= P_pv_forecast]; % 燃气轮机及储能限制 C = [C, P_gt_min <= P_gt <= P_gt_max]; C = [C, -P_es_max <= P_es <= P_es_max]; %% 目标函数 objective = sum(c_buy .* P_buy + c_gas * (P_gt ./ eta_gt) + ...); ops = sdpsettings('solver', 'cplex', 'verbose', 2, 'showprogress', 1); optimize(C, objective, ops);注意代码里有个细节:E_tes变量定义成了1×(T+1),这样递推式 E_tes(t+1) 在t从1到24时刚好覆盖所有时段。如果定义成1×T,最后一段会索引越界。
3.3 Cplex求解参数设置与求解策略
用Yalmip调用Cplex时,几个求解参数对结果和效率影响很大。第一个是MIP的相对间隙容忍度,sdpsettings里对应参数是'cplex.mip.tolerances.mipgap',默认是0.0001,如果模型规模大跑不动,可以放宽到0.01或0.05,工程上能接受1%~5%的优化偏差来换求解速度。
第二个是求解时间上限'cplex.timelimit'。我习惯设一个240秒或600秒的上限,避免偶然遇到难解的MILP模型时程序卡死。几百个连续变量加上几十个整数变量的问题,Cplex通常几秒内能解决,但一旦引入了机组启停0-1变量,求解时间会指数增长,这个时候控制时间上限就很重要。
第三个是显示级别'verbose',调试时设2能看分支定界进度,最终运行设1或0避免刷屏。另外一个很实用的参数是'cplex.threads',默认会使用所有CPU核,但如果是多任务并行测试,建议限制为4或8,防止求解器占满机器资源。
求解之前还有一个加分项:如果系统里没有0-1变量(比如不建模机组启停),这个优化问题本质是LP,纯LP问题Cplex几乎秒解。所以排查问题时建议先跑一个不含整数变量的简化版,确认模型物理上没矛盾,再逐步引入整数变量升级难度。
4. 典型日运行结果分析与可行性验证
4.1 场景设计与基础数据准备
为了让结果可复现,我给一个典型夏季场景的数据模板。电负荷曲线取典型的双峰曲线,热负荷夏季取夜间低谷期稍高。风电预测出力取夜间大、白天小,光伏出力取午间尖峰,DNI取晴天的钟形曲线,午间最高约达到800到900 W/m²。
光热电站关键参数可以按塔式50MW级设计:镜场面积约60万平方米,综合效率0.48,发电模块额定功率50MW,储热罐容量500MWh,最大充放热功率80MW,充放热效率取0.95。这个参数量级参考了实际塔式电站的设计数据,但不属于任何特定项目,方便计算。
我运行下来的典型结果呈以下特征:白天10点到14点是DNI高峰,也是光伏大发时段,此时电价低、光伏出力大,光热电站的调度策略倾向于把集热场产热量存入储热罐而不是直接发电;到18点到21点晚高峰,电价高、光伏为零,储热罐释放热量驱动发电机满发,配合燃气轮机共同支撑负荷。这种“正午储热、傍晚放热”的行为就是光热电站时间搬移能力的直观体现。
4.2 结果解读:储热罐如何起到削峰填谷作用
把储热罐24小时的SOC曲线拉出来看,会发现一个非常典型的“浅V型”变化:白天集热初期储热量逐步爬升,到下午15点左右达到最高点,之后随着放热发电快速下降,到24点回到初始储能量附近。整个曲线的平滑程度取决于储热容量和DNI波动的配合情况。
具体到运行成本层面的收益,主要体现在两个方面:一是晚高峰购电量的下降,原本需要高价从电网买的电,被储热罐里的热能替代了,这部分节省最直观;二是弃风弃光率下降,虽然光热本身不消耗风电,但它在负荷高峰提供的出力减轻了其他机组的压力,让燃气轮机不用满负荷爬坡,系统调节余量更大,间接增加了新能源的消纳空间。
一个值得注意的反直觉现象是:储热罐并不是越大越好。储热容量超过某个阈值后,因为集热场产热总量有限,储热罐在白天根本“吃不饱”,大容量反而让储热量长期处在低水平,弃热风险上升,投资收益比下降。这个拐点要通过敏感性分析去找,也正是这类优化项目里最有意思的工程结论。
4.3 对比分析:无储热光热与带储热光热的差别
我在同样场景下做了一个对照组:光热电站强制取消储热系统(发电功率必须等于集热场瞬时产热量),其他参数保持不变,对比两轮优化结果。
对比结论非常清晰:带储热的光热系统总运行成本下降了约11%到15%,晚高峰购电峰值功率下降更明显。这个降幅的区间会随DNI和负荷曲线变化,但方向是稳定的。原因在于储热打破了光热电站“只能即时发电”的约束,让光热电站在一天中从“被动跟随太阳”变成“主动匹配电价和负荷”,灵活性价值直接体现在经济运行收益上。
另外一个有意思的对比指标是弃热率。带储热方案里,DHI波动导致的瞬时弃热可以被储热罐吸收,弃热率明显低于无储热方案。这两个指标放一起,就能比较完整地论证光热储热模块在综合能源系统中的必要性。
5. 常见问题与调优经验实录
5.1 模型求解慢或者内存爆掉怎么办
遇到求解慢先别急着怪Cplex,绝大多数情况是模型本身“问得太难”。排查顺序我一般是这样:先数数有多少0-1整数变量,一般来说,整数变量数量几百个以内,对Cplex都不是难事,但如果你引入了机组启停和分段线性化组合,整数变量可能膨胀到上千甚至数千个,这时候就要考虑简化模型结构。
常用的降复杂度手段包括:一是把机组的启停状态从逐时段0-1变量改为最少运行/停机时间聚合约束,减少变量数量;二是把分段线性化的段数从4段降到3段;三是把24小时等间隔时段改成非均匀时段(比如高峰间隔1小时、夜间间隔2小时),但这会丢失精度,不建议在论文里这么做。还有一个“伪降复杂度”的招是给每个0-1变量加上本地边界,比如已知某时段不可能停机就直接把对应状态固定为1。
如果模型不是整数规模问题而是纯LP也慢,那大概率是数值问题,见下一节。
5.2 数值病态与收敛性问题的处理
做这类模型,最容易被忽视的是量纲不一致。功率用kW,储热量用kWh,而热量接力用的却是MW和MWh,混在一起以后,约束矩阵里有些系数是几百、有些是0.001,求解器的数值稳定性会变得很差,出现“初始可行解找不到”或“对偶间隙不收敛”这类诡异现象。
我的建议是建模一开始就统一单位制,最常见的是全部折算成MW和MWh。具体做法是:电功率用MW,储热罐储热量用MWh,热功率也用MW,所有能量递推公式里功率乘以时间步长才是能量,在24小时模型里时间步长是1小时,所以数值上刚好1MW×1h=1MWh,不会出现数量级混乱。
另外,给所有变量设置显式的上下界也非常重要。Yalmip虽然能自动推断很多变量的边界,但像E_tes这种递推变量,不给定界会让求解器在分支定界时做大量无用搜索。每一条约束都不妨问自己一句:这个变量的物理上下限是什么?写上去,求解速度往往能提升一个档次。
5.3 Cplex安装与多求解器共存的问题
搜索热词里很多人问“MATLAB关联安装了Gurobi,还能装Cplex吗”,答案是肯定的。Cplex和Gurobi在MATLAB里通过Yalmip调用时没有直接冲突,它们以不同的求解器文件存在,Yalmip会按路径识别。关键是在安装后确认addpath路径和求解器优先级设置没有问题。
在Yalmip中,你可以用sdpsettings('solver', 'cplex')显式指定使用Cplex,也可以用solverString变量切换求解器。如果设了默认求解器之后发现Cplex没被识别,先在MATLAB里运行yalmiptest,这个命令会列出所有能被Yalmip识别的求解器状态。Cplex显示“found”但调用报License错误,通常是许可文件的路径或者环境变量没配对,重新执行安装向导里的许可证配置步骤就能解决。
还有一个常见坑是Cplex版本与MATLAB版本不匹配。高版本MATLAB对于非常老的Cplex动态链接库不兼容,建议使用近两三年发布的Cplex版本配合当前主流MATLAB版本,这种搭配问题最少。真遇到“无法加载动态链接库”的报错,优先检查电脑有没有安装对应的MATLAB运行时组件,而不是急着重装Cplex。
按照这套流程操作下来,我自己的经验是:绝大多数“代码跑不通”“结果不合理”的问题都能在建模和数值层面找到原因。整套模型的核心价值在于,光热电站这个带储热环节的特殊单元,为综合能源系统提供了一个在时间维度上自由搬移热能的可调节支点,把这个支点建模建清楚,后续无论扩展碳交易、需求响应还是多场景随机优化,都只是在现有骨架上加约束和变量的问题。