sg_flow_lim

SG_FLOW_LIM  main function to run transmission line capacity assignment
   [LINE_CAPACITY, ZPR, XX, RR, REF, BR_LINK] = ...
       SG_FLOW_LIM(LINK_IDS, A, ZPR, PL_SETTING, PG_SETTING, ...
                   BR_OVERLOAD, REFSYS_STAT)

   Input:
       link_ids - matrix (m x 2) indicating the branch terminal buses
       A - line admittance matrix
       Zpr - generated line impedances
       PL_setting - matrix (Nl x 2) indicating the load bus numbers and
           load setting at each load bus
       Pg_setting - matrix (Ng x 2) indicating the generation bus numbers
           and  generation dispatch at each generation bus
       br_overload - branch overload option (see SG_OPTIONS)
       refsys_stat - statistics of reference system

   Output:
       Line_capacity - matrix (M x 2) indicating the line capacities
           of each branch

0001 function [Line_capacity, Zpr, xx, rr, ref, br_link] = ...
0002     sg_flow_lim(link_ids, A, Zpr, PL_setting, Pg_setting, br_overload, refsys_stat)
0003 %SG_FLOW_LIM  main function to run transmission line capacity assignment
0004 %   [LINE_CAPACITY, ZPR, XX, RR, REF, BR_LINK] = ...
0005 %       SG_FLOW_LIM(LINK_IDS, A, ZPR, PL_SETTING, PG_SETTING, ...
0006 %                   BR_OVERLOAD, REFSYS_STAT)
0007 %
0008 %   Input:
0009 %       link_ids - matrix (m x 2) indicating the branch terminal buses
0010 %       A - line admittance matrix
0011 %       Zpr - generated line impedances
0012 %       PL_setting - matrix (Nl x 2) indicating the load bus numbers and
0013 %           load setting at each load bus
0014 %       Pg_setting - matrix (Ng x 2) indicating the generation bus numbers
0015 %           and  generation dispatch at each generation bus
0016 %       br_overload - branch overload option (see SG_OPTIONS)
0017 %       refsys_stat - statistics of reference system
0018 %
0019 %   Output:
0020 %       Line_capacity - matrix (M x 2) indicating the line capacities
0021 %           of each branch
0022 
0023 %   SynGrid
0024 %   Copyright (c) 2017-2018, Electric Power and Energy Systems (EPES) Research Lab
0025 %   by Hamidreza Sadeghian and Zhifang Wang, Virginia Commonwealth University
0026 %
0027 %   This file is part of SynGrid.
0028 %   Covered by the 3-clause BSD License (see LICENSE file for details).
0029 
0030 a_stab = refsys_stat.stab(1) ; b_stab = refsys_stat.stab(2) ; g_stab = refsys_stat.stab(3) ; d_stab = refsys_stat.stab(4);
0031 Overload_b = refsys_stat.Overload_b ; 
0032 mu_beta = refsys_stat.mu_beta ;
0033 Tab_2D_FlBeta = refsys_stat.Tab_2D_FlBeta ;
0034 M = length(link_ids);
0035 N = length (A(1,:));
0036 %% generate phi for calculation of X
0037 % phi=random('stable',a_stab,b_stab,g_stab,d_stab,M,1);
0038 phi = sg_random_stable(a_stab,b_stab,g_stab,d_stab,M,1);
0039 phi_pos1 = find(phi >90);
0040 phi_pos2 = find(phi <0);
0041 phi_pos = [phi_pos1;phi_pos2];
0042 if ~isempty(phi_pos)
0043     LL=1;
0044     while LL==1
0045         S_phi = sg_random_stable(a_stab,b_stab,g_stab,d_stab,length(phi_pos),1);
0046         phi(phi_pos) = S_phi;
0047         phi_pos1 = find(phi >90);
0048         phi_pos2 = find(phi <0);
0049         phi_pos = [phi_pos1;phi_pos2];
0050         if isempty(phi_pos)
0051             LL=0;
0052         else LL=1;
0053         end
0054     end
0055 end
0056 X_first = Zpr.*sind(phi);
0057 R = Zpr.*cosd(phi);
0058 %% DC power flow
0059 if N<300
0060     nl = 1;
0061 else
0062     nl = 2;
0063 end
0064 X_uni = ones(length(Zpr),1);
0065 [I M_ref] = max(Pg_setting(:,2));
0066 ref = Pg_setting(M_ref,1);
0067 %% Impedance
0068 bus = zeros(N, 17);
0069 bus(:, 1) = (1:N)';
0070 bus(:, 2) = 1;
0071 bus(ref, 2) = 3;
0072 bus(PL_setting(:,1), 3) = PL_setting(:,2);
0073 bus(:, 8) = 1;
0074 bus(:, 9) = 0;
0075 gen = zeros(length(Pg_setting), 25);
0076 gen(:, [1 2]) = Pg_setting;
0077 gen(:, 8) = Pg_setting(:,2)>0;
0078 gen(:, 6) = 1;
0079 branch = zeros(length(link_ids), 21);
0080 branch(:,[1 2]) = link_ids;
0081 branch(:, [3 4]) = [R,X_uni];
0082 branch(:, 11) = 1;
0083 branch(:, 14) = 0;
0084 branch(:, 16) = 0;
0085 %% form MATPOWER case struct
0086 mpt.version = '2';
0087 mpt.baseMVA = 100;
0088 mpt.bus = bus;
0089 mpt.gen = gen;
0090 mpt.branch = branch;
0091 mpopt = mpoption('out.all', 0, 'verbose',0);
0092 mpt = rundcpf(mpt,mpopt);
0093 
0094 flow = mpt.branch(:,14);
0095 theta_deg = mpt.bus(:,9);
0096 
0097 for kk = 1:nl
0098     Z_BrF = (1:length(link_ids))';
0099     Z_BrF (:,[2,3]) = link_ids;
0100     Z_BrF (:,[4,5]) = [Zpr , abs(flow)];
0101     Z_BrF_sort = sortrows(Z_BrF,-5);
0102     Zpr_sort = sort(Zpr); % sort in ascending order
0103     Zpr_asgn = zeros(length(Zpr),2);
0104     for j= 1:length(Zpr)
0105         Zpr_asgn (Z_BrF_sort(j,1),1) = Z_BrF_sort(j,1);
0106         Zpr_asgn (Z_BrF_sort(j,1),2) = Zpr_sort(j,1);
0107     end
0108     Zpr_new_sort = sortrows(Zpr_asgn,2);
0109     if N>1200
0110         as = 0.3;
0111         an = 0.8;
0112     else
0113         as = 0.2;
0114         an = 0.2;
0115     end
0116     
0117     N_swap = round(as.*M); % control variable for the number of swaps
0118     N_neib = round(an.*M);% control variable for the range of neigboring
0119     for k = 1:N_swap;
0120         xch_idx1 = floor(M.*rand);
0121         if xch_idx1 ==0
0122             xch_idx1 =1;
0123         end
0124         xch_disf = floor(N_neib.*rand);
0125         xch_dis = min( xch_disf , (M -xch_idx1));
0126         xch_idx2 = xch_idx1 + xch_dis;
0127         Zpr_new_sort ([xch_idx1 xch_idx2] , 2) = Zpr_new_sort ([xch_idx2 xch_idx1] , 2); % swap elements of matrix
0128     end
0129     Zpr_new = sortrows(Zpr_new_sort, 1);
0130     X = Zpr_new(:,2).*sind(phi);
0131     R = Zpr_new(:,2).*cosd(phi);
0132     mpt.branch(:,4) = X;
0133     mpopt = mpoption('out.all', 0, 'verbose',0);
0134     mpt = rundcpf(mpt,mpopt);
0135     flow = mpt.branch(:,14);
0136     theta_deg = mpt.bus(:,9);
0137     Zpr = Zpr_new(:,2);
0138 end
0139 mpopt = mpoption('out.all', 0, 'verbose',0);
0140 mpt = rundcpf(mpt,mpopt);
0141 flow = mpt.branch(:,14);
0142 theta_deg = mpt.bus(:,9);
0143 
0144 %% Line connection
0145 bus = zeros(N, 17);
0146 bus(:, 1) = (1:N)';
0147 bus(:, 2) = 1;
0148 bus(ref, 2) = 3;
0149 bus(PL_setting(:,1), 3) = PL_setting(:,2);
0150 bus(:, 8) = 1;
0151 bus(:, 9) = theta_deg;
0152 gen = zeros(length(Pg_setting), 25);
0153 gen(:, [1 2]) = Pg_setting;
0154 gen(:, 8) = Pg_setting(:,2)>0;
0155 gen(:, 6) = 1;
0156 branch = zeros(length(link_ids), 21);
0157 branch(:,[1 2]) = link_ids;
0158 branch(:, [3 4]) = [R,X];
0159 branch(:, 11) = 1;
0160 branch(:, 14) = flow;
0161 branch(:, 16) = -flow;
0162 %% form MATPOWER case struct
0163 mpp.version = '2';
0164 mpp.baseMVA = 100;
0165 mpp.bus = bus;
0166 mpp.gen = gen;
0167 mpp.branch = branch;
0168 mpopt = mpoption('out.all', 0, 'verbose',0);
0169 mpp = rundcpf(mpp,mpopt);
0170 theta = mpp.bus(:,9);
0171 d = repmat(theta,1,length(theta));
0172 dd = repmat(theta',length(theta),1);
0173 Bus_theta = dd-d;
0174 DD = [];
0175 DD = Bus_theta(:);
0176 Branch = link_ids;
0177 AA = A;
0178 XX = X;
0179 TT_target = 10.^(0.3196.*log10(N)+0.8324);
0180 cc =0;
0181 if max(abs(DD))> TT_target + 2
0182     while max(abs(DD)) > TT_target + 2
0183         cc = cc+1;
0184         theta = mpp.bus(:,9);
0185         d = repmat(theta,1,length(theta));
0186         dd = repmat(theta',length(theta),1);
0187         Bus_theta = dd-d;
0188         DD = [];
0189         DD = Bus_theta(:);
0190         Bus_delta = abs(triu(Bus_theta,0));
0191         [Rc,Cc] = ndgrid(1:size(Bus_delta,1),1:size(Bus_delta,2));
0192         [b,idx] = sort(Bus_delta(:),'descend');
0193         Bus_idx = [Rc(idx),Cc(idx)];
0194         Bus_phase_diff  = Bus_idx;
0195         Bus_phase_diff (:,3) = b;
0196         New_row = length(mpp.branch)+1;
0197         mpp.branch(New_row,1:2) = Bus_phase_diff(1,1:2);
0198         mpp.branch(New_row,3)  = 0.001 + rand*0.001;
0199         mpp.branch(New_row,4)  = 0.002 + rand*0.003;
0200         mpp.branch(New_row,11)  = 1;
0201         
0202         %% remove the branch
0203         Bus = mpp.bus;
0204         Branch = mpp.branch;
0205         Bus_degs=zeros(length(Bus(:,1)),1);
0206         for Bii=1:length(Bus(:,1))
0207             count=0;
0208             for Li=1:length(Branch)
0209                 if Branch(Li,1)==Bus(Bii,1)
0210                     count=count+1;
0211                 elseif Branch(Li,2)==Bus(Bii,1)
0212                     count=count+1;
0213                 end
0214             end
0215             Bus_degs(Bii)=count;
0216         end
0217         Br_SL = 0.8; % parameter to select the range for branch removing
0218         Branch_deg_data = (1:length(mpp.branch))';
0219         Branch_deg_data (:,[2,3,4]) = mpp.branch(:,[1,2,4]); % collect branch 'from' 'to' and 'X'
0220         Branch_deg_data (:,5) = Bus_degs(Branch_deg_data(:,2)); % 'from' bus degree
0221         Branch_deg_data (:,6) = Bus_degs(Branch_deg_data(:,3)); % 'to' bus degree
0222         Branch_deg_data = sortrows(Branch_deg_data,4); % sort rows based on the X
0223         Branch_select = Branch_deg_data (round(Br_SL*length(mpp.branch)):end,:); % select candidate branches
0224         if N<40
0225             select_node = 2;
0226         else
0227             select_node = 3;
0228         end
0229         Branch_select(Branch_select(:,5)<select_node,:) = []; % remove branches with 'from' node degree < select_node, (2 or 3)
0230         Branch_select(Branch_select(:,6)<select_node,:) = []; % remove branches with 'to' node degree < select_node
0231         Row_num = ceil(rand*length(Branch_select(:,1)));
0232         if Row_num > 0 && Row_num <= length(Branch_select(:,1))
0233             Remove_idx = Branch_select(Row_num,1); % select the branch number to remove by generating random number in rang of 'Branch_select' length
0234             Rem_branch = mpp.branch;
0235             mpp.branch(Remove_idx,:) = [];
0236         end
0237         mpopt = mpoption('out.all', 0, 'verbose',0);
0238         mpp = rundcpf(mpp,mpopt);
0239         if mpp.success == 0
0240             mpp.branch = Rem_branch;
0241             mpp = rundcpf(mpp,mpopt);
0242         end
0243         xx = mpp.branch(:,4);
0244         rr = mpp.branch(:,3);
0245         br_link = mpp.branch(:,[1,2]);
0246     end
0247     xx = mpp.branch(:,4);
0248     rr = mpp.branch(:,3);
0249 end
0250 cc;
0251 mpp = rundcpf(mpp,mpopt);
0252 xx = mpp.branch(:,4);
0253 rr = mpp.branch(:,3);
0254 Zpr = (rr.^2+xx.^2).^(0.5);
0255 br_link = mpp.branch(:,[1,2]);
0256 flow = mpp.branch(:,14);
0257 Flow=abs(mpp.branch(:,14));
0258 meanflowww = mean(Flow);
0259 thetadifff = max(mpp.bus(:,9)) - min(mpp.bus(:,9));
0260 line_indx = colon(1,length(mpp.branch))';
0261 Max_flow = max(Flow);
0262 M = length(mpp.branch);
0263 Norm_Flow = [line_indx,Flow./Max_flow];
0264 %% generate line capacity factor
0265 Overload_bet_n = round(M*Overload_b);
0266 beta = sg_exprnd(mu_beta,M,1);
0267 L_pos = find(beta >1);
0268 if ~isempty(L_pos)
0269     LL=1;
0270     while LL==1
0271         S_beta = sg_exprnd(mu_beta,length(L_pos),1);
0272         beta(L_pos) = S_beta;
0273         L_pos = find(beta >1);
0274         if isempty(L_pos)
0275             LL=0;
0276         else LL=1;
0277         end
0278     end
0279 end
0280 if br_overload && Overload_bet_n > 0
0281     Overload_bet = 1 + 0.2*rand(Overload_bet_n,1);
0282     ind = randi([1,M],Overload_bet_n,1);
0283     beta(ind) = Overload_bet;
0284 end
0285 beta = sort(beta); % to avoid generation of line capacity factor equal to zero
0286 beta_0 = find(beta == 0);
0287 beta_0_rand = 0.01 + 0.099.*rand(length(beta_0),1);
0288 beta(1:beta_0) = beta_0_rand;
0289 [Norm_Flow_setting,Corr_FL] = Assignment_Fl(Norm_Flow,beta,Tab_2D_FlBeta);
0290 
0291 % Line_capacity = [link_ids, Norm_Flow_setting(:,4).*Max_flow];
0292 Line_capacity = [br_link, Norm_Flow_setting(:,4).*Max_flow];
0293 Total_Trans_Cap = 10 ^(0.2875*(log10(N))^2-1.457*log10(N)+7.734);
0294 Line_capacity (Line_capacity(:,3)<=2,3) = 5 +100*rand;
0295 % Total_Trans_Cap1 = 11790*(N).^0.617;
0296 % check the total backbone transmission capacity
0297 % sometimes the total transmission capacity scaling cause to decrease the
0298 % capaceties, increasing the over rating operation of lines
0299 % Over_lines = find(abs(mpp.branch(:,14)) > Line_capacity(:,3)+0.001);
0300 % Over_lines_beta = abs(mpp.branch(Over_lines,14))./ Line_capacity(Over_lines,3);
0301 % k = 0;
0302 % while k ==0
0303 % if (sum(Line_capacity(:,3))>1.05*Total_Trans_Cap ||  sum(Line_capacity(:,3))<.95*Total_Trans_Cap)
0304 %     Line_capacity(:,3) = Line_capacity(:,3).*(Total_Trans_Cap./sum(Line_capacity(:,3)));
0305 %     k =0;
0306 % else
0307 %     k =1;
0308 % end
0309 % F_find = find(abs(mpp.branch(:,14)) >= Line_capacity(:,3));
0310 % if ~isempty(F_find)
0311 % Line_capacity(F_find,3) = abs(mpp.branch(F_find,14))+1;
0312 % end
0313 % end
0314 % Line_capacity(Over_lines,3) = abs(mpp.branch(Over_lines,14))./Over_lines_beta;
0315 mpp.branch(:,6) = Line_capacity(:,3);
0316 mpp = rundcpf(mpp,mpopt);
0317 Line_capacity(:,3);
0318 
0319 %%
0320 function [flow,theta_diff,ref,diff,theta_deg] = DFPF(N,A,X,link_ids,PL_setting,Pg_setting)
0321 
0322 bus=zeros(N,2);
0323 bus(PL_setting(:,1),1)=PL_setting(:,2);
0324 bus(Pg_setting(:,1),2)=Pg_setting(:,2);
0325 [~, ind_max]=max(bus(:,2));
0326 ref=ind_max(1);
0327 Pbus=(bus(:, 2) - bus(:, 1))./100;
0328 B=A'*diag(1./X)*A;
0329 pv=Pg_setting(:,1);
0330 pq=(1:N)';
0331 pq(pv) = [];
0332 pv(pv==ref)=[];
0333 Va0=zeros(N,1);
0334 % Eliminating the corresponding row & col. i.e. row & col. of reference bus
0335 [Va, ~] = dcpf(B, Pbus, Va0, ref, pv, pq);
0336 theta_deg=Va.*(180/pi);
0337 flow=diag(1./X)*A*Va;
0338 flow = flow.*100; % actual value (Base MVA 100)
0339 diff = max(theta_deg)-min(theta_deg);
0340 for i=1:length(flow)
0341     theta_diff(i)=theta_deg(link_ids(i,1))-theta_deg(link_ids(i,2));
0342 end
0343 
0344 %%
0345 function [Norm_Flow_setting,Corr_FL]=Assignment_Fl(Norm_Flow,beta,Tab_2D_FlBeta)
0346 
0347 M=length(beta);
0348 % calculate actual nnumber of each cell in 2D table nased on the total
0349 % numbers
0350 Tab_2D_M=round(Tab_2D_FlBeta.*M);
0351 % check the mismatch comes from "round". Add or subtract from maximum
0352 % number in the matrix
0353 if sum(sum(Tab_2D_M))<M
0354     [~,ind_max_tab]=max(Tab_2D_M(:));
0355     [I_row, I_col] = ind2sub(size(Tab_2D_M),ind_max_tab);
0356     Tab_2D_M(I_row,I_col)=Tab_2D_M(I_row,I_col)+(M-sum(sum(Tab_2D_M)));
0357 elseif sum(sum(Tab_2D_M))>M
0358     [~,ind_max_tab]=max(Tab_2D_M(:));
0359     [I_row, I_col] = ind2sub(size(Tab_2D_M),ind_max_tab);
0360     Tab_2D_M(I_row,I_col)=Tab_2D_M(I_row,I_col)-(sum(sum(Tab_2D_M))-M);
0361 end
0362 %% calculate target number of buses based on power flows and 2D table
0363 N_F=zeros(1,14); % matrix for total number of buses in each node degree category
0364 N_F = sum(Tab_2D_M,1);
0365 sum(sum(Tab_2D_M,1));
0366 sum(sum(Tab_2D_M,2));
0367 %% sort power flows
0368 Sort_norm_flow = sortrows(Norm_Flow,2);
0369 %% assign branches to the power flow categories
0370 Fl1=Sort_norm_flow(1:N_F(1,1),:);
0371 Sort_norm_flow(1:N_F(1,1),:)=[];
0372 
0373 Fl2=Sort_norm_flow(1:N_F(1,2),:);
0374 Sort_norm_flow(1:N_F(1,2),:)=[];
0375 
0376 Fl3=Sort_norm_flow(1:N_F(1,3),:);
0377 Sort_norm_flow(1:N_F(1,3),:)=[];
0378 
0379 Fl4=Sort_norm_flow(1:N_F(1,4),:);
0380 Sort_norm_flow(1:N_F(1,4),:)=[];
0381 
0382 Fl5=Sort_norm_flow(1:N_F(1,5),:);
0383 Sort_norm_flow(1:N_F(1,5),:)=[];
0384 
0385 Fl6=Sort_norm_flow(1:N_F(1,6),:);
0386 Sort_norm_flow(1:N_F(1,6),:)=[];
0387 
0388 Fl7=Sort_norm_flow(1:N_F(1,7),:);
0389 Sort_norm_flow(1:N_F(1,7),:)=[];
0390 
0391 Fl8=Sort_norm_flow(1:N_F(1,8),:);
0392 Sort_norm_flow(1:N_F(1,8),:)=[];
0393 
0394 Fl9=Sort_norm_flow(1:N_F(1,9),:);
0395 Sort_norm_flow(1:N_F(1,9),:)=[];
0396 
0397 Fl10=Sort_norm_flow(1:N_F(1,10),:);
0398 Sort_norm_flow(1:N_F(1,10),:)=[];
0399 N_F(1,11);
0400 Sort_norm_flow(1:N_F(1,11),:);
0401 Sort_norm_flow;
0402 Fl11=Sort_norm_flow(1:N_F(1,11),:);
0403 Sort_norm_flow(1:N_F(1,11),:)=[];
0404 
0405 Fl12=Sort_norm_flow(1:N_F(1,12),:);
0406 Sort_norm_flow(1:N_F(1,12),:)=[];
0407 
0408 Fl13=Sort_norm_flow(1:N_F(1,13),:);
0409 Sort_norm_flow(1:N_F(1,13),:)=[];
0410 
0411 Fl14=Sort_norm_flow(1:N_F(1,14),:);
0412 Sort_norm_flow(1:N_F(1,14),:)=[];
0413 
0414 
0415 %% calculate target number of line capacity factors
0416 N_B=zeros(1,14); % matrix for total number of buses in each node degree category
0417 N_B = sum(Tab_2D_M,2)';
0418 
0419 %% sort line capacity factors a
0420 Sort_beta=sort(beta);
0421 
0422 B1=Sort_beta(1:N_B(1,1),:);
0423 Sort_beta(1:N_B(1,1),:)=[];
0424 
0425 B2=Sort_beta(1:N_B(1,2),:);
0426 Sort_beta(1:N_B(1,2),:)=[];
0427 
0428 B3=Sort_beta(1:N_B(1,3),:);
0429 Sort_beta(1:N_B(1,3),:)=[];
0430 
0431 B4=Sort_beta(1:N_B(1,4),:);
0432 Sort_beta(1:N_B(1,4),:)=[];
0433 
0434 B5=Sort_beta(1:N_B(1,5),:);
0435 Sort_beta(1:N_B(1,5),:)=[];
0436 
0437 B6=Sort_beta(1:N_B(1,6),:);
0438 Sort_beta(1:N_B(1,6),:)=[];
0439 
0440 B7=Sort_beta(1:N_B(1,7),:);
0441 Sort_beta(1:N_B(1,7),:)=[];
0442 
0443 B8=Sort_beta(1:N_B(1,8),:);
0444 Sort_beta(1:N_B(1,8),:)=[];
0445 
0446 B9=Sort_beta(1:N_B(1,9),:);
0447 Sort_beta(1:N_B(1,9),:)=[];
0448 
0449 B10=Sort_beta(1:N_B(1,10),:);
0450 Sort_beta(1:N_B(1,10),:)=[];
0451 
0452 B11=Sort_beta(1:N_B(1,11),:);
0453 Sort_beta(1:N_B(1,11),:)=[];
0454 
0455 B12=Sort_beta(1:N_B(1,12),:);
0456 Sort_beta(1:N_B(1,12),:)=[];
0457 
0458 B13=Sort_beta(1:N_B(1,13),:);
0459 Sort_beta(1:N_B(1,13),:)=[];
0460 
0461 B14=Sort_beta(1:N_B(1,14),:);
0462 Sort_beta(1:N_B(1,14),:)=[];
0463 
0464 B15=Sort_beta(1:N_B(1,15),:);
0465 Sort_beta(1:N_B(1,15),:)=[];
0466 
0467 B16=Sort_beta(1:N_B(1,16),:);
0468 Sort_beta(1:N_B(1,16),:)=[];
0469 %% assign grouped power flows to the 2D table cells
0470 Fl_cell=cell(1,14);
0471 Fl_cell{1}=Fl1; Fl_cell{2}=Fl2; Fl_cell{3}=Fl3; Fl_cell{4}=Fl4; Fl_cell{5}=Fl5; Fl_cell{6}=Fl6; Fl_cell{7}=Fl7;
0472 Fl_cell{8}=Fl8; Fl_cell{9}=Fl9; Fl_cell{10}=Fl10; Fl_cell{11}=Fl11; Fl_cell{12}=Fl12; Fl_cell{13}=Fl13; Fl_cell{14}=Fl14;
0473 
0474 B_cell=cell(1,14);
0475 B_cell{1}=B1; B_cell{2}=B2; B_cell{3}=B3; B_cell{4}=B4; B_cell{5}=B5; B_cell{6}=B6; B_cell{7}=B7;
0476 B_cell{8}=B8; B_cell{9}=B9; B_cell{10}=B10; B_cell{11}=B11; B_cell{12}=B12; B_cell{13}=B13; B_cell{14}=B14; B_cell{15}=B15; B_cell{16}=B16;
0477 
0478 for ff=14:-1:1
0479     Fl_num=1;
0480     for bb=16:-1:1
0481         if Tab_2D_M(bb,ff)>0
0482             [samp_B,ind_B]=sg_datasample(B_cell{bb},Tab_2D_M(bb,ff),'replace',false);
0483             Fl_cell{ff}(Fl_num:(Fl_num + Tab_2D_M(bb,ff)-1),3)=samp_B;
0484             Fl_num=Fl_num+Tab_2D_M(bb,ff);
0485             B_cell{bb}(ind_B)=[];
0486         else
0487             Fl_cell{ff}(Fl_num:(Fl_num + Tab_2D_M(bb,ff)-1),3)=B_cell{bb}(1:Tab_2D_M(bb,ff));
0488         end
0489     end
0490 end
0491 % Flow_st is a martrix (M by 4) inclouds line numbers, normalized
0492 % power flow and line capacity factor
0493 Flow_st=[Fl_cell{1};Fl_cell{2};Fl_cell{3};Fl_cell{4};Fl_cell{5};Fl_cell{6};Fl_cell{7};Fl_cell{8};Fl_cell{9};Fl_cell{10};Fl_cell{11};Fl_cell{12};Fl_cell{13};Fl_cell{14}];
0494 Norm_Flow_setting = sortrows(Flow_st,1);
0495 Norm_Flow_setting(:,4) = Norm_Flow_setting(:,2)./ Norm_Flow_setting(:,3);
0496 Corr_FL = corr(Norm_Flow_setting(:,2), Norm_Flow_setting(:,3));