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

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SOURCE CODE
