MATLAB光伏功率曲线K-means聚类：从数据预处理到典型日提取全流程

光伏时间序列聚类这件事我最早接触是在做某地分布式光伏消纳评估的时候。当时手里攥着一年多的小时级功率数据一共几千条曲线每条都长得不太一样有的中午平得像一条直线有的像锯齿一样上下乱跳还有的早上九点就飙到满发又迅速掉下来。甲方说要“摸清不同天气类型下的出力特征”可没人告诉我该分几类、怎么分、分完怎么用。后来用MATLAB把K-means聚类跑通又对着结果琢磨了很久才算是把这条链路彻底理顺。这篇就把我从数据预处理到聚类结果解释的完整过程写下来尤其是那些文档里不会写、但实操中一定会踩的坑。先说明白这篇文章适合谁手里有光伏电站的功率或辐照度历史数据想用MATLAB做场景划分、典型日提取、或给预测模型做训练集标签的工程师和研究生。默认你会基础的MATLAB操作知道表格怎么读、数组怎么切片但不需要你懂太多统计背景——我会把K-means里跟时间序列相关的原理掰开讲清楚。1. 为什么给光伏功率曲线做聚类从业务需求到数学问题你先想一个问题一条光伏功率曲线本质上是气象条件经过物理转换后的一张“快照”。晴天是一条光滑的单峰曲线多云天是高频抖动的波浪线阴雨天则是一条贴近横轴的矮线。不同曲线背后对应的是不同的辐照度变化模式、不同的温度区间、甚至不同的空气污染程度。如果我们能把上千条曲线自动分成几个典型的“形状组”每一组对应一类天气过程那后面很多事情都好办。1.1 聚类的三类典型用途场景划分、典型日提取与预测建模以我实际做过的项目为例聚类结果至少有三种直接用途。第一是电网侧的调度场景划分。调度员不可能逐一查看未来一天的功率预测曲线但如果把历史出力归成“晴、多云、阴、雨”四五个典型模式再给每个模式配上负荷特性和调节策略调度预案就变成查表操作了。这里聚类的输入是全天24点功率曲线输出是模式标签。第二是光伏电站的典型日提取。做电站后评估或发电量折减分析时我们需要从一整年的数据里挑出几个“代表日”来算PR值Performance Ratio。直接取平均值会抹掉峰谷特性而聚类后取每个簇的质心曲线作为典型日既能保留形状特征又能代表该场景的中等水平。第三是预测模型的训练集划分。如果你在用BP神经网络或LSTM做超短期功率预测把历史样本按聚类标签分开训练多个子模型往往比一个模型硬扛所有天气要好——我在多个数据集上验证过分场景建模的RMSE通常能下降8%到15%。这其实就是“分而治之”的朴素思想。1.2 为什么选K-means而不是其他聚类算法做时间序列聚类可选算法不少层次聚类、DBSCAN、高斯混合模型、甚至深度嵌入聚类。但我在工程实践中绝大多数场景还是会先用K-means原因很实在。K-means的计算复杂度是O(n·k·d·iter)对几千条24维曲线来说在普通台式机上跑完也就几秒钟。层次聚类需要维护一个n×n的距离矩阵样本上万时就非常吃内存。DBSCAN的ε和minPts参数对密度变化敏感而光伏曲线的簇密度在不同天气类型下差异很大调参会调到怀疑人生。K-means只有一个K值要定而且定K值有比较成熟的启发式方法后面会讲这对工程交付来说非常重要。还有一点很多人忽略K-means的“均值”特性对光伏曲线来说反而是优点。簇质心是所有样本逐点平均的结果天然就是一条平滑的典型曲线而层次聚类选出的代表样本可能正好是一条带测量噪声的异常曲线。你要给甲方看“典型晴天出力曲线”一条平滑的质心显然比一条毛糙的实测曲线更有说服力。当然K-means也有硬伤它假设簇是凸形的、大小相近对异常值敏感。这在光伏曲线上表现为——如果某天数据里混入了限电段出力被硬生生削平它会被单独拉成一个簇或者严重扭曲相邻簇的质心。所以后面我专门用一节讲数据清洗这一步偷懒的话前面所有功夫都白费。2. 数据预处理决定聚类上限的隐形战场很多教程会把数据预处理一笔带过直接读数据、调kmeans函数、画图完事。但以我的经验预处理对光伏聚类结果的影响远大于算法本身。你可以把K-means想象成一个雕刻师预处理就是给他一块什么成色的石料——石料里有裂缝杂质雕刻师再厉害也雕不出好东西。2.1 原始数据清洗的三道关口第一关是时间戳对齐。光伏数据来自数据采集器SCADA或独立气象站经常出现缺失、重复、时间不同步。我的做法是先统一生成一个完整的时间向量比如从2023年1月1日0点到2024年1月1日0点逐小时再用retime或outerjoin把原始数据映射上去缺失的时段先标记为NaN绝不急着插值。第二关是异常值判别与限电段识别。光伏功率的物理上限是装机容量乘以一个略大于1的系数考虑超发超过这个上限的值基本是采集错误。下限是负数夜间负值通常是逆变器自消耗或测量偏置直接截断为0。但最隐蔽的是限电段功率长时间平贴在某个固定值比如装机容量的60%且波动极小同时辐照度数据正常——这种“削顶”曲线如果不剔除聚类会分出一个没有物理意义的“限电平顶簇”。第三关是夜间时段取舍。这直接影响你聚类的输入维度。光伏出力在夜间恒为0如果保留夜间点所有曲线在横轴上重合K-means会天然把这些维度视为“相似”变相强化了白天的形状差异。但夜间段没有任何区分度还会稀释白天曲线的差异。我实测后发现把每天的数据截取到“日出前半小时到日落后半小时”比如春夏秋三季取6点到19点冬季取7点到17点聚类结果明显更清晰。下面这段代码展示如何实现时段截取并剔除无效天% 假设 data 是 timetable包含变量 PowerkW和 Time % 设定每天的分析时段这里用 [6, 19] 作为默认窗口可根据季节调整 window_start 6; window_end 19; % 提取每天的日期标签 day_labels dateshift(data.Time, start, day); [days, ~, day_id] unique(day_labels); n_days length(days); % 遍历每一天截取时段并判断有效性 power_daily NaN(n_days, window_end - window_start 1); for i 1:n_days idx (day_id i); t_day data.Time(idx); p_day data.Power(idx); % 筛选在窗口内的小时 hour_vals hour(t_day); win_idx (hour_vals window_start) (hour_vals window_end); if sum(win_idx) (window_end - window_start 1) % 完全覆盖窗口才保留 power_daily(i, :) p_day(win_idx); end end % 剔除含 NaN 的行缺测以及最大值小于 20% 装机容量的天夜间/无效天 % 假设 installed_capacity 为已知装机容量 valid_mask ~any(isnan(power_daily), 2) (max(power_daily, [], 2) 0.2 * installed_capacity); power_daily power_daily(valid_mask, :);这里有个容易忽略的细节有效性判断里我用了“最大值大于20%装机容量”而非“全天总和大于某个阈值”因为多云天的全天发电量可能很低但仍有分析价值而全天都趴在零附近的无效天才是我们要剔除的。2.2 特征缩放要不要标准化以及如何标准化K-means基于欧氏距离特征的量纲直接影响距离计算。光伏功率曲线的量纲是kW如果某条曲线峰值是5000 kW另一条是300 kW两者在“峰值高度”这一维度上的差异会完全淹没形状差异。所以标准化是必须的。但对时间序列来说标准化不是简单的zscore而是有策略的。我常用的有两种方式一是逐日最大值归一化每条曲线除以当天的最大值把所有曲线压到[0, 1]区间。这样保留的是“形状信息”丢掉的是“幅度信息”。适合你想按天气过程分组的场景因为晴天和多云天可能峰值都很高区别主要在波动形态。二是全样本统一标准化先计算所有数据的均值和标准差再统一变换。这样保留幅度差异适合你想把“辐照资源好坏”也作为分类依据的场景比如把强辐照晴天和弱辐照晴天分开。实际项目中我强烈建议两种都跑一遍然后对比轮廓系数和业务可解释性选更合理的一个。代码如下% 方式一逐日最大值归一化 pow_norm power_daily ./ max(power_daily, [], 2); % 方式二全局 z-score 标准化 mu mean(power_daily, all); sigma std(power_daily, 0, all); pow_zscore (power_daily - mu) / sigma;补充一个容易踩的坑光伏电站出力跟装机容量强相关如果你在做多个电站的对比聚类必须先除以各自的装机容量得到“归一化出力per-unit”否则大电站会把小电站的曲线全部吞掉聚出个寂寞。再补充一个进阶技巧如果数据是15分钟级的高频采样曲线维度会达到96点甚至更高直接聚类会产生“维度灾难”——欧氏距离在高维空间趋于平均化。这时可以先做主成分分析降到5~10维或者提取手工特征峰值时刻、波动率、爬坡率、总发电量等再聚类效果通常会更好。3. K-means的数学机制与MATLAB实现细节K-means不是简单调个函数就行尤其在时间序列上有几个机制层面的细节直接决定聚类质量。下面结合MATLAB的kmeans函数参数逐一说清楚。3.1 距离度量为什么默认选平方欧氏距离kmeans默认使用平方欧氏距离这并非没有原因。K-means的目标函数是簇内误差平方和SSESum of Squared Errors它等价于在欧氏距离下的硬划分。在时间序列场景里平方欧氏距离对“逐点差异”求和恰好放大了波形错位和幅值差异物理意义明确两个时刻点相差越大对聚类目标的贡献是二次增长。正因如此除非曲线做过特殊处理比如DTW距离否则我一般不改默认度量。kmeans函数签名很长这里给出最常用的参数组合rng(42); % 固定随机种子保证结果可复现 k 4; % 先假定4类稍后用肘部法则验证 [idx, C, sumd] kmeans(pow_norm, k, ... Distance, sqeuclidean, ... Replicates, 10, ... MaxIter, 500, ... Display, final);参数逐一说一下ReplicatesK-means对初始质心敏感不同初始化可能收敛到不同的局部最优。设10次重复、每次随机选初始质心MATLAB会返回SSE最小的一次结果。对几千条样本来说10次重复的计算成本微乎其微但能显著提升聚类稳定性。MaxIter默认100次迭代有时候不够尤其在高维曲线上。500次能在更多情况下收敛且不会明显增加耗时。Display设为final可以输出每次重复的SSE方便我们观察各次重复间的差异判断是否碰到了不理想的局部最优。3.2 时间序列特有陷阱相位对齐问题直接用原始功率曲线做K-means隐含一个假设所有曲线在时间轴上已对齐且“中午12点的特征在每条曲线上都是中午12点”。光伏数据大多满足这一点因为太阳运行轨迹有规律但也正是这个“有规律”会掩盖一个重要问题——不同月份、不同经度的光伏曲线峰值出现时刻其实在漂移。夏天峰值一般在正午12点到13点冬天在11点到12点东西跨度大的地区甚至会差出1~2个小时。这个漂移是物理规律但你做聚类时如果目标是识别“晴天vs多云vs阴天”的天气模式相位漂移会产生一个副作用同一个天气模式下夏季曲线和冬季曲线因为“峰位不同”被分到不同的簇。解决办法有两个。一是把曲线按峰值时刻对齐每天找最大值对应的小时把所有曲线向左或向右平移让峰值都落在同一时刻。这样聚类聚焦在波形形状而非相位。二是不解决相位问题反而把它当成分类依据如果你压根就想区分“资源型”场景不同光伏利用小时数那保留相位差异可能是对的。关键是你得清楚自己在做什么——这在聚类里比算法本身更重要。峰值对齐的实现思路比较简单function aligned align_peak(x) [~, peak_idx] max(x); center_idx ceil(length(x) / 2); shift center_idx - peak_idx; aligned circshift(x, shift); % 注意circshift 是循环移位会引入边界伪影 % 更稳妥的办法是用 NaN 填充边界再插值这里为了展示简化处理 end需要提醒的是circshift循环移位会把曲线末尾的数据挪到开头这在光伏曲线上会产生“夜间0值突然出现在清晨”的伪影影响聚类质量。工程上更推荐的做法是截取峰值前后固定窗口比如峰值前4小时、后6小时作为聚类输入宁可减少维度也要保证特征的物理一致性。3.3 距离度量在时间序列上的替代方案除了平方欧氏距离偶尔也会用相关性距离或DTW动态时间规整。相关性距离1减去皮尔逊相关系数对幅值不敏感只关心形状趋势——听起来很适合光伏但我用下来的经验是它会把“两个形状相似但一个容量差异大”的电站聚在一起这在某些业务场景比如电网调度要看实际可调出力反而不利。DTW可以处理时间轴伸缩问题对日照时长随季节变化的场景有天然优势。但DTW没有原生的kmeans实现你得自己写距离矩阵再丢给kmedoids而不是kmeans每一次迭代都要重新计算质心的DTW中心计算量翻好几倍。我的建议是先用欧氏距离跑一版如果簇内曲线在业务上明显牵强比如同一个簇里混入了峰位相差2小时的晴天再考虑DTW或前述的峰值对齐方案。4. K值选择与聚类评估怎么确定“分几类”最合理K值的确定是K-means实操里最让人纠结的一步因为不存在一个完美的“正确答案”。有统计指标可以帮忙但最终判据一定是业务含义是否自洽。4.1 肘部法在光伏场景下的实际表现传统的肘部法画SSE随K的变化曲线找“拐点”。但光伏曲线的SSE曲线通常不像教科书里那么理想——它往往是平滑下降的没有明确的肘部。我处理过多个电站数据SSE曲线大多形态类似K2到3下降很快K4到6趋缓K7以后基本是线性缓降。这种曲线的“肘”一般出现在K3或K4的位置。这里有个实操技巧不要只看SSE绝对值而是看SSE下降率也就是相邻两个K之间SSE的边际下降百分比。当边际下降率首次跌破某个阈值比如15%时对应的K就是一个合理选择。这样可以避免依赖肉眼判断的模糊性。K_max 10; sse_vals zeros(1, K_max); for k 1:K_max [~, ~, sumd] kmeans(pow_norm, k, Replicates, 5, MaxIter, 500); sse_vals(k) sum(sumd); end % 计算边际下降率 marginal_reduction -diff(sse_vals) ./ sse_vals(1:end-1); % 找第一个边际下降率低于 0.15 的 K k_elbow find(marginal_reduction 0.15, 1, first); fprintf(肘部法建议 K %d\n, k_elbow);4.2 轮廓系数量化“类内紧、类间松”的辅助判据轮廓系数的思想很直观对每个样本计算它到同簇其他样本的平均距离a(i)再计算它到最近的其他簇所有样本的平均距离b(i)轮廓系数s(i) (b(i)-a(i)) / max(a(i), b(i))。s(i)接近1说明样本分类非常清晰接近-1说明它可能分错了簇。全样本的平均轮廓系数越大聚类质量越好。MATLAB自带的silhouette函数画图一目了然figure; silhouette(pow_norm, idx); title(K4 时的轮廓系数图);我实际跑过的一个电站数据K3时平均轮廓系数0.65K4是0.62K5是0.58。如果只看轮廓系数选K3更优但画出来的簇中心显示K3的分法里第三簇把“多云”和“阴天”强行合并了业务上无法解释。这就是我强调的——指标是参考业务含义是最终裁判。4.3 业务判据三种天气类型的可解释性检验聚类结果必须接受一个终极提问每个簇能不能用一句话说清楚它是什么场景我的经验是光伏曲线聚类后通常有几种稳定的业务标签晴天簇质心曲线平滑峰值高且接近正午早晚迅速攀升和下降日发电量占比高。多云/晴间多云簇质心曲线有明显波动峰值次高可能有多个小峰或平台期。阴天簇质心曲线低矮平滑全天变化平缓峰值一般低于装机容量的40%。雨天/恶劣天气簇如果样本足够多质心曲线极低且伴有长时段波动。把簇质心画成同一张图如果四个簇之间“长得明显不一样”K4就是一个好的选择。如果发现两个簇质心几乎重合说明K选大了或数据本身可分性弱需要回到预处理阶段找原因。5. 数据输入格式、维度选择与kmeans批处理注意事项很多人在第一步就卡住了MATLAB的kmeans要求输入是一个矩阵每行一个观测、每列一个特征。对光伏时间序列来说矩阵的建立方式决定了整个分析的成败。5.1 如何把timetable转化为kmeans吃的矩阵如果你用readtimetable读取了CSV得到的是一个包含多天数据的timetable要做聚类必须先把“连续时间序列”转换为“每日一条曲线”的二维矩阵。这里给出一个完整的转换示例% 读取数据 data readtimetable(pv_power.csv); % 假设数据每小时一条连续记录一年 % 提取日期和小时 d dateshift(data.Time, start, day); h hour(data.Time); % 确定每天的小时范围 hours_of_day unique(h); n_hours length(hours_of_day); % 理论上24 days_unique unique(d); n_days length(days_unique); % 初始化矩阵每行一天每列一个时刻 X NaN(n_days, n_hours); for i 1:n_days day_mask (d days_unique(i)); for j 1:n_hours hour_mask (h hours_of_day(j)); vals data.Power(day_mask hour_mask); if ~isempty(vals) X(i, j) mean(vals); % 若同一小时有多条记录取平均 end end end % 清洗删除全为NaN或全为0的行 valid_rows ~all(isnan(X), 2) any(X 0.05 * installed_capacity, 2); X X(valid_rows, :);这段代码看起来繁琐但它把所有边界情况都处理了缺失时段置为NaN而不是00会污染聚类重复记录取平均无效天直接剔除。它比用reshape快进快出要稳重得多——因为reshape要求数据严格均匀一旦中间有缺失整个数组就错位了。5.2 维度全部拿进来还是只拿“日间时段”前面说了夜间时段全是0会稀释聚类效果。除了时间窗口选择还有一个维度层面的技巧不要用全天的24维数据而是截取活跃时段比如6点到19点共14维。截取后维度变低一方面计算更快另一方面有效规避了“维数灾难”。还有个附加好处画图解读时14维曲线比24维曲线更清爽给甲方汇报也更直观。5.3 大规模数据下的内存与计算优化如果你手头不是一年数据而是五年多电站数据比如几万条日曲线kmeans直接跑也还能扛但有几个优化技巧值得掌握用Options, statset(UseParallel, true)开启并行计算多核加速明显。对超大数据集先用一个子集比如1万条做聚类得到质心再用Start, C0把质心作为初始值对全量数据跑一次迭代收敛。这样可以避免每次随机初始化的不确定性也大幅缩短总计算时间。如果数据量超过内存能舒服处理的范围比如超过50万条先做datasample随机抽样聚类再用knnsearch把剩余样本分配到最近的簇。这在工程上是合法的近似做法。% 子集聚类确定初始质心 sub_idx datasample(1:size(X, 1), 10000, Replace, false); C0 kmeans(X(sub_idx, :), k, Replicates, 3); % 全量数据用 C0 初始化继续聚类 [idx, C] kmeans(X, k, Start, C0, MaxIter, 300);6. 聚类结果的可视化与解释簇质心、样本分布与时间热力图聚类不是为了得到一个标签数组而是为了从数据里读出业务故事。可视化的核心是让“每个簇代表什么场景”一目了然。6.1 簇质心曲线与样本分布图基础操作是把每类的质心曲线画在一起同时用半透明细线叠印每类的若干条样本曲线展示类内离散程度figure; hold on; colors lines(k); for j 1:k cluster_curves pow_norm(idx j, :); % 画样本细线最多画50条避免太乱 n_show min(50, size(cluster_curves, 1)); plot(6:19, cluster_curves(1:n_show, :), -, Color, [colors(j, :), 0.15]); end % 画质心粗线 for j 1:k plot(6:19, C(j, :), LineWidth, 2.5, Color, colors(j, :), DisplayName, [簇 num2str(j)]); end xlabel(时刻 (h)); ylabel(归一化功率); legend show; grid on;这张图能直观回答四个簇是否足够区分如果两簇的粗线基本重合说明业务上合并更合理如果某簇的细线满天飞说明该类内部还有细分空间可能需要增大K或引入更多特征。6.2 簇样本的时间分布查看季节性模式聚类标签结合时间信息能发现很有价值的规律。比如某个簇集中在夏季某个簇在冬季这本身就是对光伏资源特征的侧面描述。我常用histogram画簇标签的月份分布% 假设 month_labels 是从原始日期提取的月份 figure; histogram(month_labels(idx 1), BinEdges, 0.5:12.5, FaceColor, colors(1, :)); hold on; histogram(month_labels(idx 2), BinEdges, 0.5:12.5, FaceColor, colors(2, :)); % ... 依次画完各个簇 xlabel(月份); ylabel(样本数); legend(簇1, 簇2, ...);如果某个簇的样本在月份上高度集中比如簇1只出现在5到8月那这个簇很可能是“夏季高辐照晴天模式”与簇3的“冬季晴天模式”峰位偏早、峰值偏低形成呼应。这种发现是PPT汇报的亮点也是你对数据理解深度的体现。6.3 一个实用的“按簇着色的日曲线堆叠图”还有一个可视化技巧把所有样本曲线按日期顺序纵向排列每条曲线用颜色映射功率大小并按聚类标签分组排列形成类似热力图的展示。这样整个一年的出力模式“一眼望穿”甲方看了直呼专业。figure; % 先把样本按簇和日期排序 [~, order] sort(idx); imagesc(pow_norm(order, :)); % 在Y轴刻度上标记簇边界 colormap(parula); xlabel(时刻); ylabel(样本序号); title(按聚类结果排序的日功率曲线热力图);排序后同一簇的曲线在图像上呈水平带状排列颜色模式高度一致簇与簇之间有肉眼可见的色差边界——这种全局视图比任何统计指标都有说服力。7. 进阶主题聚类结果如何反哺光伏功率预测与调度聚类的终点不是画图而是下场应用。我在实际项目中验证过的三个应用路径如下供你参考。7.1 分场景预测模型的训练集划分做法非常简单用聚类标签把原始数据集分割成几个子集分别训练预测模型可以是BP、LSTM、或者简单的线性回归。预测时先用一个分类器比如基于清晨辐照度或数值天气预报判断当天属于哪个场景再调用对应的子模型。具体代码示意% 假设 train_data 有特征矩阵 Feat 和功率标签 Power % cluster_labels 是聚类结果 for j 1:k idx_j (cluster_labels j); model_j fitrnet(Feat(idx_j, :), Power(idx_j), LayerSizes, [10, 5]); models{j} model_j; end % 预测阶段先用简单规则或浅层分类器判断场景 % scene_pred predict_scene(nwp_features); % power_pred models{scene_pred}.predict(Feat_new);这里的核心逻辑是不同天气场景下功率与气象特征的映射关系不同。与其让一个模型去拟合所有场景的复杂非线性映射不如拆成几个较简单的问题。7.2 典型场景曲线在储能容量配置中的应用做光储联合系统的容量配置时需要一组“代表性日曲线”作为仿真的输入。如果直接用全年8760小时数据仿真时长不可接受如果只取平均日曲线又会抹掉波动特征。用聚类选出的每簇质心曲线作为典型场景再按每个簇的样本占比作为场景权重这样既压缩了计算量又保持了原有出力分布的多样性。这在HOMER、PyPSO等工具里都是标准的场景削减思路。我的一个实际案例用K4个典型日替代全年365条曲线做储能容量优化得到的最优容量配置与用全数据做滚动优化的结果相差在3%以内但仿真时间从4小时缩短到20分钟。这就是聚类的工程价值——用精度换速度而精度损失在可接受范围。7.3 聚类的动态更新模型漂移问题最后提醒一个容易被忽略的问题光伏电站的数据分布不是一成不变的。组件衰减、清洗周期、周边新建建筑遮挡、逆变器更换都会让功率曲线的形态逐年漂移。如果你在2023年的数据上做聚类得到的情景模型直接丢到2026年用很可能失效。所以聚类模型不能“一次定终身”。我的习惯是每季度或每半年重新聚类一次并且对比新旧标签分布是否发生了显著变化。如果某簇的样本占比突然从30%掉到15%那一定发生了什么值得关注的事情——这种可解释的变化往往是运维诊断的线索也是你做数据分析师价值的体现。如果你打算把聚类做成一个周期性运行的分析流程建议把建模脚本打包成函数或App输入新数据、输出更新后的标签和质心并自动生成对比图表。这一步自动化虽然前期要投入时间但长期来看省下的精力是几何级的。我个人的建议是不要停留在“聚类能跑出图”的层面而是要跟业务人员反复确认每个簇对应什么现场运行状态把数学簇翻译成业务语言。只有完成了这一步K-means才真正从学术玩具变成工程工具。