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

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SOURCE CODE
