sgvm_IndClass

0001 classdef sgvm_IndClass < handle
0002 %SGVM_INDCLASS individual placement in the variations method
0003 %   IND = SGVM_INDCLASS(MPC, TYPE, GEN, OPT) evolves individuals based on
0004 %   MPC and the TYPE of permuation (NODE, BRANCH, INIT) that is specified.
0005 %
0006 %   The SGVM_INDCLASS contains a MATPOWER case and methods to modify it
0007 %   towards a more desireable solution by permuting properties as well as
0008 %   adding shunt elements.
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     properties
0018         mpc
0019         id
0020         pid
0021         ng
0022         nb
0023         nl
0024         gen
0025         scale_s
0026         overload_frac
0027         call
0028         deltaf
0029     end
0030     methods
0031         function obj = sgvm_IndClass(mpc, type, gen, opt)
0032             if isa(mpc, 'sgvm_IndClass')
0033               parent = mpc;
0034               mpc    = parent.mpc;
0035               obj.scale_s = parent.scale_s;
0036               obj.overload_frac = parent.overload_frac;
0037               obj.pid = parent.id;
0038             else
0039               obj.scale_s = opt.vm.nodeperm.scale_s;
0040               obj.overload_frac = 1;
0041               obj.pid = '--N/A--';
0042             end
0043             obj.ng = size(mpc.gen, 1);
0044             obj.nb = size(mpc.bus, 1);
0045             obj.nl = size(mpc.branch, 1);
0046             obj.gen = gen;
0047 
0048             switch type
0049               case 'init'
0050                 obj.call = 'init';
0051                 obj.mpc = mpc;
0052                 obj.strip_result();
0053                 obj.mpc = sgvm_branch_perm(obj.mpc);
0054                 obj.mpc = sgvm_nodeperm_init(obj.mpc);
0055                 obj.assign_bustypes();
0056                 obj.setid();
0057                 obj.solve(opt);
0058             %         obj = sgvm_IndClass(obj.mpc, 'branch', obj.gen, opt);
0059               case 'branch'
0060                 obj.call = 'branch';
0061                 [~, nviolations] = obj.checkoverloads(mpc);
0062                 if nviolations < 15
0063                   obj.overload_frac = max(0.80, opt.vm.branchperm.overload_frac_factor*obj.overload_frac);
0064                 else
0065                   obj.overload_frac = 1;
0066                 end
0067                 opt.vm.branchperm.overload_frac = obj.overload_frac;
0068                 if opt.vm.branchperm.verbose > 0
0069                   fprintf('\t(child of %s) overload_frac = %0.3f\n', obj.pid(1:6), obj.overload_frac);
0070                 end
0071                 [obj.mpc, flag] = sgvm_branch_perm(mpc, opt);
0072                 obj.setid();
0073                 if flag
0074                   obj.clear_shunts();
0075                   obj.solve(opt);
0076                 elseif opt.vm.branchperm.verbose > 0
0077                   fprintf('\t(child of %s) Identity permutation returned (i.e no improvement).\n', obj.pid(1:6));
0078                 end
0079               case 'node'
0080                 obj.call = 'node';
0081                 if opt.vm.nodeperm.verbose > 0 && ~opt.vm.parallel.use
0082                   fprintf('\t(child of %s)---Time stats:\n', obj.pid(1:6));
0083                 end
0084                 [~, nviolations] = obj.checkoverloads(mpc);
0085                 if nviolations < 10
0086                   opt.vm.nodeperm.usedv = true;
0087                   obj.scale_s = max(0.75, opt.vm.nodeperm.scale_s_factor*obj.scale_s);
0088                 else
0089                   obj.scale_s = min(1, obj.scale_s/opt.vm.nodeperm.scale_s_factor);
0090                 end
0091                 opt.vm.nodeperm.scale_s = obj.scale_s;
0092                 tcalc = tic;
0093                 [Pbus,Qbus] = sgvm_calc_injection_delta(mpc, opt);
0094                 if opt.vm.nodeperm.verbose > 0
0095                   fprintf('\t(child of %s) injection calc: scale_s = %0.3f, time %0.3f\n', obj.pid(1:6), obj.scale_s, toc(tcalc));
0096                 end
0097                 tperm = tic;
0098                 perm = sgvm_deltainjection2perm(Pbus, Qbus, opt);
0099                 if opt.vm.nodeperm.verbose > 0
0100                   fprintf('\t(child of %s) sgvm_deltainjection2perm: time %0.3f\n', obj.pid(1:6), toc(tperm));
0101                 end
0102                 if all(perm == (1:size(mpc.bus,1)).')
0103                   if opt.vm.nodeperm.verbose > 0
0104                      fprintf('\t(child of %s) Identity permutation returned (i.e no improvement).\n', obj.pid(1:6));
0105                   end
0106                   obj.mpc = mpc;
0107                 else
0108                   obj.mpc = sgvm_perform_permute(mpc, perm, 'perm');
0109                   obj.assign_bustypes();
0110                   obj.setid();
0111                   tsolve = tic;
0112                   obj.clear_shunts();
0113                   obj.solve(opt);
0114                   if opt.vm.nodeperm.verbose > 0
0115                     fprintf('\tInd %s (child of %s) solve time: %0.3f\n', obj.id(1:6), obj.pid(1:6), toc(tsolve));
0116                   end
0117                 end
0118             end
0119             obj.setid();
0120             if isfield(mpc,'f') && obj.solcheck()
0121                 obj.deltaf = obj.mpc.f - mpc.f;
0122             else
0123                 obj.deltaf = NaN;
0124             end
0125         end
0126 
0127         function setid(obj)
0128             % ID is a hash of some of the key data
0129             % hashing is largely copied from:
0130             % https://www.mathworks.com/matlabcentral/answers/45323-how-to-calculate-hash-sum-of-a-string-using-java
0131 
0132             define_constants;
0133             data = [obj.mpc.bus(:,PD); obj.mpc.bus(:,QD);...
0134                 obj.mpc.branch(:,BR_R); obj.mpc.branch(:,BR_X); obj.mpc.branch(:,BR_B);...
0135                 obj.mpc.branch(:,TAP); obj.mpc.branch(:,SHIFT); obj.mpc.gen(GEN_BUS)];
0136             str = mat2str(data,8);
0137             if have_fcn('octave')
0138                 obj.id = hash('SHA256', str);
0139             else
0140                 md = java.security.MessageDigest.getInstance('SHA-256');
0141                 obj.id = sprintf('%2.2x', typecast(md.digest(uint8(str)), 'uint8'));
0142             end
0143         end
0144 
0145         function obj = solve(obj, opt)
0146             % ensure softlimits are turned on and solve opf
0147 
0148             define_constants;
0149             obj.mpc.softlims = opt.vm.softlims;
0150             if ~toggle_softlims(obj.mpc, 'status')
0151               obj.mpc = toggle_softlims(obj.mpc, 'on');
0152             end
0153             obj.toggle_var_support('on'); % add fictitious generators for var support
0154             mpopt = opt.mpopt;
0155             if ~isempty(opt.vm.opflogpath) && strcmp(mpopt.opf.ac.solver, 'IPOPT')
0156               fname = fullfile(opt.vm.opflogpath,sprintf('Ind%s_%s.out', obj.id(1:6), obj.call));
0157               mpopt.ipopt.opts.output_file =fname;
0158             end
0159             if strcmp(obj.call, 'init')
0160               maxit = sgvm_get_max_iter(mpopt);
0161               if maxit < 1000
0162                   mpopt = sgvm_set_max_iter(mpopt, 1000);
0163               end
0164             elseif obj.gen >= opt.vm.ea.generations
0165               obj.mpc.softlims.PMAX.hl_mod = 'none';
0166               obj.mpc.softlims.PMIN.hl_mod = 'none';
0167             end
0168             if (obj.gen > 0) && obj.gen <= opt.vm.ea.generations
0169               % progressively make reactive support less expensive and line violations more expensive
0170               obj.mpc.softlims.QMAX.cost = opt.vm.softlims.RATE_A.cost + (opt.vm.ea.generations - obj.gen)/opt.vm.ea.generations * (opt.vm.softlims.QMAX.cost - opt.vm.softlims.RATE_A.cost);
0171               obj.mpc.softlims.QMIN.cost = obj.mpc.softlims.QMAX.cost;
0172 
0173               obj.mpc.softlims.RATE_A.cost = opt.vm.softlims.QMAX.cost - (opt.vm.ea.generations - obj.gen)/opt.vm.ea.generations * (opt.vm.softlims.QMAX.cost - opt.vm.softlims.RATE_A.cost);
0174 
0175               % progressively shrink tolerable voltage range
0176               obj.mpc.softlims.VMAX.hl_val = 1.1 + (opt.vm.ea.generations - obj.gen)/opt.vm.ea.generations * (opt.vm.softlims.VMAX.hl_val - 1.1);
0177               obj.mpc.softlims.VMIN.hl_val = 0.9 - (opt.vm.ea.generations - obj.gen)/opt.vm.ea.generations * (0.9 - opt.vm.softlims.VMIN.hl_val);
0178             end
0179             while true
0180               if opt.verbose > 1
0181                 maxit = sgvm_get_max_iter(mpopt);
0182                 fprintf('  Ind %s (child of %s): solving with %d max iterations\n', obj.id(1:6), obj.pid(1:6), maxit);
0183               end
0184               obj.mpc = rundcopf(obj.mpc, mpopt);
0185               mpopt.opf.start = 2;
0186               r = runopf(obj.mpc, mpopt);
0187               if opt.verbose > 1
0188                 fprintf('  Ind %s (child of %s): runopf complete with status %d\n', obj.id(1:6), obj.pid(1:6), r.success);
0189               end
0190               if ~r.success
0191                 % if initial solution fails try using a powerflow initialization
0192                 tmp = struct('baseMVA', r.baseMVA, 'bus', r.bus,...
0193                     'branch', r.branch, 'gen', r.gen, ...
0194                     'gencost', r.gencost);
0195                 rpf = runpf(tmp, opt.mpopt);
0196                 if opt.verbose > 1
0197                   fprintf('  Ind %s (child of %s): runpf complete with status %d\n', obj.id(1:6), obj.pid(1:6), rpf.success);
0198                 end
0199                 if rpf.success
0200                   obj.mpc_from_pf(rpf, r.softlims);
0201                   break
0202                 else
0203                   if strcmp(obj.mpc.softlims.PMAX.hl_mod, 'none')
0204                       obj.mpc.softlims.PMAX.hl_mod = 'remove';
0205                       obj.mpc.softlims.PMIN.hl_mod = 'remove';
0206                   %end
0207                   %switch mpopt.opf.ac.solver
0208                   %    case {'MIPS', 'IPOPT'}
0209                   %        maxit = sgvm_get_max_iter(mpopt, 1000);
0210                   %        if maxit < 1000
0211                   %            warning('sgvm_IndClass/solve: increasing Ind %s (child of %s) maximum iteration to %d.', obj.id(1:6), obj.pid(1:6), 2*maxit)
0212                   %            mpopt = sgvm_set_max_iter(mpopt, 2*maxit);
0213                   %        else
0214                   %            obj.mpc = r;
0215                   %            warning('sgvm_IndClass/solve: Unable to solve opf (Ind %s, child of %s).', obj.id(1:6), obj.pid(1:6))
0216                   %            break
0217                   %        end
0218                   %    otherwise
0219                   else
0220                     obj.mpc = r;
0221                     warning('sgvm_IndClass/solve: Unable to solve opf (Ind %s, child of %s).', obj.id(1:6), obj.pid(1:6))
0222                     break
0223                   end
0224                 end
0225               else
0226                 obj.mpc = r;
0227                 break
0228               end
0229             end
0230             obj.toggle_var_support('off'); % remove fictitious generators for var support
0231             obj.clean_opf_fields();
0232         end
0233 
0234         function toggle_var_support(obj,on_off)
0235             % toggle variable shunts at every node on or off
0236 
0237             define_constants;
0238             if strcmp(on_off, 'on')
0239               new_gen = zeros(obj.nb, size(obj.mpc.gen,2));
0240               new_gen(:,[GEN_BUS, VG, MBASE, GEN_STATUS]) = ...
0241                   [(1:obj.nb)', ones(obj.nb,1), obj.mpc.baseMVA*ones(obj.nb,1), ones(obj.nb,1)];
0242               obj.mpc.gen = [obj.mpc.gen; new_gen];
0243 
0244               cst = 0;%obj.mpc.softlims.QMAX.cost(1);
0245               new_gencost = zeros(obj.nb, size(obj.mpc.gencost,2));
0246               new_gencost(:,[MODEL, NCOST, COST]) = [2*ones(obj.nb,1), 2*ones(obj.nb,1), cst*ones(obj.nb,1)];
0247 
0248               obj.mpc.gencost = [obj.mpc.gencost; new_gencost];
0249 
0250               obj.mpc.softlims.QMAX.idx   = (obj.ng+1:obj.ng+obj.nb)';
0251               obj.mpc.softlims.QMIN.idx   = (obj.ng+1:obj.ng+obj.nb)';
0252               obj.mpc.softlims.PMAX.idx   = (obj.ng+1:obj.ng+obj.nb)';
0253               obj.mpc.softlims.PMIN.idx   = (obj.ng+1:obj.ng+obj.nb)';
0254             elseif strcmp(on_off, 'off')
0255               obj.mpc.gen     = obj.mpc.gen(1:obj.ng, :);
0256               obj.mpc.gencost = obj.mpc.gencost(1:obj.ng, :);
0257 
0258               %% Convert to generators to shunt elements
0259               if isfield(obj.mpc.softlims.QMAX, 'overload')
0260                 bsh = (obj.mpc.softlims.QMAX.overload(obj.ng+1:end) - obj.mpc.softlims.QMIN.overload(obj.ng+1:end))./obj.mpc.bus(:,VM).^2;
0261               else
0262                 bsh = zeros(obj.nb,1);
0263               end
0264               if isfield(obj.mpc.softlims.PMIN, 'overload')
0265                 % IMPORTANT: this could produce NEGATIVE real shunt elements. That is CLEARLY not realistic.
0266                 % These elements will be removed at a later point, they are simply a way of indicating the added injections for a future
0267                 % child object.
0268                 gsh = (obj.mpc.softlims.PMIN.overload(obj.ng+1:end) - obj.mpc.softlims.PMAX.overload(obj.ng+1:end))./obj.mpc.bus(:,VM).^2;
0269               else
0270                 gsh = zeros(obj.nb,1);
0271               end
0272               obj.mpc.bus(:, [GS, BS]) = [gsh, bsh];
0273             else
0274               error('sgvm_IndClass/toggle_var_support: on_off must be either ''on'' or ''off''.')
0275             end
0276         end
0277 
0278         function clear_shunts(obj)
0279             % remove shunt elements from mpc object
0280             define_constants;
0281             obj.mpc.bus(:, [GS, BS]) = 0;
0282         end
0283 
0284         function clean_opf_fields(obj)
0285             % remove opf fields that will not be accessed and are therefore taking up a lot of space
0286             fields = {'om', 'x', 'mu', 'var', 'lin', 'qdc', 'raw', 'nle', 'nli'};
0287             for f = fields
0288               if isfield(obj.mpc, f{1})
0289                 obj.mpc = rmfield(obj.mpc, f{1});
0290               end
0291             end
0292         end
0293 
0294         function obj = mpc_from_pf(obj, r, softlims)
0295             % take a powerflow solved case, r, and convert it to what appears like
0296             % a solved opf case
0297 
0298             define_constants;
0299             % indicate that r solved but not exactly in the expected way.
0300             % by using 2 a test like (if r.success) still passes since r.success is
0301             % not 0.
0302             r.success = 2;
0303             cost = sum(totcost(r.gencost, r.gen(:, PG)));
0304             r.softlims = softlims;
0305             for prop = fieldnames(r.softlims).'
0306               if strcmp(r.softlims.(prop{:}).hl_mod, 'none')
0307                 continue
0308               end
0309               switch prop{:}
0310                 case 'RATE_A'
0311                   Sf = sqrt(r.branch(:,PF).^2 + r.branch(:,QF).^2);
0312                   St = sqrt(r.branch(:,PT).^2 + r.branch(:,QT).^2);
0313                   %maximum apparent flow on branch
0314                   S  = max(Sf,St);
0315                   idx = r.softlims.RATE_A.idx;
0316                   r.softlims.RATE_A.overload = zeros(size(r.branch, 1), 1);
0317                   r.softlims.RATE_A.overload(idx) = max(0, S(idx) - r.branch(idx, RATE_A));
0318                 case {'PMAX', 'QMAX'}
0319                   idx = (1:size(r.gen,1)).'; %r.softlims.(prop{:}).idx;
0320                   r.softlims.(prop{:}).overload = zeros(size(r.gen,1),1);
0321                   r.softlims.(prop{:}).overload(idx) = max(0, ...
0322                       r.gen(idx, eval([prop{:}(1), 'G'])) - r.gen(idx, eval(prop{:})));
0323                 case 'VMAX'
0324                   idx = r.softlims.VMAX.idx;
0325                   r.softlims.VMAX.overload = zeros(size(r.bus, 1), 1);
0326                   r.softlims.VMAX.overload(idx) = max(0, r.bus(idx, VM) - r.bus(idx, VMAX));
0327                 case 'VMIN'
0328                   idx = r.softlims.VMIN.idx;
0329                   r.softlims.VMIN.overload = zeros(size(r.bus, 1), 1);
0330                   r.softlims.VMIN.overload(idx) = max(0, r.bus(idx, VMIN) - r.bus(idx, VM));
0331                 case 'ANGMAX'
0332                   delta = calc_branch_angle(r);
0333                   idx = r.softlims.ANGMAX.idx;
0334                   r.softlims.ANGMAX.overload = zeros(size(r.branch, 1), 1);
0335                   r.softlims.ANGMAX.overload(idx) = max(0, delta(idx) - r.branch(idx,ANGMAX));
0336                 case 'ANGMIN'
0337                   delta = calc_branch_angle(r);
0338                   idx = r.softlims.ANGMIN.idx;
0339                   r.softlims.ANGMIN.overload = zeros(size(r.branch, 1), 1);
0340                   r.softlims.ANGMIN.overload(idx) = max(0, delta(idx) - r.branch(idx,ANGMIN));
0341                 case {'PMIN', 'QMIN'}
0342                   idx = (1:size(r.gen,1)).'; %r.softlims.(prop{:}).idx;
0343                   r.softlims.(prop{:}).overload = zeros(size(r.gen,1),1);
0344                   r.softlims.(prop{:}).overload(idx) = max(0, ...
0345                      r.gen(idx, eval(prop{:})) - r.gen(idx, eval([prop{:}(1), 'G'])));
0346               end
0347               cost = cost + r.softlims.(prop{:}).cost(1) * sum(r.softlims.(prop{:}).overload);
0348             end
0349             r.f = cost; % total cost at current solution
0350             obj.mpc = r;
0351         end
0352 
0353         function bool = solcheck(obj)
0354             % checks whether mpc has been solved and if so, whether it converged
0355             bool = 0;
0356             if isfield(obj.mpc, 'success')
0357               if obj.mpc.success
0358                 bool = 1;
0359               end
0360             end
0361         end
0362 
0363         function ind = iscomplete(obj, level)
0364             % checks if mpc satisfies all the necessary constraint.
0365             if nargin < 2
0366               level = 1;
0367             end
0368             ind = obj.checkoverloads(obj.mpc);
0369             if level > 1
0370               define_constants;
0371               % check voltage constraints
0372               vtest = obj.mpc.bus(:,VM) <= (obj.mpc.bus(:,VMAX) + 1e-4);
0373               vtest = vtest & (obj.mpc.bus(:,VM) >= (obj.mpc.bus(:,VMIN) - 1e-4));
0374               ind = ind + all(vtest);
0375             end
0376         end
0377 
0378         function bool = comp_perm(obj, x)
0379             % check whether obj.mpc is the sampe as x.mpc by comparing their hash
0380             bool = isequal(obj.id, x.id);
0381         end
0382 
0383         function obj = strip_result(obj)
0384             % keep only bus, branch, gen and gencost, and baseMVA fields of mpc
0385             % and remove powerflow results
0386             % change all buses to PQ buses
0387             tmp.bus     = obj.mpc.bus(:,1:13);
0388             tmp.branch  = obj.mpc.branch(:,1:13);
0389             tmp.gen     = obj.mpc.gen(:, 1:21);
0390             tmp.gencost = obj.mpc.gencost;
0391             tmp.baseMVA = obj.mpc.baseMVA;
0392             obj.mpc = tmp;
0393         end
0394 
0395         function obj = assign_bustypes(obj)
0396             % all generator buses are PV, largest generator bus is REF,
0397             % and all others are PQ
0398             define_constants;
0399             obj.mpc.bus(:,BUS_TYPE) = PQ;
0400             % all gen buses are pv buses
0401             obj.mpc.bus(obj.mpc.gen(:, GEN_BUS), BUS_TYPE) = 2;
0402 
0403             % reference bus is bus with largest gen
0404             [~,tmpidx] = max(obj.mpc.gen(:, PMAX));
0405             refbus = obj.mpc.gen(tmpidx, GEN_BUS);
0406             obj.mpc.bus(refbus, BUS_TYPE) = 3;
0407         end
0408 
0409         function obj = add_shunts(obj, opt)
0410             % add shunts to the case to meet the voltage requirements and try to
0411             % avoid negative LMPS. The algorithm can fail in 2 ways:
0412             % 1) the add_shunts routine fails to solve the opf
0413             % 2) the could be branch flow violations at the end of the procedure
0414             % for case 1, the routine restores a backup of the original object
0415             % and returns.
0416             % for case 2 the overloads should be indicated by the softlims field
0417             % and handled elswhere.
0418 
0419             define_constants;
0420             bkup = obj;
0421             bkupflag = 1;
0422             cnt = 0;
0423             while true
0424                 cnt = cnt + 1;
0425                 obj.clear_shunts();
0426                 if ~isempty(opt.vm.opflogpath) && strcmp(mpopt.opf.ac.solver, 'IPOPT')
0427                     fname = fullfile(opt.vm.opflogpath,sprintf('Ind%s_shunts.out', obj.id(1:6)));
0428                     mpopt.ipopt.opts.output_file =fname;
0429                 end
0430                 opt.mpopt.opf.start = 2;
0431                 if opt.verbose > 1
0432                     fprintf('  Ind %s (child of %s): solving.\n', obj.id(1:6), obj.pid(1:6) );
0433                 end
0434                 [r, ~] = sgvm_add_shunts(obj.mpc, opt.mpopt, opt.vm.shunts);
0435                 if ~r.success
0436                     obj = bkup;
0437                     obj.clean_opf_fields();
0438                     return
0439                 end
0440                 if opt.vm.shunts.verbose > 0
0441                     fprintf('  Ind %s (child of %s): sgvm_add_shunts complete.\n', obj.id(1:6), obj.pid(1:6));
0442                 end
0443                 if (obj.branch_violations(r) == 0) || (cnt > 5)
0444                     if isfield(r, 'softlims')
0445                         r = rmfield(r, 'softlims');
0446                     end
0447                     if toggle_softlims(r, 'status')
0448                         r = toggle_softlims(r, 'off');
0449                     end
0450 
0451                     % solve WITHOUT softlims
0452                     tmp = runopf(r, opt.mpopt);
0453                     if opt.vm.shunts.verbose > 1
0454                         fprintf('result with No softlimits:\n')
0455                         printpf(tmp)
0456                     end
0457                     % update the saved backup
0458                     if bkupflag
0459                       bkup.mpc = tmp;
0460                       bkupflag = 0;
0461                     elseif  tmp.success
0462                       bkup.mpc = tmp;
0463                     end
0464                     if (min(tmp.bus(:, LAM_P)) > 0 )
0465                         if opt.vm.shunts.verbose > 0
0466                             fprintf('  Ind %s (child of %s): successful completion.\n', obj.id(1:6), obj.pid(1:6));
0467                         end
0468                         break
0469                     elseif cnt > 5
0470                         if opt.vm.shunts.verbose > 0
0471                             fprintf('  Ind %s (child of %s): shunt iteration threshold exceeded.\n', obj.id(1:6), obj.pid(1:6));
0472                         end
0473                         break
0474                     elseif opt.vm.shunts.shift_in > 0.04
0475                         if opt.vm.shunts.verbose > 0
0476                             fprintf('  Ind %s (child of %s): Maximum shift_in of 0.5 reached.\n', obj.id(1:6), obj.pid(1:6));
0477                         end
0478                         break
0479                     else
0480                         opt.vm.shunts.shift_in = opt.vm.shunts.shift_in + 0.005;
0481                         if opt.vm.shunts.verbose > 0
0482                             fprintf('  Ind %s (child of %s): negative LMP detected, tightning opt.vm.shunts.shift_in to %0.3f.\n', obj.id(1:6), obj.pid(1:6), opt.vm.shunts.shift_in);
0483                         end
0484                         continue
0485                     end
0486                 end
0487                 obj.mpc = r;
0488                 permcnt = 0;
0489                 while ~obj.iscomplete() && (permcnt < 5)
0490                     obj.scale_s = 1;
0491                     obj.overload_frac = 1;
0492                     permcnt = permcnt + 1;
0493                     [nviolations, idx] = obj.branch_violations(obj.mpc);
0494                     if strcmp(obj.call, 'branch')
0495                         type = 'node';
0496                     elseif strcmp(obj.call, 'node')
0497                         type = 'branch';
0498                     elseif strcmp(obj.call, 'init')
0499                         type = 'branch';
0500                     end
0501                     if opt.vm.shunts.verbose > 0
0502                         fprintf('  Ind %s (child of %s): %d overloaded branches detected. Performing %s permutation.\n', obj.id(1:6), obj.pid(1:6), nviolations, type)
0503                         if nviolations < 10
0504                             obj.print_flow_stats(idx)
0505                         end
0506                     end
0507                     opttmp = opt;
0508                     opttmp.vm.softlims.RATE_A.cost = 2*opt.vm.softlims.QMAX.cost;
0509                     opttmp.vm.softlims.QMAX.cost   = 0.5*opt.vm.softlims.RATE_A.cost;
0510                     opttmp.vm.softlims.QMIN.cost   = opttmp.vm.softlims.QMAX.cost;
0511                     opttmp.vm.softlims.VMAX.cost   = opttmp.vm.softlims.QMAX.cost;
0512                     opttmp.vm.softlims.VMIN.cost   = opttmp.vm.softlims.QMAX.cost;
0513                     opttmp.vm.softlims.VMAX.hl_val = 1.1;
0514                     opttmp.vm.softlims.VMIN.hl_val = 0.85;
0515                     opttmp.vm.nodeperm.scale_s_factor = 1;
0516                     opttmp.vm.branchperm.overload_frac_factor = 1;
0517                     obj = sgvm_IndClass(obj, type, opt.vm.ea.generations + 1, opttmp);
0518                 end
0519             end
0520 
0521             %if isfield(r, 'softlims')
0522             %   r = rmfield(r, 'softlims');
0523             %end
0524             %if toggle_softlims(r, 'status')
0525             %   r = toggle_softlims(r, 'off');
0526             %end
0527             %
0528             %% solve WITHOUT softlims
0529             %tmp = runopf(r, opt.mpopt);
0530                 %fprintf('result with No softlimits:\n')
0531                 %printpf(tmp)
0532 
0533             obj.mpc = tmp;
0534             obj.clean_opf_fields();
0535         end
0536 
0537         function print_flow_stats(obj, idx)
0538             if nargin == 1
0539                 idx = (1:size(obj.mpc.branch, 1)).';
0540             end
0541             define_constants;
0542             sf = sqrt(obj.mpc.branch(idx, PF).^2 + obj.mpc.branch(idx, QF).^2);
0543             st = sqrt(obj.mpc.branch(idx, PT).^2 + obj.mpc.branch(idx, QT).^2);
0544             ra = obj.mpc.branch(idx, RATE_A);
0545             fprintf('----------------------------------------------------\n')
0546             fprintf('   flow statistic\n')
0547             fprintf('----------------------------------------------------\n')
0548             fprintf(' branch id  |  SF (MVA)  |  SF (MVA)  | RATE (MVA) |\n')
0549             fprintf('------------|------------|------------|------------|\n')
0550             fprintf('%11d |%11.4f |%11.4f |%11.4f |\n', [idx, sf, st, ra].')
0551         end
0552 
0553         %% ------  private methods  -----------
0554         function [bool, nviolations] = checkoverloads(obj, mpc)
0555             % checks if mpc satisfies all the necessary constraint.
0556 
0557             if isfield(mpc, 'softlims')
0558                 if isfield(mpc.softlims, 'RATE_A')
0559                   nviolations = sum(mpc.softlims.RATE_A.overload > 1e-4);
0560                 else
0561                   nviolations = obj.branch_violations(mpc);
0562                 end
0563                 if isfield(mpc.softlims, 'PMAX') && ~strcmp(mpc.softlims.PMAX.hl_mod, 'none')
0564                   genoverload = sum(mpc.softlims.PMAX.overload > 1e-4);
0565                 else
0566                   genoverload = obj.gen_violations(mpc);
0567                 end
0568             else
0569                 nviolations = obj.branch_violations(mpc);
0570                 genoverload = obj.gen_violations(mpc);
0571             end
0572             bool =  (nviolations == 0) && (genoverload == 0);
0573             bool = bool && mpc.success;
0574         end
0575 
0576         function [nviolations, idx] = branch_violations(obj, mpc)
0577 
0578             define_constants;
0579             Sf = sqrt(mpc.branch(:,PF).^2 + mpc.branch(:,QF).^2);
0580             St = sqrt(mpc.branch(:,PT).^2 + mpc.branch(:,QT).^2);
0581             %maximum apparent flow on branch in per-unit
0582             S  = max(Sf,St) / mpc.baseMVA;
0583             % line ratings in per-unit
0584             ra  = mpc.branch(:,RATE_A) / mpc.baseMVA;
0585             ra(ra == 0) = Inf;
0586 
0587             nviolations = sum(S > (ra + 1e-4));
0588             idx = find(S > (ra + 1e-4));
0589         end
0590 
0591         function genoverload = gen_violations(obj, mpc)
0592           define_constants;
0593           genoverload = sum(mpc.gen(:,PG) > mpc.gen(:,PMAX) + 1e-4);
0594         end    end
0595 end