(文章复现)基于变异粒子群算法的主动配电网故障恢复策略matlab代码

%% 孤岛划分的程序 %% 步骤1：相关参数 mpc = IEEE33; % 配电网数据 Nb = 33; % 节点数 Nt = 24; % 时段数 Y_DG = [6, 8, 12, 23, 30]; % 光储系统接入位置 wf = 0.01*ones(1,Nb); % 负荷等级系数 wf([6, 12, 23, 30]) = 1; % 一级负荷 wf([1,5,8,14,18,21,22,27,29,32]) = 0.1; % 二级负荷 Ct = zeros(Nb, Nt); % 负荷时变需求系数 Ct([6, 12, 23, 30], :) = 10; % 一级负荷的需求度恒为10 Ct(setdiff(1:Nb, [6, 12, 23, 30]), :) = ones(Nb - 4 , 1)*10*[7.77960493860530,5.32481765803490,5.04294379929200,5.19839957470577,2.98814682221510,3.42926158320840,7.20423714015130,10.6852586809857,12.5476776005102,15.5211669393227,24.8696812348447,33.2302134035133,24.2862651332590,24.7988707426448,36.3317495052460,27.4948240683687,17.8808252142146,14.0969740403730,18.9045593211401,20.4888032159577,22.2151328923961,19.7627571226032,21.4368039620930,17.6810254063144]/36.3317; Ft = Ct.*(wf'*ones(1, Nt)); % 优先恢复系数 % 光伏出力 P_PV = [425, 200, 480, 120, 600; 700, 500, 720, 620, 720]; PL = [3305 , 4870]; % 负荷总量 P_load = mpc.bus(:,3)/3.715*PL; % 各节点负荷 a_cut = zeros(1,Nb); % 削减系数 C_load1 = [3,13,20,24]; % 完全可控的负荷集合 a_cut(C_load1) = 1; % 完全可控 C_load2 = [2,7,10]; % 50%可控的负荷集合 a_cut(C_load2) = 0.5; % 50%可控 %% 步骤2：读取故障支路 fault_branch = [5, 25, 20]; % 故障支路 %% 步骤3：确定故障时刻 time0 = [7 , 15]; time = time0(index); % 当前时刻 if index ~= 1 && index ~= 2 error('输入有误请重新输入！！！') end %% 步骤4：确定优先恢复集合 F_Rload = {[6, 12, 15, 23, 28, 30], [6, 12, 23, 30, 32]}; %% 步骤5：采用BFS确定孤岛包含的负荷 island_node = cell(1, length(Y_DG)); % 各个孤岛包含的负荷节点 island_node_all = []; for k = 1:length(Y_DG) [island_node{k} , ~] = BFS(Y_DG(k) , P_PV(index,k)); island_node_all = [island_node_all, setdiff(island_node{k}, island_node_all)]; end % 将交叉的孤岛融合 index_island = zeros(length(Y_DG)); % 孤岛是否交叉的标志，例如index_island(1,2)=1，表示孤岛1和2有交叉 for k = 1:length(Y_DG) for kk = k + 1 : length(Y_DG) if ~isempty(intersect(island_node{k} , island_node{kk})) index_island(k , kk) = 1; % 如果孤岛出现交叉情况就把标志为置为1 end end end island_node_cross = cell(1, length(Y_DG) - sum(index_island(:))); % 将交叉的孤岛融合后的孤岛划分方案 island_cross_DG = cell(1, length(Y_DG) - sum(index_island(:))); % 融合后孤岛中包含的光伏节点 island_node_cross{1} = island_node{1}; island_cross_DG{1} = 1; kk = 2; for k = 2:length(Y_DG) if sum(index_island(k , :)) == 0 && sum(index_island(: , k)) == 0 island_node_cross{kk} = island_node{k}; island_cross_DG{kk} = k; kk = kk + 1; else cross_num = [find(index_island(k , :)) , find(index_island(: , k))]; for kkk = 1:length(cross_num) island_node_cross{min(cross_num)} = union(island_node_cross{min(cross_num)} , island_node{k}); island_cross_DG{min(cross_num)} = union(island_cross_DG{min(cross_num)} , k); end end end %% 步骤6：判断优先恢复集中的负荷是否被恢复，以未恢复的负荷所在节点作为根节点，利用 DFS 法将其接入最近的光伏电站 island_num = length(island_node_cross); % 融合后的孤岛数量 for k = 1:length(F_Rload{index}) if ~ismember(F_Rload{index}(k), island_node_all) for kk = 1:length(Y_DG) [island_node0{kk}, distance0(kk)] = dfs_shortest_paths(F_Rload{index}(k), Y_DG(kk)); end [~, min_index] = min(distance0); for kk = 1:island_num if ismember(min_index , island_cross_DG{kk}) island_node_cross{kk} = union(island_node0{min_index}, island_node_cross{kk}); end end end end %% 步骤7：削减可控负荷 island_power = zeros(1, island_num); % 孤岛内光储系统总出力 island_load = zeros(1, island_num); % 孤岛内负荷总量 for k = 1:island_num island_power(k) = sum(P_PV(index, island_cross_DG{k})); % 孤岛内光储系统总出力 island_load(k) = sum(P_load(island_node_cross{k} , index)); % 孤岛内负荷总量 if island_power(k) < island_load(k) % 如果光储系统无法满足负荷需求，需要削减负荷 for kk = 1:length(island_node_cross{k}) if ismember(island_node_cross{k}(kk), C_load1) island_load(k) = island_load(k) - P_load(island_node_cross{k}(kk) , index); island_node_cross{k} = setdiff(island_node_cross{k} , island_node_cross{k}(kk)); if island_power(k) >= island_load(k) break end elseif ismember(island_node_cross{k}(kk), C_load2) island_load(k) = island_load(k) - 0.5*P_load(island_node_cross{k}(kk) , index); if island_power(k) >= island_load(k) break end end end end if island_power(k) < island_load(k) % 如果削减了可控负荷后还是无法满足需求，则按照重要度等级进行切负荷操作 [~ , Ft_index] = sort(Ft(:, time)); % 把负荷按重要度等级从大到小排列 for kk = 1:Nb if ismember(kk, island_node_cross{k}) island_load(k) = island_load(k) - P_load(kk , index); island_node_cross{k} = setdiff(island_node_cross{k} , kk); if island_power(k) >= island_load(k) break end end end end end %% 步骤8：输出孤岛划分方案 for k = 1:island_num disp(['**********************************第', num2str(k), '个孤岛**********************************']) disp(['包含光储系统' , num2str(island_cross_DG{k}) , '，光储系统总功率为：' , num2str(island_power(k)) , 'kW']) disp(['包含节点：' , num2str(island_node_cross{k}) , '，孤岛内负荷总需求为：' , num2str(island_load(k)) , 'kW']) end