PSM怎么使用probit回归

转载

是大魔术师 2024-08-01 08:31:26

文章标签 PSM怎么使用probit回归算法回归 matlab 数据 文章分类 机器学习人工智能

1 内容介绍

风电功率预测为电网规划提供重要的依据,研究风电功率预测方法对确保电网在安全稳定运行下接纳更多的风电具有重要的意义.针对极限学习机(ELM)回归模型预测结果受输入参数影响的问题,现将粒子群优化算法(PSO)应用于ELM中,提出了一种基于粒子群优化极限学习机的风功率预测方法.该方法首先将数值天气预报信息(NWP)数据进行数据预处理,并构建出训练样本集,随后建立ELM模型,利用粒子群算法优化ELM中的输入权值和阈值,从而建立起基于NWP和PSO-ELM风功率预测模型.对华东地区3个不同装机容量的风场NWP数据进行实验.结果表明:该方法的预测精度高且稳定性能好,能够为风电场功率预测以及风电并网安全可靠性提供科学有效的参考依据.

风的随机性和波动性导致了风电功率的不平稳性，风电场的接入影响电网的稳定运行环境［１］。对风电场的输出功率进行精准预测，有助于电网人员根据风电场输出功率的变化调整发电计划，减少电网的备用容量以节约能源的消耗。因此，有效的风电功率预测方法可保障电网的稳定经济运行。目前，国内外研究人员做了大量的工作，研究方法主要包括物理方法、统计方法以及神经网络方法［２］。其中物理方法要求风机相关物理信息，建模复杂且精度不稳定。统计方法模型简单、数据单一，但预测的精度受时间的限制，预测时间越长精度越低，通常应用于超短期预测。神经网络方法泛化能力强，能够处理回归问题，适用于风功率预测。极限学习机（ＥｘｔｒｅｍｅＬｅａｒｎｉｎｇＭａｃｈｉｎｅ，ＥＬＭ）是一类基于前馈神经网络的算法［３］。ＥＬＭ因结构简单、学习效率高，被用以解决回归、聚类、分类等问题。ＥＬＭ算法的优化问题被很多学者研究，王浩等［４］利用遗传算法优化极限学习机，并将风电系统参数模糊化，从而提高预测模型的精度。龙浩［５］提出了加权极限学习机方法，旨在解决样本数据不均衡的问题。王文锦等［６］利用蜂群算法优化ＥＬＭ，旨在提高模型的稳定性能。王宏刚等［７］将蜻蜓算法（ＤｒａｇｏｎｆｌｙＡｌｇｏｒｉｔｈｍ，ＤＡ）分布式应用于ＥＬＭ，优化初始化输入权重和阈值的影响，有效提高了电能质量扰动识别率。风电功率的预测主要基于当日的数值天气预报信息（ＮｕｍｅｒｉｃａｌＷｅａｔｈｅｒＰｒｅｄｉｃｔｉｏｎ，ＮＷＰ）。ＮＷＰ包括风速、风向、温度、湿度以及气压［８］。实测ＮＷＰ数据在同一风功率下存在奇异值以及波动的问题。关于风功率ＮＷＰ数据的预处理问题，符杨等［９］针对ＮＷＰ数据不准确、爬坡事件频发等原因将ＮＷＰ数据分类，并对功率波动进行预测。杨茂等［１０］利用经验模态分解对风功率数据进行分解去噪重构，一定程度上减少了噪点对预测结果的影响。杨家然等［１１］利用模糊聚类的方法将原始信号进行分类，采用不同的模型组合预测。本文将粒子群优化算法（ＰａｒｔｉｃｌｅＳｗａｒｍＯｐｔｉｍｉｚａｔｉｏｎ，ＰＳＯ）和ＥＬＭ结合，欲提高传统ＥＬＭ预测模型的精度和稳定性能，从而为风电功率预测技术提供新的方法。

基本理论

１．１ＥＬＭ算法

ＥＬＭ因学习效率高被普遍应用于风电功率预测、变压器故障诊断、风机故障诊断等方面。该算法任意赋值输入层权重和阈值，训练过程中无需改变模型参数，仅设定隐层神经元个数，便可通过最小二乘法获得输出层权重。

PSM怎么使用probit回归_算法

PSM怎么使用probit回归_matlab_02

1．２ＰＳＯ－ＥＬＭ预测模型如前文所述，ＥＬＭ的初始输入权值和阈值是随机确定的，其训练的效果会受初始值影响。因此，采用ＰＳＯ优化ＥＬＭ的输入权重和阈值，可避免盲目性训练ＥＬＭ模型。ＰＳＯ算法优化ＥＬＭ的步骤中，先初始化ＰＳＯ参数，包括粒子群的规模、空间维度、惯性参数ｗ、学习因子ｃ１和ｃ２、迭代次数和最大速度ｖｍａｘ等。ＰＳＯ－ＥＬＭ预测模型是将每个粒子对应的输入权值和阈值代入ＥＬＭ预测模型中，将ＥＬＭ学习样本输出与实际输出的均方误差（ＭｅａｎＳｑｕａｒｅｄＥｒｒｏｒ，ＭＳＥ）作为ＰＳＯ的适应度。将粒子的当前适应度与最优适应度做对比，若比最优适应度小，说明当前输入权值和阈值所建立的ＥＬＭ模型进行预测产生的均方误差较小，则将当前适应度更新为最优适应度，将当前位置更新为Ｐｂ，否则保持最优适应度不变。同理比较适应度和全局适应度，更新Ｇｂ。当迭代次数达到最大值或适应度达到设定值时停止算法。ＰＳＯ优化得到的最优输入权值Ｗ和阈值ｂ后，代入ＥＬＭ模型中进行预测。具体步骤如图１所示。

PSM怎么使用probit回归_回归_03

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_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)';
% 随机生成训练集、测试集
 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 flagfigure;
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
 endfunction 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
 endfunction 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