sgvm_data2mpc

SGVM_DATA2MPC create seed MATPOWER case by sampling DATA
   MPC = SGVM_DATA2MPC(DATA, TOPO, OPT)

0001 function mpc = sgvm_data2mpc(data, topo, opt)
0002 %SGVM_DATA2MPC create seed MATPOWER case by sampling DATA
0003 %   MPC = SGVM_DATA2MPC(DATA, TOPO, OPT)
0004 
0005 %   SGVM_DATA2MPC takes parameter data structure DATA and topology, TOPO,
0006 %   and returns a basic matpower case, mpc.
0007 
0008 %   SynGrid
0009 %   Copyright (c) 2018, Power Systems Engineering Research Center (PSERC)
0010 %   by Eran Schweitzer, Arizona State University
0011 %
0012 %   This file is part of SynGrid.
0013 %   Covered by the 3-clause BSD License (see LICENSE file for details).
0014 
0015 %% check inputs and get constants
0016 if nargin < 3
0017     smpl_opt = opt_default(struct());
0018 else
0019     smpl_opt = opt_default(opt);
0020 end
0021 
0022 nb = max(topo(:));
0023 if ~all(sort(unique(topo(:))) == (1:nb).')
0024     error('sgvm_data2mpc: topology must use consecutive buses 1:nb')
0025 end
0026 nl = size(topo,1);
0027 define_constants;
0028 
0029 %% power base
0030 if isfield(data, 'baseMVA')
0031     baseMVA = data.baseMVA;
0032 else
0033     baseMVA = smpl_opt.baseMVA_default;
0034 end
0035 
0036 %% branch matrix
0037 branch = branch_initialize(nl, topo);
0038 if isstruct(data.branch)
0039     tmp = 0;
0040     for f = {'BR_R', 'BR_X', 'BR_B', 'RATE_A', 'TAP', 'SHIFT'}
0041         if ~isfield(data.branch, f{:})
0042             error('sgvm_data2mpc: when provided as a structure, data.branch must have field %s', f{:})
0043         end
0044         if tmp == 0
0045             tmp = length(data.branch.(f{:}));
0046         else
0047             if tmp ~= length(data.branch.(f{:}))
0048                 error('sgvm_data2mpc: when provided as a structure, each field of data.branch must be the same length')
0049             end
0050         end
0051     end
0052     tmpbranch = [data.branch.BR_R, data.branch.BR_X, data.branch.BR_B,...
0053                             data.branch.RATE_A, data.branch.TAP, data.branch.SHIFT];
0054 else
0055     tmpbranch = data.branch;
0056 end
0057 
0058 switch smpl_opt.branch
0059 case {0, 'none'}
0060     % use samples as is
0061     if (size(data.branch,1) == nl) && (size(data.branch,2) == 6) % r, x, b, rate, tap, shift
0062         branch(:,[BR_R, BR_X, BR_B, RATE_A, TAP, SHIFT]) = tmpbranch;
0063     else
0064         error(['sgvm_data2mpc: when using sample method 0:''none'' the provided data matrix must be of length nl.',...
0065         'nl = %d, # of samples = %d'], nl, size(data.branch,1))
0066     end
0067 case {1 , 'direct'}
0068     % sample nl times uniformly at random from data
0069     idx = randi(size(tmpbranch,1), nl, 1);
0070     branch(:,[BR_R, BR_X, BR_B, RATE_A, TAP, SHIFT]) = tmpbranch(idx,:);
0071 case {2 , 'kde'}
0072     % fit a gaussian kde to the data and sample
0073     xfrm_mask = tmpbranch(:,5) ~= 0;
0074     nxfrm = round(nl*sum(xfrm_mask)/size(tmpbranch,1));
0075     nline = nl - nxfrm;
0076     % first sample lines
0077     tmp = tmpbranch(~xfrm_mask,[1:4]);
0078     cnstmask = all(tmp == tmp(1,:), 1);
0079     kde = sgvm_gaussian_kde(tmp(:,~cnstmask));
0080     maxval = max(tmp(:,~cnstmask), [], 1);
0081     minval = min(tmp(:,~cnstmask), [], 1);
0082     rvs = sgvm_sample_dist(nline, kde, minval, maxval);
0083     idx = [BR_R, BR_X, BR_B, RATE_A];
0084     branch(1:nline, idx(~cnstmask)) = rvs;
0085     branch(1:nline, idx(cnstmask))  = repmat(tmp(1,cnstmask), nline, 1);
0086     % next sample xfrms
0087     tmp = tmpbranch(xfrm_mask,:);
0088     cnstmask = all(tmp == tmp(1,:), 1);
0089     kde = sgvm_gaussian_kde(tmp(:,~cnstmask));
0090     maxval = max(tmp(:,~cnstmask), [], 1);
0091     minval = min(tmp(:,~cnstmask), [], 1);
0092     rvs = sgvm_sample_dist(nxfrm, kde, minval, maxval);
0093     idx = [BR_R, BR_X, BR_B, RATE_A, TAP, SHIFT];
0094     branch(nline+1:end, idx(~cnstmask)) = rvs;
0095     branch(nline+1:end, idx(cnstmask))  = repmat(tmp(1,cnstmask), nxfrm, 1);
0096 end
0097 branch(:,RATE_A) = round(branch(:,RATE_A));
0098 
0099 if all(branch(:,RATE_A) == 0) || all(isinf(branch(:,RATE_A)))
0100   warning('sgvm_data2mpc: There are no line ratings in the data. Using default value of %d MVA for all lines', smpl_opt.rate_a_default)
0101   branch(:,RATE_A) = smpl_opt.rate_a_default;
0102 end
0103 %% node properties
0104 bus = bus_initialize(nb);
0105 
0106 %% setup load data in a matrix [PD, QD]
0107 if isstruct(data.load)
0108     tmp = 0;
0109     for f = {'PD','QD'}
0110         if ~isfield(data.load, f{:})
0111             error('sgvm_data2mpc: when provided as a structure, data.load must have field %s', f{:})
0112         end
0113         if tmp == 0
0114             tmp = length(data.load.(f{:}));
0115         else
0116             if tmp ~= length(data.load.(f{:}))
0117                 error('sgvm_data2mpc: when provided as a structure, each field of data.load must be the same length.')
0118             end
0119         end
0120     end
0121     tmpload = [data.load.PD, data.load.QD];
0122 else
0123     tmpload = data.load;
0124 end
0125 
0126 %% setup gen data in a matrix [QMAX, QMIN, PMAX, PMIN]
0127 % if GEN_BUS is provided it is saved in genbus, otherwise, a vector of length 1:length(gen data) is created
0128 % genbusflag indicates whether GEN_BUS was provided or not
0129 if isstruct(data.gen)
0130     tmp = 0;
0131     for f = {'QMAX', 'QMIN', 'PMAX', 'PMIN'}
0132         if ~isfield(data.gen, f{:})
0133             error('sgvm_data2mpc: when provided as a structure, data.gen must have field %s', f{:})
0134         end
0135         if tmp == 0
0136             tmp = length(data.gen.(f{:}));
0137         else
0138             if tmp ~= length(data.gen.(f{:}))
0139                 error('sgvm_data2mpc: when provided as a structure, each field of data.gen must be the same length.')
0140             end
0141         end
0142     end
0143     if isfield(data.gen, 'GEN_BUS')
0144         genbus = data.gen.GEN_BUS;
0145         if length(genbus) ~= tmp
0146                 error('sgvm_data2mpc: when provided as a structure, each field of data.gen must be the same length.')
0147         end
0148         genbusflag = 1;
0149     else
0150         genbusflag = 0;
0151         genbus= (1:tmp).';
0152     end
0153     tmpgen = [data.gen.QMAX, data.gen.QMIN, data.gen.PMAX, data.gen.PMIN];
0154 else
0155     switch size(data.gen,2)
0156     case 4
0157         tmpgen = data.gen;
0158         genbusflag = 0;
0159         genbus = (1:size(tempgen,1)).';
0160     case 5
0161         tmpgen = data.gen(:,2:end);
0162         genbus = data.gen(:,1);
0163         genbusflag = 1;
0164     otherwise
0165         error('sgvm_data2mpc: when providing data.gen as a matrix it must either have 4 or 5 columns.')
0166     end
0167 end
0168 
0169 %% setup gencost data
0170 if isfield(data, 'gencost')
0171     % generator cost provided
0172     if size(data.gencost, 1) ~= size(tmpgen, 1)
0173         error('sgvm_data2mpc: when providing data.gencost it must have the same number of entries as data.gen.')
0174     end
0175     tmpgencost = data.gencost;
0176 else
0177     % no gencost provided create uniform cost
0178     tmp = size(tmpgen, 1);
0179     %             MODEL = 2     STARTUP/SHUTDOWN=0, NCOST = 2,   LINCOST = smpl_opt.lincost,   fixed cost = 0
0180     tmpgencost = [2*ones(tmp,1), zeros(tmp,2),    2*ones(tmp,1), smpl_opt.lincost*ones(tmp,1), zeros(tmp,1)];
0181 end
0182 
0183 %% sample node properties
0184 switch smpl_opt.node
0185 case {0 , 'none'}
0186     % use samples as is
0187     ng = size(tmpgen,1);
0188     gen = gen_initialize(ng);
0189     gencost = tmpgencost;
0190     if size(tmpload,1) ~= nb
0191         error(['sgvm_data2mpc: when using sample method 0: ''none'' the provided data must have size of nb, ',...
0192         'the number of buses. nb = %d, # of load samples = %d'], nb, size(tmpload,1))
0193     end
0194     bus(:, [PD, QD]) = tmpload;
0195     gen(:, [QMAX, QMIN, PMAX, PMIN]) = tmpgen;
0196     if genbusflag
0197         gen(:,GEN_BUS) = genbus;
0198     else
0199         error('sgvm_data2mpc: when using sample method 0: ''none'' generation data must include bus number')
0200     end
0201 case {1 , 'direct'}
0202     idx = randi(size(tmpload,1), nb, 1);
0203     bus(:, [PD, QD]) = tmpload(idx,:);
0204     % directly sampe the load data (gen data is always directly sampled)
0205     if genbusflag && smpl_opt.usegenbus
0206         % if generation bus is specified simply sample buses nb times
0207         %nbus x 1 vector with number of generators at each bus of input data
0208         genperbus = full(sparse(genbus, 1, 1, size(tmpload,1), 1));
0209         ng = sum(genperbus(idx)); % total number of generators
0210         gen = gen_initialize(ng);
0211         gencost = gencost_initialize(ng, size(tmpgencost,2));
0212         ptr = 0;
0213         for bidx = 1:nb
0214             k = idx(bidx); %bus bidx in output is bus k of input data
0215             gidx = find(genbus == k); %rows of gen data connected to bus k
0216             if isempty(gidx)
0217                 continue
0218             end
0219             gen(ptr+1:ptr+genperbus(k), [GEN_BUS, QMAX, QMIN, PMAX, PMIN]) = [ones(genperbus(k),1) * bidx, tmpgen(gidx,:)];
0220             gencost(ptr+1:ptr+genperbus(k),:) = tmpgencost(gidx,:);
0221             ptr = ptr + genperbus(k);
0222         end
0223         if ptr ~= ng
0224             error('sgvm_data2mpc: error while generating the gen matrix %d generators were expected but %d were created', ng, ptr)
0225         end
0226     else
0227         [gen, gencost] = gen_sample(bus, tmpload, tmpgen, tmpgencost, genbus, smpl_opt);
0228     end
0229 case {2 , 'kde'}
0230     % sample load with kde (gen data is always directly sampled)
0231     % first determine number of no-load buses
0232     noloadmask  = (tmpload(:,1) == 0) & (tmpload(:,2) == 0);
0233     noloadbuses = round(nb*sum(noloadmask)/size(tmpload,1));
0234     
0235     % Negative PD
0236     negpdmask   = (tmpload(:,1) < 0);
0237     if sum(negpdmask) < 3
0238         negpdbuses = 0;
0239     else
0240         negpdbuses  = round(nb*sum(negpdmask)/size(tmpload,1));
0241     end
0242     % fit kde to positive load buses
0243     kde = sgvm_gaussian_kde(tmpload(~noloadmask & ~negpdmask,:));
0244     maxval = max(tmpload(~noloadmask & ~negpdmask, :), [], 1);
0245     minval = min(tmpload(~noloadmask & ~negpdmask, :), [], 1);
0246     pospdbuses = nb-noloadbuses-negpdbuses;
0247     rvs = sgvm_sample_dist(pospdbuses, kde, minval, maxval);
0248     bus(1:pospdbuses, [PD, QD]) = rvs;
0249     % fit kde to negative load buses (if any)
0250     if negpdbuses > 0
0251         kde = sgvm_gaussian_kde(tmpload(negpdmask,:));
0252         maxval = max(tmpload(negpdmask, :), [], 1);
0253         minval = min(tmpload(negpdmask, :), [], 1);
0254         rvs = sgvm_sample_dist(negpdbuses, kde, minval, maxval);
0255         bus(pospdbuses+1:pospdbuses + negpdbuses, [PD, QD]) = rvs;
0256     end
0257     % threshold small values to 0 or to minimum value (plus some noise)
0258     thresholdpd = min( abs(tmpload(tmpload(:,1) ~= 0, 1)) );
0259     thresholdqd = min( abs(tmpload(tmpload(:,2) ~= 0, 2)) );
0260     pdbuses  = find((abs(bus(:, PD)) < thresholdpd) & (bus(:,PD) ~= 0));
0261     pd2zero = (thresholdpd - abs(bus(pdbuses, PD)) ) <= 0.5*thresholdpd;
0262     bus(pdbuses(pd2zero), PD)  = 0;
0263     bus(pdbuses(~pd2zero), PD) = sign(bus(pdbuses(~pd2zero), PD)).*thresholdpd.*(1 + rand(sum(~pd2zero),1));
0264     
0265     qdbuses  = find((abs(bus(:, QD)) < thresholdqd) & (bus(:,QD) ~= 0));
0266     qd2zero = (thresholdqd - abs(bus(qdbuses, QD)) ) <= 0.5*thresholdqd;
0267     bus(qdbuses(qd2zero), QD)  = 0;
0268     bus(qdbuses(~qd2zero), QD) = sign(bus(qdbuses(~qd2zero), QD)).*thresholdqd.*(1 + rand(sum(~qd2zero),1));
0269     % sample generators
0270     [gen, gencost] = gen_sample(bus, tmpload, tmpgen, tmpgencost, genbus, smpl_opt);
0271 end
0272 
0273 %% assemble mpc
0274 mpc = struct('version', 2, 'baseMVA', baseMVA, 'bus', bus, 'branch', branch, 'gen', gen, 'gencost', gencost);
0275 
0276 %% adjust gen to load ratio
0277 if ~any(strcmp(smpl_opt.node, {0, 'none'}))
0278     mpc = gen2load_adjust(mpc, tmpload, tmpgen, tmpgencost, genbus, genbusflag, smpl_opt);
0279 end
0280 %% utility functions
0281 
0282 function [gen, gencost] = gen_sample(bus, tmpload, tmpgen, tmpgencost, genbus, opt)
0283 % given a bus matrix alread populated with values this function
0284 % samples from the tmpgen matrix and creates the gen matrix
0285 
0286 define_constants;
0287 nb = size(bus,1);
0288 % buses with no load
0289 noloadbuses  = find( (bus(:,PD) == 0) & (bus(:,QD) == 0));
0290 % buses with load
0291 loadbuses    = find( (bus(:,PD) ~= 0) | (bus(:,QD) ~= 0));
0292 
0293 in_gbusnums = unique(genbus);       %unique *input data* generator bus numbers
0294 multgens = sparse(genbus, 1, 1); % multgens(i) = number of generators at input bus i
0295 
0296 % number buses with generation
0297 if opt.ngbuses > 0
0298     % if specified explicitely by user input
0299     ngbuses = opt.ngbuses;
0300 else
0301     % scale # of input generator buses to total input number of buses to the number of buses nb
0302     ngbuses = round(nb*length(in_gbusnums)/size(tmpload,1));
0303 end
0304 
0305 % sample fraction of no load and no gen buses from [0.6,0.95]
0306 noload_nogen = round(length(noloadbuses)*(0.6 + 0.35*rand()));
0307 noload_gen   = length(noloadbuses) - noload_nogen; % # of noload buses *with* generation
0308 % get # of buses with both load and generation
0309 gen_load     = max(ngbuses - noload_gen, 0); %max(x, 0) is in case noload_gen > ngbuses
0310 
0311 usedbuses = [];
0312 % gen indices to generation to be attached to no load and load buses
0313 noload_gidx = in_gbusnums(randi(length(in_gbusnums),noload_gen,1));
0314 load_gidx   = in_gbusnums(randi(length(in_gbusnums),gen_load, 1));
0315 
0316 ng = full(sum(multgens(noload_gidx)) + sum(multgens(load_gidx))); % number of generators
0317 gen = gen_initialize(ng);
0318 gencost = gencost_initialize(ng, size(tmpgencost,2));
0319 ptr = 0;
0320 % go through the no load buses and create entries in gen matrix
0321 for gidx = sgvm_ensure_col_vect(noload_gidx).'
0322     mask = genbus == gidx; % mask of all generators in input data connected to bus gidx
0323     gbus = noloadbuses(randi(length(noloadbuses))); %bus to attach generators
0324     while ismember(gbus, usedbuses)
0325         gbus = noloadbuses(randi(length(noloadbuses)));
0326     end
0327     gen(ptr+1:ptr+sum(mask), [GEN_BUS, QMAX, QMIN, PMAX, PMIN]) = [ones(sum(mask),1)*gbus, tmpgen(mask,:)];
0328     gencost(ptr+1:ptr+sum(mask),:) = tmpgencost(mask,:);
0329     usedbuses = [usedbuses, gbus];
0330     ptr = ptr + sum(mask);
0331 end
0332 % go through the load buses and create entries in gen matrix
0333 for gidx = sgvm_ensure_col_vect(load_gidx).'
0334     mask = genbus == gidx;
0335     gbus = loadbuses(randi(length(loadbuses))); %bus to attach generators
0336     while ismember(gbus, usedbuses)
0337         gbus = loadbuses(randi(length(loadbuses)));
0338     end
0339     gen(ptr+1:ptr+sum(mask), [GEN_BUS, QMAX, QMIN, PMAX, PMIN]) = [ones(sum(mask),1)*gbus, tmpgen(mask,:)];
0340     gencost(ptr+1:ptr+sum(mask),:) = tmpgencost(mask,:);
0341     usedbuses = [usedbuses, gbus];
0342     ptr = ptr + sum(mask);
0343 end
0344 if ptr ~= ng
0345     error('sgvm_data2mpc: error while generating the gen matrix %d generators were expected but %d were created', ng, ptr)
0346 end
0347 
0348 function mpc = gen2load_adjust(mpc, tmpload, tmpgen, tmpgencost, genbus, genbusflag, smpl_opt)
0349 % resample generators (and possibly load) to move target the desired
0350 % ratio of generation capacity to load.
0351 % If the actual ratio is SMALLER than desired samples are replaced
0352 % until the ratio crosses the threshold. Samples are only kept if the
0353 % move the ratio in the "correct" direction.
0354 
0355 define_constants;
0356 if smpl_opt.usegen2load
0357     % use ratio in data
0358     gen2load = sum(tmpgen(:,3))/sum(tmpload(:,1));
0359 else
0360     % otherwise randomly select a ratio between 1.3 and 1.6 (based on
0361     % Birchfield Metric paper)
0362     gen2load = 1.3 + 0.3*rand();
0363 end
0364 
0365 ratiotest = gen2load - sum(mpc.gen(:,PMAX))/sum(mpc.bus(:,PD));
0366 changedir = sign(ratiotest);
0367 while ratiotest*changedir > 0.05
0368     % pick a random generator bus
0369     gbus0 = mpc.gen(randi(size(mpc.gen,1)), GEN_BUS);
0370     gidx0 = mpc.gen(:,GEN_BUS) == gbus0;
0371     % pick a new entry in the generation data
0372     gbus1 = genbus(randi(size(tmpgen,1)));
0373     gidx1 = genbus == gbus1;
0374     
0375     loadnew = sum(mpc.bus(:,PD));
0376     if any(strcmp(smpl_opt.node,{1, 'direct'})) && genbusflag && smpl_opt.usegenbus
0377         % load and generation are sampled together
0378         loadnew = loadnew - mpc.bus(gbus0,PD) + tmpload(gbus1,1);
0379     end
0380     gennew  = sum(mpc.gen(:,PMAX)) - sum(mpc.gen(gidx0,PMAX)) + sum(tmpgen(gidx1, 3));
0381     if ( (changedir > 0) && ((gennew/loadnew) > sum(mpc.gen(:,PMAX))/sum(mpc.bus(:,PD)))) || ...
0382        ( (changedir < 0) && ((gennew/loadnew) < sum(mpc.gen(:,PMAX))/sum(mpc.bus(:,PD))))
0383         % if more generation needed and new ratio shows MORE generation OR
0384         % if less generation needed and new ratio shows LESS generation
0385         % make change
0386         if any(strcmp(smpl_opt.node,{1, 'direct'})) && genbusflag && smpl_opt.usegenbus
0387             mpc.bus(gbus0,[PD,QD]) = tmpload(gbus1, :);
0388         end
0389         mpc.gen(gidx0, :) = [];
0390         mpc.gencost(gidx0, :) = [];
0391         gen     = gen_initialize(sum(gidx1));
0392         gen(:, [GEN_BUS, QMAX, QMIN, PMAX, PMIN]) = [ones(sum(gidx1),1)*gbus0, tmpgen(gidx1,:)];
0393         gencost = tmpgencost(gidx1,:);
0394         mpc.gen = [mpc.gen; gen];
0395         mpc.gencost = [mpc.gencost; gencost];
0396         
0397         ratiotest = gen2load - sum(mpc.gen(:,PMAX))/sum(mpc.bus(:,PD));
0398     end
0399 end
0400 
0401 
0402 function branch = branch_initialize(nl, topo)
0403 % initialize branch matrix with topology and all branches active
0404 
0405 branch = zeros(nl,13);   %inialize branch matrix
0406 branch(:,[1, 2]) = topo; % 1,2 = F_BUS, T_BUS
0407 branch(:, 11)    = 1;    % 11 = BR_STATUS
0408 
0409 function bus = bus_initialize(nb)
0410 % initialize bus matrix with consecutive bus numbers
0411 % all buses in one area and loss zone and set as PQ buses
0412 
0413 define_constants;
0414 bus = zeros(nb, 13);
0415 bus(:,BUS_I) = (1:nb).';
0416 bus(:,BUS_TYPE) = 1;
0417 bus(:,BUS_AREA) = 1;
0418 bus(:,VM)       = 1;
0419 bus(:,ZONE)     = 1;
0420 bus(:,VMAX)     = 1.05;
0421 bus(:,VMIN)     = 0.95;
0422 
0423 function gen = gen_initialize(ng)
0424 % initial gen matrix
0425 
0426 gen = zeros(ng, 21);
0427 gen(:, 8) = 1; %8 is GEN_STATUS
0428 
0429 function gencost = gencost_initialize(ng, gencostcols)
0430 % initialize the gencost matrix
0431 
0432 gencost = zeros(ng, gencostcols);
0433 
0434 function smpl_opt = opt_default(smpl_opt)
0435 % Create opt data structure with defaults.
0436 % Go through the data structure and update user inputs where necessary
0437 opt = sgvm_smplopts();
0438 smpl_opt = nested_struct_copy(opt, smpl_opt);