sgvm_branch_perm

SGVM_BRANCH_PERM permute the branch properties of the mpc case
   [MPC, FLAG] = SGVM_BRANCH_PERM(MPC, OPT)

   Permuting the branch parameters simply means perumuting
   the rows of the mpc.branch matrix, while keeping
   columns F_BUS and T_BUS (1,2) fixed.

0001 function [mpc, flag] = sgvm_branch_perm(mpc, opt)
0002 %SGVM_BRANCH_PERM permute the branch properties of the mpc case
0003 %   [MPC, FLAG] = SGVM_BRANCH_PERM(MPC, OPT)
0004 %
0005 %   Permuting the branch parameters simply means perumuting
0006 %   the rows of the mpc.branch matrix, while keeping
0007 %   columns F_BUS and T_BUS (1,2) fixed.
0008 
0009 %   SynGrid
0010 %   Copyright (c) 2018, Power Systems Engineering Research Center (PSERC)
0011 %   by Eran Schweitzer, Arizona State University
0012 %
0013 %   This file is part of SynGrid.
0014 %   Covered by the 3-clause BSD License (see LICENSE file for details).
0015 
0016 %% some constants
0017 [F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, RATE_C, ...
0018     TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST, ...
0019     ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch;
0020 nl = size(mpc.branch,1);
0021 
0022 %% random permutation
0023 if QT > size(mpc.branch,2)
0024     % if no power flow results, create a random permutation
0025     perm = randperm(nl);
0026     mpc.branch(:,BR_R:ANGMAX) = mpc.branch(perm, BR_R:ANGMAX);
0027     return
0028 end
0029 
0030 %%
0031 if ~isfield(opt.vm.branchperm, 'overload_frac')
0032     overload_frac = 1;
0033 else
0034     overload_frac = opt.vm.branchperm.overload_frac;
0035 end
0036 Sf = sqrt(mpc.branch(:,PF).^2 + mpc.branch(:,QF).^2);
0037 St = sqrt(mpc.branch(:,PT).^2 + mpc.branch(:,QT).^2);
0038 %maximum apparent flow on branch in per-unit
0039 S  = max(Sf,St) / mpc.baseMVA;
0040 % line ratings in per-unit
0041 r  = mpc.branch(:,RATE_A) / mpc.baseMVA;
0042 
0043 R = mpc.branch(:, BR_R); % branch resistance
0044 X = mpc.branch(:, BR_X); % branch reactance
0045 
0046 perm = (1:nl).'; % initialize identity permutation
0047 %% find overloaded branches
0048 % indices of overloaded branches
0049 overloaded = find(S./r > overload_frac);
0050 % overloaded = find(mpc.softlims.RATE_A.overload > 1e-4);
0051 
0052 % sort overloades from LARGEST to SMALLEST
0053 % [~,tmp] = sort(mpc.softlims.RATE_A.overload(overloaded), 'descend');
0054 [~,tmp] = sort(S(overloaded) - r(overloaded), 'descend');
0055 overloaded = overloaded(tmp);
0056 
0057 %% swap loop
0058 % loop over overload branches and find possible elements to swap with.
0059 % Normal case: possible branches, B, have rating GREATER than current flow,
0060 %              their current flow is SMALLER than the overloded branche's rating, and
0061 %              they have not been used in a swap yet. From these, the branch whose
0062 %              impedance most closesly matches the overloaded branch in the L2 sense is chosen.
0063 % Alternative: The requirement that rating is GREATER than overload flow is dropped.
0064 %              In this case branches are chosen based on the closest
0065 %              rating to the flow.
0066 % Otherwise:   Skip this branch, no good swaps are possible.
0067 available  = true(nl,1);
0068 for k = 1:length(overloaded)
0069     if ~available(overloaded(k))
0070         % if branch overloaded(k) was already part of a swap
0071         continue
0072     end
0073     B = find( (r*overload_frac >= S(overloaded(k))) & (S <= r(overloaded(k))*overload_frac) & available );
0074     if isempty(B)
0075         if S(overloaded(k)) <= r(overloaded(k))
0076            % if branch is just heavily loaded by not overloaded: ignore
0077            continue
0078         end
0079         B = find( S <= r(overloaded(k)) & available ); % allow overload
0080         if isempty(B)
0081             % still no possiblities: ignore
0082             continue
0083         end
0084         % in this case, ignore impedance and just pick the best rating
0085         [~, tmp] = min(S(overloaded(k)) - r(B));
0086         idx = B(tmp);
0087         end
0088     %else
0089     % minimize the distance between the impedances as well as the
0090     % rating from the actual flow.
0091 %     [~, tmp1] = min((R(overloaded(k)) - R(B)).^2 + (X(overloaded(k)) - X(B)).^2); % impedance
0092     [~, tmp1] = sort((R(overloaded(k)) - R(B)).^2 + (X(overloaded(k)) - X(B)).^2); % impedance
0093     [~, tmp2] = sort( S(overloaded(k)) - r(B) ); % best rating for branch k
0094     [~, tmp3] = sort( S(B) - r(overloaded(k)) ); % best flow to recieve branch k's rating
0095     
0096     % for each entry in B sum up its location in the three tests. Then
0097     % take the minimum to pick the final swap candidate
0098     test_stat = full(sparse([tmp1; tmp2; tmp3], 1, repmat(1:length(B), 1, 3), length(B), 1));
0099     [~, tmp]  = min(test_stat);
0100     idx = B(tmp);
0101     %end
0102     perm(overloaded(k)) = idx; % overloaded branch gets branch idx's property
0103     perm(idx)           = overloaded(k); %branch idx gets overloaded(k)'s property
0104     available([overloaded(k), idx]) = false;
0105 end
0106 
0107 %% apply permutation
0108 if ~all(perm == (1:nl).')
0109     flag = 1;
0110     mpc.branch(:,BR_R:ANGMAX) = mpc.branch(perm, BR_R:ANGMAX);
0111 else
0112     % no permutation performed
0113     flag = 0;
0114 end
0115 %% end of function
0116 return
0117 
0118 %% optimization formulation
0119 % the optimization formulation currently does not work well due to the need
0120 % to minimize a non PSD quadratic term. It is kept here for potential
0121 % implementation with a local solver (e.g. IPOPT)
0122 %% create partitions for large cases
0123 vars = struct();
0124 max_block_size = 1000;
0125 nblock = ceil(nl/max_block_size);
0126 avgnum = nl/nblock;
0127 available_ids = true(nl,1);
0128 ptr = 0;
0129 for k = 1:nblock
0130     if k == nblock
0131         ids = find(available_ids);
0132         nelem = length(ids);
0133     elseif mod(k,2) == 0
0134         nelem = floor(avgnum);
0135         ids = find(available_ids);
0136         ids = ids(randperm(length(ids), nelem)); % sample WITHOUT replacement
0137         ids = sort(ids); %sort in ascending order
0138         available_ids(ids) = false;
0139     else
0140         nelem = ceil(avgnum);
0141         ids = find(available_ids);
0142         ids = ids(randperm(length(ids), nelem)); % sample WITHOUT replacement
0143         ids = sort(ids); %sort in ascending order
0144         available_ids(ids) = false;
0145     end
0146     vars.(['Pi', num2str(k)]) = struct('first', ptr+1, 'last', ptr + nelem^2, 'ids', ids);
0147     ptr = ptr + nelem^2;
0148 end
0149 %% problem setup
0150 prob = struct();
0151 
0152 vars.t = struct('first', ptr+1, 'last', ptr+nl);
0153 ptr    = ptr + nl;
0154 
0155 vars.tr = struct('first', ptr+1, 'last', ptr+nl);
0156 ptr    = ptr + nl;
0157 vars.tx = struct('first', ptr+1, 'last', ptr+nl);
0158 ptr    = ptr + nl;
0159 % vars = struct('Pi', struct('first',1, 'last', nl*nl),...
0160 %                           't', struct('first',nl*nl + 1, 'last', nl*(nl + 1)));
0161 
0162 % the indices in Pi are ROW-wise. For exmample if nl=3:
0163 % 1 2 3
0164 % 4 5 6
0165 % 7 8 9
0166 
0167 prob.A  = [];
0168 prob.l = [];
0169 prob.u = [];
0170 
0171 total_vars = ptr; %nl + nl*nl;
0172 
0173 Sf = sqrt(mpc.branch(:,PF).^2 + mpc.branch(:,QF).^2);
0174 St = sqrt(mpc.branch(:,PT).^2 + mpc.branch(:,QT).^2);
0175 %maximum apparent flow on branch in per-unit
0176 S  = max(Sf,St) / mpc.baseMVA;
0177 % line ratings in per-unit
0178 r  = mpc.branch(:,RATE_A) / mpc.baseMVA;
0179 
0180 R = mpc.branch(:, BR_R);
0181 X = mpc.branch(:, BR_X);
0182 %% Constraints
0183 % go one branch at a time
0184 % for k = 0:nl-1
0185 % for k = 1:nl
0186 %     for kk = 1:nblock
0187 %         v = ['Pi', num2str(kk)];
0188 %         ridx = find(vars.(v).ids == k);
0189 %         if ~isempty(ridx)
0190 %             nelem = length(vars.(v).ids);
0191 %             break
0192 %         end
0193 %     end
0194 %     % variable indices of row k in Pi_v
0195 %     Pirowidx = vars.(v).first + (nelem*(ridx-1):nelem*ridx-1);
0196 %     % variable indices of column k in Pi_v
0197 %     Picolidx = vars.(v).first + (ridx-1:nelem:nelem*nelem-1);
0198 % %     Pirowidx = (nl*k+1):(nl*(k+1)); %variable ids of row k+1 in Pi
0199 % %     Picolidx = k+1:nl:nl*nl; % variable ids for column k+1 in Pi
0200 %
0201 %   % -------- overload -----------
0202 %   % t_i + (Pi(i,:) * r - S_i) >= 0 ---> t_i + Pi(i,:)*r >= S_i
0203 % %     A  = sparse(1, [Pirowidx, vars.t.first+k], [r.', 1], 1, total_vars);
0204 %     A  = sparse(1, [Pirowidx, vars.t.first+k-1], [r(vars.(v).ids).', 1], 1, total_vars);
0205 % %     lb = S(k+1);
0206 %     lb = S(k);
0207 %   ub = Inf;
0208 %   prob = update_Albub(prob, A, lb, ub);
0209 %
0210 %   % ----- stochastic rows of Pi------
0211 %   % sum(Pi(k+1,:)) = 1
0212 %     % sum(Pi_v(ridx,:)) = 1
0213 %   A = sparse(1, Pirowidx, 1, 1, total_vars);
0214 %   lb = 1;
0215 %   ub = 1;
0216 %   prob = update_Albub(prob, A, lb, ub);
0217 %
0218 %   % ----- stochastic columns of Pi------
0219 %   % sum(Pi(:,k+1)) = 1
0220 %     % sum(Pi_v(:,k)) = 1
0221 %   A = sparse(1, Picolidx, 1, 1, total_vars);
0222 %   lb = 1;
0223 %   ub = 1;
0224 %   prob = update_Albub(prob, A, lb, ub);
0225 % end
0226 rows = cell(nblock, 1);
0227 for k = 1:nblock
0228     v = ['Pi', num2str(k)];
0229     nelem = length(vars.(v).ids);
0230     
0231     % -------- overload -----------
0232     % t_i + (Pi(i,:) * r - S_i) >= 0 ---> t_i + Pi(i,:)*r >= S_i
0233     tidx = sgvm_ensure_col_vect(vars.t.first - 1 + vars.(v).ids).';
0234     A  = sparse([1:nelem, reshape(ones(nelem,1)*(1:nelem), 1, nelem*nelem) ],...
0235          [tidx, vars.(v).first:vars.(v).last], ...
0236          [ones(1,nelem), repmat(r(vars.(v).ids).', 1, nelem)], nelem, total_vars);
0237     lb = S(vars.(v).ids);
0238     ub = Inf(nelem,1);
0239     prob = update_Albub(prob, A, lb, ub);
0240     
0241     % -------- branch resistance ----------
0242     % -tr_i <= Pi(i, :)*R - R <= tr_i leads to
0243     % Pi(i, :)*R + tr_i >= R
0244     % Pi(i, :)*R - tr_i <= R
0245     tidx = sgvm_ensure_col_vect(vars.tr.first - 1 + vars.(v).ids).';
0246     A = sparse([1:nelem, reshape(ones(nelem,1)*(1:nelem), 1, nelem*nelem) ],...
0247         [tidx, vars.(v).first:vars.(v).last],...
0248         [ones(1,nelem), repmat(R(vars.(v).ids).', 1, nelem)], nelem, total_vars);
0249     lb = R(vars.(v).ids);
0250     ub = Inf(nelem,1);
0251     prob = update_Albub(prob, A, lb, ub);
0252     
0253     tidx = sgvm_ensure_col_vect(vars.tr.first - 1 + vars.(v).ids).';
0254     A = sparse([1:nelem, reshape(ones(nelem,1)*(1:nelem), 1, nelem*nelem) ],...
0255         [tidx, vars.(v).first:vars.(v).last],...
0256         [-ones(1,nelem), repmat(R(vars.(v).ids).', 1, nelem)], nelem, total_vars);
0257     lb = -Inf(nelem,1);
0258     ub = R(vars.(v).ids);
0259     prob = update_Albub(prob, A, lb, ub);
0260     
0261     % -------- branch reactance ----------
0262     % -tx_i <= Pi(i, :)*X - X_i <= tx_i leads to
0263     % Pi(i, :)*X + tx_i >= X_i
0264     % Pi(i, :)*X - tx_i <= X_i
0265     tidx = sgvm_ensure_col_vect(vars.tx.first - 1 + vars.(v).ids).';
0266     A = sparse([1:nelem, reshape(ones(nelem,1)*(1:nelem), 1, nelem*nelem) ],...
0267         [tidx, vars.(v).first:vars.(v).last],...
0268         [ones(1,nelem), repmat(X(vars.(v).ids).', 1, nelem)], nelem, total_vars);
0269     lb = X(vars.(v).ids);
0270     ub = Inf(nelem,1);
0271     prob = update_Albub(prob, A, lb, ub);
0272     
0273     tidx = sgvm_ensure_col_vect(vars.tx.first - 1 + vars.(v).ids).';
0274     A = sparse([1:nelem, reshape(ones(nelem,1)*(1:nelem), 1, nelem*nelem) ],...
0275         [tidx, vars.(v).first:vars.(v).last],...
0276         [-ones(1,nelem), repmat(X(vars.(v).ids).', 1, nelem)], nelem, total_vars);
0277     lb = -Inf(nelem,1);
0278     ub = X(vars.(v).ids);
0279     prob = update_Albub(prob, A, lb, ub);
0280     
0281     % ----- stochastic rows of Pi------
0282     % sum(Pi_v(ridx,:)) = 1
0283     A = sparse(reshape(ones(nelem,1)*(1:nelem), 1, nelem*nelem),...
0284                vars.(v).first:vars.(v).last, 1, nelem, total_vars);
0285     lb = ones(nelem,1);
0286     ub = ones(nelem,1);
0287     prob = update_Albub(prob, A, lb, ub);
0288     
0289     % ----- stochastic columns of Pi------
0290     % sum(Pi_v(:,k)) = 1
0291     A  = sparse(reshape((1:nelem).'*ones(1,nelem), 1, nelem*nelem),...
0292                 vars.(v).first:vars.(v).last, 1, nelem, total_vars);
0293     lb = ones(nelem,1);
0294     ub = ones(nelem,1);
0295     prob = update_Albub(prob, A, lb, ub);
0296     
0297     % sum(Pi_v*r_v) = sum(r_v)
0298     % note all the first column of Pi_v gets r_v(1), second column r_v(2), etc.
0299     A = sparse(1, vars.(v).first:vars.(v).last, ...
0300         repmat(r(vars.(v).ids).', 1, nelem), 1, total_vars);
0301     lb = sum(r(vars.(v).ids).');
0302     ub = sum(r(vars.(v).ids).');
0303     prob = update_Albub(prob, A, lb, ub);
0304     
0305 %     % sum(Pi_v^T*r_v) = sum(r_v)
0306 %     % note all the first column of Pi_v gets r_v(1), second column r_v(2), etc.
0307 %     A = sparse(1, vars.(v).first:vars.(v).last, ...
0308 %         reshape(ones(nelem,1)*r(vars.(v).ids).', 1, nelem*nelem), 1, total_vars);
0309 %     lb = sum(r(vars.(v).ids).');
0310 %     ub = sum(r(vars.(v).ids).');
0311 %     prob = update_Albub(prob, A, lb, ub);
0312     
0313     if k == 1
0314         rows{1} = [1, size(prob.A,1)];
0315     else
0316         rows{k} = [rows{k-1}(2)+1, size(prob.A, 1)];
0317     end
0318 end
0319 A = [];
0320 % % sum(Pi*r) = sum(r)
0321 % A  = sparse(1, vars.Pi.first:vars.Pi.last, repmat(r.',1,nl),1,total_vars);
0322 % lb = sum(r);
0323 % ub = sum(r);
0324 % prob = update_Albub(prob, A, lb, ub);
0325 
0326 %% variable limits
0327 prob.xmin = zeros(total_vars,1);
0328 prob.xmax = Inf(total_vars,1);
0329 % prob.xmax(vars.Pi.first:vars.Pi.last) = 1;
0330 prob.xmax(vars.Pi1.first:vars.(['Pi', num2str(nblock)]).last) = 1;
0331 %% linear cost
0332 
0333 % sum(t) - mu*sum(Pi_ii);
0334 % mu = 1e-4/mpc.baseMVA;
0335 % prob.f_fcn = @(x)obj_fnc(x, lambda, nblock, vars);
0336 % prob.x0 = zeros(total_vars, 1);
0337 % the indices of diagonal elements of Pi are 1, nl+2, 2*nl+3 etc.
0338 % Pidiagidx = (nl*(0:nl-1)) + 1:nl;
0339 
0340 % the diagonal indices of each Pi_v are first+ (0, nelem-1+2, 2*nelem -
0341 % 1+3) etc. = first + nelem*(0:nelem-1) + 0:nelem-1
0342 
0343 % Pidiagidx = zeros(1,nl);
0344 % ptr = 0;
0345 % for k = 1:nblock
0346 %     v = ['Pi', num2str(k)];
0347 %     nelem = length(vars.(v).ids);
0348 %     Pidiagidx(ptr+1:ptr+nelem) = vars.(v).first + nelem*(0:nelem-1) + (0:nelem-1);
0349 %     ptr = ptr + nelem;
0350 % end
0351 % prob.x0(Pidiagidx) = 1;
0352 % prob.c = sparse([vars.t.first:vars.t.last, Pidiagidx], 1, ...
0353 % [ceil(avgnum)*ones(1,nl), -mu*ones(1,nl)], total_vars, 1);
0354 
0355 prob.c = sparse([vars.t.first:vars.t.last,...
0356                  vars.tr.first:vars.t.last,...
0357                  vars.tx.first:vars.t.last], 1, 1, total_vars, 1);
0358 
0359 %% quadratic cost
0360 % sum of squares of elements in Pi_v weighted by lambda
0361 % prob.H = sparse(vars.Pi1.first:vars.(['Pi', num2str(nblock)]).last, ...
0362 %                 vars.Pi1.first:vars.(['Pi', num2str(nblock)]).last, ...
0363 %                 -lambda, total_vars,total_vars);
0364 
0365 %% solve
0366 prob.opt.verbose = opt.vm.branchperm.verbose;
0367 % for k = 1:nblock
0368 %     v = ['Pi', num2str(k)];
0369 %     nelem = length(vars.(v).ids);
0370 %     tidx = sgvm_ensure_col_vect(vars.t.first - 1 + vars.(v).ids).';
0371 %     idxtmp = [vars.(v).first:vars.(v).last, tidx];
0372 %     varstmp = struct('Pi1', struct('first', 1, 'last', nelem*nelem, 'ids', vars.(v).ids),...
0373 %         't', struct('first', nelem*nelem+1, 'last', nelem*nelem + nelem));
0374 %     probtmp = struct('A', prob.A(rows{k}(1):rows{k}(2),idxtmp), ...
0375 %         'l', prob.l(rows{k}(1):rows{k}(2)), 'u', prob.u(rows{k}(1):rows{k}(2)),...
0376 %         'x0', prob.x0(idxtmp), 'f_fcn', @(x)obj_fnc(x, lambda, 1, varstmp),...
0377 %         'opt', prob.opt);
0378 %     [x, f, exitflag, output, lambda] = mips(probtmp);
0379 % end
0380 [x, f, eflag, output, lambda] = qps_matpower(prob);
0381 
0382 if eflag
0383     % for each row find the largest entry
0384     perm = zeros(nl,1);
0385     for k = 1:nblock
0386         v = ['Pi', num2str(k)];
0387         nelem = length(vars.(v).ids);
0388         ptmp  = sgvm_extract_perm(reshape(x(vars.(v).first:vars.(v).last), nelem, nelem));
0389         perm(vars.(v).ids) = vars.(v).ids(ptmp);
0390     end
0391 %   for k = 0:nl-1
0392 %         for kk = 1:nblock
0393 %             v = ['Pi', num2str(kk)];
0394 %             ridx = find(vars.(v).ids == k+1); % row in Pi_v corresponding to branch k+1
0395 %             if ~isempty(ridx)
0396 %                 nelem = length(vars.(v).ids);
0397 %                 break
0398 %             end
0399 %         end
0400 % %         Pirowidx = (nl*k+1):(nl*(k+1)); %variable ids of row k+1 in Pi
0401 % %         [~,perm(k+1)] = max(x(Pirowidx));
0402 %         Pirowidx = vars.(v).first + (ridx-1)*nelem + (0:nelem-1); % variable indices of row ridx  in Pi_v
0403 %         [~, tmp] = max(x(Pirowidx)); % maximum value index along this row
0404 %         perm(k+1) = vars.(v).ids(tmp); % maps property at branch with maximum value of Pi to branch k+1.
0405 %   end
0406     % double check that all integers are found. this should be
0407     % the case as long as no row has uniform values 1/nl
0408     if length(unique(perm)) ~= nl
0409         error(['sgvm_branch_perm: binarzation of Pi by selecting maximum entry in each row ',...
0410         'failed to return a valid permutation.'])
0411     end
0412 else
0413     error('sgvm_branch_perm: optimization failed to converge.')
0414 end
0415 
0416 %% apply permutation
0417 mpc.branch(:,BR_R:ANGMAX) = mpc.branch(perm, BR_R:ANGMAX);
0418 
0419 function prob = update_Albub(prob, A, lb, ub)
0420 % simple utility function for updating the constraint and bounds arrays.
0421 
0422 prob.A  = [prob.A; A];
0423 prob.l  = [prob.l; lb];
0424 prob.u  = [prob.u; ub];
0425 
0426 function [f, df, d2f] = obj_fnc(x, lambda, nblock, vars)
0427 
0428 v = ['Pi', num2str(nblock)];
0429 f = sum(x(vars.t.first:vars.t.last)) + ...
0430     -lambda*sum(x(vars.Pi1.first:vars.(v).last).^2);
0431 if nargout > 1          %% gradient is required
0432     df = [ ones(vars.t.last-vars.t.first + 1, 1); ...
0433            -2*lambda*x(vars.Pi1.first:vars.(v).last)];
0434     if nargout > 2    %% Hessian is required
0435         d2f = sparse(vars.Pi1.first:vars.(v).last,...
0436             vars.Pi1.first:vars.(v).last, -2*lambda, length(x), length(x));
0437     end
0438 end