1 内容介绍
风电功率预测为电网规划提供重要的依据,研究风电功率预测方法对确保电网在安全稳定运行下接纳更多的风电具有重要的意义.针对极限学习机(ELM)回归模型预测结果受输入参数影响的问题,现将粒子群优化算法(PSO)应用于ELM中,提出了一种基于粒子群优化极限学习机的风功率预测方法.该方法首先将数值天气预报信息(NWP)数据进行数据预处理,并构建出训练样本集,随后建立ELM模型,利用粒子群算法优化ELM中的输入权值和阈值,从而建立起基于NWP和PSO-ELM风功率预测模型.对华东地区3个不同装机容量的风场NWP数据进行实验.结果表明:该方法的预测精度高且稳定性能好,能够为风电场功率预测以及风电并网安全可靠性提供科学有效的参考依据.
风的随机性和波动性导致了风电功率的不平 稳性，风电场的接入影响电网的稳定运行环境［１］。对风电场的输出功率进行精准预测，有助于电网人 员根据风电场输出功率的变化调整发电计划，减少 电网的备用容量以节约能源的消耗。因此，有效的 风电功率预测方法可保障电网的稳定经济运行。目前，国内外研究人员做了大量的工 作，研 究 方法主要包括物理方法、统计方法以及神经网络方 法［２］。其中物理方法要求风机相关物理信息，建模 复杂且 精 度 不 稳 定。统 计 方 法 模 型 简 单、数 据 单 一，但预测的精度受时间的限制，预测时间越长精 度越低，通常应用于超短期预测。神经网络方法泛 化能力强，能够处理回归问题，适用于风功率预测。极 限 学 习 机 （Ｅｘｔｒｅｍｅ Ｌｅａｒｎｉｎｇ Ｍａｃｈｉｎｅ， ＥＬＭ）是一类基于前馈神经网络的算法［３］。ＥＬＭ 因结构简单、学习效率高，被用以解决回归、聚类、 分类等问题。ＥＬＭ 算法的优化问题被很多学者研 究，王浩等［４］利用遗传算法优化极限学习机，并将 风电系统参数模糊化，从而提高预测模型的精度。龙浩［５］提出了加权极限学习机方法，旨在解决样本 数据不均衡的问题。王文锦等［６］利用蜂群算 法 优 化 ＥＬＭ，旨在提高模型的稳定性 能。王 宏 刚 等［７］ 将蜻蜓算法（ＤｒａｇｏｎｆｌｙＡｌｇｏｒｉｔｈｍ，ＤＡ）分布式应 用于 ＥＬＭ，优化初始化输入权重和阈值的影响，有 效提高了电能质量扰动识别率。风电功率的预测主要基于当日的数值天气预 报信 息 （Ｎｕｍｅｒｉｃａｌ ＷｅａｔｈｅｒＰｒｅｄｉｃｔｉｏｎ，ＮＷＰ）。ＮＷＰ包括风速、风向、温度、湿度以及气压［８］。实 测 ＮＷＰ数据在同一风功率下存在奇异值以及波 动的问题。关于风功率 ＮＷＰ数据的预处理问题， 符杨等［９］针对 ＮＷＰ数据不准确、爬坡事件频发等 原因将 ＮＷＰ 数据分类，并对功率波动进 行 预 测。杨茂等［１０］利用经验模态分解对风功率数据进行分 解去噪重构，一定程度上减少了噪点对预测结果的 影响。杨家然等［１１］利用模糊聚类的方法将原始信 号进行分类，采用不同的模型组合预测。本 文 将 粒 子 群 优 化 算 法 （Ｐａｒｔｉｃｌｅ Ｓｗａｒｍ Ｏｐｔｉｍｉｚａｔｉｏｎ，ＰＳＯ）和 ＥＬＭ 结 合，欲 提 高 传 统 ＥＬＭ 预测模型的精度和稳定性能，从而 为 风 电 功 率预测技术提供新的方法。
基本理论
１．１ ＥＬＭ 算法
ＥＬＭ 因学习效率高被普遍应用于风电功率预测、变压器故障诊断、风机故障诊断等方面。该 算法任意赋值输入层权重和阈值，训练过程中无需改变模型参数，仅设定隐层神经元个数，便可通过最小二乘法获得输出层权重。
1．２ ＰＳＯ－ＥＬＭ 预测模型 如前文所述，ＥＬＭ 的 初 始 输 入 权 值 和 阈 值 是 随机确 定 的，其 训 练 的 效 果 会 受 初 始 值 影 响。因 此，采用 ＰＳＯ 优化 ＥＬＭ 的输入权重和阈值，可避 免盲目 性 训 练 ＥＬＭ 模 型。ＰＳＯ 算 法 优 化 ＥＬＭ 的步骤中，先 初 始 化 ＰＳＯ 参 数，包括粒子群的规 模、空间维度、惯性参数ｗ、学习因子ｃ１ 和ｃ２、迭代 次数和最 大 速 度ｖｍａｘ等。ＰＳＯ－ＥＬＭ 预 测 模 型 是 将每个粒子对应的输入权值和阈值代入 ＥＬＭ 预 测模型中，将 ＥＬＭ 学习样本输出与实际输出的均 方误差（ＭｅａｎＳｑｕａｒｅｄＥｒｒｏｒ，ＭＳＥ）作为 ＰＳＯ 的 适应度。将粒子的当前适应度与最优适应度做对 比，若比最优适应度小，说明当前输入权值和阈值 所建立的 ＥＬＭ 模型进行预测产生的均方误差较 小，则将当前适应度更新为最优适应度，将当前位 置更新为Ｐｂ，否则保持最优适应度不变。同 理 比 较适应度和全局适应 度，更 新 Ｇｂ。当 迭 代 次 数 达 到最大值或适应度达到设定值时停止算法。ＰＳＯ 优化得到的最优输入权值 Ｗ 和 阈 值ｂ 后，代 入 ＥＬＭ 模型中进行预测。具体步骤如图１所示。
2 仿真代码
%% 导入数据 clear all;clc load data % 随机生成训练集、测试集 k = randperm(size(input,1)); % 训练集——1900个样本 P_train=input(k(1:1900),:)'; T_train=output(k(1:1900)); % 测试集——100个样本 P_test=input(k(1901:2000),:)'; T_test=output(k(1901:2000)); %% 归一化%% 极限学习机在回归拟合问题中的应用研究
[Pn_train,inputps] = mapminmax(P_train,-1,1); Pn_test = mapminmax('apply',P_test,inputps); % 测试集 [Tn_train,outputps] = mapminmax(T_train,-1,1); Tn_test = mapminmax('apply',T_test,outputps); tic %% 节点个数 inputnum=2; hiddennum=20; outputnum=1; %% 参数初始化 %粒子群算法中的两个参数 c1 = 1.49445; c2 = 1.49445; maxgen=100; % 进化次数 sizepop=30; % 种群规模 Vmax=10^(2); % 最大速度 Vmin=-Vmax; % 最小速度 function Y = elmpredict(P,IW,B,LW) Q = size(P,2); BiasMatrix = repmat(B,1,Q); tempH = IW * P + BiasMatrix; H = 1 ./ (1 + exp(-tempH)); % 计算预测输出% 训练集
k = randperm(size(input,1)); % 训练集——1900个样本 P_train=input(k(1:1900),:)'; T_train=output(k(1:1900)); % 测试集——100个样本 P_test=input(k(1901:2000),:)'; T_test=output(k(1901:2000)); %% 归一化 % 训练集 [Pn_train,inputps] = mapminmax(P_train,-1,1); Pn_test = mapminmax('apply',P_test,inputps); % 测试集 [Tn_train,outputps] = mapminmax(T_train,-1,1); Tn_test = mapminmax('apply',T_test,outputps); Y = (H' * LW)'; disp(zbest) x=zbest;%% 优化结束后种群中的最好个体 IW=x(1:inputnum*hiddennum); B=x(inputnum*hiddennum+1:inputnum*hiddennum+hiddennum); %赋给网络权值和阈值 IW=reshape(IW,hiddennum,inputnum); B=reshape(B,hiddennum,1); %% ELM创建/训练 LW = elmtrain(IW,B,Pn_train,Tn_train,'sig'); % 随机生成训练集、测试集 k = randperm(size(input,1)); % 训练集——1900个样本 P_train=input(k(1:1900),:)'; T_train=output(k(1:1900)); % 测试集——100个样本 P_test=input(k(1901:2000),:)'; T_test=output(k(1901:2000)); %% 归一化 % 训练集% 随机生成训练集、测试集
k = randperm(size(input,1)); % 训练集——1900个样本 P_train=input(k(1:1900),:)'; T_train=output(k(1:1900)); % 测试集——100个样本 P_test=input(k(1901:2000),:)'; T_test=output(k(1901:2000)); %% 归一化% 随机生成训练集、测试集
[Pn_train,inputps] = mapminmax(P_train,-1,1); Pn_test = mapminmax('apply',P_test,inputps); % 测试集 [Tn_train,outputps] = mapminmax(T_train,-1,1); Tn_test = mapminmax('apply',T_test,outputps); tic %% 节点个数 inputnum=2; hiddennum=20; outputnum=1; %% 参数初始化 %粒子群算法中的两个参数 c1 = 1.49445; c2 = 1.49445; maxgen=100; % 进化次数 sizepop=30; % 种群规模 Vmax=10^(2); % 最大速度 Vmin=-Vmax; % 最小速度 function Y = elmpredict(P,IW,B,LW) Q = size(P,2); BiasMatrix = repmat(B,1,Q); tempH = IW * P + BiasMatrix; H = 1 ./ (1 + exp(-tempH)); % 计算预测输出% 训练集
k = randperm(size(input,1)); % 训练集——1900个样本 P_train=input(k(1:1900),:)'; T_train=output(k(1:1900)); % 测试集——100个样本 P_test=input(k(1901:2000),:)'; T_test=output(k(1901:2000)); %% 归一化 % 训练集% 随机生成训练集、测试集
Pn_test = mapminmax('apply',P_test,inputps); % 测试集 [Tn_train,outputps] = mapminmax(T_train,-1,1); Tn_test = mapminmax('apply',T_test,outputps); %% 结果对比 result = [T_test' T_sim']; % 均方误差 E = mse(T_sim - T_test) % 决定系数 N = length(T_test); R2 = (N*sum(T_sim.*T_test)-sum(T_sim)*sum(T_test))^2/((N*sum((T_sim).^2)-(sum(T_sim))^2)*(N*sum((T_test).^2)-(sum(T_test))^2)) %% 绘图 figure plot(1:length(T_test),T_test,'r*') hold on plot(1:length(T_sim),T_sim,'b:o') xlabel('测试集样本编号') ylabel('测试集输出') title('ELM测试集输出') legend('期望输出','预测输出')[Pn_train,inputps] = mapminmax(P_train,-1,1);
plot(1:length(T_test),T_test-T_sim,'r-*') xlabel('测试集样本编号') ylabel('绝对误差') title('ELM测试集预测误差')figure
Q = size(P,2); BiasMatrix = repmat(B,1,Q); tempH = IW * P + BiasMatrix; H = 1 ./ (1 + exp(-tempH)); % 计算预测输出 Y = (H' * LW)';function Y = elmpredict(P,IW,B,LW)
k = randperm(size(input,1)); % 训练集——1900个样本 P_train=input(k(1:1900),:)'; T_train=output(k(1:1900)); % 测试集——100个样本 P_test=input(k(1901:2000),:)'; T_test=output(k(1901:2000)); %% 归一化% 随机生成训练集、测试集
[Pn_train,inputps] = mapminmax(P_train,-1,1); Pn_test = mapminmax('apply',P_test,inputps); % 测试集 [Tn_train,outputps] = mapminmax(T_train,-1,1); Tn_test = mapminmax('apply',T_test,outputps); tic %% 节点个数 inputnum=2; hiddennum=20; outputnum=1; %% 参数初始化 %粒子群算法中的两个参数 c1 = 1.49445; c2 = 1.49445; maxgen=100; % 进化次数 sizepop=30; % 种群规模 Vmax=10^(2); % 最大速度 Vmin=-Vmax; % 最小速度 function Y = elmpredict(P,IW,B,LW) Q = size(P,2); BiasMatrix = repmat(B,1,Q); tempH = IW * P + BiasMatrix; H = 1 ./ (1 + exp(-tempH)); % 计算预测输出% 训练集
%% SOTracking - Single Object Tracking
% data format
% *data format:* 1-X, 2-Y, 3-Z, 4-RANGE, 5-AZIMUTH, 6-ELEVATION, 7-DOPPLER,
% 8-POWER, 9-POWER_VALUE, 10-TIMESTAMP_MS
%% env init
clear, clc, close
addpath(genpath('./utils'));
%% param
% path and data
result_dir = './result/';
data_dir = './data/mmWave_radar_data/';
data_item = 'SOT/';
start_frame = 1;
end_frame = 10000;
traj_dim = 2; % 2d/3d trajectory
% denoise
param_denoise.dpl_thr = 0.15;
param_denoise.loc_thr = [-40, 40, 0, 40, -5, 40];
% cluster
epsilon = 5;
MinPts = 20;
obj_count = 1;
% Kalman filter
motion_type = 'ConstantVelocity'; % 'ConstantVelocity' | 'ConstantAcceleration'
param_kf = getDefaultKFParameters(motion_type);
% param.initialEstimateError = 1E5 * ones(1, 2);
% param.motionNoise= [25, 10];
param.measurementNoise = 100;
% show
% axis_range = [-5, 5, 0, 20, -2, 5];
axis_range = [-10, 10, 0, 20, -1, 5];
show_delay = 0.0;
%% denoise, cluster, KF_tracking
% ---- file info ----
datas = dir([data_dir data_item '*.txt']);
data_names = {datas.name};
data_num = length(data_names);
end_frame = min(data_num, end_frame);
if start_frame>end_frame
error("start frame over range")
end
% ---- init ----
KF = []; % KF handle
det_loc = []; % detected location
meas_traj = NaN(start_frame-1,traj_dim); % trajectory points
kf_traj = NaN(start_frame-1,traj_dim); % KF corrected trajectory points
bounding_box = NaN(start_frame-1,traj_dim*2); % bounding box
isDetected = false; % detected flag
figure;
for k = start_frame:end_frame
% ---- load data
frame=importdata([data_dir data_item data_names{k}]);
% ---- denoise ----
frame_clean = point_cloud_denoise(frame, param_denoise);
disp(['effective points num: ' num2str(size(frame_clean,1))])
% [ToDo] TBD
if size(frame_clean, 1) < 4
isDetected = false;
end
idx = DBSCAN(frame_clean(:,[1,2]),epsilon,MinPts); % DBSCAN Cluster
% delete noise points cluster(idx==0)
frame_clean(idx==0,:) = [];
idx(idx==0,:) = [];
% [idx,C] = kmeans(frame_doppler(:,[1,2]),2); % K-Means Cluster
[idx, Dg] = cluster_idx_arranege(frame_clean(:,[1,2]), idx);
disp(['cluster count:' num2str(numel(unique(idx)))])
if isempty(idx)
isDetected = false;
else
isDetected = true;
end
if isDetected
frame_obj = frame_clean(idx<=obj_count,:);
subplot(121)
gscatter3(frame_obj(:,1),frame_obj(:,2),frame_obj(:,3),idx,[],[],10,'on')
axis(axis_range)
% calc bounding box
rect_min = min(frame_obj(:,1:3),[],1);
rect_max = max(frame_obj(:,1:3),[],1);
rect_size = rect_max - rect_min;
rect_center = calcCentroid(frame_obj(:,1:3));
det_loc = rect_center(1:traj_dim);
% show bounding box
plotBoundingbox(rect_min, rect_size, [0 0 1], 'obj1', k, axis_range)
else
det_loc = NaN(1,traj_dim);
end
% Kalman Filter
[kf_loc, KF, states] = KF_step(det_loc, KF, param_kf);
if isempty(kf_loc)
kf_loc = NaN(1,traj_dim);
end
meas_traj(k,:) = det_loc;
kf_traj(k,:) = kf_loc;
% show trajectory
subplot(122)
% cmpTraj(meas_traj, kf_traj, 'plot', 'xlim', axis_range(1:2), 'ylim', axis_range(3:4));
plotTraj(kf_traj, k, axis_range)
figtitle(data_item(1:end-1),'color','blue','linewidth',4,'fontsize',15);
drawnow
pause(show_delay)
end
%% save data
data_save_dir = [result_dir data_item 'ResData/'];
if ~exist(data_save_dir,'dir')
mkdir(data_save_dir)
end
save([data_save_dir 'traj.mat'], 'meas_traj', 'kf_traj')
disp(['result data saved to: ' data_save_dir])
%% -------------------------------------------------------
%% sub functions
% get KF default parameters
function param = getDefaultKFParameters(motion_type)
if nargin<1
motion_type = 'ConstantVelocity';
end
param.motionModel = motion_type;
param.initialLocation = 'Same as first detection';
if strcmp(motion_type, 'ConstantAcceleration')
param.initialEstimateError = 1E5 * ones(1, 3);
param.motionNoise= [25, 10, 1];
param.measurementNoise = 25;
elseif strcmp(motion_type, 'ConstantVelocity')
param.initialEstimateError = 1E5 * ones(1, 2);
param.motionNoise= [25, 10];
param.measurementNoise = 25;
else
error(['No assigned motion type - ' motion_type])
end
end
function plotBoundingbox(rect_p, rect_size, clr, lgd, frame_idx, axis_range)
plotcube(rect_size, rect_p, 0.1, clr,lgd)
title(['Frame #' num2str(frame_idx) ' - 3D detection']);
xlabel('X'), ylabel('Y'), zlabel('Z');
axis(axis_range);
view(3);
grid on
end
function plotTraj(traj, frame_idx, axis_range)
if size(traj, 2)==2
plot(traj(:,1),traj(:,2),'r-o','MarkerSize',4,'LineWidth',1.5)
elseif size(traj, 2)==3
plot3(traj(:,1),traj(:,2),traj(:,3),'r-o','MarkerSize',4,'LineWidth',1.5)
end
title(['Frame #' num2str(frame_idx) ' - XY trajectory']);
xlabel('X'), ylabel('Y'), zlabel('Z');
legend('KF Traj.')
axis(axis_range);
view(2);
grid on
end
% data format % *data format:* 1-X, 2-Y, 3-Z, 4-RANGE, 5-AZIMUTH, 6-ELEVATION, 7-DOPPLER, % 8-POWER, 9-POWER_VALUE, 10-TIMESTAMP_MS%% MOTracking - Multiple Object Tracking
clear, clc, close addpath(genpath('./utils'));%% env init
% path and data result_dir = './result/'; data_dir = './data/mmWave_radar_data/'; data_item = 'MOT/';%% param
end_frame = 10000; traj_dim = 2; % 2d/3d trajectorystart_frame = 1;
param_denoise.dpl_thr = 0.1;% denoise
param_det.minObjPoints = 30; param_det.DBSCAN_epsilon = 0.3; param_det.DBSCAN_MinPts = 30; % param_det.max_obj_count = 2;% detection
motion_type = 'ConstantVelocity'; % 'ConstantVelocity' | 'ConstantAcceleration' param_kf = getDefaultKFParameters(motion_type); param.initialEstimateError = [200 50]; param.motionNoise= [100 25]; param.measurementNoise = 600;% Kalman filter
axis_range = [-5, 5, 0, 20, -2, 5]; % axis_range = [-20, 20, 0, 20, -10, 10];% show
% ---- file info ---- datas = dir([data_dir data_item '*.txt']); data_names = {datas.name}; data_num = length(data_names); end_frame = min(data_num, end_frame); if start_frame>end_frame error("start frame over range") end%% denoise, cluster, KF_tracking
KF = []; % KF handle tracks = initializeTracks(); % object tracks nextId = 1; % next track id meas_traj = NaN(start_frame-1,traj_dim); % trajectory points % isDetected = false; % detected flag% ---- init ----
figure;
% ---- load data frame = importdata([data_dir data_item data_names{k}]); % ---- denoise ---- frame_clean = point_cloud_denoise(frame, param_denoise); disp(['clean points num: ' num2str(size(frame_clean,1))]) % ---- detect ---- [centroids, bboxes, obj_frame, obj_idx, obj_features] = getDetections(frame_clean,param_det); disp(['cluster count:' num2str(size(centroids,1))]) % ---- track ---- % predict new locations of last location (for cost calculation) tracks = predictNewLocationsOfTracks(tracks);for k = start_frame:end_frame
[assignments, unassignedTracks, unassignedDetections] = ... detectionToTrackAssignment(tracks, centroids, obj_frame, obj_idx, obj_features);% determine assignment of detection to tracks
tracks = updateAssignedTracks(tracks, assignments, centroids, bboxes); % update unassigned tracks; tracks = updateUnassignedTracks(tracks, unassignedTracks); % update track states tracks = updateTrackStates(tracks); % create new tracks(tracks); [tracks,nextId] = createNewTracks(tracks, unassignedDetections, ... centroids, bboxes, obj_features, param_kf, nextId, k);% undate assigned tracks
showTrackingResults(obj_frame, obj_idx, tracks, k, axis_range, data_item(1:end-1))% display track results
end
data_save_dir = [result_dir data_item 'ResData/']; if ~exist(data_save_dir,'dir') mkdir(data_save_dir) end save([data_save_dir 'track.mat'], 'meas_traj', 'tracks') disp(['result data saved to: ' data_save_dir])%% save data
%% sub functions % get KF default parameters function param = getDefaultKFParameters(motion_type) if nargin<1 motion_type = 'ConstantVelocity'; end param.motionModel = motion_type; param.initialLocation = 'Same as first detection'; if strcmp(motion_type, 'ConstantAcceleration') param.initialEstimateError = 1E5 * ones(1, 3); param.motionNoise= [25, 10, 1]; param.measurementNoise = 25; elseif strcmp(motion_type, 'ConstantVelocity') param.initialEstimateError = 1E5 * ones(1, 2); param.motionNoise= [25, 10]; param.measurementNoise = 25; else error(['No assigned motion type - ' motion_type]) end end%% -------------------------------------------------------
function showTrackingResults(obj_frame, obj_idx, tracks, frame_idx, axis_range, fig_name) % init % default bbox color clr = [1 0 0; 0 1 0; 0 0 1; 0 1 1; 1 0 1; 1 1 0]; show_delay = 0.0; % show 3d condition minVisibleCount = 5; % minimal consecutive appearing frame count maxInvisibleCount = 5; % maximal consecutive disappearing frame count % get normal tracks & effect tracks normal_track_ind = ... [tracks(:).totalVisibleCount] > minVisibleCount &... [tracks(:).consecutiveInvisibleCount] < maxInvisibleCount &... strcmp([tracks(:).state],"normal"); normalTracks = tracks(normal_track_ind);% show tracking results
[tracks(:).totalVisibleCount] > minVisibleCount &... ~strcmp([tracks(:).state],"noise"); effectTracks = tracks(effect_track_ind); clf(gcf) % clear figure before new display if ~ isempty(effectTracks) && ~isempty(obj_idx) % show 3d subplot(121) % scatter3(obj_frame(:,1),obj_frame(:,2),obj_frame(:,3),10,'filled') gscatter3(obj_frame(:,1),obj_frame(:,2),obj_frame(:,3),obj_idx,0.3*ones(10,3),[],10,'off')effect_track_ind = ...
plotBoundingbox(normalTracks(m).bbox(1:3), normalTracks(m).bbox(4:6), clr(normalTracks(m).id,:), ['obj' num2str(normalTracks(m).id)], axis_range) end legend % view(2) % 2D view for debugfor m = 1:length(normalTracks)
subplot(122) hold on for m = 1:length(effectTracks) plotTraj(effectTracks(m).traj_rec(:,1:2), axis_range, ['obj' num2str(effectTracks(m).id)], clr(effectTracks(m).id,:)) end hold off legend end % fig info figtitle([fig_name ' - Frame #' num2str(frame_idx)],'color','blue','linewidth',4,'fontsize',15); drawnow pause(show_delay)% show 2D
end
3 运行结果
4 参考文献
[1]赵睿智, and 丁云飞. "基于粒子群优化极限学习机的风功率预测." 上海电机学院学报 22.4():6.
[2]郭城, 刘新忠, and 苗宇. "基于粒子群优化极限学习机的轧制力预测." 冶金自动化 45.S01():4.
博主简介：擅长智能优化算法、神经网络预测、信号处理、元胞自动机、图像处理、路径规划、无人机等多种领域的Matlab仿真，相关matlab代码问题可私信交流。
部分理论引用网络文献，若有侵权联系博主删除。
