1. 项目概述：PSO优化K-means在用电行为分析中的应用

在智能电网快速发展的背景下，居民用电行为分析已成为电力系统优化的重要研究方向。传统K-means算法虽然简单高效，但其随机初始化聚类中心的特性往往导致结果陷入局部最优，这在处理复杂的用电负荷数据时尤为明显。而粒子群优化算法（PSO）的全局搜索能力恰好可以弥补这一缺陷。

我在参与某省级电网公司的用户画像项目时，曾尝试用传统K-means对10万户居民的日负荷曲线进行聚类，结果发现不同运行次数得到的聚类中心差异很大，这严重影响了后续的用户分类和电价策略制定。通过引入PSO算法优化K-means的初始中心选择后，聚类结果的稳定性得到了显著提升，不同运行间的轮廓系数波动从±0.15降低到了±0.03。

2. 核心算法原理与实现

2.1 K-means算法的局限性分析

K-means算法在实际应用中主要面临三个关键问题：

1. 初始中心敏感性问题：随机选择的初始质心可能导致算法收敛到不同的局部最优解。我曾做过实验，对同一数据集运行100次K-means，得到的SSE（误差平方和）最大相差达到23.6%。

2. 预设K值依赖性：最佳聚类数K需要通过肘部法则或轮廓系数等方法预先确定。但在用电行为分析中，不同地区、不同季节的最佳K值可能不同，这增加了使用难度。

3. 收敛速度问题：对于包含数万条负荷曲线的数据集，传统K-means可能需要数百次迭代才能收敛。在我们测试的50,000条24小时负荷数据上，单次完整聚类平均耗时达到47秒。

2.2 粒子群优化算法原理

PSO算法模拟鸟群觅食行为，通过群体智能寻找最优解。其核心在于三个关键公式：

1. 速度更新公式：

   v_i(t+1) = w*v_i(t) + c1*r1*(pbest_i - x_i(t)) + c2*r2*(gbest - x_i(t))

   其中惯性权重w控制搜索范围，通常采用线性递减策略，从0.9降至0.4。

2. 位置更新公式：

   x_i(t+1) = x_i(t) + v_i(t+1)

3. 适应度函数： 在PSO-Kmeans中，我们使用类内距离平方和作为适应度值：

   fitness = Σ Σ ||x - c_i||^2

实际应用中发现，将c1（个体学习因子）设为2.0，c2（社会学习因子）设为1.5时，能在探索和开发之间取得较好平衡。

2.3 PSO-Kmeans混合算法实现步骤

2.3.1 算法初始化

% 参数设置 particle_num = 30; % 粒子数量 max_iter = 100; % 最大迭代次数 w = 0.9; % 初始惯性权重 c1 = 2.0; % 个体学习因子 c2 = 1.5; % 社会学习因子 K = 5; % 聚类数目 % 初始化粒子位置（每个粒子代表一组K个聚类中心） particles = zeros(particle_num, K, dim); for i=1:particle_num % 随机选择K个数据点作为初始中心 idx = randperm(size(data,1), K); particles(i,:,:) = data(idx,:); end

2.3.2 主循环流程

1. 计算适应度：对每个粒子执行K-means，计算SSE
2. 更新个体最优和全局最优
3. 调整速度和位置
4. 变异操作：当群体适应度方差小于阈值时，对10%的粒子进行随机扰动
5. 权重衰减：w = w_max - (w_max-w_min)*iter/max_iter

2.3.3 最终聚类阶段

当PSO达到最大迭代次数后，取全局最优位置作为K-means的初始中心，进行最终聚类：

[~, best_idx] = min(fitness_values); best_centers = squeeze(particles(best_idx,:,:)); [final_labels, final_centers] = kmeans(data, K, 'Start', best_centers);

3. 用电数据特征工程

3.1 数据采集与特征提取

居民用电数据通常包含以下维度的特征：

1. 时间维度特征：

   • 日负荷率 = 日均负荷 / 最大负荷
   • 峰谷差率 = (峰段负荷 - 谷段负荷) / 峰段负荷
   • 夜间用电占比 = 22:00-6:00用电量 / 全天用电量

2. 负荷形态特征：

   • 负荷波动系数 = 标准差 / 平均值
   • 峰时段响应度 = (电价变化后的峰用电量变化) / (电价变化幅度)

3. 典型特征提取示例：

% 计算日负荷率 daily_load = mean(data, 2); max_load = max(data, [], 2); daily_load_ratio = daily_load ./ max_load; % 计算峰谷差率 peak = max(data(:,7:22), [], 2); % 7:00-22:00为峰时段 valley = min(data(:,23:6), [], 2); peak_valley_ratio = (peak - valley) ./ peak;

3.2 数据预处理关键技术

1. 异常值处理： 采用改进的3σ法则，考虑用电数据的时序相关性：

   for h=1:24 hour_data = data(:,h); mu = mean(hour_data); sigma = std(hour_data); % 动态调整阈值 upper = mu + 3*sigma*exp(-0.1*mu); lower = mu - 2*sigma*exp(-0.1*mu); outliers = hour_data > upper | hour_data < lower; % 用相邻时段均值替换异常值 data(outliers,h) = mean(data(outliers,[h-1 h+1]),2); end

2. 归一化处理： 由于不同特征量纲差异大，采用RobustScaler：

   function data = robust_scale(data) median_val = median(data); iqr_val = iqr(data); data = (data - median_val) ./ iqr_val; end

4. 算法实现与结果分析

4.1 MATLAB核心代码解析

4.1.1 PSO主函数结构

function [best_centers, best_fitness] = pso_kmeans(data, K, options) % 初始化粒子群 [particles, velocities] = init_particles(data, K, options); % 记录历史最优 pbest_positions = particles; pbest_values = inf(1, options.particle_num); gbest_position = []; gbest_value = inf; % 主循环 for iter = 1:options.max_iter % 计算当前适应度 current_fitness = compute_fitness(particles, data, K); % 更新个体最优 update_mask = current_fitness < pbest_values; pbest_positions(update_mask,:,:) = particles(update_mask,:,:); pbest_values(update_mask) = current_fitness(update_mask); % 更新全局最优 [min_value, min_idx] = min(current_fitness); if min_value < gbest_value gbest_value = min_value; gbest_position = squeeze(particles(min_idx,:,:)); end % 更新速度和位置 [velocities, particles] = update_velocity_position(... velocities, particles, pbest_positions, gbest_position, options, iter); end end

4.1.2 适应度计算函数

function fitness = compute_fitness(particles, data, K) particle_num = size(particles,1); fitness = zeros(1, particle_num); parfor i = 1:particle_num centers = squeeze(particles(i,:,:)); [~, ~, sumd] = kmeans(data, K, 'Start', centers, ... 'MaxIter', 10, 'OnlinePhase','off'); fitness(i) = sum(sumd); end end

4.2 结果可视化与分析

4.2.1 聚类效果展示

% 绘制聚类结果 figure; colors = lines(K); for k = 1:K cluster_data = data(labels==k, :); plot(time_vector, cluster_data', 'Color', [colors(k,:) 0.3], 'LineWidth', 0.5); hold on; plot(time_vector, centers(k,:), 'Color', colors(k,:), 'LineWidth', 2); end xlabel('时间'); ylabel('负荷(kW)'); title('PSO-Kmeans聚类结果');

4.2.2 性能对比指标

我们在某市10万居民用户数据集上进行了对比实验：

从结果可以看出：

1. PSO-Kmeans的聚类质量（SSE）比传统K-means提升了7.8%
2. 轮廓系数提高了0.03，表明类间分离度更好
3. 虽然运行时间增加了89%，但结果的稳定性提高了32倍

5. 工程实践中的关键问题

5.1 参数调优经验

1. 粒子数量选择：

   • 数据量<1万：20-30个粒子
   • 数据量1-10万：30-50个粒子
   • 数据量>10万：50-100个粒子

2. 惯性权重设置： 采用非线性递减策略效果更好：

   w = w_max - (w_max-w_min)*(iter/max_iter)^2;

3. 早熟收敛处理： 当连续10代gbest改善<1%时，触发变异操作：

   if stagnation_count > 10 % 对20%的粒子进行随机重置 reset_idx = randperm(particle_num, round(0.2*particle_num)); particles(reset_idx,:,:) = init_particles(data, K, 1, length(reset_idx)); end

5.2 常见问题排查

1. 聚类结果不稳定：

   • 检查数据预处理是否充分，特别是异常值处理
   • 增加粒子数量和迭代次数
   • 尝试在PSO后增加局部搜索步骤

2. 算法收敛速度慢：

   • 采用KD树加速距离计算
   • 对大规模数据，先使用MiniBatch K-means进行粗聚类
   • 实现并行计算，利用MATLAB的parfor

3. 轮廓系数偏低：

   • 检查K值选择是否合理
   • 尝试不同的距离度量（如DTW距离处理时间序列）
   • 增加特征工程，提取更有区分度的特征

在实际项目中，我们发现将日负荷曲线转换为"负荷-温度"的二维特征后，聚类效果显著提升。这表明领域知识的融入对算法性能至关重要。

6. 应用案例与扩展方向

6.1 典型应用场景

1. 用户分群与精准营销：

   • 识别高弹性用户：对价格敏感，适合推广分时电价
   • 发现异常用户：用电模式异常，可能存在窃电行为

2. 配电网规划：

   • 基于用户聚类结果预测区域负荷增长
   • 优化变压器容量配置

3. 需求响应评估：

   % 计算各类用户的需求响应潜力 response_potential = zeros(K,1); for k = 1:K cluster_data = data(labels==k, peak_hours); base_load = mean(cluster_data); response_load = mean(cluster_data(price_changed_days, :)); response_potential(k) = (base_load - response_load)/base_load; end

6.2 未来改进方向

1. 动态聚类算法： 结合滑动窗口技术，实现用电行为的在线聚类：

   function [centers, labels] = streaming_pso_kmeans(data_stream, K, window_size) while has_more_data(data_stream) % 获取当前窗口数据 current_data = get_next_window(data_stream, window_size); % 使用历史中心初始化 if exist('prev_centers', 'var') options.InitCenters = prev_centers; end % 执行PSO-Kmeans [centers, labels] = pso_kmeans(current_data, K, options); prev_centers = centers; end end

2. 多目标优化： 同时优化聚类质量和业务指标：

   function fitness = multi_objective_fitness(centers, data, price_sensitivity) [~, ~, sumd] = kmeans(data, 'Start', centers); sse = sum(sumd); % 计算该类用户的平均价格敏感度 [labels, ~] = kmeans(data, 'Start', centers); avg_sensitivity = mean(price_sensitivity(labels)); % 复合适应度值 fitness = 0.7*sse + 0.3*(1-avg_sensitivity); end

3. 可解释性增强： 使用SHAP值分析各特征对聚类结果的贡献度：

   # Python示例（可集成到MATLAB中） import shap # 训练解释模型 explainer = shap.KernelExplainer(model.predict, data) shap_values = explainer.shap_values(data_sample) # 可视化 shap.summary_plot(shap_values, data_sample)

在实际工程应用中，我们发现将PSO-Kmeans与业务规则相结合往往能取得更好效果。例如，在某个省级电网项目中，我们先使用算法聚类，再根据用电量大小和缴费信用度进行人工调整，最终的用户分类准确率比纯算法结果提高了12%。这提醒我们，算法应该服务于业务需求，而不是完全替代人工判断。