项目标题里最核心的四个词，“源荷随机特征”“热电联供微网”“优化”“Matlab实现”，分开看都不算冷门，但组合在一起就是典型的综合能源系统方向的研究课题了。这两年做微电网优化的朋友应该都能感觉到，纯电微网的套路已经非常成熟，真正容易出问题也容易出成果的，恰恰是这种带热负荷、带热电联供机组、还要把风电光伏和负荷随机性一起塞进优化模型的场景。这篇文章我就围绕这个题目，把从模型拆解到Matlab代码落地的整个思路捋一遍，也会把我自己调试过程中踩过的坑一并写出来，给正在做相关课题或者准备复现类似工作的同学一个参考。
这个研究到底解决什么问题？一句话来说：在风电、光伏出力以及电/热负荷都存在不确定性的情况下，怎么安排热电联供机组、燃气锅炉、储能等设备的出力计划，让系统在满足供电和供热需求的同时，运行成本最低、排放尽量小，而且还能扛得住源荷波动。适合正在做微电网优化调度、综合能源系统规划或运行优化的同学看，也适合刚接手Matlab+Yalmip这类工具、想快速跑通一个完整算例的入门者。
1. 项目解读：热电联供微网优化到底在做什么
1.1 源荷随机特征这句话怎么拆解
“源”指的是风电机组和光伏阵列这种出力受天气影响的分布式电源，出力曲线不是人为可控的，来一阵云、起一阵风，出力就变了。“荷”指的是用户侧的电负荷和热负荷，居民用电有早晚高峰，工厂用电跟着生产班次走，冬天的热负荷还跟室外温度强相关。这些源和荷都有一个共同特点：用确定性曲线来描述它们是不充分的，实际运行中它们一直在随机波动。
做热电联供微网优化研究时，如果无视这种随机性，把光伏出力设成一个固定的预测曲线，然后对着这个曲线去做优化，最后得到的结果通常是非常乐观的。遇到实际运行时风光出力比预测低，或者负荷比预测高，原先的调度方案就会出现功率缺额，供热温度也可能掉出允许区间。所以课题名里专门强调“考虑源荷随机特征”，本质上就是要让优化模型在面对不确定性时依然有足够的鲁棒性和经济性。
随机特征具体怎么刻画，一般有两层意思。第一层是概率分布，比如风速近似服从Weibull分布，太阳辐照度可以用Beta分布描述，负荷预测误差通常假设为正态分布。第二层是时间相关性，源荷的预测误差在相邻时段之间不是相互独立的，而是有自相关的特性。这两层信息最终都可以转化成优化模型中的随机场景，用多个场景来代表可能发生的源荷情况。
1.2 为什么需要联合考虑源和电负荷、热负荷的随机性
很多初接触这个课题的朋友会有一个疑问：既然热负荷相对平稳，供热需求的变化也比电负荷缓慢，那能不能只考虑电负荷和风光出力的随机性，热负荷直接用确定值？我在实际做课题时发现，这种简化在纯供热系统里问题不大，但在热电联供微网里会带来明显偏差。原因在于热电联供机组存在“以热定电”的运行约束，也就是机组的发电功率和供热功率是耦合在一起的。热负荷一旦出现波动，机组的电出力上下限就会跟着变，电功率平衡也会受影响。
举个直观的例子，一个严寒天的早晨，热负荷预测值假设是10MW，但实际因为气温骤降变成了11.5MW。机组为了多供1.5MW热量，其电出力就会被抬高或者压低（具体要看机组运行区间和运行模式）。如果优化模型里没有把这种热负荷的随机波动考虑进去，系统中其他电源的计划就可能完全偏离实际需求，储能充放电计划也会错位。因此，至少要把电、热两类负荷以及风光出力都纳入随机场景的构建，才算真正做到了“源荷随机特征”。
联合考虑还有一个实际收益，就是可以利用电热负荷之间的时序互补性。比如夜间风电出力大、电负荷低，这时候可以通过电锅炉或者热泵把多余电能转化为热能存储起来，白天再放热。但如果电负荷和热负荷的随机性没有被描述在同一个框架里，这种协调优化的机会就很难被模型捕捉到，做出来的调度方案自然不够经济。
2. 系统建模与随机性处理
2.1 热电联供微网的典型结构与设备模型
要搭建优化模型，首先要把微网里有哪些设备、怎么连接、能量怎么流动，用数学语言描述出来。一个典型的热电联供微网包括以下模块：
- 热电联产机组（CHP）：这是系统的核心供能设备，可以同时输出电功率和热功率。它的技术特点是热电出力比在一定范围内可调，通常用多边形的可行运行域来描述。简化建模时也可以用固定热电比，但精度稍差。
- 燃气锅炉：作为补充热源，在CHP供热不足时启动，特点是响应快、热效率高，但消耗燃气。
- 风电机组和光伏阵列：可再生能源电源，出力具有随机性，通常作为负的负荷或者可控程度有限的电源。预测值取自气象数据，随机性在后面用场景描述。
- 电储能（电池）和蓄热罐：实现电、热两种能量在时间维度上的平移。电池有充放电功率限制和容量约束，蓄热罐有储/放热功率限制和热容量上限。
- 电负荷和热负荷：用户需求，分为预测值和实际值。预测值用于场景生成，实际值用于仿真验证。
设备模型的核心是输入-输出关系。以CHP机组为例，常见建模方式是：
[ C_{chp,t} = a \cdot P_{chp,t} + b \cdot H_{chp,t} + c ]
其中(C_{chp,t})是燃料成本，(P_{chp,t})是电功率，(H_{chp,t})是热功率，(a)、(b)、(c)是效率系数。这个二次函数在求解时往往线性化成分段线性函数。电池模型需要考虑充放电效率，蓄热罐模型需要考虑热损失系数。把这些数学模型写清楚，优化模型的基础才算打牢。
2.2 风光出力与负荷的随机特征建模方法
源荷随机特征的建模，本质上是把不确定性从“不可描述”变成“可以用场景描述”。我在实际项目中用过的比较可行的方法有这么几种：
第一种是参数化概率分布法。风速用两参数Weibull分布拟合，太阳辐照度用Beta分布拟合，负荷预测误差假设为均值为0的正态分布。这种方法的优点是简单，缺点是无法表达复杂的时空相关性。
第二种是场景生成法。基于历史数据，通过蒙特卡洛抽样或者拉丁超立方抽样，生成大量包含时序相关性的源荷场景。注意这里不能直接按独立分布抽样，否则生成的负荷曲线会非常毛糙，和真实负荷曲线的平滑特性不符。一般做法是先对预测误差的时间序列用ARMA模型或者ARIMA模型建模，再从这个模型里抽样得到误差场景，叠加到预测曲线上。
第三种是数据驱动的场景生成，比如用生成对抗网络去拟合源荷联合分布。这种方法适合数据量比较大的场合，我在论文复现里用得不多，因为代码复杂度和训练成本都比较高，对课题组没有专门算力的情况不太友好。
实际工程中，最稳的组合是“预测曲线 + ARMA误差抽样 + 拉丁超立方采样”。既能保证时序相关性，又能用较少样本覆盖大概率空间，后续的场景削减也方便。
2.3 场景生成与削减：让随机问题变成可计算问题
有了随机场景之后，优化问题就变成多场景下的随机优化问题。理论上场景数越大，不确定性描述得越精准，但计算量也越大。一个24小时调度问题，如果场景数取1000，每个场景有24个时段，决策变量就会膨胀到几十万个，尤其加入二进制变量之后，求解难度会急剧上升。所以实践中几乎都会做场景削减。
比较常用的削减方法是快速前向选择法或者同步回代消除法。核心思想是：在保留原始场景集合概率特征的前提下，删除那些与其余场景距离较近、代表性弱的场景，把1000个场景削减到50个或者20个。每个保留场景赋予一个概率权重，所有保留场景的权重之和等于1。
用Matlab实现场景削减时，我偏向直接用现成的函数。如果不想依赖额外工具箱，可以自己写一个简单的贪心削减：每次找一对欧式距离最小的场景，删除其中一个，把被删场景的概率加到保留场景上。虽然不如同步回代消除法那么精确，但对大多数课题研究来说已经足够了。
场景削减后的一个细节是，削减后的集合所对应的期望成本会和原始场景集合的期望成本有偏差。做结果分析时要留意这个偏差的规模，一般控制在5%以内就可以接受。如果偏差太大，说明初始场景数不够或者削减比例过高，需要回调参数。
3. 优化模型与求解策略
3.1 目标函数：经济性、环保性与供能可靠性的取舍
热电联供微网优化的目标函数通常不是单一的。最经典的是运行成本最小，成本项包括：
- 燃气成本：CHP机组和燃气锅炉消耗的燃料费用。
- 购电成本：微网与配电网交换功率的费用。如果允许向配电网卖电，还要加入售电收益项。
- 设备启停成本：CHP机组和锅炉的启动成本，这个虽然单次不大，但对机组频繁启停有一定抑制作用。
环保性作为第二个目标时，可以引入碳排放成本，将CO2排放量按碳价折算到目标函数中。也可以用多目标优化的形式同时优化经济性和排放量，但多目标问题的求解复杂度比单目标高不少。研究性课题里，把排放折算成经济成本放入单目标函数，处理起来简洁得多。
供能可靠性作为第三个目标时，最常见的处理方式是引入失负荷惩罚。在随机优化框架下，如果一个场景下实际发生的负荷超过了可用出力，那么缺额部分就会被切掉，切负荷量乘以惩罚系数进入目标函数。这个系数如果设得太低，优化结果会保持很小的备用，经济性很好但风险大；设得太高，备用过大，成本上升。我自己的经验是，惩罚系数至少取单位购电成本的10到20倍，才能有效抑制切负荷的发生。
3.2 约束条件：功率平衡、爬坡、储能、热负荷平衡
约束条件是整个模型里最容易出bug的地方。我列一下必须覆盖的几类，并说明它们各自在模型中起什么作用。
电功率平衡约束要求每个时段内，所有电源出力加上储能放电、加上从配电网购入的电功率，等于电负荷加储能充电加向配电网售出的功率。如果模型允许切负荷，右侧还需要加上切负荷变量。这条约束如果写得不对，求解出的结果要么明显违反物理规律，要么干脆无解。
热功率平衡约束类似，但要注意的是CHP机组的高温出水经过换热站后供给热负荷，蓄热罐可以吸收或释放热量。蓄热罐的蓄热状态量还必须满足时序递推方程，并受容量的上下限约束。
机组爬坡约束描述的是CHP机组在相邻时段之间电出力的变化速率限制。如果不加这条约束，优化模型可能会给出一个锯齿状出力曲线，现实中机组完全跟不上。同类约束也适用于燃气锅炉的出力变化。储能约束则包括充放电功率上限、容量上下限，以及同一时段不能同时充放电的逻辑约束，这个逻辑约束通常用二进制变量实现。
还有一个容易被忽略的约束是CHP运行区间约束。实际机组的电出力和热出力不是任意组合都可以的，存在一个可行域，需要用一组线性不等式去逼近。不少初学的同学只在模型里写了电出力和热出力的上下限，结果得到的调度点落在了机组的不可运行区域，仿真时才发现不对。
3.3 求解思路：混合整数线性规划与场景法的结合
优化模型经过线性化处理后，因为里面含有二进制变量，所以是一个混合整数线性规划（MILP）问题。场景法把随机性问题转化为一个大规模确定性MILP问题，然后用成熟的商业求解器求解。
在Matlab环境里，用Yalmip作为建模工具，底层调Cplex、Gurobi或者开源的Cbc都可以。我个人最常用的是Yalmip+Gurobi的组合，建模代码简洁，求解效率也高。Gurobi对MILP问题的剪枝效果在同类求解器里数一数二，对场景规模较大的算例，能明显感觉到时间优势。没有Gurobi许可证的话，Cplex也可以，Cbc免费但求解速度要慢不少，适合模型比较小的验证场景。
求解思路的细节是：把所有场景下的决策变量统一放进一个大的优化问题里，场景之间的变量通过“非预期性约束”耦合。最简单的耦合方式是让第一个时段的决策变量在所有场景下保持一致，因为实际运行时在第一个时段还没有获取未来信息，调度决策不能因为场景不同而变化。如果忽略这一点，就会得到“事前知道未来”的非预期性解，结果过于乐观。这一点在做随机优化研究时非常关键。
4. Matlab代码实现要点
4.1 代码整体结构与数据准备
拿到这个课题之后，我建议先把代码的结构想清楚，不要一上来就写优化模型。一个可维护的Matlab代码通常分成四个部分：数据输入、场景生成、模型构建、求解与结果输出。
数据输入部分需要准备的基础数据包括：24小时的电负荷预测曲线和热负荷预测曲线，风速和光照强度的历史序列或预测序列，CHP机组、燃气锅炉、电池、蓄热罐的技术参数，天然气价格、购电电价等经济参数。这些参数用结构体或者表格存起来就行，我用得比较习惯的是把所有参数放在一个结构体参数表params里，后续模型里引用起来非常方便。
场景生成部分把源荷不确定性转化成场景集合。这里要留意的是，用来做随机优化的风电、光伏和负荷场景，需要来自同一个随机源，不能各编各的。比如用ARMA模型生成误差序列时，得把风、光、电负荷、热负荷的误差序列一起模拟，保留它们之间的相关性。一种简便做法是假设各随机变量之间的相关系数矩阵，然后用经Cholesky分解处理后的随机数序列来抽样。
模型构建部分，Yalmip的语法我个人觉得是能提高效率的。每个变量用sdpvar定义，二进制变量额外指定类型，约束用循环填入一个cell数组，目标函数直接写表达式，最后调用optimize。代码可读性比直接用Cplex的底层API要高很多。
求解与结果输出部分，除了打印出目标函数值和各设备出力，还应把电平衡、热平衡的残差输出来，方便定位约束写错的地方。仿真结果如果画成图，建议把预测场景和实际场景画在同一条时间轴上，能直观看出随机优化的价值。
4.2 核心函数与求解器调用（Yalmip/Gurobi举例）
下面给一个Yalmip建模的核心代码片段，覆盖CHP和储能的基本约束。这个是我调试过很多次的模板，写给大家做个参考。
% T为时段数，Ns为场景数，x可能包含三维变量（时段、场景、设备） P_chp = sdpvar(T, Ns, 'full'); % CHP电出力 H_chp = sdpvar(T, Ns, 'full'); % CHP热出力 P_bat = sdpvar(T, Ns, 'full'); % 电池出力，正为放电，负为充电 S_bat = sdpvar(T+1, Ns, 'full'); % 电池荷电状态 u_chp = binvar(T, Ns); % CHP启停状态 Constraints = []; % CHP运行区间约束（简化） Constraints = [Constraints, P_chp <= 5000*u_chp, P_chp >= 0]; Constraints = [Constraints, H_chp <= 3000*u_chp, H_chp >= 0]; % 电池充放电约束 Constraints = [Constraints, S_bat(:,1+1) == S_bat(:,1) + ...]; % 实际写的时候按“上一时刻SOC + 充电功率*效率 - 放电功率/效率”递推 Constraints = [Constraints, -200 <= P_bat <= 200]; Constraints = [Constraints, 0 <= S_bat <= 1000]; % 电功率平衡 Constraints = [Constraints, P_chp + P_wind_scene - P_load_scene + P_bat + P_grid == 0]; % 目标函数 Objective = sum(sum(C_gas * ... )) + sum(sum(C_grid .* P_grid));这段代码虽然简化，核心意思很清楚：所有设备变量都按“时段 x 场景”两个维度展开，约束逐条添加。注意binvar声明二进制变量的位置，如果场景数很大，二进制变量会非常多，这时候尽量在Yalmip里用维度矩阵方式声明，比用循环一件一件定义要快得多。
调用Gurobi之前，不要忘了在代码开头设置求解器参数：
options = sdpsettings('solver','gurobi','gurobi.MIPGap',0.001,'verbose',2); optimize(Constraints, Objective, options);MIPGap设为0.001表示允许相对最优间隙在千分之一，既保证精度又不至于卡死。如果是调试阶段，可以放宽到0.01，会快很多。
4.3 参数标定与代码调试心得
参数标定这件事，想认真做起来比写优化模型还费时间。设备参数一般在手册上都能查个大概，但实际使用时要根据自己系统的规模调整。比如CHP的燃料系数a、b、c，如果直接把某区域电网的平均供电煤耗套进去，成本和设备实际消耗的对应关系就会失真。我的习惯是先做一次确定性优化的基准运行，把结果和历史运行数据对照一下，如果各设备出力比例和实际相差很大，就反推是哪个效率参数偏了。
代码调试的第一个重灾区是量纲一致性。Matlab本身不自动检查单位，最容易出现的是电功率用了kW，但热功率还停留在MW，导致热平衡约束完全失效。我建议在数据准备阶段就统一单位，比如全部用kW和kWth，天然气消耗量用kW或者m³单独转换确认。第二个重灾区内置变量名冲突，比如循环变量i和表达式里面的虚数单位i冲突，这个在复数运算多的脚本里特别容易出问题。第三个是Yalmip的x、y等函数容易和变量名撞车，所以自定义变量名尽量用带语义的字符串。
调试的时候建议先跑一个单场景的确定性优化，确认代码逻辑没有报错、结果趋势合理，再扩展成多场景随机优化。这样定位问题会快很多，不然几十个场景一起跑，一旦出问题很难判断是模型问题还是场景数据问题。
5. 仿真案例与结果分析
5.1 典型日场景设置
为了说明随机优化的价值，我做了一个典型日仿真，系统包括1台CHP机组、1台燃气锅炉、一组电储能和一组蓄热罐。风电和光伏各配置了一定容量，并给它们设置了比较强的波动场景。
电负荷预测曲线设成了传统的双峰曲线，早高峰在8点到11点，晚高峰在18点到22点。热负荷则是夜间和清晨偏高，白天略低，这是北方供暖期典型特征。风电场景设成了夜间大发、白天减少的模式，光伏则相反，集中在中午出力。这个时间是典型案例：晚上风电多但电负荷相对低，正好可以用来蓄热或者给电池充电；白天光伏大发，但晚高峰时光伏出力已经衰减，需要靠储能和机组顶上。
场景生成阶段，我用ARMA(1,1)模型生成了误差序列，初始场景数取500，削减后保留20个。这样模型规模不至于太大，又能比较充分地代表不确定性。为了对比，还设置了一个“确定性优化”对照组，即以预测曲线为唯一场景进行优化。
5.2 对比实验：确定性优化 vs 随机优化
两种方式下，优化得到的机组调度计划有显著差异。确定性优化给出的方案在一些时段非常激进，比如它会在夜间风电预测出力较大时尽量压低CHP电出力，把多余的电用来给电池充电或蓄热罐蓄热。这个计划在预测曲线成立的条件下成本很低，但一旦实际风电出力达不到预测值，系统的电功率就会出现缺口。多场景模拟回测时，确定性方案出现切负荷和切热负荷的时段数明显高于随机优化方案。
随机优化因为同时考虑了20个场景，会在调度时主动留出一部分备用空间。例如在风电出力预测较大但波动也大的时段，它不会把所有风电都视为可靠出力，而是把储能的充电功率稍微压低，或者让CHP保持在一个更稳妥的电出力水平。代价是随机优化期望成本比确定性优化高出4%到6%，但最恶劣场景下的失负荷量大幅降低，而且供热可靠性指标明显更好。
这个“成本升高但风险下降”的权衡，正是随机优化研究的核心价值所在。做结果分析时，不能只看成本绝对值，还要把可靠性指标和风险指标放到一起看。我习惯用条件风险价值CVaR来衡量风险，把最差5%场景的成本平均值作为风险指标，这样优化时就可以在目标中加入CVaR项，得到更稳健的调度方案。当然这个属于扩展内容，基础版本只做期望成本最小化也完全可以。
5.3 关键结果解读与能效指标
仿真结果输出后，需要用几个关键指标来衡量系统的运行水平。第一个是CHP平均热电比，用来判断机组工作点是否落在较优区间。第二个是电储能和蓄热罐的日循环次数，如果一天内频繁满充满放，说明储能被过度调用，实际运行寿命会受影响。第三个是综合能效比，即系统输出电热之和与输入燃料热值之比。热电联供系统的优越性就体现在这个综合能效比上，一般可以做到0.8以上。
我在结果里还专门检查了热平衡约束的残差。随机优化后，由于蓄热罐在其中起了缓冲作用，大多数时段热平衡都是严格满足的。但在极端场景下，出现过蓄热罐达到容量上限导致多余热量无法吸收的时段，这时候系统不得不通过降低CHP热出力来配合，热负荷若仍不足，就必须由燃气锅炉顶上。这种边界情况如果不做仿真回测，单看优化结果是察觉不到的。
做结果分析时，另一个值得关注的细节是：随机优化中不同场景下同一时段的CHP出力会有一定散布范围，这个散布范围其实是系统应对不确定性的“策略空间”。可以把它画成概率带图，和确定性优化的单点曲线作对比，视觉冲击力很强，也便于在论文或汇报中说明随机优化的必要性。
6. 常见问题与排查技巧实录
6.1 场景数量怎么选才合理
场景数量的选择没有绝对标准，但有一条经验区间：如果采用蒙特卡洛抽样，初始场景数取200到1000之间，削减到10到30个，效果通常都不错。太少，比如只留5个场景，随机性表达不够，优化结果会和确定性方案差不多；太多，比如留100个，计算时间成倍增长，而且贡献有限。
我自己试过初始场景500、削减到50和削减到20两种情况，期望成本差异不超过2%，但求解时间后者快了一倍多。所以如果只是做课程设计或者论文中的一个算例，我建议削减到15到25个即可。如果是更复杂的热电联供系统，需要加入机组组合变量，建议场景数控制在20以内，否则MILP的整数变量会撑爆内存。
6.2 求解慢、内存占用高怎么办
场景法最麻烦的问题就是求解时间太长。慢的时候，一个24小时、20场景、含二进制变量的MILP问题可能也要跑几分钟。如果模型规模大，跑几个小时也正常。我常用的加速手段有三个。
第一，尽量用连续变量近似二进制变量，这需要仔细分析哪些逻辑约束是必须的。比如电池充放电互斥约束可以先用一个很小的惩罚项代替，让模型允许极小概率下同时充放电，求解速度会有提升，但要在结果检查时确认不会有实际问题。第二，在Yalmip里定义变量时用维度矩阵，不要用cell矩阵，后者的索引开销很大。第三，给Gurobi设置MIPGap，比如0.005，而不是默认的0，这样求解器可以早早返回一个可接受的精化解，能省下大量时间。
一个容易被忽略的点是：场景削减后的场景权重总和是否严格等于1。如果削减代码里有数值误差导致权重和不等于1，目标函数中的成本项就会整体失真，求解器会出现一些奇怪的解。我建议在调用optimize之前，加一行断言检查权重和与1的偏差。
6.3 热负荷与电负荷的时间尺度错配问题
热电联供微网里最隐蔽的问题，是电负荷和热负荷的波动特性完全不同。电负荷响应速度是秒级到分钟级，热负荷因为有建筑热惯性的缓冲，变化相对平缓，而且供热管网的传输延迟也让热负荷的“实际感受”滞后于热源出力。
如果优化模型用同样的时间分辨率去处理电和热，比如都用一小时一个时段，热负荷延迟效应就会被完全忽略。这样得到的蓄热罐放热策略，可能在真实系统中过早地把热量放完了，等热负荷真正升高的时候，蓄热罐里已经没有余量。解决这个问题有两种常见手段：一是把热负荷加上一段滞后模型，比如把当前时段的热负荷与上一时段的平均供热输出绑定；二是在热平衡约束中引入一个热惯性系数，模拟建筑允许的温度波动范围。这些措施会增加模型复杂度，但能明显提高调度结果在实际运行中的可行性。
我还有一个经验：如果条件允许，热负荷模型的时段可以设成电负荷的两倍，即热力部分的调度周期拉长到2小时，这样能自然缓解热惯性问题。当然这要看具体系统的实际情况，不是所有场景都适用。
除了上面这个热惯性问题，还有一类典型坑来自蓄热罐的离散状态。有些同学把蓄热罐的热损失系数设成常数，但实际运行时蓄热罐顶部和底部的温度分层不一样，热损失不是一个固定比例。做研究时用常数近似可以接受，但如果做工程应用，建议蓄热罐模型至少分两层温度节点，否则长期运行的累计误差会非常可观。
最后再分享一点我的实际体会
一个研究课题的价值，不完全体现在算法有多高深，更多在于能不能把问题背后的物理特性看清楚。源荷随机特征这件事，看起来只是数学建模，但真正理解了热电联供机组“以热定电”的运行边界，理解了蓄热罐在时间维度上的缓冲作用，再回头看那些随机场景、多场景优化，就都有了落地的抓手。
如果想把这个课题向下延伸，比较值得做的是把模型扩展到多微网协同，或者把运行优化和容量配置放到同一个框架里做两阶段优化。前者需要面对更复杂的荷侧不确定性，后者则要处理设备投资与运行成本的耦合。无论哪条路，Matlab加Yalmip加成熟求解器的这套技术栈都还能继续用，模型改起来也比较顺手。最后给一个小建议：拿到任何优化模型的时候，务必先花半小时检查约束的物理意义，再去盯求解器的设置，很多看起来是“求解失败”的问题，本质上都是约束里藏着一两个量纲或者索引错误。
