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,v 1.5 2010/12/03 15:33:33 cvs Exp $
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  = mpopt(10);        %% PF_DC        : 1 = DC OPF, 0 = AC OPF
0042 alg = mpopt(11);        %% OPF_ALG
0043 verbose = mpopt(31);    %% VERBOSE
0044 
0045 %% define named indices into data matrices
0046 [PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ...
0047     VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus;
0048 [GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, GEN_STATUS, PMAX, PMIN, ...
0049     MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN, PC1, PC2, QC1MIN, QC1MAX, ...
0050     QC2MIN, QC2MAX, RAMP_AGC, RAMP_10, RAMP_30, RAMP_Q, APF] = idx_gen;
0051 [F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, RATE_C, ...
0052     TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST, ...
0053     ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch;
0054 [PW_LINEAR, POLYNOMIAL, MODEL, STARTUP, SHUTDOWN, NCOST, COST] = idx_cost;
0055 
0056 %% data dimensions
0057 nb   = size(mpc.bus, 1);    %% number of buses
0058 nl   = size(mpc.branch, 1); %% number of branches
0059 ng   = size(mpc.gen, 1);    %% number of dispatchable injections
0060 if isfield(mpc, 'A')
0061   nusr = size(mpc.A, 1);    %% number of linear user constraints
0062 else
0063   nusr = 0;
0064 end
0065 if isfield(mpc, 'N')
0066   nw = size(mpc.N, 1);      %% number of general cost vars, w
0067 else
0068   nw = 0;
0069 end
0070 
0071 if dc
0072   %% ignore reactive costs for DC
0073   mpc.gencost = pqcost(mpc.gencost, ng);
0074 
0075   %% reduce A and/or N from AC dimensions to DC dimensions, if needed
0076   if nusr || nw
0077     acc = [nb+(1:nb) 2*nb+ng+(1:ng)];   %% Vm and Qg columns
0078     if nusr && size(mpc.A, 2) >= 2*nb + 2*ng
0079       %% make sure there aren't any constraints on Vm or Qg
0080       if any(any(mpc.A(:, acc)))
0081         error('opf_setup: attempting to solve DC OPF with user constraints on Vm or Qg');
0082       end
0083       mpc.A(:, acc) = [];               %% delete Vm and Qg columns
0084     end
0085     if nw && size(mpc.N, 2) >= 2*nb + 2*ng
0086       %% make sure there aren't any costs on Vm or Qg
0087       if any(any(mpc.N(:, acc)))
0088         [ii, jj] = find(mpc.N(:, acc));
0089         ii = unique(ii);    %% indices of w with potential non-zero cost terms from Vm or Qg
0090         if any(mpc.Cw(ii)) || (isfield(mpc, 'H') && ~isempty(mpc.H) && ...
0091                 any(any(mpc.H(:, ii))))
0092           error('opf_setup: attempting to solve DC OPF with user costs on Vm or Qg');
0093         end
0094       end
0095       mpc.N(:, acc) = [];               %% delete Vm and Qg columns
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 && 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 Vm(gen(:, GEN_BUS)) = gen(:, VG);   %% buses with gens, init Vm from gen data
0133 Pg   = gen(:, PG) / baseMVA;
0134 Qg   = gen(:, QG) / baseMVA;
0135 Pmin = gen(:, PMIN) / baseMVA;
0136 Pmax = gen(:, PMAX) / baseMVA;
0137 Qmin = gen(:, QMIN) / baseMVA;
0138 Qmax = gen(:, QMAX) / baseMVA;
0139 
0140 if dc               %% DC model
0141   %% more problem dimensions
0142   nv    = 0;            %% number of voltage magnitude vars
0143   nq    = 0;            %% number of Qg vars
0144   q1    = [];           %% index of 1st Qg column in Ay
0145 
0146   %% power mismatch constraints
0147   [B, Bf, Pbusinj, Pfinj] = makeBdc(baseMVA, bus, branch);
0148   neg_Cg = sparse(gen(:, GEN_BUS), 1:ng, -1, nb, ng);   %% Pbus w.r.t. Pg
0149   Amis = [B neg_Cg];
0150   bmis = -(bus(:, PD) + bus(:, GS)) / baseMVA - Pbusinj;
0151 
0152   %% branch flow constraints
0153   il = find(branch(:, RATE_A) ~= 0 & branch(:, RATE_A) < 1e10);
0154   nl2 = length(il);         %% number of constrained lines
0155   lpf = -Inf * ones(nl2, 1);
0156   upf = branch(il, RATE_A) / baseMVA - Pfinj(il);
0157   upt = branch(il, RATE_A) / baseMVA + Pfinj(il);
0158 
0159   user_vars = {'Va', 'Pg'};
0160   ycon_vars = {'Pg', 'y'};
0161 else                %% AC model
0162   %% more problem dimensions
0163   nv    = nb;           %% number of voltage magnitude vars
0164   nq    = ng;           %% number of Qg vars
0165   q1    = 1+ng;         %% index of 1st Qg column in Ay
0166 
0167   %% dispatchable load, constant power factor constraints
0168   [Avl, lvl, uvl]  = makeAvl(baseMVA, gen);
0169   
0170   %% generator PQ capability curve constraints
0171   [Apqh, ubpqh, Apql, ubpql, Apqdata] = makeApq(baseMVA, gen);
0172 
0173   user_vars = {'Va', 'Vm', 'Pg', 'Qg'};
0174   ycon_vars = {'Pg', 'Qg', 'y'};
0175 end
0176 
0177 %% voltage angle reference constraints
0178 Vau = Inf * ones(nb, 1);
0179 Val = -Vau;
0180 Vau(refs) = Va(refs);
0181 Val(refs) = Va(refs);
0182 
0183 %% branch voltage angle difference limits
0184 [Aang, lang, uang, iang]  = makeAang(baseMVA, branch, nb, mpopt);
0185 
0186 %% basin constraints for piece-wise linear gen cost variables
0187 if alg == 545 || alg == 550     %% SC-PDIPM or TRALM, no CCV cost vars
0188   ny = 0;
0189   Ay = sparse(0, ng+nq);
0190   by =[];
0191 else
0192   ipwl = find(gencost(:, MODEL) == PW_LINEAR);  %% piece-wise linear costs
0193   ny = size(ipwl, 1);   %% number of piece-wise linear cost vars
0194   [Ay, by] = makeAy(baseMVA, ng, gencost, 1, q1, 1+ng+nq);
0195 end
0196 if any(gencost(:, MODEL) ~= POLYNOMIAL & gencost(:, MODEL) ~= PW_LINEAR)
0197     error('opf_setup: some generator cost rows have invalid MODEL value');
0198 end
0199 
0200 
0201 %% more problem dimensions
0202 nx    = nb+nv + ng+nq;  %% number of standard OPF control variables
0203 if nusr
0204   nz = size(mpc.A, 2) - nx; %% number of user z variables
0205   if nz < 0
0206     error('opf_setup: user supplied A matrix must have at least %d columns.', nx);
0207   end
0208 else
0209   nz = 0;               %% number of user z variables
0210   if nw                 %% still need to check number of columns of N
0211     if size(mpc.N, 2) ~= nx;
0212       error('opf_setup: user supplied N matrix must have %d columns.', nx);
0213     end
0214   end
0215 end
0216 
0217 %% construct OPF model object
0218 om = opf_model(mpc);
0219 if ~isempty(pwl1)
0220   om = userdata(om, 'pwl1', pwl1);
0221 end
0222 if dc
0223   om = userdata(om, 'Bf', Bf);
0224   om = userdata(om, 'Pfinj', Pfinj);
0225   om = userdata(om, 'iang', iang);
0226   om = add_vars(om, 'Va', nb, Va, Val, Vau);
0227   om = add_vars(om, 'Pg', ng, Pg, Pmin, Pmax);
0228   om = add_constraints(om, 'Pmis', Amis, bmis, bmis, {'Va', 'Pg'}); %% nb
0229   om = add_constraints(om, 'Pf',  Bf(il,:), lpf, upf, {'Va'});      %% nl
0230   om = add_constraints(om, 'Pt', -Bf(il,:), lpf, upt, {'Va'});      %% nl
0231   om = add_constraints(om, 'ang', Aang, lang, uang, {'Va'});        %% nang
0232 else
0233   om = userdata(om, 'Apqdata', Apqdata);
0234   om = userdata(om, 'iang', iang);
0235   om = add_vars(om, 'Va', nb, Va, Val, Vau);
0236   om = add_vars(om, 'Vm', nb, Vm, bus(:, VMIN), bus(:, VMAX));
0237   om = add_vars(om, 'Pg', ng, Pg, Pmin, Pmax);
0238   om = add_vars(om, 'Qg', ng, Qg, Qmin, Qmax);
0239   om = add_constraints(om, 'Pmis', nb, 'nonlinear');
0240   om = add_constraints(om, 'Qmis', nb, 'nonlinear');
0241   om = add_constraints(om, 'Sf', nl, 'nonlinear');
0242   om = add_constraints(om, 'St', nl, 'nonlinear');
0243   om = add_constraints(om, 'PQh', Apqh, [], ubpqh, {'Pg', 'Qg'});   %% npqh
0244   om = add_constraints(om, 'PQl', Apql, [], ubpql, {'Pg', 'Qg'});   %% npql
0245   om = add_constraints(om, 'vl',  Avl, lvl, uvl,   {'Pg', 'Qg'});   %% nvl
0246   om = add_constraints(om, 'ang', Aang, lang, uang, {'Va'});        %% nang
0247 end
0248 
0249 %% y vars, constraints for piece-wise linear gen costs
0250 if ny > 0
0251   om = add_vars(om, 'y', ny);
0252   om = add_constraints(om, 'ycon', Ay, [], by, ycon_vars);          %% ncony
0253 end
0254 
0255 %% add user vars, constraints and costs (as specified via A, ..., N, ...)
0256 if nz > 0
0257   om = add_vars(om, 'z', nz, z0, zl, zu);
0258   user_vars{end+1} = 'z';
0259 end
0260 if nusr
0261   om = add_constraints(om, 'usr', mpc.A, lbu, ubu, user_vars);      %% nusr
0262 end
0263 if nw
0264   user_cost.N = mpc.N;
0265   user_cost.Cw = Cw;
0266   if ~isempty(fparm)
0267     user_cost.dd = fparm(:, 1);
0268     user_cost.rh = fparm(:, 2);
0269     user_cost.kk = fparm(:, 3);
0270     user_cost.mm = fparm(:, 4);
0271   end
0272   if ~isempty(H)
0273     user_cost.H = H;
0274   end
0275   om = add_costs(om, 'usr', user_cost, user_vars);
0276 end
0277 
0278 %% execute userfcn callbacks for 'formulation' stage
0279 om = run_userfcn(userfcn, 'formulation', om);
opf_setup

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SOURCE CODE
