Description of opf

0001 function om = opf_setup(mpc, mpopt)
0002 %OPF  Constructs an OPF model object from a MATPOWER case struct.
0003 %   OM = OPF_SETUP(MPC, MPOPT)
0004 %
0005 %   Assumes that MPC is a MATPOWER case struct with internal indexing,
0006 %   all equipment in-service, etc.
0007 %
0008 %   See also OPF, EXT2INT, OPF_EXECUTE.
0009 
0010 %   MATPOWER
0011 %   Copyright (c) 1996-2016, Power Systems Engineering Research Center (PSERC)
0012 %   by Ray Zimmerman, PSERC Cornell
0013 %   and Carlos E. Murillo-Sanchez, PSERC Cornell & Universidad Nacional de Colombia
0014 %
0015 %   This file is part of MATPOWER.
0016 %   Covered by the 3-clause BSD License (see LICENSE file for details).
0017 %   See http://www.pserc.cornell.edu/matpower/ for more info.
0018 
0019 %% options
0020 dc  = strcmp(upper(mpopt.model), 'DC');
0021 alg = upper(mpopt.opf.ac.solver);
0022 use_vg = mpopt.opf.use_vg;
0023 
0024 %% define named indices into data matrices
0025 [PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ...
0026     VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus;
0027 [GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, GEN_STATUS, PMAX, PMIN, ...
0028     MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN, PC1, PC2, QC1MIN, QC1MAX, ...
0029     QC2MIN, QC2MAX, RAMP_AGC, RAMP_10, RAMP_30, RAMP_Q, APF] = idx_gen;
0030 [F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, RATE_C, ...
0031     TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST, ...
0032     ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch;
0033 [PW_LINEAR, POLYNOMIAL, MODEL, STARTUP, SHUTDOWN, NCOST, COST] = idx_cost;
0034 
0035 %% data dimensions
0036 nb   = size(mpc.bus, 1);    %% number of buses
0037 nl   = size(mpc.branch, 1); %% number of branches
0038 ng   = size(mpc.gen, 1);    %% number of dispatchable injections
0039 if isfield(mpc, 'A')
0040   nusr = size(mpc.A, 1);    %% number of linear user constraints
0041 else
0042   nusr = 0;
0043 end
0044 if isfield(mpc, 'N')
0045   nw = size(mpc.N, 1);      %% number of general cost vars, w
0046 else
0047   nw = 0;
0048 end
0049 
0050 if dc
0051   %% ignore reactive costs for DC
0052   mpc.gencost = pqcost(mpc.gencost, ng);
0053 
0054   %% reduce A and/or N from AC dimensions to DC dimensions, if needed
0055   if nusr || nw
0056     acc = [nb+(1:nb) 2*nb+ng+(1:ng)];   %% Vm and Qg columns
0057     if nusr && size(mpc.A, 2) >= 2*nb + 2*ng
0058       %% make sure there aren't any constraints on Vm or Qg
0059       if any(any(mpc.A(:, acc)))
0060         error('opf_setup: attempting to solve DC OPF with user constraints on Vm or Qg');
0061       end
0062       mpc.A(:, acc) = [];               %% delete Vm and Qg columns
0063     end
0064     if nw && size(mpc.N, 2) >= 2*nb + 2*ng
0065       %% make sure there aren't any costs on Vm or Qg
0066       if any(any(mpc.N(:, acc)))
0067         [ii, jj] = find(mpc.N(:, acc));
0068         ii = unique(ii);    %% indices of w with potential non-zero cost terms from Vm or Qg
0069         if any(mpc.Cw(ii)) || (isfield(mpc, 'H') && ~isempty(mpc.H) && ...
0070                 any(any(mpc.H(:, ii))))
0071           error('opf_setup: attempting to solve DC OPF with user costs on Vm or Qg');
0072         end
0073       end
0074       mpc.N(:, acc) = [];               %% delete Vm and Qg columns
0075     end
0076   end
0077 else    %% AC
0078   if use_vg     %% adjust bus voltage limits based on generator Vg setpoint
0079     %% gen connection matrix, element i, j is 1 if, generator j at bus i is ON
0080     Cg = sparse(mpc.gen(:, GEN_BUS), (1:ng)', mpc.gen(:, GEN_STATUS) > 0, nb, ng);
0081     Vbg = Cg * sparse(1:ng, 1:ng, mpc.gen(:, VG), ng, ng);
0082     Vmax = max(Vbg, [], 2); %% zero for non-gen buses, else max Vg of gens @ bus
0083     ib = find(Vmax);                %% buses with online gens
0084     Vmin = max(2*Cg - Vbg, [], 2);  %% same as Vmax, except min Vg of gens @ bus
0085     Vmin(ib) = 2 - Vmin(ib);
0086 
0087     if use_vg == 1      %% use Vg setpoint directly
0088         mpc.bus(ib, VMAX) = Vmax(ib);   %% max set by max Vg @ bus
0089         mpc.bus(ib, VMIN) = Vmin(ib);   %% min set by min Vg @ bus
0090     elseif use_vg > 0 && use_vg < 1     %% fractional value
0091         %% use weighted avg between original Vmin/Vmax limits and Vg
0092         mpc.bus(ib, VMAX) = (1-use_vg) * mpc.bus(ib, VMAX) + use_vg * Vmax(ib);
0093         mpc.bus(ib, VMIN) = (1-use_vg) * mpc.bus(ib, VMIN) + use_vg * Vmin(ib);
0094     else
0095         error('opf_setup: option ''opf.use_vg'' (= %g) cannot be negative or greater than 1', use_vg);
0096     end
0097   end
0098 end
0099 
0100 %% convert single-block piecewise-linear costs into linear polynomial cost
0101 pwl1 = find(mpc.gencost(:, MODEL) == PW_LINEAR & mpc.gencost(:, NCOST) == 2);
0102 % p1 = [];
0103 if ~isempty(pwl1)
0104   x0 = mpc.gencost(pwl1, COST);
0105   y0 = mpc.gencost(pwl1, COST+1);
0106   x1 = mpc.gencost(pwl1, COST+2);
0107   y1 = mpc.gencost(pwl1, COST+3);
0108   m = (y1 - y0) ./ (x1 - x0);
0109   b = y0 - m .* x0;
0110   mpc.gencost(pwl1, MODEL) = POLYNOMIAL;
0111   mpc.gencost(pwl1, NCOST) = 2;
0112   mpc.gencost(pwl1, COST:COST+1) = [m b];
0113 end
0114 
0115 %% create (read-only) copies of individual fields for convenience
0116 [baseMVA, bus, gen, branch, gencost, Au, lbu, ubu, mpopt, ...
0117     N, fparm, H, Cw, z0, zl, zu, userfcn] = opf_args(mpc, mpopt);
0118 
0119 %% warn if there is more than one reference bus
0120 refs = find(bus(:, BUS_TYPE) == REF);
0121 if length(refs) > 1 && mpopt.verbose > 0
0122   errstr = ['\nopf_setup: Warning: Multiple reference buses.\n', ...
0123               '           For a system with islands, a reference bus in each island\n', ...
0124               '           may help convergence, but in a fully connected system such\n', ...
0125               '           a situation is probably not reasonable.\n\n' ];
0126   fprintf(errstr);
0127 end
0128 
0129 %% set up initial variables and bounds
0130 Va   = bus(:, VA) * (pi/180);
0131 Vm   = bus(:, VM);
0132 Pg   = gen(:, PG) / baseMVA;
0133 Qg   = gen(:, QG) / baseMVA;
0134 Pmin = gen(:, PMIN) / baseMVA;
0135 Pmax = gen(:, PMAX) / baseMVA;
0136 Qmin = gen(:, QMIN) / baseMVA;
0137 Qmax = gen(:, QMAX) / baseMVA;
0138 
0139 if dc               %% DC model
0140   %% more problem dimensions
0141   nv    = 0;            %% number of voltage magnitude vars
0142   nq    = 0;            %% number of Qg vars
0143   q1    = [];           %% index of 1st Qg column in Ay
0144 
0145   %% power mismatch constraints
0146   [B, Bf, Pbusinj, Pfinj] = makeBdc(baseMVA, bus, branch);
0147   neg_Cg = sparse(gen(:, GEN_BUS), 1:ng, -1, nb, ng);   %% Pbus w.r.t. Pg
0148   Amis = [B neg_Cg];
0149   bmis = -(bus(:, PD) + bus(:, GS)) / baseMVA - Pbusinj;
0150 
0151   %% branch flow constraints
0152   il = find(branch(:, RATE_A) ~= 0 & branch(:, RATE_A) < 1e10);
0153   nl2 = length(il);         %% number of constrained lines
0154   if nl2
0155     upf = branch(il, RATE_A) / baseMVA - Pfinj(il);
0156     upt = branch(il, RATE_A) / baseMVA + Pfinj(il);
0157   else
0158     upf = [];
0159     upt = [];
0160   end
0161 
0162   user_vars = {'Va', 'Pg'};
0163   ycon_vars = {'Pg', 'y'};
0164 else                %% AC model
0165   %% more problem dimensions
0166   nv    = nb;           %% number of voltage magnitude vars
0167   nq    = ng;           %% number of Qg vars
0168   q1    = 1+ng;         %% index of 1st Qg column in Ay
0169 
0170   %% dispatchable load, constant power factor constraints
0171   [Avl, lvl, uvl]  = makeAvl(baseMVA, gen);
0172   
0173   %% generator PQ capability curve constraints
0174   [Apqh, ubpqh, Apql, ubpql, Apqdata] = makeApq(baseMVA, gen);
0175 
0176   user_vars = {'Va', 'Vm', 'Pg', 'Qg'};
0177   ycon_vars = {'Pg', 'Qg', 'y'};
0178 end
0179 
0180 %% voltage angle reference constraints
0181 Vau = Inf(nb, 1);
0182 Val = -Vau;
0183 Vau(refs) = Va(refs);
0184 Val(refs) = Va(refs);
0185 
0186 %% branch voltage angle difference limits
0187 [Aang, lang, uang, iang]  = makeAang(baseMVA, branch, nb, mpopt);
0188 
0189 %% basin constraints for piece-wise linear gen cost variables
0190 if (strcmp(alg, 'PDIPM') && mpopt.pdipm.step_control) || strcmp(alg, 'TRALM')
0191   %% SC-PDIPM or TRALM, no CCV cost vars
0192   ny = 0;
0193   Ay = sparse(0, ng+nq);
0194   by =[];
0195 else
0196   ipwl = find(gencost(:, MODEL) == PW_LINEAR);  %% piece-wise linear costs
0197   ny = size(ipwl, 1);   %% number of piece-wise linear cost vars
0198   [Ay, by] = makeAy(baseMVA, ng, gencost, 1, q1, 1+ng+nq);
0199 end
0200 if any(gencost(:, MODEL) ~= POLYNOMIAL & gencost(:, MODEL) ~= PW_LINEAR)
0201     error('opf_setup: some generator cost rows have invalid MODEL value');
0202 end
0203 
0204 
0205 %% more problem dimensions
0206 nx    = nb+nv + ng+nq;  %% number of standard OPF control variables
0207 if nusr
0208   nz = size(mpc.A, 2) - nx; %% number of user z variables
0209   if nz < 0
0210     error('opf_setup: user supplied A matrix must have at least %d columns.', nx);
0211   end
0212 else
0213   nz = 0;               %% number of user z variables
0214   if nw                 %% still need to check number of columns of N
0215     if size(mpc.N, 2) ~= nx;
0216       error('opf_setup: user supplied N matrix must have %d columns.', nx);
0217     end
0218   end
0219 end
0220 
0221 %% construct OPF model object
0222 om = opf_model(mpc);
0223 if ~isempty(pwl1)
0224   om = userdata(om, 'pwl1', pwl1);
0225 end
0226 if dc
0227   om = userdata(om, 'Bf', Bf);
0228   om = userdata(om, 'Pfinj', Pfinj);
0229   om = userdata(om, 'iang', iang);
0230   om = add_vars(om, 'Va', nb, Va, Val, Vau);
0231   om = add_vars(om, 'Pg', ng, Pg, Pmin, Pmax);
0232   om = add_constraints(om, 'Pmis', Amis, bmis, bmis, {'Va', 'Pg'}); %% nb
0233   om = add_constraints(om, 'Pf',  Bf(il,:), -upt, upf, {'Va'});     %% nl2
0234   om = add_constraints(om, 'ang', Aang, lang, uang, {'Va'});        %% nang
0235 else
0236   om = userdata(om, 'Apqdata', Apqdata);
0237   om = userdata(om, 'iang', iang);
0238   om = add_vars(om, 'Va', nb, Va, Val, Vau);
0239   om = add_vars(om, 'Vm', nb, Vm, bus(:, VMIN), bus(:, VMAX));
0240   om = add_vars(om, 'Pg', ng, Pg, Pmin, Pmax);
0241   om = add_vars(om, 'Qg', ng, Qg, Qmin, Qmax);
0242   om = add_constraints(om, 'Pmis', nb, 'nonlinear');
0243   om = add_constraints(om, 'Qmis', nb, 'nonlinear');
0244   om = add_constraints(om, 'Sf', nl, 'nonlinear');
0245   om = add_constraints(om, 'St', nl, 'nonlinear');
0246   om = add_constraints(om, 'PQh', Apqh, [], ubpqh, {'Pg', 'Qg'});   %% npqh
0247   om = add_constraints(om, 'PQl', Apql, [], ubpql, {'Pg', 'Qg'});   %% npql
0248   om = add_constraints(om, 'vl',  Avl, lvl, uvl,   {'Pg', 'Qg'});   %% nvl
0249   om = add_constraints(om, 'ang', Aang, lang, uang, {'Va'});        %% nang
0250 end
0251 
0252 %% y vars, constraints for piece-wise linear gen costs
0253 if ny > 0
0254   om = add_vars(om, 'y', ny);
0255   om = add_constraints(om, 'ycon', Ay, [], by, ycon_vars);          %% ncony
0256 end
0257 
0258 %% add user vars, constraints and costs (as specified via A, ..., N, ...)
0259 if nz > 0
0260   om = add_vars(om, 'z', nz, z0, zl, zu);
0261   user_vars{end+1} = 'z';
0262 end
0263 if nusr
0264   om = add_constraints(om, 'usr', mpc.A, lbu, ubu, user_vars);      %% nusr
0265 end
0266 if nw
0267   user_cost.N = mpc.N;
0268   user_cost.Cw = Cw;
0269   if ~isempty(fparm)
0270     user_cost.dd = fparm(:, 1);
0271     user_cost.rh = fparm(:, 2);
0272     user_cost.kk = fparm(:, 3);
0273     user_cost.mm = fparm(:, 4);
0274   end
0275   if ~isempty(H)
0276     user_cost.H = H;
0277   end
0278   om = add_costs(om, 'usr', user_cost, user_vars);
0279 end
0280 
0281 %% execute userfcn callbacks for 'formulation' stage
0282 om = run_userfcn(userfcn, 'formulation', om);
opf_setup

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SOURCE CODE
