sgvm_mpc_perm

SGVM_MPC_PERM main VARIATIONS MODE script
   [R, STATUS] = SGVM_MPC_PERM(MPC, OPT)

   Permutes the branch and node properties in MPC in an
   attempt to generate a better powerflow solution.

   Inputs
       MPC - MATPOWER case to be permuted. Should have been "sanitized"
           by some input process. For example, all buses should be
             consecutive and (at present) no shunt elements
             should exists.
       OPT - SYNGRID options structure

   Outputs
     R      -  array of matpower casese
     STATUS -  0,1,2. 2 means success, no line violations and no voltage
            violations. 1 means success, no line violations but
            voltage violations. 0 either means mpc did not solve OR
            the mpc solves but only with softlims since there are line
            violations.

0001 function [r, status] = sgvm_mpc_perm(mpc, opt)
0002 %SGVM_MPC_PERM main VARIATIONS MODE script
0003 %   [R, STATUS] = SGVM_MPC_PERM(MPC, OPT)
0004 %
0005 %   Permutes the branch and node properties in MPC in an
0006 %   attempt to generate a better powerflow solution.
0007 %
0008 %   Inputs
0009 %       MPC - MATPOWER case to be permuted. Should have been "sanitized"
0010 %           by some input process. For example, all buses should be
0011 %             consecutive and (at present) no shunt elements
0012 %             should exists.
0013 %       OPT - SYNGRID options structure
0014 %
0015 %   Outputs
0016 %     R      -  array of matpower casese
0017 %     STATUS -  0,1,2. 2 means success, no line violations and no voltage
0018 %            violations. 1 means success, no line violations but
0019 %            voltage violations. 0 either means mpc did not solve OR
0020 %            the mpc solves but only with softlims since there are line
0021 %            violations.
0022 
0023 %   SynGrid
0024 %   Copyright (c) 2018, Power Systems Engineering Research Center (PSERC)
0025 %   by Eran Schweitzer, Arizona State University
0026 %
0027 %   This file is part of SynGrid.
0028 %   Covered by the 3-clause BSD License (see LICENSE file for details).
0029 
0030 %% handle inputs
0031 define_constants;
0032 msg = 'sgvm_mpc_perm: error with input mpc:\n';
0033 tests = false(2,1);
0034 tests(1) = ~all(mpc.bus(:,BUS_I) == (1:size(mpc.bus,1)).');
0035 if tests(1)
0036     msg = [msg, '\tbuses need to be consecutive starting at 1.', ...
0037     'Consider using sgvm_data2mpc(sgvm_mpc2data(mpc)).\n'];
0038 end
0039 tests(2) = any(mpc.bus(:, GS) | mpc.bus(:, BS));
0040 if tests(2)
0041     msg = [msg,'\tInput mpc should not have any shunt elements ',...
0042     'GS or BS.\n'];
0043 end
0044 if any(tests)
0045     error(msg)
0046 end
0047 if ~isfield(opt, 'vm')
0048     error('sgvm_mpc_perm: options structure must have a field ''vm'' with the necessary subfields, see sg_options and the SynGrid Manual')
0049 end
0050 
0051 if ~strcmp(opt.mpopt.opf.ac.solver, 'IPOPT')
0052   errorstr = sprintf(['sgvm_mpc_perm: currently sgvm_mpc_perm should be used with IPOPT as the ac-opf solver.\n',...
0053          'Unfortunately results are too unstable with other interior-points solvers.\n',...
0054          'The IPOPT binaries via the PARDISO project for MATPOWER can be found at:\n',...
0055          '\thttps://pardiso-project.org/.']);
0056   warning(errorstr);
0057 end
0058 
0059 %% add softlimits
0060 % the softlimit settings are (at least for now) not changeable by
0061 % the user.
0062 % Softlims are only enabled for line rating and voltage magnitude
0063 max_gen_cost = max(margcost(mpc.gencost, mpc.gen(:, PMAX)));
0064 default_cost = 1000;
0065 
0066 % The following cost calculation effectively equates 1 MVA flow violation
0067 % with 0.1 pu voltage violation.
0068 % This is intentionally fairly generous w.r.t voltage violation since
0069 % We are more ok with voltage violation than rating.
0070 rateacost = max(default_cost, 2*max_gen_cost);
0071 vcost     = max(default_cost*10, 2*10*max_gen_cost);
0072 gencost   = max(default_cost*10, 2*10*max_gen_cost);
0073 opt.vm.softlims = struct('RATE_A', struct('hl_mod', 'remove', 'cost', rateacost),...
0074                       'VMAX', struct('hl_mod', 'replace', 'hl_val', 1.25, 'cost', vcost),...
0075                       'VMIN', struct('hl_mod', 'replace', 'hl_val', 0.8, 'cost', vcost),...
0076                       'QMAX', struct('hl_mod', 'remove', 'cost', gencost), ...
0077                       'QMIN', struct('hl_mod', 'remove', 'cost', gencost),...
0078                       'PMAX', struct('hl_mod', 'remove', 'cost', 2*gencost), ...
0079                       'PMIN', struct('hl_mod', 'remove', 'cost', 2*gencost));
0080 
0081 %% (optional) parallel pool
0082 if opt.vm.parallel.numcores == 0
0083   opt.vm.parallel.use = 0;
0084 elseif opt.vm.parallel.use && ~exist('parfor','builtin')
0085   warning('sgvm_mpc_perm: parallelization requested but it appears the Parallel Computing Toolbox is not installed. Setting parallel.use=0')
0086   opt.vm.parallel.use = 0;
0087 end
0088 if opt.vm.parallel.use
0089   if ~isempty(gcp('nocreate'))
0090     p = gcp;
0091   else
0092     if opt.vm.parallel.numcores > 0
0093       p = parpool(opt.vm.parallel.numcores);
0094     else
0095       p = parpool('local');
0096     end
0097   end
0098 else
0099     p = [];
0100 end
0101 %% start time
0102 tstart = tic;
0103 %% initialize first generation
0104 if opt.verbose > 0
0105   fprintf('Initializing first generation of %d cases...\n', opt.vm.ea.inds)
0106 end
0107 psi = sgvm_GenerationClass(mpc, opt.vm.ea.inds, opt);
0108 if opt.verbose > 0
0109   fprintf('\tinitialization complete (%0.3f sec).\n', toc(tstart))
0110   if opt.verbose > 1
0111     fprintf('Current stats:\n')
0112     psi.stats()
0113   end
0114 end
0115 %% generations loop
0116 for g = 1:opt.vm.ea.generations
0117   tgen = tic;
0118   if opt.verbose > 0
0119     fprintf('-------------------------------------------------------\n')
0120     fprintf('---- Generation Loop %d\n', g)
0121     fprintf('-------------------------------------------------------\n')
0122   end
0123   %% update gen count
0124   psi.update_gen();
0125   %% branch permutation
0126   tbperm = tic;
0127   if opt.verbose > 0
0128     fprintf('Performing branch permutation loop (%d cases)...\n', length(psi.inds))
0129   end
0130     psi.perm_loop('branch', opt);
0131   if opt.verbose > 0
0132     fprintf('\tbranch permutation loop completed (%0.3f sec).\n', toc(tbperm))
0133     if opt.verbose > 1 && ~psi.exitflag
0134       fprintf('Current stats:\n')
0135       psi.stats()
0136       fprintf('Stash stats:\n')
0137       psi.stash.stats()
0138     end
0139   end
0140   %% check solutions
0141   if length(psi.soln) >= opt.vm.ea.select
0142     psi.picksoln(opt.vm.ea.select);
0143     if psi.exitflag
0144       if opt.verbose > 0
0145         fprintf('---- Solution exit flag satisfied (%0.3f sec) \n', g, toc(tgen))
0146         fprintf('\tObjective Range: %0.5g -- %0.5g\n', psi.soln{1}.mpc.f, psi.soln{end}.mpc.f)
0147         fprintf('-------------------------------------------------------\n')
0148       end
0149       psi.inds = {};
0150       break
0151     end
0152   end
0153   %% midway cleanup
0154   % keep initial number of inds or half of current number, whichever is greater
0155   psi.select(max(opt.vm.ea.inds, floor(length(psi.inds)/2)), 0, opt)
0156   %% node permutation
0157   tnperm = tic;
0158   if opt.verbose > 0
0159     fprintf('Performing node permutation loop (%d cases)...\n', length(psi.inds))
0160   end
0161     psi.perm_loop('node', opt);
0162   if opt.verbose > 0
0163     fprintf('\tnode permutation loop completed (%0.3f sec).\n', toc(tnperm))
0164     if opt.verbose > 1
0165       fprintf('Current stats:\n')
0166       psi.stats()
0167       fprintf('Stash stats:\n')
0168       psi.stash.stats()
0169     end
0170   end
0171   %% check solutions
0172   if length(psi.soln) >= opt.vm.ea.select
0173     psi.picksoln(opt.vm.ea.select);
0174     if psi.exitflag
0175       if opt.verbose > 0
0176         fprintf('---- Solution exit flag satisfied (%0.3f sec) \n', g, toc(tgen))
0177         fprintf('\tObjective Range: %0.5g -- %0.5g\n', psi.soln{1}.mpc.f, psi.soln{end}.mpc.f)
0178       end
0179       psi.inds = {};
0180       break
0181     end
0182   end
0183   %% selection
0184   if g < opt.vm.ea.generations
0185     psi.select(opt.vm.ea.inds - opt.vm.ea.randnew, opt.vm.ea.randnew, opt);
0186     if opt.verbose > 0
0187       fprintf('---- Generation Loop %d Complete (%0.3f sec) \n', g, toc(tgen))
0188       fprintf('\tObjective Range: %0.5g -- %0.5g\n', psi.inds{1}.mpc.f, psi.inds{end}.mpc.f)
0189       fprintf('\t Total possible solutions found: %d\n', length(psi.solnlist))
0190     end
0191   else
0192     if length(psi.soln) < opt.vm.ea.select
0193       % select best non-solution individuals
0194       psi.select(opt.vm.ea.select - length(psi.soln), 0);
0195       % merge in solved individuals from stash
0196       psi.merge_stash();
0197       psi.inds_remove_soln();
0198       % combine for solution
0199       psi.soln = horzcat(psi.soln, psi.inds);
0200     end
0201       psi.picksoln(opt.vm.ea.select);
0202       psi.inds = {};
0203       if opt.verbose > 0
0204         fprintf('---- Final Generation Loop %d Complete (%0.3f sec) \n', g, toc(tgen))
0205         fprintf('\tObjective Range: %0.5g -- %0.5g\n', psi.soln{1}.mpc.f, psi.soln{end}.mpc.f)
0206         fprintf('\t Total possible solutions found: %d\n', length(psi.solnlist))
0207       end
0208   end
0209 end
0210 %% Reactive planning
0211 if opt.verbose > 0
0212   tshunts = tic;
0213   fprintf('-------------------------------------------------------\n')
0214   fprintf('---- Reactive Planning Stage \n')
0215   fprintf('\tAdding shunts to %d cases...\n', length(psi.soln))
0216 end
0217 psi.reactive_planning(opt);
0218 if opt.verbose > 0
0219   fprintf('\tReactive planning completed (%0.3f sec).\n', toc(tshunts))
0220 end
0221 %%
0222 if opt.verbose > 0
0223   fprintf('-------------------------------------------------------\n')
0224   fprintf('---- Completed (%0.3f sec).\n', toc(tstart))
0225   fprintf('\tObjective Range: %0.5g -- %0.5g\n', psi.soln{1}.mpc.f, psi.soln{end}.mpc.f)
0226   psi.stats('soln')
0227 end
0228 
0229 %% delete parallel pool
0230 if ~isempty(p)
0231   delete(p)
0232 end
0233 
0234 %% collect result
0235 [r, status] = psi.mpc_export() ;