Home > matpower7.1 > extras > syngrid > lib > sgvm_deltainjection2perm.m

sgvm_deltainjection2perm

PURPOSE ^

SGVM_DELTAINJECTION2PERM convert change in injection to a permutation.

SYNOPSIS ^

function perm = sgvm_deltainjection2perm(Pbus, Qbus, opt)

DESCRIPTION ^

SGVM_DELTAINJECTION2PERM convert change in injection to a permutation.
   PERM = SGVM_DELTAINJECTION2PERM(PBUS, QBUS, OPT)

   Determine a permutation of the injection vector to best achieve the
   new desired vector. Pbus and Qbus are structures with fields 'old'
   and 'new'. The goal is to permute 'old' so it looks as close to 'new'
   as possible.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 function perm = sgvm_deltainjection2perm(Pbus, Qbus, opt)
0002 %SGVM_DELTAINJECTION2PERM convert change in injection to a permutation.
0003 %   PERM = SGVM_DELTAINJECTION2PERM(PBUS, QBUS, OPT)
0004 %
0005 %   Determine a permutation of the injection vector to best achieve the
0006 %   new desired vector. Pbus and Qbus are structures with fields 'old'
0007 %   and 'new'. The goal is to permute 'old' so it looks as close to 'new'
0008 %   as possible.
0009 
0010 %   SynGrid
0011 %   Copyright (c) 2018, Power Systems Engineering Research Center (PSERC)
0012 %   by Eran Schweitzer, Arizona State University
0013 %
0014 %   This file is part of SynGrid.
0015 %   Covered by the 3-clause BSD License (see LICENSE file for details).
0016 
0017 nb = length(Pbus.old);
0018 % loopsolve = opt.vm.nodeperm.loopsolve;
0019 % some thresholding
0020 for fld = {'new', 'old'}
0021     Pbus.(fld{:})(abs(Pbus.(fld{:})) < 1e-6) = 0;
0022     Qbus.(fld{:})(abs(Qbus.(fld{:})) < 1e-6) = 0;
0023 end
0024 perm = zeros(nb,1); %initialize permutation vector
0025 %% non optimization approach
0026 err = sqrt( (Pbus.new - Pbus.old).^2 + (Qbus.new - Qbus.old).^2 );
0027 [~, idx] = sort(err, 'descend');
0028 available = (1:nb).';
0029 for k = idx.'
0030     [~, tmp] = min(abs(Pbus.new(k) - Pbus.old(available)) + abs(Qbus.new(k) - Qbus.old(available)));
0031     perm(k) = available(tmp);
0032     available(tmp) = [];
0033 end
0034 %% verify permutation
0035 if ~all(sort(perm) == (1:nb).')
0036     error('deltainjection2per: failed to return a valid permutation.')
0037 end
0038 %% end of function
0039 return
0040 %% create partitions for large cases
0041 vars = struct();
0042 max_block_size = 1000;
0043 nblock = ceil(nb/max_block_size);
0044 avgnum = nb/nblock;
0045 available_ids = true(nb,1);
0046 ptr = 0;
0047 for k = 1:nblock
0048     if k == nblock
0049         ids = find(available_ids); % automatically sorted
0050         nelem = length(ids);
0051     elseif mod(k,2) == 0
0052         nelem = floor(avgnum);
0053         ids = find(available_ids);
0054         ids = ids(randperm(length(ids), nelem)); % sample WITHOUT replacement
0055         ids = sort(ids); %sort in ascending order
0056         available_ids(ids) = false;
0057     else
0058         nelem = ceil(avgnum);
0059         ids = find(available_ids);
0060         ids = ids(randperm(length(ids), nelem)); % sample WITHOUT replacement
0061         ids = sort(ids); %sort in ascending order
0062         available_ids(ids) = false;
0063     end
0064     vars.(['Pi', num2str(k)]) = struct('first', ptr+1, 'last', ptr + nelem^2, 'ids', ids);
0065     ptr = ptr + nelem^2;
0066 end
0067 %% problem setup
0068 
0069 if ~loopsolve
0070     vars.tp = struct('first', ptr+1, 'last', ptr+nb);
0071     ptr = ptr + nb;
0072     vars.tq = struct('first', ptr+1, 'last', ptr+nb);
0073     ptr = ptr + nb;
0074     total_vars = ptr;
0075     prob = prob_init();
0076     % initial point
0077     prob.x0 = zeros(total_vars,1);
0078     prob.x0(vars.tp.first:vars.tp.last) = abs(Pbus.new - Pbus.old);
0079     prob.x0(vars.tq.first:vars.tq.last) = abs(Qbus.new - Qbus.old);
0080 end
0081 % vars = struct('Pi',struct('first',1, 'last', nb*nb),...
0082 %               'tp', struct('first', nb*nb+1, 'last', nb*(nb+1)),...
0083 %             'tq', struct('first',nb*(nb+1) + 1,'last', nb*(nb +2)));
0084 
0085 % the indices in Pi are ROW-wise. For exmample if nb=3:
0086 % 1 2 3
0087 % 4 5 6
0088 % 7 8 9
0089 
0090 % total_vars = 2*nb + nb*nb;
0091 
0092 %% Constraints
0093 
0094 for k = 1:nblock
0095     v = ['Pi', num2str(k)];
0096     nelem = length(vars.(v).ids);
0097     if loopsolve
0098         prob = prob_init();
0099         varsloc = struct('Pi', struct('first', 1, 'last', nelem*nelem),...
0100                   'tp', struct('first',nelem*nelem+1, 'last', nelem*(nelem + 1)),...
0101                   'tq', struct('first',nelem*(nelem + 1)+1, 'last', nelem*(nelem + 2)));
0102         total_vars_loc = nelem*(nelem + 2);
0103     end
0104     % --real power P ---------------
0105     % -tp <= Pnew - Pi*Pold <= tp becomes
0106     % -Inf <= Pi*Pold - tp <= Pnew
0107     % Pnew <= Pi*Pold + tp <= Inf
0108     if ~loopsolve
0109         tidx = sgvm_ensure_col_vect(vars.tp.first - 1 + vars.(v).ids).';
0110         A = sparse([1:nelem, reshape(ones(nelem,1)*(1:nelem), 1, nelem*nelem) ],...
0111             [tidx, vars.(v).first:vars.(v).last], ...
0112             [-ones(1,nelem), repmat(Pbus.old(vars.(v).ids).', 1, nelem)], ...
0113             nelem, total_vars);
0114     else
0115         tidx = varsloc.tp.first:varsloc.tp.last;
0116         A = sparse([1:nelem, reshape(ones(nelem,1)*(1:nelem), 1, nelem*nelem) ],...
0117             [tidx, varsloc.Pi.first:varsloc.Pi.last], ...
0118             [-ones(1,nelem), repmat(Pbus.old(vars.(v).ids).', 1, nelem)], ...
0119             nelem, total_vars_loc);
0120     end
0121     lb = -Inf(nelem,1);
0122     ub = Pbus.new(vars.(v).ids);
0123     prob = update_Albub(prob, A, lb, ub);
0124 
0125     if ~loopsolve
0126         A = sparse([1:nelem, reshape(ones(nelem,1)*(1:nelem), 1, nelem*nelem) ],...
0127             [tidx, vars.(v).first:vars.(v).last], ...
0128             [ones(1,nelem), repmat(Pbus.old(vars.(v).ids).', 1, nelem)], ...
0129             nelem, total_vars);
0130     else
0131         A = sparse([1:nelem, reshape(ones(nelem,1)*(1:nelem), 1, nelem*nelem) ],...
0132             [tidx, varsloc.Pi.first:varsloc.Pi.last], ...
0133             [ones(1,nelem), repmat(Pbus.old(vars.(v).ids).', 1, nelem)], ...
0134             nelem, total_vars_loc);
0135     end
0136     lb = Pbus.new(vars.(v).ids);
0137     ub = Inf(nelem,1);
0138     prob = update_Albub(prob, A, lb, ub);
0139 
0140     % --reactive power Q ---------------
0141     % -tq <= Qnew - Pi*Qold <= tq becomes
0142     % -Inf <= Pi*Qold - tq <= Qnew
0143     % Qnew <= Pi*Qold + tq <= Inf
0144     if ~loopsolve
0145         tidx = sgvm_ensure_col_vect(vars.tq.first - 1 + vars.(v).ids).';
0146         A = sparse([1:nelem, reshape(ones(nelem,1)*(1:nelem), 1, nelem*nelem) ],...
0147             [tidx, vars.(v).first:vars.(v).last], ...
0148             [-ones(1,nelem), repmat(Qbus.old(vars.(v).ids).', 1, nelem)], ...
0149             nelem, total_vars);
0150     else
0151         tidx = varsloc.tq.first:varsloc.tq.last;
0152         A = sparse([1:nelem, reshape(ones(nelem,1)*(1:nelem), 1, nelem*nelem) ],...
0153             [tidx, varsloc.Pi.first:varsloc.Pi.last], ...
0154             [-ones(1,nelem), repmat(Qbus.old(vars.(v).ids).', 1, nelem)], ...
0155             nelem, total_vars_loc);
0156     end
0157     lb = -Inf(nelem,1);
0158     ub = Qbus.new(vars.(v).ids);
0159     prob = update_Albub(prob, A, lb, ub);
0160     
0161     if ~loopsolve
0162         A = sparse([1:nelem, reshape(ones(nelem,1)*(1:nelem), 1, nelem*nelem) ],...
0163             [tidx, vars.(v).first:vars.(v).last], ...
0164             [ones(1,nelem), repmat(Qbus.old(vars.(v).ids).', 1, nelem)], ...
0165             nelem, total_vars);
0166     else
0167         A = sparse([1:nelem, reshape(ones(nelem,1)*(1:nelem), 1, nelem*nelem) ],...
0168             [tidx, varsloc.Pi.first:varsloc.Pi.last], ...
0169             [ones(1,nelem), repmat(Qbus.old(vars.(v).ids).', 1, nelem)], ...
0170             nelem, total_vars_loc);
0171     end
0172     lb = Qbus.new(vars.(v).ids);
0173     ub = Inf(nelem,1);
0174     prob = update_Albub(prob, A, lb, ub);
0175     
0176     % ----- stochastic rows of Pi------
0177     % sum(Pi_v(ridx,:)) = 1
0178     if ~loopsolve
0179         A = sparse(reshape(ones(nelem,1)*(1:nelem), 1, nelem*nelem),...
0180                    vars.(v).first:vars.(v).last, 1, nelem, total_vars);
0181     else
0182         A = sparse(reshape(ones(nelem,1)*(1:nelem), 1, nelem*nelem),...
0183                    varsloc.Pi.first:varsloc.Pi.last, 1, nelem, total_vars_loc);
0184     end
0185     lb = ones(nelem,1);
0186     ub = ones(nelem,1);
0187     prob = update_Albub(prob, A, lb, ub);
0188     
0189     % ----- stochastic columns of Pi------
0190     % sum(Pi_v(:,k)) = 1
0191     if ~loopsolve
0192         A  = sparse(reshape((1:nelem).'*ones(1,nelem), 1, nelem*nelem),...
0193                     vars.(v).first:vars.(v).last, 1, nelem, total_vars);
0194     else
0195         A  = sparse(reshape((1:nelem).'*ones(1,nelem), 1, nelem*nelem),...
0196                     varsloc.Pi.first:varsloc.Pi.last, 1, nelem, total_vars_loc);
0197     end
0198     lb = ones(nelem,1);
0199     ub = ones(nelem,1);
0200     prob = update_Albub(prob, A, lb, ub);
0201     
0202     % ------ sum of permutation = sum of initial vector ------
0203     % sum(Pi_v*Pold_v) = sum(Pold_v)
0204     % note all the first column of Pi_v gets Pold_v(1), second column Pold_v(2), etc.
0205     if ~loopsolve
0206         A = sparse(1, vars.(v).first:vars.(v).last, ...
0207             repmat(Pbus.old(vars.(v).ids).', 1, nelem), 1, total_vars);
0208     else
0209         A = sparse(1, varsloc.Pi.first:varsloc.Pi.last, ...
0210             repmat(Pbus.old(vars.(v).ids).', 1, nelem), 1, total_vars_loc);
0211     end
0212     lb = sum(Pbus.old(vars.(v).ids).');
0213     ub = sum(Pbus.old(vars.(v).ids).');
0214     prob = update_Albub(prob, A, lb, ub);
0215     
0216     % sum(Qi_v*Qold_v) = sum(Qold_v)
0217     if ~loopsolve
0218         A = sparse(1, vars.(v).first:vars.(v).last, ...
0219             repmat(Qbus.old(vars.(v).ids).', 1, nelem), 1, total_vars);
0220     else
0221         A = sparse(1, varsloc.Pi.first:varsloc.Pi.last, ...
0222             repmat(Qbus.old(vars.(v).ids).', 1, nelem), 1, total_vars_loc);
0223     end
0224     lb = sum(Qbus.old(vars.(v).ids).');
0225     ub = sum(Qbus.old(vars.(v).ids).');
0226     prob = update_Albub(prob, A, lb, ub);
0227     
0228     clear A;
0229     
0230     if ~loopsolve
0231         prob.x0(vars.(v).first + (0:nelem-1) + nelem*(0:nelem-1)) = 1;
0232     else
0233         %% variable limits
0234         prob.xmin = zeros(total_vars_loc,1);
0235         prob.xmax = Inf(total_vars_loc,1);
0236         prob.xmax(varsloc.Pi.first:varsloc.Pi.last) = 1;
0237         %% linear cost
0238         % sum(tp) + sum(tq)
0239         prob.c = sparse([varsloc.tp.first:varsloc.tp.last, varsloc.tq.first:varsloc.tq.last], ...
0240             1, 1, total_vars_loc, 1);
0241         %% initial point
0242         prob.x0 = zeros(total_vars_loc,1);
0243         prob.x0(varsloc.Pi.first + (0:nelem-1) + nelem*(0:nelem-1)) = 1;
0244         prob.x0(varsloc.tp.first:varsloc.tp.last) = abs(Pbus.new(vars.(v).ids) - Pbus.old(vars.(v).ids));
0245         prob.x0(varsloc.tq.first:varsloc.tq.last) = abs(Qbus.new(vars.(v).ids) - Qbus.old(vars.(v).ids));
0246         %% solve
0247         prob.opt.verbose = opt.vm.nodeperm.verbose;
0248         [x, f, eflag, output, lambda] = qps_matpower(prob);
0249         if eflag
0250             ptmp  = sgvm_extract_perm(reshape(x(varsloc.Pi.first:varsloc.Pi.last), nelem, nelem).', opt.vm.nodeperm);
0251             perm(vars.(v).ids) = vars.(v).ids(ptmp);
0252         else
0253             error('sgvm_deltainjection2perm: optimization failed to converge.')
0254         end
0255     end
0256 end
0257 
0258 if loopsolve
0259     %% verify permutation
0260     if length(unique(perm)) ~= nb
0261         error(['deltainjection2per: binarzation of Pi by selecting maximum entry in each row ',...
0262         'failed to return a valid permutation.'])
0263     end
0264 
0265 else
0266     %% variable limits
0267     prob.xmin = zeros(total_vars,1);
0268     prob.xmax = Inf(total_vars,1);
0269     prob.xmax(vars.Pi1.first:vars.(['Pi', num2str(nblock)]).last) = 1;
0270     %% linear cost
0271     % sum(tp) + sum(tq)
0272     prob.c = sparse([vars.tp.first:vars.tp.last, vars.tq.first:vars.tq.last], ...
0273         1, 1, total_vars, 1);
0274     %% solve
0275     prob.opt.verbose = opt.vm.nodeperm.verbose;
0276     [x, f, eflag, output, lambda] = qps_matpower(prob);
0277 
0278     if eflag
0279         for k = 1:nblock
0280             v = ['Pi', num2str(k)];
0281             nelem = length(vars.(v).ids);
0282             ptmp  = sgvm_extract_perm(reshape(x(vars.(v).first:vars.(v).last), nelem, nelem).', opt.vm.nodeperm);
0283             perm(vars.(v).ids) = vars.(v).ids(ptmp);
0284         end
0285         %% verify permutation
0286         if length(unique(perm)) ~= nb
0287             error(['deltainjection2per: binarzation of Pi by selecting maximum entry in each row ',...
0288             'failed to return a valid permutation.'])
0289         end
0290     else
0291         error('sgvm_deltainjection2perm: optimization failed to converge.')
0292     end
0293 
0294 end
0295 
0296 %% helper functions
0297 function prob = update_Albub(prob, A, lb, ub)
0298 % simple utility function for updating the constraint and bounds arrays.
0299 
0300 prob.A  = [prob.A; A];
0301 prob.l  = [prob.l; lb];
0302 prob.u  = [prob.u; ub];
0303 
0304 function prob = prob_init()
0305 prob = struct();
0306 prob.A  = [];
0307 prob.l  = [];
0308 prob.u  = [];
Generated on Fri 09-Oct-2020 11:21:31 by m2html © 2005