0001 function perm = sgvm_deltainjection2perm(Pbus, Qbus, opt)
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017 nb = length(Pbus.old);
0018
0019
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);
0025
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
0035 if ~all(sort(perm) == (1:nb).')
0036 error('deltainjection2per: failed to return a valid permutation.')
0037 end
0038
0039 return
0040
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);
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));
0055 ids = sort(ids);
0056 available_ids(ids) = false;
0057 else
0058 nelem = ceil(avgnum);
0059 ids = find(available_ids);
0060 ids = ids(randperm(length(ids), nelem));
0061 ids = sort(ids);
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
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
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
0082
0083
0084
0085
0086
0087
0088
0089
0090
0091
0092
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
0105
0106
0107
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
0141
0142
0143
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
0177
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
0190
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
0203
0204
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
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
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
0238
0239 prob.c = sparse([varsloc.tp.first:varsloc.tp.last, varsloc.tq.first:varsloc.tq.last], ...
0240 1, 1, total_vars_loc, 1);
0241
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
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
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
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
0271
0272 prob.c = sparse([vars.tp.first:vars.tp.last, vars.tq.first:vars.tq.last], ...
0273 1, 1, total_vars, 1);
0274
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
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
0297 function prob = update_Albub(prob, A, lb, ub)
0298
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 = [];