光伏出力曲线K-means聚类：MATLAB完整流程与实战踩坑

做光伏功率预测或者储能容量配置的人应该都见过这种数据同一个电站晴天的出力曲线是一根干净的单峰多云天是锯齿状上下跳动阴雨天直接贴着零轴。用MATLAB对这些历史曲线做K-means聚类把它们分成几类典型形态是光伏曲线聚类研究里最常见、也最实用的一步。下面是我在实际项目里跑通的一套完整流程从数据清洗、归一化、K值选择到结果验证和踩坑记录都会讲到代码可以直接复制修改适合正在做新能源数据分析、功率预测或者并网规划的研究生和工程师参考。1. 为什么光伏曲线先聚类再分析形态差异带来的连锁问题1.1 三种典型形态对下游分析的直接影响光伏出力本质上由气象条件驱动而气象条件的变化不是连续的渐变而是相对稳定的几种模态。晴天时出力曲线近似光滑单峰从早上爬起来到中午前后到达峰值再落下去多云天气里云层反复遮挡曲线出现大量锯齿状波动阴雨天整体出力很低曲线贴着横轴缓慢起伏。除此之外季节、纬度、空气污染程度还会改变峰值的高度和出现时间。如果不去管这些差异直接把几千天历史曲线揉在一起算平均值结果就是一条两头都不靠的平庸曲线。做功率预测时模型用这种均值曲线打底晴天严重低估、阴天严重高估做储能配置时如果用全年的平均出力来设计系统碰到连续阴雨天很容易把容量估计不足做并网消纳分析时不同天气模式下光伏出力对电网调峰的压力也完全不同。这些场景都需要把典型日先分出来再针对每种典型日分别制定策略。这正是聚类的价值把历史上形态相似的日期归为一组每组得到一条可以代表该组特性的典型出力曲线同时给出每个组出现的概率。这个过程在电力系统研究里也叫典型场景提取或者场景缩减是预测、调度、规划类课题里的常见前置步骤。1.2 K-means为什么适合这类曲线聚类K-means的原理不复杂先随机选K个初始中心把每条曲线分给距离最近的中心然后重新计算每组的中心反复迭代直到中心不再变化。优化的目标函数是簇内平方和WCSS也就是所有曲线到各自所属中心的欧氏距离平方之和。选择K-means的根本原因是它和这类数据的问题结构匹配。光伏曲线经过时间对齐、缺失值处理后每条曲线都落在同一个长度、同一个采样间隔的向量空间里本身就是一个规整的矩阵。K-means在这种矩阵上运行速度很快几千条曲线几十次迭代就能收敛。而且MATLAB内置的kmeans函数接口稳定老版本到新版本几乎没有变化统计工具箱自带不需要额外装任何东西。当然K-means有它的短板用欧氏距离意味着它对曲线的绝对数值和尺度很敏感这一点后面归一化时要特别小心它假设簇是凸的、近似球形的如果某两类曲线形态非常纠缠效果就不好。但对光伏曲线聚类这个任务来说天气模态之间的差异足够大工程上完全够用。相比之下层次聚类虽然不需要预设K但样本量大时树状图很难解释高斯混合模型GMM能给出概率归属但计算更重作为第一步的探索分析反而显得复杂了。2. 预处理和数据约定聚类结果好坏的前置因素2.1 清洗规则哪些日期可以直接剔除很多人拿到数据第一件事就是跑聚类结果发现簇乱七八糟回过头来才意识到是脏数据惹的祸。光伏数据来自SCADA系统或者电站采集终端通讯中断、逆变器停机、电网限电、传感器故障都会造成记录异常常见的表现有三种全天出力恒为0、全天数值极低且毫无波动、或者某几个时段突然变成0再跳回正常。这些日子不是真实的天气模式混进聚类里只会产生孤立簇或噪声簇。我的做法是先对每条曲线做两次检查。第一次是粗过滤剔除全天最大出力低于阈值比如装机容量的0.5%的日期这类日子基本是停机或数据故障第二次是逐点检查统计全天出力小于等于0的点数占比如果超过三成说明这条曲线存在明显的数据断裂也直接剔除。需要注意的是如果原始数据没有截断到白天时段深夜的零值是真实的物理出力而不是异常不能用这个规则过滤这一点要和下文的时段截断配合起来。2.2 时间对齐与活跃时段截断同一电站的记录通常间隔是固定的比如15分钟或30分钟一个点但有些老旧电站会出现少采、多采的情况尤其是雷雨天气前后通讯多次中断会导致某些日期只有七八成的点。这时候必须做时间对齐把每条曲线重采样到统一的网格上。MATLAB里用retime同步时间表或者简单一点按照时间索引做interp1线性插值插值之前先确保每天的起止时间一致。对齐之后还有一件很多文章不会强调的事把活跃时段截出来。一天96个点15分钟间隔里大概有40多个点是夜间零出力这些零值点没有任何天气信息但它们在欧氏距离计算里贡献了一个巨大的、几乎恒定的分量。换句话说所有曲线的夜间部分都长得一样这部分会稀释白天形态的差异导致聚类结果对白天的细节不敏感。我自己第一次跑的时候用了完整96点平均轮廓系数只有0.31截取到7:00到18:00之后跳到0.46提升非常明显。如果电站坐标已知更讲究一点的做法是按照每天的实际日出日落时间动态截取比如只保留日照时段前后各扩展15分钟。这样能进一步去掉清晨和傍晚低出力段对距离计算的干扰。2.3 归一化路线形状优先还是量级优先归一化是这一步里最需要想清楚的决策它直接决定聚类分组的逻辑。K-means用的是欧氏距离曲线数值大不大、峰值高不高都会直接影响距离和分组。如果按全局最大值归一化也就是把所有曲线除以历史最大值或装机容量那么曲线之间的绝对出力水平被保留。这种做法的聚类结果会同时考虑形状和量级高功率日和低功率日自然会被分开哪怕它们形状几乎一样。如果按日最大值归一化每一天的曲线除以自身当天的最大值所有曲线都落在0到1区间聚类的逻辑变成只看形状同样是晴天夏天高功率的日子和冬天低功率的日子会因为形状相似被归到一类。怎么选完全取决于业务目标。做天气类型识别和预测模型的训练集划分时通常用日最大值归一化因为模型关心的是形态对应的气象条件而不是当天发了多少电做储能容量配置和消纳分析时量级本身就是核心信息建议用装机容量归一化保留发电量差异。这个选择没有标准答案但有标准检验方法看聚类结果和实际天气记录的对应关系哪个更符合业务解释就用哪个。2.4 要不要先做特征工程再聚类对曲线直接聚类是把全部形态信息都交给算法好处是自动化程度高不需要人为设计特征坏处是采样点越多噪声波动的影响越大距离计算被大量的局部锯齿主导。另一个常见路线是先压缩特征再聚类比如每天提取日发电量、峰值时刻、峰值大小、出力波动方差这几个标量组成特征向量后聚类。特征化的好处是结果更稳健对局部的分钟级波动不敏感坏处是丢了曲线局部的形状细节比如同样是上午多云下午晴的模式日发电量和峰值可能差不多但形状完全不同。我的建议是先用曲线直接聚类跑通流程把簇数和天气对应关系想清楚如果发现分类结果有难以解释的浑浊地带再尝试特征化聚类作为对比。两条路线在MATLAB里改动很小特征化只是把矩阵换成一个特征矩阵kmeans函数不需要任何变化。3. K值选择从肘部法则、轮廓系数到业务语义交叉验证3.1 肘部法则与轮廓系数的配合逻辑K-means需要用户指定K这是它最不方便的地方也是经验成分最大的地方。K太小晴天、多云、阴天被硬并到一起簇内差异巨大K太大每个簇都分得很碎失去典型场景的意义。常用的量化判断有两个。肘部法则看WCSS随K变化的曲线随着K增大WCSS单调下降但下降速度会在某个K之后明显变缓形成肘部这个点被认为是最佳K。原理不难理解簇数从2增加到3时如果恰好把原本混杂的一大簇劈成了两个明显不同的子簇簇内距离会大幅下降再继续增加K时劈开的只是本来就比较紧凑的簇降幅自然变小。轮廓系数是另一个角度对每个点计算a到同簇其他点的平均距离和b到最近其他簇所有点的平均距离系数s(b-a)/max(a,b)范围从-1到1。s接近1说明该点离自己簇很近、离最近的邻居簇很远聚类效果好接近0说明在边界上接近-1说明可能被分错了。全体点的平均轮廓系数就是常用的综合指标。单独用任何一个指标都有风险。肘部法在曲线平滑时拐点不明显轮廓系数则偶尔会在K偏大的时候给出虚高值。我习惯两个指标一起看再配合每个簇的样本量做判断。3.2 循环评估K的MATLAB代码这一段代码是我每次做聚类都会先跑的先把不同K的指标都算出来再决定取哪个K。数据用上一章处理好的data_norm已经是按日最大值归一化并截取过活跃时段的矩阵。% data_norm: nDays x nPoints逐日归一化的出力矩阵 K_range 2:8; wcss zeros(length(K_range), 1); avg_silh zeros(length(K_range), 1); min_cluster_size zeros(length(K_range), 1); for i 1:length(K_range) K K_range(i); rng(42); % 固定随机种子保证结果可复现 [idx_tmp, ~, sumd] kmeans(data_norm, K, Replicates, 20, MaxIter, 500); wcss(i) sum(sumd); avg_silh(i) mean(silhouette(data_norm, idx_tmp)); % 记录最小的簇样本数防止出现只有一个样本的空簇 counts histcounts(idx_tmp, 1:(K1)); min_cluster_size(i) min(counts); end % 输出对照表 T table(K_range, wcss, avg_silh, min_cluster_size, ... VariableNames, {K, WCSS, AvgSilhouette, MinClusterSize}); disp(T); % 画图 figure; subplot(1,2,1); plot(K_range, wcss, o-, LineWidth, 1.5); xlabel(K); ylabel(WCSS); title(肘部法则曲线); grid on; subplot(1,2,2); plot(K_range, avg_silh, o-, LineWidth, 1.5); xlabel(K); ylabel(平均轮廓系数); title(平均轮廓系数曲线); grid on;代码里有两个细节值得说明。第一是rng(42)固定随机种子这一点很关键后面会专门讲第二是每个K都跑了20次Replicates避免单次初始化陷入局部最优得到稳定的WCSS值。silhouette函数会逐个样本计算系数样本量几千条时速度稍慢但完全可接受。我拿真实数据跑过一次结果大致如下KWCSS平均轮廓系数最小簇样本数2152.30.38126398.70.4385472.10.4657558.90.429651.20.391K从3到4提升明显从4到5时平均轮廓系数不升反降而且最小簇样本数骤降到9天说明某些簇被过度拆分了。K6时甚至出现只有一个样本的孤立簇完全失去了典型场景的意义。所以我最终选了K4。3.3 4类还是5类业务语义说了算量化指标只是参考聚类的最终归属是业务。我做了几次之后对K的语义有了直观认知K3对应晴、多云、阴雨三个大气候K4通常把阴雨拆成阴天和雨天或者把晴天拆成夏季强晴天和春秋弱晴天K5以上就要小心了多出来的簇往往是某些特殊极端日的残留比如强台风过境后的异常低出力日或者某天下午短暂放晴的混合形态。判断一个K值是否合理可以把每个簇的质心曲线画出来再配合每个簇对应的真实日期去查天气记录。如果第4簇的日期几乎全部落在降雨日第3簇全部落在阴天日说明这个分组和气象语义吻合良好如果某个簇同时混着雨天和强晴天说明K多了或少了需要回去调整。4. MATLAB主程序从原始矩阵到聚类结果的一站式实现4.1 主程序代码与运行逻辑这一节给出完整的聚类主程序。假设你已经把原始数据整理成了data矩阵和dates日期数组data的每一行是一条日曲线dates是与行对应的日期列表。代码从读取数据开始到保存结果结束中间包含清洗、截取、归一化、聚类的全部环节。% 1. 数据读取与基础信息 load(pv_data.mat, data, dates); nDays size(data, 1); nPoints size(data, 2); % 假设15分钟采样一天96点只取7:00-18:00索引29到73 time_idx 29:73; data_day data(:, time_idx); time_axis (7:0.25:18); % 与time_idx对应的时间刻度 % 2. 异常日清洗 max_vals max(data_day, [], 2); % 剔除全天最大出力过低的日子以及零值点占比超过30%的日子 valid max_vals 0.005 sum(data_day 0, 2) length(time_idx) * 0.3; data_day data_day(valid, :); dates dates(valid, :); max_vals max_vals(valid); % 3. 归一化按日最大值归一到[0,1] data_norm data_day ./ max_vals; % 4. K值评估快速版只跑2到8 K_range 2:8; avg_silh zeros(length(K_range), 1); for i 1:length(K_range) rng(42); [idx_tmp, ~, ~] kmeans(data_norm, K_range(i), Replicates, 20); avg_silh(i) mean(silhouette(data_norm, idx_tmp)); end % 5. 选定K并执行最终聚类 K_final 4; rng(42); [idx, C, sumd] kmeans(data_norm, K_final, Replicates, 20, MaxIter, 500); % 6. 结果整理与保存 clusterResult struct(); for k 1:K_final clusterResult.dates{k} dates(idx k); clusterResult.prop(k) mean(idx k); avg_max_k mean(max_vals(idx k)); % 质心从归一化空间还原到实际功率空间 clusterResult.C_real(k, :) C(k, :) * avg_max_k; end clusterResult.K K_final; clusterResult.time_axis time_axis; save(cluster_result.mat, clusterResult, idx, C, sumd); disp(聚类完成结果已保存到 cluster_result.mat);关键的参数是Replicates和MaxIter。Replicates表示算法从不同初始点开始重复运行的次数每次都会得到一个局部最优解最后返回其中WCSS最小的那次。设置20次在几千条数据上图个稳如果数据量特别大可以降到10。MaxIter控制单次迭代上限默认100在数据规整时够用但我习惯设500遇到波动剧烈的曲线也不必担心提前停止。rng(42)这一行不是可选项后面讲复现性时你会看到它的必要性。4.2 质心曲线还原从归一化空间回到实际功率聚类是在归一化空间里做的得到的质心C也是归一化的形状。要呈现给业务看必须把质心还原回有功功率的单位。这一步要认真做不能想当然地用一个全局最大值乘回去。正确做法是按簇还原每个簇有其自身的日最大功率平均水平因为日最大值归一化时除以的是每天的峰值某簇还原后的质心应当是该簇内曲线平均形状乘以该簇实际功率的中位数或均值。主程序里我用的是avg_max_k mean(max_vals(idx k))对这个簇的日最大出力做平均然后乘到归一化质心C(k,:)上。得出的C_real(k,:)就是该簇的典型有功功率曲线单位与原始数据一致。如果当初选的是装机容量归一化路线质心还原就更简单直接乘以装机容量就行。但无论哪种路线都要在结果输出前把单位换算清楚否则后续拿聚类结果去做预测或优化时数值对不上会很痛苦。4.3 输出内容簇归属、占比和典型日库聚类完成后至少有四类信息需要整理出来每个簇的质心曲线、每个簇包含的日期列表、每个簇的样本占比、每个样本的簇归属标签idx。日期列表和占比要组合起来看。样本占比直接把簇的权重量化了比如分析某地光伏全年发电特性时强晴天占比42%、多云占比31%、阴天占比19%、雨天占比8%这个描述比任何统计指标都直观。日期列表则方便后续回溯如果想核对某个簇是否与某段时间的天气异常吻合直接从clusterResult.dates{k}里取日期去查气象记录即可。我还会把这四类信息打包成结构体保存后续做功率预测分模型训练或者储能优化时直接load进来用不用重新聚类。保持输出的一致性是工程实践里很重要但很容易被忽视的一点。5. 结果验证轮廓图判读与天气记录对号入座5.1 轮廓图怎么看零下点多不多是关键单独看平均轮廓系数可能被平均骗过去最好的办法是把每个样本的轮廓系数画出来也就是MATLAB的silhouette函数直接出图。figure; silhouette(data_norm, idx); title(sprintf(K%d 的轮廓图, K_final));一张好的轮廓图有比较明显的特征每个簇的样本轮廓值从高到低排列成刀把形绝大部分样本的轮廓值大于0.2只有少数样本落到负值区。负值意味着这个样本距离相邻簇的中心比距离自己簇的中心更近也就是被分错了。负值样本占比超过10%时聚类质量就有问题。如果某几个簇的轮廓值普遍偏低分布像一片高原而不是刀刃状通常说明这些簇之间边界模糊数据本身没有明显的分组结构或者K取多了把本来连续的形态硬生生切开。反过来如果所有簇的轮廓值都很高平均超过0.5也要警惕可能是数据中存在大量重复或近似重复的记录比如同一电站连续多日的限电记录被当成真实天气样本收录了。5.2 和天气记录对号入座一张交叉表说明问题聚类是纯无监督的算法不知道哪组对应晴天哪组对应雨天它只是按形态分组。因此最有力的验证方法是用外部标签——当地气象站的历史天气记录——做交叉对照。我在项目里翻过气象记录后把每一天标注成晴、多云、阴、雨四类再和聚类标签做交叉统计。得到的结果大概长这样簇标签晴多云阴雨簇142300簇252862簇304308簇400537这个表一出来K4的合理性显而易见簇1几乎只包含晴天簇4几乎只包含雨天簇2和簇3则以多云和阴天为主体夹杂少量过渡天气这完全符合晴天—多云—阴天—雨天的气象渐变关系。如果交叉表呈现严重的交叉比如某个簇里晴天和雨天各占一半那就要警惕了可能是归一化方向不对把量级差异和形状差异混在一起可能是数据里有限电、断数等噪声干扰了形态也可能是K没选对某个真正独立的模态被遗漏了。5.3 聚类效果不好时的排查顺序遇到聚类效果不满意时我建议按照下面的顺序排查不要一上来就换算法先看原始数据有没有问题把每个簇里的日期列出来手动抽几天画出曲线看看是否存在明显的坏数据。再看归一化方向是不是和业务目标一致如果两种归一化路线的结果差异巨大说明数据里量级信息和形状信息在打架。重新审视K值回头看WCSS曲线和轮廓系数曲线特别关注最小簇样本数出现单独一两天的簇通常说明K偏大。确认活跃时段有没有截取正确如果夜间的零值点仍然大量参与距离计算轮廓系数会整体偏低且各个簇之间差异不清晰。最后才是考虑换距离度量或者换算法这一点放到后面扩展部分讲。这五步走完大多数情况下问题都出在前两步反而是算法本身背了锅。6. 实际项目里踩过的坑与后续扩展思路6.1 归一化一变分类逻辑全变了这个坑我记忆深刻。第一次做聚类时我用整体全局最大值做归一化得到的结果是夏季高功率晴天单独成一簇春秋天低功率晴天又成另一簇两簇形状几乎一样只是高度不同。从距离计算的角度完全合理但从天气类型识别的角度看毫无意义——我想要的晴天被量级差异劈成了两半。改成按日最大值归一化之后两类晴天合到了一起簇才和天气记录对上了号。这件事给我的教训是归一化不是算法前的一个例行公事它实际上在替你定义什么叫相似。是想让形状相似的曲线聚在一起还是让绝对出力相近的曲线聚在一起必须在动手前想明白否则算法给出什么结果你都无法合理使用。6.2 同一份代码跑两次结果不一样K-means的初始中心是随机选择的如果Replicates设成1两次运行的结果很可能不同某条位于两簇边界的曲线这次归了簇2下次归了簇3质心曲线也随之轻微变化。这在做研究、写报告的时候是个大问题复现不了结果。我的解决办法是在脚本开头固定随机种子也就是rng(42)或者任意一个你喜欢的整数。固定种子之后无论跑多少次只要数据和参数不变结果完全一致。如果你需要和别人对比实验结果这一步必须做到。Replicates从1提高到20也能显著降低落到很差局部最优的概率但并不能让随机性消失真正让结果可复现的还是固定种子。6.3 零点段越多距离越失真之前提到夜间零值段会稀释距离这里展开讲一下机制。欧氏距离计算的是逐点差值的平方和零点段的两个点如果都是0它们的差值贡献为0。这意味着任意两条曲线之间的相似度里有相当一部分来自大家夜间都是0这个显然的地理事实而不是来自真正需要关注的白天出力形态。极端情况下如果两条曲线只在白天完全相反一条峰值在上午、一条峰值在下午但因为夜间都很低欧氏距离依然可能很小。针对这个情况除了截取活跃时段还可以考虑用逐日动态裁剪把每天的曲线截到各自日出日落范围内再做一次插值对齐到统一长度。这样做后轮廓系数通常会进一步提高因为所有参与距离计算的点都在针对性地描述形态差异。6.4 K设太大空簇和离群日就来了当K值大于数据天然模态数时最典型的症状就是出现极小簇甚至单样本簇。这些孤立样本往往是某些极端日的形态比如一场午后雷阵雨把晴天打成了上午晴下午雨的混合曲线本来这种混合形态可以作为中间过渡归到相邻的大簇里但如果K给得足够多算法会专门为这一两天的曲线单开一个簇。对于这些离群日我的做法是先确认它不是坏数据然后把它归到最邻近的簇或者单独标记为特殊日用于后续的异常检测研究。不要轻易删掉这些样本因为它们虽然不是典型场景但在极端情况分析里恰恰有研究价值——比如光伏电站防雷策略、低出力日备用容量配置都是靠这些离群日才能把问题暴露出来的。6.5 进一步可以做点什么DTW、KShape、GMM把K-means链路做扎实之后有两条典型的进阶方向可以走。一条是针对时间轴轻微错位的问题。光伏曲线虽然对齐了时间网格但实际错位还是存在同一片云从东边飘过来和从西边飘过来出力波峰出现的时间差半小时甚至更久欧氏距离会把这些错位当成很大的差异。如果这种错位在你的数据里很明显可以尝试动态时间规整DTW作为距离度量或者直接用专门为时间序列聚类设计的KShape算法。代价是计算量成倍增加几千条曲线需要提前算距离矩阵内存和时间都要考虑。另一条是用概率模型替代硬划分。高斯混合模型和K-means的关系很直接K-means相当于GMM在协方差为单位阵时的硬判决版本。GMM的好处是每条曲线得到一个属于各个簇的概率而不是硬性标签。在后续做预测时如果一个曲线介于晴天和多云之间可以把两种天气下的预测结果按概率加权比单纯选一个标签更细腻。我做这类分析的一点个人体会是聚类的模型部分其实只占一小部分工作量真正花时间的是数据清洗、归一化决策和结果验证。K-means算法的参数就那么几个跑通代码一天就够但要把每个簇解释清楚、让聚类结果经得起业务推敲往往需要反复对照天气记录、排查异常日期。建议第一次跑通流程后把每个簇的质心曲线和对应日期的天气记录并排放到一起看一眼只要做过一次你对K的取值和归一化的方向就会有很直接的直觉。