【回归预测-ELM预测】基于粒子群算法PSO优化极限学习机预测附matlab代码

matlab科研助手

已于 2022-08-02 17:06:01 修改

阅读量1.4k

点赞数 2

分类专栏：机器学习组合优化系列matlab完整代码文章标签：算法回归 matlab

于 2022-07-27 22:18:00 首次发布

本文链接：https://blog.csdn.net/matlab_dingdang/article/details/126024381

版权

机器学习组合优化系列matlab完整代码专栏收录该内容

448 篇文章 80 订阅

订阅专栏

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); %% 结果对比 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('期望输出','预测输出')

figure plot(1:length(T_test),T_test-T_sim,'r-*') xlabel('测试集样本编号') ylabel('绝对误差') title('ELM测试集预测误差')

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)); % 计算预测输出 Y = (H' * LW)';

%% 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;

   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

%% MOTracking - Multiple 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 = 'MOT/';

start_frame = 1; end_frame = 10000; traj_dim = 2; % 2d/3d trajectory

% denoise param_denoise.dpl_thr = 0.1;

% detection param_det.minObjPoints = 30; param_det.DBSCAN_epsilon = 0.3; param_det.DBSCAN_MinPts = 30; % param_det.max_obj_count = 2;

% Kalman filter motion_type = 'ConstantVelocity'; % 'ConstantVelocity' | 'ConstantAcceleration' param_kf = getDefaultKFParameters(motion_type); param.initialEstimateError = [200 50]; param.motionNoise = [100 25]; param.measurementNoise = 600;

% show axis_range = [-5, 5, 0, 20, -2, 5]; % axis_range = [-20, 20, 0, 20, -10, 10];

%% 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 tracks = initializeTracks(); % object tracks nextId = 1; % next track id meas_traj = NaN(start_frame-1,traj_dim); % trajectory points % 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(['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);

% determine assignment of detection to tracks [assignments, unassignedTracks, unassignedDetections] = ... detectionToTrackAssignment(tracks, centroids, obj_frame, obj_idx, obj_features);

% undate assigned 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);

% display track results showTrackingResults(obj_frame, obj_idx, tracks, k, axis_range, data_item(1:end-1))

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 'track.mat'], 'meas_traj', 'tracks') 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

% show tracking results 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);

effect_track_ind = ... [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')

for m = 1:length(normalTracks) 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 debug

% show 2D 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)

end