sgvm_minqflow

SGVM_MINQFLOW minimize magintude of reactive flows by adding shunts
   OUT = SGVM_MINQFLOW(MPC, OPT)

   WARNING: NOT CURRENTLY USED

0001 function out = sgvm_minqflow(mpc, opt)
0002 %SGVM_MINQFLOW minimize magintude of reactive flows by adding shunts
0003 %   OUT = SGVM_MINQFLOW(MPC, OPT)
0004 %
0005 %   WARNING: NOT CURRENTLY USED
0006 
0007 %   SynGrid
0008 %   Copyright (c) 2018, Power Systems Engineering Research Center (PSERC)
0009 %   by Eran Schweitzer, Arizona State University
0010 %
0011 %   This file is part of SynGrid.
0012 %   Covered by the 3-clause BSD License (see LICENSE file for details).
0013 
0014 %% parse inputs
0015 if nargin < 2
0016     mpopt = mpoption('opf.dc.solver', 'MIPS', 'opf.ac.solver', 'IPOPT', 'mips.step_control', 1);
0017     mpopt.out.all = 0;
0018     mpopt.verbose = 0;
0019     opt.limqflow = [];
0020 %     dqmax         = 200/mpc.baseMVA;
0021 %     percent_pick  = 0.1;
0022 %     tmag          = 0.1; %default minimum shunt is  0.1 MW.
0023 elseif ~isfield(opt, 'limqflow')
0024     mpopt        = opt.mpopt;
0025     opt.limqflow = [];
0026 end
0027 opt.limqflow = opt_default(opt.limqflow);
0028 dqmax        = opt.limqflow.shunt_max/mpc.baseMVA;
0029 percent_lines= opt.limqflow.percent_lines;
0030 percent_pick = opt.limqflow.percent_pick;
0031 tmag         = opt.limqflow.tmag;
0032 %%
0033 [PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ...
0034     VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus;
0035 [F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, RATE_C, ...
0036     TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST, ...
0037     ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch;
0038 nb = size(mpc.bus,1);
0039 nl = size(mpc.branch,1);
0040 out = mpc;
0041 
0042 %%
0043 if QF > size(mpc.branch,2)
0044     % no solution in mpc try solving with softlimits
0045     if isfield(mpc, 'softlims')
0046         mpc = rmfield(mpc, 'softlims');
0047     end
0048     mpopt.opf.softlims.default = 0;
0049     mpc.softlims.RATE_A.hl_mod = 'remove';
0050     mpc.softlims.VMAX.hl_mod = 'remove';
0051     mpc.softlims.VMIN = struct('hl_mod', 'replace', 'hl_val', 0);
0052     if ~toggle_softlims(mpc, 'status')
0053         mpc = toggle_softlims(mpc,'on');
0054     end
0055     r = runopf(mpc,mpopt);
0056     if ~r.success
0057         error('sgvm_minqflow: could not solve OPF even after enabling soft limits.')
0058     end
0059     qold = r.branch(:,QF)/r.baseMVA;
0060     H    = sgvm_acptdf(r);
0061 else
0062     qold = mpc.branch(:,QF)/mpc.baseMVA;
0063     H    = sgvm_acptdf(mpc);
0064 end
0065 %%
0066 tmp = quantile(abs(qold), 1 - percent_lines);
0067 nluse = find(abs(qold) > tmp);
0068 % Hq   = H(nl+1:end,nb+1:end);
0069 Hq   = H(nl+nluse,nb+1:end);
0070 clear H;
0071 nl    = length(nluse);
0072 %% setup problem
0073 vars = struct('qnew', (1:nl)', 'dqinj', (nl+1:nl+nb)', 'dqabs', (nl+nb+1:nl+2*nb)');
0074 total_vars = nl + 2*nb;
0075 %% constraints
0076 % qnew - Hq*dqinj >= +qold
0077 % qnew + Hq*dqinj >= -qold
0078 % dqabs >= +dqinj ----> dqabs - dqinj >= 0
0079 % dqabs >= -dqinj ----> dqabs + dqinj >= 0
0080 
0081 %           dqnew         dqinj       dqabs
0082 prob.A  = [speye(nl)    , -Hq       , sparse(nl,nb);
0083            speye(nl)    , +Hq       , sparse(nl,nb);
0084            sparse(nb,nl), -speye(nb), speye(nb)   ;
0085            sparse(nb,nl), +speye(nb), speye(nb)];
0086 prob.l = [qold; -qold; zeros(2*nb,1)];
0087 prob.u = Inf(2*nl + 2*nb,1);
0088 
0089 %% cost
0090 % objective: sum(qnew) + sum(dqabs)
0091 prob.c = sparse([vars.qnew; vars.dqabs], 1, 1, total_vars,1);
0092 
0093 %% variable limits
0094 % qnew >= 0
0095 % -dqmax <= dqinj <= dqmax
0096 % dqabs >= 0
0097 prob.xmin = [zeros(nl,1); -dqmax*ones(nb,1); zeros(nb,1)];
0098 prob.xmax = [Inf(nl,1); dqmax*ones(nb,1); Inf(nb,1)];
0099 
0100 %% solve
0101 [x, f, eflag, output, lambda] = qps_matpower(prob);
0102 
0103 if eflag
0104     dqinj = x(vars.dqinj)*mpc.baseMVA;
0105     dqabs = x(vars.dqabs)*mpc.baseMVA;
0106     
0107     % threshold number of shunts
0108     threshold = quantile(dqabs,1-percent_pick);
0109     dqinj(dqabs <= threshold) = 0;
0110     
0111     % apply minimum magnitude threshold of tmag
0112     % Note: there shouldn't be must since we already selected the largest
0113     % 10% shunts
0114     mask = (abs(dqinj) > 0) & (abs(dqinj) < tmag);
0115     dqinj(mask) = sign(dqinj(mask))*tmag;
0116     out.bus(:,BS) = out.bus(:,BS) + dqinj;
0117 else
0118     error('sgvm_minqflow: optimization failed to converge.')
0119 end
0120 
0121 %% utility functions
0122 function opt = opt_default(opt)
0123 % setup default options
0124 optdef = struct();
0125 optdef.shunt_max     = 200;
0126 optdef.percent_lines = 0.5; %consider the upper 50% of loaded lines
0127 optdef.percent_pick  = 0.1;
0128 optdef.tmag          = 0.1; %default minimum shunt is  0.1 MW.
0129 
0130 if ~isempty(opt)
0131     opt = struct_compare(optdef, opt);
0132 else
0133     opt = optdef;
0134 end
0135 
0136 function b = struct_compare(a, b)
0137 % compares structure b to structure a.
0138 % if b lacks a field in a it is added
0139 % this is performed recursively, so if if a.x is a structure
0140 % and b has field x, the the function is called on (a.x, b.x)
0141 for f = fieldnames(a).'
0142     if ~isfield(b, f{:})
0143         b.(f{:}) = a.(f{:});
0144     elseif isstruct(a.(f{:}))
0145         b.(f{:}) = struct_compare(a.(f{:}), b.(f{:}));
0146     end
0147 end