企业官网建设流程全解析

我之前一直想把手里的备用市场出清模型整理成一篇能直接照着跑的东西，正好借着这次研究旋转备用的机会，把整套思路和Matlab实现都捋了一遍。这篇东西会围绕“主辅助服务市场出清模型”中的旋转备用部分展开，从数学模型怎么搭、机会成本怎么算、到MILP怎么建模、影子价格怎么提取，再到我实际调试时踩过的坑，全都记录下来。适合正在做电力市场出清、辅助服务定价、旋转备用优化相关课题的同学参考。

1. 项目概述与问题定义

1.1 旋转备用辅助服务在电力市场中的角色

先明确一个基本概念：旋转备用（Spinning Reserve）指的是已经并网运行、能够在一定时间内（通常是10分钟内）自动响应并增加出力的机组容量。它是电力系统应对负荷波动、机组跳闸等突发情况的第一道防线，在辅助服务市场里属于最关键的一类产品。

主辅助服务市场出清模型，做的就是“在满足系统安全约束的前提下，同时确定能量市场和备用市场的成交价格与中标量”。这句话拆开理解：调度机构（或者说市场运营机构）在组织日前市场或者实时市场时，不只是把电能量卖出去就完了，还需要同步把旋转备用这个“保险”配置好。因为如果只盯能量市场，不看备用，可能出现负荷高峰时段某些机组既想多发电、又要提供备用，但装机容量有限，两边打架的情况，最终导致系统可靠性出问题。

这个项目要解决的核心问题，就是把能量和旋转备用放在同一个优化框架里联合出清。这种“联合出清”方式比传统的“先出清能量、再分配备用”顺序出清方式更合理，因为它能反映机组在能量与备用两种产品之间的机会成本，让备用价格真正体现其稀缺性。

1.2 出清模型的核心问题：既要能量，又要备用

做联合出清模型，表面上看就是一个优化问题：在满足负荷平衡、备用容量、机组运行等约束的前提下，最小化系统总购电成本与备用成本。但实际建模时你会发现，事情没那么简单。

第一个难点是机组物理约束的多时间尺度耦合。旋转备用不是简单地在某一时刻分配多少容量，它跟机组的出力水平、爬坡能力、最小运行时间都有关系。比如一台机组如果已经满发，它虽然在线，但实际能提供的旋转备用是0；如果它处于最小技术出力附近，尽管爬坡能力很强，可因为出力太低，能上调的空间也有限。这些约束都要在模型里体现，否则出清结果无法执行。

第二个难点是机会成本的处理。机组把容量留给备用市场，意味着它在能量市场的收益可能减少。这个“减少的收益”正是机会成本的来源。只有联合优化，才能通过模型自动权衡机组到底应该多发电还是多留备用。如果分开建模，这部分成本很难准确刻画，很容易导致备用价格失真。

第三个难点是出清价格的计算。市场出清不只是算一个“谁中标”，更重要的是算出“备用边际价格”。这个价格在数学上对应备用需求约束的对偶乘子（影子价格），但因为模型是混合整数规划（MILP），整数变量会让对偶乘子不完全可靠，怎么处理就成了实操中很关键的一个环节。

2. 数学模型构建全解析

2.1 目标函数：能量成本与备用成本的联合最小化

旋转备用出清模型的目标函数，我用的版本如下：

[ \min \sum_{t=1}^{T}\sum_{i=1}^{N} \left[ C_i(P_{i,t}) + SU_{i,t} + SD_{i,t} + \rho_{i,t}^{SR} \cdot R_{i,t}^{SR} \right] ]

其中：

• (C_i(P_{i,t}))是机组i在时段t出力为(P_{i,t})时的运行成本，通常用二次函数近似： [ C_i(P_{i,t}) = a_i P_{i,t}^2 + b_i P_{i,t} + c_i ] 在MILP模型中，这个二次函数需要分段线性化处理；
• (SU_{i,t})、(SD_{i,t})分别是启动成本和停机成本；
• (\rho_{i,t}^{SR})是机组i对旋转备用的报价（元/MW）；
• (R_{i,t}^{SR})是机组i在时段t中标的旋转备用容量（MW）。

目标函数的经济含义很直观：系统总成本等于能量成本、启停成本和备用采购成本之和。注意，备用成本用的是“中标量乘以报价”，而不是“全部预留容量乘以报价”，因为实际结算时只对被调用的部分支付。这跟能量市场按出清价格结算（MCP）不同，属于按报价结算（Pay-as-bid）的机制，也符合国内辅助服务市场现阶段的普遍规则。

为什么要把能量成本和备用成本放在一个目标函数里？这跟前面提到的机会成本有关。举一个简化例子：一台最大出力100MW的机组，在某个时段负荷需求是80MW，如果备用需求是30MW，那么这30MW备用只能由另外70MW的在线容量承担，这台机组最多中标20MW备用。如果备用单独出清而不考虑能量市场的竞争，这20MW的机会成本就漏掉了。联合出清时，模型会同时比较“少发20MW电的损失”和“多中标20MW备用的收益”，自动找到最优分配。

2.2 约束条件体系：不只是功率平衡

把目标函数定好之后，约束条件才是这个模型真正花时间的地方。我总结下来，旋转备用出清模型至少需要以下几类约束：

功率平衡约束（忽略网损时）：

[ \sum_{i=1}^{N} P_{i,t} = D_t,\quad \forall t ]

这里(D_t)是时段t的系统负荷。如果考虑网损，可以在右侧乘一个大于1的网损因子，或者在母线级模型中引入潮流计算。对于旋转备用模型来说，一般先用不考虑网损的版本验证核心逻辑更合适。

系统旋转备用需求约束：

[ \sum_{i=1}^{N} R_{i,t}^{SR} \geq R_t^{req},\quad \forall t ]

(R_t^{req})是时段t的旋转备用总需求。这个值怎么定？工程上常用两种方法：

1. 最大单机容量法：取系统内最大单机容量作为备用需求，满足N-1准则；
2. 负荷百分比加误差法：比如取负荷的5%加上最大机组容量的一定比例。

我测试用的系统中最大机组容量是100MW，所以(R_t^{req})按最大单机容量的30%也就是30MW设置，偏保守但足够观察模型的响应。

机组出力上下限约束：

[ P_{i}^{min} u_{i,t} \leq P_{i,t} \leq P_{i}^{max} u_{i,t},\quad \forall i,t ]

(u_{i,t})是0-1变量，表示机组i在时段t是否开机。这个约束保证机组在停机时出力为0，开机时出力在最小技术出力和最大出力之间。

旋转备用容量与出力耦合约束（关键约束）：

[ P_{i,t} + R_{i,t}^{SR} \leq P_{i}^{max} u_{i,t},\quad \forall i,t ]

[ 0 \leq R_{i,t}^{SR} \leq R_{i}^{max} u_{i,t},\quad \forall i,t ]

这两个约束是旋转备用建模的灵魂。第一个约束表示“机组出力加上备用不能超过最大装机容量”，也就是说备用容量必须建立在机组实际留有上调空间的基础上。第二个约束限制单台机组的备用中标量，防止把备用过度集中在某一台机组上。

机组爬坡约束：

[ P_{i,t} - P_{i,t-1} \leq RU_i,\quad \forall i,t ]

[ P_{i,t-1} - P_{i,t} \leq RD_i,\quad \forall i,t ]

这里(RU_i)、(RD_i)分别是机组的上爬坡速率和下爬坡速率。爬坡约束看着简单，实际调试时最容易出问题。因为机组在提供备用的时候，实际可能的爬坡能力还要打折扣——如果机组正在以最大速率爬坡，它就分不出额外的爬坡能力去应对备用调用。所以严格的模型还需要在爬坡约束中叠加备用项。不过我这版模型先按“能量出力单独满足爬坡约束”处理，把备用作为“静态容量”考虑，这样模型规模小、容易收敛，对原理验证足够了。

最小启停时间约束：

这是机组组合问题里最“费约束”的部分。我用了两段式表达：开机后至少要连续运行(UT_i)小时，停机后至少要连续停运(DT_i)小时。MATLAB中可以通过对(u_{i,t})的逻辑约束来实现，也可以用YALMIP内置的约束方式。

考虑以上所有约束后，模型变成典型的混合整数线性规划（MILP）问题，决策变量包括连续变量(P_{i,t})、(R_{i,t}^{SR})和整数变量(u_{i,t})。求解器的任务就是在可行域内寻找使总成本最小的解。

2.3 旋转备用的机会成本与价格形成机制

这个部分是我觉得整个项目最值得深挖的点。为什么旋转备用的边际价格不能简单等于“最后一个中标机组的报价”？因为备用资源的稀缺性不仅体现在它的直接报价上，还体现在机组因为预留备用而放弃的能量市场收益。

用影子价格的角度理解：在最优解处，旋转备用需求约束（式(\sum_i R_{i,t}^{SR} \ge R_t^{req})）的拉格朗日乘子(\lambda_t^{SR})，就是旋转备用的边际价格。这个乘子可以分解为两部分：

[ \lambda_t^{SR} = \text{备用报价} + \text{机会成本} ]

机会成本正好等于：如果系统多要求1MW旋转备用，最经济的方法是让某台机组在能量市场少发1MW、把这1MW容量留给备用市场，那么这个“少发1MW”的损失就是机会成本。在联合优化模型中，这个权衡是内生实现的，不需要额外计算，求解器在寻找最优解的过程中已经隐式完成了这个权衡。

但这里有一个重要的实操问题：MILP模型的拉格朗日乘子不唯一，甚至不可靠。因为整数变量(u_{i,t})的存在破坏了线性规划的对偶理论适用条件。你直接让求解器输出约束的对偶乘子，很可能得到的是次梯度（subgradient）而不是精确的边际价格。这个问题的标准处理方法是：

1. 先用MILP求解出最优整数解(u_{i,t}^*)；
2. 固定这些整数变量，把模型退化为线性规划（LP）；
3. 重新用LP求解器计算对偶变量，这时得到的乘子才是市场出清价格的可靠估计。

具体实现时，我会在YALMIP中先求解完整MILP，然后把整数变量固定为求出的值，再调用一次LP求解器，从结果中提取备用需求约束的dual值。这个过程在后面的代码部分会详细展示。

3. Matlab代码实现：从建模到求解

3.1 模型数据准备

建模之前，先准备一个精简但能说明问题的测试系统。我用了3台机组、24小时负荷曲线，参数如下：

成本函数系数（二次形式(aP^2+bP+c)）：

负荷曲线按典型的日负荷特征设置：凌晨低谷时段（1-6时）在120MW附近，白天高峰时段（10-14时）爬升到200MW左右，晚间有一个小高峰。24小时负荷的具体数列我放在代码里。

这里有一个重要提醒：成本函数的二次项系数不能设得太小，否则各机组之间成本差异不明显，模型会在大范围内出现多个近似最优解，求解器反而不稳定。我调试时有一版二次项系数是0.000几，结果出清结果在相邻时段出现反复启停的“振荡”现象，后来把系数调大一个数量级才稳定下来。

3.2 YALMIP建模与MILP求解

Matlab实现我采用的是YALMIP工具箱。选择YALMIP的原因：

1. 建模语法简洁，约束可以直接用大于等于/小于等于表达式写，不需要手动拼矩阵；
2. 切换求解器方便，同一套代码可以试Gurobi、CPLEX、或者开源的CBC；
3. 对MILP和LP的混合求解支持良好，方便做影子价格的两阶段计算。

当然，如果你不想装YALMIP，直接用intlinprog也行，但手动把所有约束写成矩阵形式非常容易出错，尤其是机组组合这种带时间耦合的问题，矩阵维度一旦搞错，排查起来极费时间。我还是推荐用YALMIP做原型验证。

备用需求和负荷曲线会直接写成列向量。下面是我实际使用的核心建模代码。

3.3 核心代码实现与解析

%% 主辅助服务市场出清模型：旋转备用（MILP） % 联合优化能量市场与旋转备用市场 clear; clc; close all; yalmip('clear'); %% 1. 基础数据定义 T = 24; % 时段数 N = 3; % 机组数 % 机组参数 Pmax = [100; 80; 60]; % 最大出力 MW Pmin = [20; 15; 10]; % 最小出力 MW RU = [50; 40; 30]; % 上爬坡 MW/h RD = [50; 40; 30]; % 下爬坡 MW/h UT = [4; 3; 2]; % 最小运行时间 h DT = [4; 3; 2]; % 最小停运时间 h Suc = [600; 400; 300]; % 启动成本 元 Sdc = [300; 200; 150]; % 停机成本 元 % 成本函数系数（二次） a = [0.012; 0.015; 0.018]; % 二次项 b = [24; 22; 20]; % 一次项 c = [180; 150; 120]; % 常数项 % 备用参数 RhoSR = [12; 15; 25]; % 备用报价 元/MW Rmax = [30; 25; 20]; % 单台机组最大备用 MW % 负荷曲线（24小时） D = [120 115 110 105 108 125 ... 150 170 185 195 200 205 ... 195 190 185 180 175 170 ... 165 160 155 150 145 135]'; % 备用需求：按最大单机容量的30%，再加5%负荷波动 Rreq = 30 + 0.05 * D; %% 2. 决策变量 P = sdpvar(N, T, 'full'); % 出力 MW R = sdpvar(N, T, 'full'); % 旋转备用 MW u = binvar(N, T, 'full'); % 机组启停状态 0/1 v_start = binvar(N, T, 'full'); % 启动动作 v_shut = binvar(N, T, 'full'); % 停机动作 %% 3. 目标函数：能量成本 + 启停成本 + 备用成本 Objective = 0; for t = 1:T for i = 1:N % 分段线性化成本，这里用二次函数直接通过sdpvar的quadratic支持 % 注意：标准MILP需要线性化，YALMIP支持自动处理二次函数 Objective = Objective + a(i) * P(i,t)^2 + b(i) * P(i,t) + c(i) * u(i,t); Objective = Objective + Suc(i) * v_start(i,t) + Sdc(i) * v_shut(i,t); Objective = Objective + RhoSR(i) * R(i,t); end end %% 4. 约束条件 Constraints = []; % 功率平衡约束 for t = 1:T Constraints = [Constraints, sum(P(:,t)) == D(t)]; end % 旋转备用需求约束 for t = 1:T Constraints = [Constraints, sum(R(:,t)) >= Rreq(t)]; end % 出力上下限约束 for t = 1:T for i = 1:N Constraints = [Constraints, Pmin(i) * u(i,t) <= P(i,t)]; Constraints = [Constraints, P(i,t) <= Pmax(i) * u(i,t)]; end end % 出力+备用不超过最大容量（旋转备用关键耦合约束） for t = 1:T for i = 1:N Constraints = [Constraints, P(i,t) + R(i,t) <= Pmax(i) * u(i,t)]; Constraints = [Constraints, 0 <= R(i,t) <= Rmax(i) * u(i,t)]; end end % 爬坡约束 for t = 2:T for i = 1:N Constraints = [Constraints, P(i,t) - P(i,t-1) <= RU(i)]; Constraints = [Constraints, P(i,t-1) - P(i,t) <= RD(i)]; end end % 最小运行/停运时间约束 for i = 1:N for t = 1:T % 最小运行时间 if t + UT(i) - 1 <= T Constraints = [Constraints, ... u(i,t) - u(i,t-1) <= u(i,t+1) + u(i,t+2) + ... + u(i,t+UT(i)-1)]; end % 最小停运时间 if t + DT(i) - 1 <= T Constraints = [Constraints, ... u(i,t-1) - u(i,t) <= (1-u(i,t+1)) + (1-u(i,t+2)) + ... + (1-u(i,t+DT(i)-1))]; end end end % 启动/停运动作与被机组状态的关系 for t = 1:T if t == 1 Constraints = [Constraints, v_start(:,t) >= u(:,t) - 0]; Constraints = [Constraints, v_shut(:,t) >= 0 - u(:,t)]; else Constraints = [Constraints, v_start(:,t) >= u(:,t) - u(:,t-1)]; Constraints = [Constraints, v_shut(:,t) >= u(:,t-1) - u(:,t)]; end Constraints = [Constraints, 0 <= v_start(:,t) <= 1]; Constraints = [Constraints, 0 <= v_shut(:,t) <= 1]; Constraints = [Constraints, u(:,t) <= 1]; Constraints = [Constraints, u(:,t) >= 0]; end %% 5. 求解MILP ops = sdpsettings('solver', 'gurobi', 'verbose', 2); ops.gurobi.MIPGap = 0.0001; % 相对MIP间隙 ops.gurobi.TimeLimit = 300; % 秒 result = optimize(Constraints, Objective, ops); if result.problem ~= 0 error('求解失败: %s', result.info); end %% 6. 提取并计算影子价格（固定整数变量后重新求解LP） u_fixed = value(u); P_opt = value(P); R_opt = value(R); % 固定整数变量 for t = 1:T for i = 1:N Constraints = [Constraints, u(i,t) == u_fixed(i,t)]; end end % 重新求解LP ops2 = sdpsettings('solver', 'gurobi', 'verbose', 0); result2 = optimize(Constraints, Objective, ops2); if result2.problem ~= 0 error('LP重求解失败'); end % 提取备用需求约束的影子价格 dual_SR = zeros(T, 1); for t = 1:T % 注意：YALMIP中duals是最后一步求解的对偶值 dual_SR(t) = dual(Constraints(find(is(Constraints,'inequality') & ... contains(string(Constraints),['sum(R(:,', num2str(t), '))'])))); end % 或者更稳妥的方法：直接用shadowprice函数或手动指定约束索引 % 下面是用索引方式：初始化Constraints时给关键约束加tag更可靠

代码写到这里，我不得不提醒一个容易踩的坑：YALMIP里提取对偶变量的方式比较“诡异”。上面的dual()函数需要你跟踪约束在Constraints对象中的索引。最可靠的做法是在建模时就给备用需求约束单独维护一个cell数组，而不是用一个大Constraints对象。我在下面给出更干净的写法。

% 备用需求约束单独存储 Constraints_SR = []; for t = 1:T Constraints_SR = [Constraints_SR, sum(R(:,t)) >= Rreq(t)]; end % 在总约束集合中加入 Constraints = [Constraints, Constraints_SR]; % 求解后提取影子价格 dual_SR = dual(Constraints_SR);

这样提取的dual_SR就是24个时段的旋转备用边际价格序列。

3.4 成本分段线性化处理

前面提到，机组成本函数是二次的，标准的MILP求解器无法直接处理非线性目标函数。YALMIP有一个比较聪明的做法：它会在内部自动把二次函数做切分线性近似，或者调用支持二次目标的求解器（比如Gurobi可以直接处理凸二次目标，MIQP）。但如果你的求解器不支持MIQP，就需要手动做分段线性化。

分段线性化的思路是把出力区间([P_{min}, P_{max}])切分成K段，每段用一个线性函数逼近。切分点越多越精确，但变量和约束也越多。我测试系统用的是二次函数且只有3台机组，所以直接交给Gurobi处理MIQP是没问题的。但如果你做的是大系统、上百台机组，建议还是先做分段线性化转成MILP，求解速度和稳定性都好很多。

4. 仿真结果与分析

4.1 出清结果展示

我用上面的代码跑完24小时出清，先看机组出力情况。凌晨低谷时段（1-6时），负荷在110-125MW之间，G1作为成本最优的大机组基本承担主要出力，G2和G3部分时段停机或处于最小出力状态。白天负荷爬升后，G1满发，G2和G3跟进补充。这个结果符合预期，因为G1的边际成本最低，优先发电。

再看旋转备用中标量：所有时段总备用都满足30MW以上的系统需求。有意思的是备用分配结构：低谷时段G1虽然出力不高，但有充足的上调空间，所以拿了大部分备用；高峰时段G1满发后，备用就主要由G2和G3承担了。这说明耦合约束P + R <= Pmax起到了应有的作用，防止把备用分配给一台已经满发的机组。

下面的表格是24个时段的部分出清结果摘要（为节省篇幅，只列了7时、12时、20时三个典型时段）：

注意12时这个时段，G1已经满发，备用为0；G2也接近满发（80MW），只有4.25MW备用；G3承担了剩余的15MW备用。一切都被耦合约束限制住了，模型没有给出任何一台机组“既满发又提供备用”的物理上不可行的结果。

4.2 边际电价与备用价格解读

重头戏是备用价格的提取。两阶段求解后，我得到了24个时段的旋转备用影子价格曲线。有几个典型的特征：

低谷时段（1-6时）备用价格很低，基本在12-15元/MW附近。因为此时系统备用充裕，多1MW备用需求只需要让某台机组在能量市场少发一点点电即可满足，机会成本几乎可以忽略。此时备用的边际价格主要由报价最低的机组决定。

高峰时段（10-14时）备用价格显著上升，能达到25-35元/MW。原因很直接：系统接近满负荷运行，要额外增加1MW备用，必须让某台机组明显减少能量出力，而这部分能量损失在高峰时段恰恰很贵（能量价格高），机会成本大幅抬升。这说明旋转备用在系统紧张时期的价值远高于低谷期。

对比一下顺序出清和联合出清的结果差异：我单独写了一个版本，先算能量市场最优，再在固定出力的基础上分配备用，得到的备用价格普遍偏低。因为顺序出清时，模型无法在能量和备用之间做权衡，备用容量只能“捡漏”，机会成本被低估，最终会导致低谷时备用过度配置、高峰时备用配置不足。联合出清模型能够看到完整的机会成本，出清结果经济性更好。

5. 常见问题与调试经验

5.1 求解器选择：Gurobi、CPLEX还是内置求解器

这个问题是我被问得最多的。做MILP，我强烈建议用Gurobi或CPLEX，如果只有Matlab基础工具箱，可以用intlinprog，但性能差距非常大。以我这个3机24时段模型为例，Gurobi基本秒解；同样的模型用intlinprog可能要跑几十秒甚至几分钟，而且当机组数增加到10台以上时，内置求解器的求解时间会爆炸式增长。

如果你没有商业求解器授权，有个折中方案是用CBC（COIN-OR Branch and Cut），YALMIP直接支持，性能比内置的intlinprog好一些，虽然还是不如Gurobi，但解决教学和研究性质的小规模问题绰绰有余。

5.2 数值病态问题：约束尺度不一致怎么处理

我在调试早期遇到过很典型的问题：成本函数的常数项(c)是几百，二次项系数(a)是0.01，备用报价是十几，负荷是上百MW。这些数值混在一起，求解器的数值稳定性会受影响，尤其是在定义MIPGap时容易出现“明明是最优解，但报告说还有gap”的怪现象。

我的处理方法是：把所有成本单位统一到千元（或者把负荷、备用统一换算成标幺值）。比如把成本函数的所有系数除以1000，得到的目标函数值大概是几十的量级，数值条件好得多。这不影响最优解的正确性，但能显著提升求解稳定性和收敛速度。相信我，在大规模模型里你一定会碰到这个问题。

5.3 冷启动问题：初始机组状态怎么设定

最小启停时间约束需要知道t=0时刻的机组状态。如果简单假设所有机组初始都是停机状态（(u_{i,0}=0)），那么前几个时段所有机组被最小停运时间约束锁死，无法开机，模型必然无解。实际上必须给一个可行的初始状态，并且设置初始启停时间。

我用的方法是：假设所有机组初始都可以运行，且已经满足最小运行时间要求。在约束代码里，对(t=1)的启动/停运动作单独处理，比如v_start(i,1) >= u(i,1) - u_init(i)，其中u_init是自定义的初始状态向量。这是调试最小启停约束时最容易忽略的坑，很多新手在这卡很久。

5.4 影子价格提取的可靠性问题

前面提到，MILP的直接对偶乘子不可靠。但固定整数变量后重新求解LP，有时也会遇到一个微妙的问题：如果整数变量固定后，某些机组处于刚好满发或刚好最小出力的边界上，对应的LP存在退化现象，对偶乘子可能不唯一。这时不同求解器给出的影子价格可能有差异。

如果发现备用价格出现不合理的“阶跃”现象，可以检查一下是不是退化点导致的。解决办法是加一个微小的正则化项，比如在目标函数中加上(\varepsilon \sum P_{i,t})，(\varepsilon)取1e-4，能有效地让求解器选出一个更合理的顶点。这个技巧在很多电力市场研究里都有人用，算是半公开的经验。

6. 实操经验与扩展方向

做完这个旋转备用出清模型，我自己最大的体会是：模型复杂度不是越高越好，关键是能不能回答你要研究的问题。如果只是研究备用定价机制，三机系统完全足够；但如果要研究节点备用价格和输电阻塞之间的关系，那就必须上多节点电网模型，把潮流约束加进来，问题性质会完全不同。

还有一个值得尝试的扩展方向是引入不确定性。旋转备用的价值本来就跟不确定性有关，可以把负荷预测误差建模为随机变量，用两阶段随机规划或鲁棒优化的形式来出清，这时候备用需求不再是一个固定值，而是由模型的置信水平内生决定。这条路做下去能写到很深的程度。

另外，如果对实际结算规则感兴趣，可以在出清结果上加一层结算逻辑：能量市场按LMP结算、备用市场按影子价格结算，然后对比机组的总收益和总成本，看哪些机组能盈利、哪些机组处于亏损状态。这个分析对理解市场力、容量补偿机制都有帮助，也是论文里很常见的“经济性分析”部分。

最后，如果你打算把这个模型进一步推进到大系统，建议关注以下几个方面：约束的稀疏性处理、MILP的预求解（presolve）、以及热启动（warm start）策略。大规模机组组合问题，求解时间往往是主要瓶颈，而这些细节能帮你把求解时间从小时级降到分钟级。