从零到一：基于MATLAB/YALMIP/CPLEX的电力系统机组组合优化实战

1. 初识电力系统机组组合优化

第一次接触电力系统机组组合优化问题时，我正为一个省级电网的调度项目头疼。当时手头有6台发电机组的运行数据，需要制定24小时的最优启停计划。这个看似简单的任务，实际上涉及到复杂的数学建模和优化计算。后来我发现，MATLAB+YALMIP+CPLEX这个黄金组合，能帮我们轻松搞定这类问题。

机组组合优化（Unit Commitment）的核心目标很简单：在满足电力需求的前提下，合理安排发电机组的启停和出力，使总运行成本最低。但实际操作中会遇到各种约束条件，比如机组最小运行时间、爬坡速率限制、网络潮流安全等。这就好比你要安排一个家庭聚餐，既要考虑每个人的时间安排（约束条件），又要尽量节省总体开销（目标函数）。

MATLAB作为技术计算领域的瑞士军刀，提供了强大的数值计算和可视化能力。而YALMIP就像是个聪明的翻译官，能把复杂的数学模型转换成CPLEX能理解的语言。CPLEX则是IBM开发的商用优化求解器，特别擅长处理混合整数规划问题。三者配合使用，就像组成了一个高效的问题解决团队。

2. 环境搭建与工具准备

2.1 软件安装指南

工欲善其事，必先利其器。首先需要安装MATLAB基础环境，我推荐使用R2020b或更新版本。安装完成后，还需要获取两个关键工具箱：

1. YALMIP工具箱：这是一个免费的MATLAB建模语言，可以从官网直接下载
2. CPLEX优化器：IBM提供的商业软件，需要申请学术许可证或购买商业版

安装YALMIP时，只需将下载的压缩包解压到MATLAB工具箱目录，然后在MATLAB命令行运行：

addpath(genpath('yalmip目录路径')) savepath

CPLEX的安装稍微复杂些。Windows用户可以直接运行安装程序，Linux用户则需要配置环境变量。安装完成后，在MATLAB中测试是否成功：

try cplex = Cplex('test'); disp('CPLEX安装成功！'); catch disp('请检查CPLEX安装和路径配置'); end

2.2 数据准备技巧

在实际项目中，我习惯将数据整理成三个Excel工作表：

• 机组参数：包括最小/最大出力、成本系数、爬坡速率等
• 网络参数：线路阻抗、容量限制等
• 负荷曲线：24小时系统总需求

读取数据的MATLAB代码可以这样写：

% 读取机组参数 gen_data = xlsread('UC_data.xlsx', '机组参数'); min_output = gen_data(:,3); % 最小出力(MW) max_output = gen_data(:,4); % 最大出力(MW) cost_a = gen_data(:,5); % 成本系数a cost_b = gen_data(:,6); % 成本系数b cost_c = gen_data(:,7); % 成本系数c % 读取24小时负荷数据 load_data = xlsread('UC_data.xlsx', '负荷曲线'); hourly_load = load_data(end,2:25); % 取最后一行作为总负荷

3. 数学模型构建详解

3.1 问题分析与变量定义

机组组合问题需要考虑两类决策变量：

1. 二元变量u(i,t)：表示机组i在时段t的状态（0停机，1运行）
2. 连续变量p(i,t)：表示机组i在时段t的出力(MW)

在MATLAB中使用YALMIP定义这些变量：

num_gen = 6; % 机组数量 num_hours = 24; % 时段数量 % 定义决策变量 u = binvar(num_gen, num_hours, 'full'); % 机组状态 p = sdpvar(num_gen, num_hours, 'full'); % 机组出力

3.2 目标函数构建

目标是最小化总成本，包括：

1. 发电成本：通常用二次函数表示
2. 启停成本：机组启动和停机产生的固定成本

对应的数学表达式为：

min Σ [c_i(p_it) + H_i*(u_it - u_i(t-1))^+ + J_i*(u_i(t-1) - u_it)^+]

其中c_i(p_it) = a_ip_it² + b_ip_it + c_i

YALMIP实现代码：

% 计算发电成本 generation_cost = 0; for i = 1:num_gen for t = 1:num_hours generation_cost = generation_cost + cost_a(i)*p(i,t)^2 + cost_b(i)*p(i,t) + cost_c(i); end end % 计算启停成本 startup_cost = 0; shutdown_cost = 0; for i = 1:num_gen for t = 2:num_hours startup_cost = startup_cost + H(i) * max(u(i,t) - u(i,t-1), 0); shutdown_cost = shutdown_cost + J(i) * max(u(i,t-1) - u(i,t), 0); end end total_cost = generation_cost + startup_cost + shutdown_cost;

3.3 约束条件设置

机组组合问题需要满足多种约束条件：

1. 功率平衡约束：

constraints = []; for t = 1:num_hours constraints = [constraints, sum(p(:,t)) == hourly_load(t)]; end

2. 机组出力限制：

for i = 1:num_gen for t = 1:num_hours constraints = [constraints, p(i,t) >= min_output(i) * u(i,t), p(i,t) <= max_output(i) * u(i,t)]; end end

3. 爬坡速率限制：

for i = 1:num_gen for t = 2:num_hours constraints = [constraints, p(i,t) - p(i,t-1) <= ramp_up(i) * u(i,t-1), p(i,t-1) - p(i,t) <= ramp_down(i) * u(i,t)]; end end

4. 模型求解与结果分析

4.1 求解器配置与参数调优

CPLEX求解器有很多参数可以调整，对于机组组合问题，我通常会设置：

options = sdpsettings('solver','cplex',... 'cplex.timelimit',3600,... % 时间限制1小时 'cplex.mip.tolerances.mipgap',0.01,... % 允许1%的gap 'cplex.mip.strategy.heuristicfreq',100,... % 启发式频率 'verbose',1); % 显示求解过程

第一次运行时，建议先用小规模测试案例，设置较短的时间限制，快速验证模型是否正确。我曾经犯过一个错误，把爬坡速率的单位弄错了（MW/h vs MW/15min），导致求解结果完全不合理。

4.2 结果可视化技巧

求解完成后，我们可以提取并可视化结果：

% 提取最优解 u_opt = value(u); p_opt = value(p); total_cost_opt = value(total_cost); % 绘制机组启停状态 figure; imagesc(u_opt); colormap([1 1 1; 0 0.5 0]); % 白色表示停机，绿色表示运行 xlabel('时段'); ylabel('机组'); title('机组启停计划'); % 绘制机组出力曲线 figure; hold on; for i = 1:num_gen plot(1:num_hours, p_opt(i,:), 'LineWidth',2); end plot(1:num_hours, hourly_load, 'k--', 'LineWidth',2); xlabel('时段'); ylabel('出力(MW)'); legend('机组1','机组2','机组3','机组4','机组5','机组6','系统负荷');

4.3 常见问题排查

在实际应用中，可能会遇到各种问题：

1. 模型不可行：检查约束条件是否冲突，特别是爬坡速率和最小运行时间约束
2. 求解时间过长：尝试线性化目标函数，或增加MIP gap容忍度
3. 结果不合理：检查单位是否一致，参数输入是否正确

我曾经遇到一个案例，因为忽略了机组最小运行时间约束，导致求解器给出了频繁启停的计划，这在现实中是完全不可行的。后来添加了以下约束后问题得到解决：

% 最小运行时间约束 for i = 1:num_gen for t = 2:num_hours if u_opt(i,t) > u_opt(i,t-1) % 如果机组启动 constraints = [constraints, sum(u_opt(i,t:min(t+min_up_time(i)-1,num_hours))) >= min_up_time(i)*(u_opt(i,t)-u_opt(i,t-1))]; end end end

5. 实战技巧与性能优化

5.1 模型线性化方法

二次成本函数会导致求解时间大幅增加。我们可以采用分段线性化来近似：

% 分段线性化参数 num_segments = 4; % 分段数 segment_length = (max_output - min_output)/num_segments; % 定义分段出力变量 p_seg = sdpvar(num_gen, num_hours, num_segments, 'full'); % 修改目标函数 linear_cost = 0; for i = 1:num_gen for t = 1:num_hours for s = 1:num_segments linear_cost = linear_cost + marginal_cost(i,s) * p_seg(i,t,s); end end end % 添加分段约束 for i = 1:num_gen for t = 1:num_hours constraints = [constraints, p(i,t) == min_output(i)*u(i,t) + sum(p_seg(i,t,:)), p_seg(i,t,:) >= 0, p_seg(i,t,:) <= segment_length(i)]; end end

5.2 并行计算加速

对于大规模问题，可以使用MATLAB的并行计算功能：

% 开启并行池 if isempty(gcp('nocreate')) parpool('local',4); % 使用4个worker end % 在求解设置中启用并行 options.cplex.parallel = 1; options.cplex.threads = 4;

5.3 实用调试技巧

调试优化模型时，我常用的几个技巧：

1. 先求解松弛问题（忽略整数约束）：

relaxed_constraints = constraints; for i = 1:num_gen for t = 1:num_hours relaxed_constraints = [relaxed_constraints, 0 <= u(i,t) <= 1]; end end optimize(relaxed_constraints, total_cost, options);

2. 固定部分变量，分阶段求解：

% 先求解前12小时 partial_constraints = constraints(:,1:12); optimize(partial_constraints, total_cost_partial, options); % 固定前12小时结果，再求解剩余时段 fixed_constraints = [constraints(:,13:24), u(:,1:12) == value(u(:,1:12))]; optimize(fixed_constraints, total_cost_remaining, options);

3. 使用回调函数监控求解过程：

function stop = mycallback(~,~) best = cplex.getBestObjValue(); current = cplex.getIncumbentObjValue(); fprintf('当前解: %.2f, 最优边界: %.2f, Gap: %.2f%%\n',... current, best, 100*(current-best)/current); stop = false; end options.cplex.mip.callback = @mycallback;