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 https://matpower.org 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 vcart = ~dc && mpopt.opf.v_cartesian;
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 %% define flag to indicate whether we are tied to legacy formulation
0037 %% implemented by legacy MINOS and PDIPM solvers (e.g. with hard-coded
0038 %% costs and constrain
0039 if strcmp(alg, 'MINOPF') || strcmp(alg, 'PDIPM') || ...
0040         strcmp(alg, 'TRALM') || strcmp(alg, 'SDPOPF')
0041     legacy_formulation = 1;
0042     if vcart
0043         error('Option ''opf.v_cartesian'' = 1 is not compatible with ''opf.solver.ac''=''%s''.', alg);
0044     end
0045     if mpopt.opf.current_balance
0046         error('Option ''opf.current_balance'' = 1 is not compatible with ''opf.solver.ac''=''%s''.', alg);
0047     end
0048 else
0049     legacy_formulation = 0;
0050 end
0051 if ~dc && ( ~isempty(mpopt.exp.sys_wide_zip_loads.pw) && ...
0052                     ~isequal(mpopt.exp.sys_wide_zip_loads.pw, [1 0 0]) || ...
0053             ~isempty(mpopt.exp.sys_wide_zip_loads.qw) && ...
0054                     ~isequal(mpopt.exp.sys_wide_zip_loads.qw, [1 0 0]) )
0055     if vcart
0056         warning('Voltage dependent loads are not supported with option ''opf.v_cartesian'' = 1. Reverting to constant power load model.');
0057     end
0058     if mpopt.opf.current_balance
0059         warning('Voltage dependent loads are not supported with option ''opf.current_balance'' = 1. Reverting to constant power load model.');
0060     end
0061 end
0062 
0063 %% data dimensions
0064 nb   = size(mpc.bus, 1);    %% number of buses
0065 nl   = size(mpc.branch, 1); %% number of branches
0066 ng   = size(mpc.gen, 1);    %% number of dispatchable injections
0067 nnle = 0;                   %% number of nonlinear user-defined equality cons
0068 nnli = 0;                   %% number of nonlinear user-defined inequality cons
0069 if isfield(mpc, 'A')
0070   nlin = size(mpc.A, 1);    %% number of linear user constraints
0071 else
0072   nlin = 0;
0073 end
0074 if isfield(mpc, 'N')
0075   nw = size(mpc.N, 1);      %% number of general cost vars, w
0076 else
0077   nw = 0;
0078 end
0079 
0080 if dc
0081   %% ignore reactive costs for DC
0082   mpc.gencost = pqcost(mpc.gencost, ng);
0083 
0084   %% reduce A and/or N from AC dimensions to DC dimensions, if needed
0085   if nlin || nw
0086     acc = [nb+(1:nb) 2*nb+ng+(1:ng)];   %% Vm and Qg columns
0087     if nlin && size(mpc.A, 2) >= 2*nb + 2*ng
0088       %% make sure there aren't any constraints on Vm or Qg
0089       if any(any(mpc.A(:, acc)))
0090         error('opf_setup: attempting to solve DC OPF with user constraints on Vm or Qg');
0091       end
0092       mpc.A(:, acc) = [];               %% delete Vm and Qg columns
0093     end
0094     if nw && size(mpc.N, 2) >= 2*nb + 2*ng
0095       %% make sure there aren't any costs on Vm or Qg
0096       if any(any(mpc.N(:, acc)))
0097         [ii, jj] = find(mpc.N(:, acc));
0098         ii = unique(ii);    %% indices of w with potential non-zero cost terms from Vm or Qg
0099         if any(mpc.Cw(ii)) || (isfield(mpc, 'H') && ~isempty(mpc.H) && ...
0100                 any(any(mpc.H(:, ii))))
0101           error('opf_setup: attempting to solve DC OPF with user costs on Vm or Qg');
0102         end
0103       end
0104       mpc.N(:, acc) = [];               %% delete Vm and Qg columns
0105     end
0106   end
0107 else    %% AC
0108   if use_vg     %% adjust bus voltage limits based on generator Vg setpoint
0109     %% gen connection matrix, element i, j is 1 if, generator j at bus i is ON
0110     Cg = sparse(mpc.gen(:, GEN_BUS), (1:ng)', mpc.gen(:, GEN_STATUS) > 0, nb, ng);
0111     Vbg = Cg * sparse(1:ng, 1:ng, mpc.gen(:, VG), ng, ng);
0112     Vmax = max(Vbg, [], 2); %% zero for non-gen buses, else max Vg of gens @ bus
0113     ib = find(Vmax);                %% buses with online gens
0114     Vmin = max(2*Cg - Vbg, [], 2);  %% same as Vmax, except min Vg of gens @ bus
0115     Vmin(ib) = 2 - Vmin(ib);
0116 
0117     if use_vg == 1      %% use Vg setpoint directly
0118         mpc.bus(ib, VMAX) = Vmax(ib);   %% max set by max Vg @ bus
0119         mpc.bus(ib, VMIN) = Vmin(ib);   %% min set by min Vg @ bus
0120         mpc.bus(ib, VM) = mpc.bus(ib, VMAX);
0121     elseif use_vg > 0 && use_vg < 1     %% fractional value
0122         %% use weighted avg between original Vmin/Vmax limits and Vg
0123         mpc.bus(ib, VMAX) = (1-use_vg) * mpc.bus(ib, VMAX) + use_vg * Vmax(ib);
0124         mpc.bus(ib, VMIN) = (1-use_vg) * mpc.bus(ib, VMIN) + use_vg * Vmin(ib);
0125     else
0126         error('opf_setup: option ''opf.use_vg'' (= %g) cannot be negative or greater than 1', use_vg);
0127     end
0128   end
0129   if isfield(mpc, 'user_constraints')
0130     if isfield(mpc.user_constraints, 'nle')
0131       for k = 1:length(mpc.user_constraints.nle)
0132         nnle = nnle + mpc.user_constraints.nle{k}{2};
0133       end
0134     end
0135     if isfield(mpc.user_constraints, 'nli')
0136       for k = 1:length(mpc.user_constraints.nli)
0137         nnli = nnli + mpc.user_constraints.nli{k}{2};
0138       end
0139     end
0140   end
0141 end
0142 
0143 %% convert single-block piecewise-linear costs into linear polynomial cost
0144 pwl1 = find(mpc.gencost(:, MODEL) == PW_LINEAR & mpc.gencost(:, NCOST) == 2);
0145 % p1 = [];
0146 if ~isempty(pwl1)
0147   x0 = mpc.gencost(pwl1, COST);
0148   y0 = mpc.gencost(pwl1, COST+1);
0149   x1 = mpc.gencost(pwl1, COST+2);
0150   y1 = mpc.gencost(pwl1, COST+3);
0151   m = (y1 - y0) ./ (x1 - x0);
0152   b = y0 - m .* x0;
0153   mpc.gencost(pwl1, MODEL) = POLYNOMIAL;
0154   mpc.gencost(pwl1, NCOST) = 2;
0155   mpc.gencost(pwl1, COST:COST+1) = [m b];
0156 end
0157 
0158 %% create (read-only) copies of individual fields for convenience
0159 [baseMVA, bus, gen, branch, gencost, Au, lbu, ubu, mpopt, ...
0160     N, fparm, H, Cw, z0, zl, zu, userfcn] = opf_args(mpc, mpopt);
0161 
0162 %% warn if there is more than one reference bus
0163 refs = find(bus(:, BUS_TYPE) == REF);
0164 if length(refs) > 1 && mpopt.verbose > 0
0165   errstr = ['\nopf_setup: Warning: Multiple reference buses.\n', ...
0166               '           For a system with islands, a reference bus in each island\n', ...
0167               '           may help convergence, but in a fully connected system such\n', ...
0168               '           a situation is probably not reasonable.\n\n' ];
0169   fprintf(errstr);
0170 end
0171 
0172 %% set up initial variables and bounds
0173 Va   = bus(:, VA) * (pi/180);
0174 Vau = Inf(nb, 1);       %% voltage angle limits
0175 Val = -Vau;
0176 Vau(refs) = Va(refs);   %% voltage angle reference constraints
0177 Val(refs) = Va(refs);
0178 Pg   = gen(:, PG) / baseMVA;
0179 Pmin = gen(:, PMIN) / baseMVA;
0180 Pmax = gen(:, PMAX) / baseMVA;
0181 if ~dc
0182   Vm   = bus(:, VM);
0183   Qg   = gen(:, QG) / baseMVA;
0184   Qmin = gen(:, QMIN) / baseMVA;
0185   Qmax = gen(:, QMAX) / baseMVA;
0186   if vcart
0187     V = Vm .* exp(1j*Va);
0188     Vr = real(V);
0189     Vi = imag(V);
0190   end
0191 end
0192 
0193 %% find/prepare polynomial generator costs
0194 cpg = [];
0195 cqg = [];
0196 [pcost qcost] = pqcost(mpc.gencost, ng);
0197 ip0 = find(pcost(:, MODEL) == POLYNOMIAL & pcost(:, NCOST) == 1);   %% constant
0198 ip1 = find(pcost(:, MODEL) == POLYNOMIAL & pcost(:, NCOST) == 2);   %% linear
0199 ip2 = find(pcost(:, MODEL) == POLYNOMIAL & pcost(:, NCOST) == 3);   %% quadratic
0200 ip3 = find(pcost(:, MODEL) == POLYNOMIAL & pcost(:, NCOST) > 3);    %% cubic or greater
0201 if ~isempty(ip2) || ~isempty(ip1) || ~isempty(ip0)
0202     kpg = zeros(ng, 1);
0203     cpg = zeros(ng, 1);
0204     if ~isempty(ip2)
0205         Qpg = zeros(ng, 1);
0206         Qpg(ip2) = 2 * pcost(ip2, COST) * baseMVA^2;
0207         cpg(ip2) = cpg(ip2) + pcost(ip2, COST+1) * baseMVA;
0208         kpg(ip2) = kpg(ip2) + pcost(ip2, COST+2);
0209     else
0210         Qpg = [];   %% no quadratic terms
0211     end
0212     if ~isempty(ip1)
0213         cpg(ip1) = cpg(ip1) + pcost(ip1, COST) * baseMVA;
0214         kpg(ip1) = kpg(ip1) + pcost(ip1, COST+1);
0215     end
0216     if ~isempty(ip0)
0217         kpg(ip0) = kpg(ip0) + pcost(ip0, COST);
0218     end
0219 end
0220 if ~isempty(qcost)
0221     iq0 = find(qcost(:, MODEL) == POLYNOMIAL & qcost(:, NCOST) == 1);   %% constant
0222     iq1 = find(qcost(:, MODEL) == POLYNOMIAL & qcost(:, NCOST) == 2);   %% linear
0223     iq2 = find(qcost(:, MODEL) == POLYNOMIAL & qcost(:, NCOST) == 3);   %% quadratic
0224     iq3 = find(qcost(:, MODEL) == POLYNOMIAL & qcost(:, NCOST) > 3);    %% cubic or greater
0225     if ~isempty(iq2) || ~isempty(iq1) || ~isempty(iq0)
0226         kqg = zeros(ng, 1);
0227         cqg = zeros(ng, 1);
0228         if ~isempty(iq2)
0229             Qqg = zeros(ng, 1);
0230             Qqg(iq2) = 2 * qcost(iq2, COST) * baseMVA^2;
0231             cqg(iq2) = cqg(iq2) + qcost(iq2, COST+1) * baseMVA;
0232             kqg(iq2) = kqg(iq2) + qcost(iq2, COST+2);
0233         else
0234             Qqg = [];   %% no quadratic terms
0235         end
0236         if ~isempty(iq1)
0237             cqg(iq1) = cqg(iq1) + qcost(iq1, COST) * baseMVA;
0238             kqg(iq1) = kqg(iq1) + qcost(iq1, COST+1);
0239         end
0240         if ~isempty(iq0)
0241             kqg(iq0) = kqg(iq0) + qcost(iq0, COST);
0242         end
0243     end
0244 end
0245 
0246 %% branch voltage angle difference limits
0247 [Aang, lang, uang, iang]  = makeAang(baseMVA, branch, nb, mpopt);
0248 nang = length(iang);
0249 
0250 if dc               %% DC model
0251   %% check generator costs
0252   if ~isempty(ip3)
0253     error('opf_setup: DC OPF cannot handle polynomial costs with higher than quadratic order.');
0254   end
0255 
0256   %% more problem dimensions
0257   nv    = 0;            %% number of voltage magnitude vars
0258   nq    = 0;            %% number of Qg vars
0259   q1    = [];           %% index of 1st Qg column in Ay
0260 
0261   %% power mismatch constraints
0262   [B, Bf, Pbusinj, Pfinj] = makeBdc(baseMVA, bus, branch);
0263   neg_Cg = sparse(gen(:, GEN_BUS), 1:ng, -1, nb, ng);   %% Pbus w.r.t. Pg
0264   Amis = [B neg_Cg];
0265   bmis = -(bus(:, PD) + bus(:, GS)) / baseMVA - Pbusinj;
0266 
0267   %% branch flow constraints
0268   il = find(branch(:, RATE_A) ~= 0 & branch(:, RATE_A) < 1e10);
0269   nl2 = length(il);         %% number of constrained lines
0270   if nl2
0271     upf = branch(il, RATE_A) / baseMVA - Pfinj(il);
0272     upt = branch(il, RATE_A) / baseMVA + Pfinj(il);
0273   else
0274     upf = [];
0275     upt = [];
0276   end
0277 
0278   user_vars = {'Va', 'Pg'};
0279   ycon_vars = {'Pg', 'y'};
0280 else                %% AC model
0281   %% more problem dimensions
0282   nv    = nb;           %% number of voltage magnitude vars
0283   nq    = ng;           %% number of Qg vars
0284   q1    = 1+ng;         %% index of 1st Qg column in Ay
0285 
0286   %% find branches with flow limits
0287   il = find(branch(:, RATE_A) ~= 0 & branch(:, RATE_A) < 1e10);
0288   nl2 = length(il);         %% number of constrained lines
0289 
0290   %% build admittance matrices
0291   [Ybus, Yf, Yt] = makeYbus(baseMVA, bus, branch);
0292 
0293   %% dispatchable load, constant power factor constraints
0294   [Avl, lvl, uvl]  = makeAvl(mpc);
0295 
0296   %% generator PQ capability curve constraints
0297   [Apqh, ubpqh, Apql, ubpql, Apqdata] = makeApq(baseMVA, gen);
0298 
0299   if vcart
0300       user_vars = {'Vr', 'Vi', 'Pg', 'Qg'};
0301       nodal_balance_vars = {'Vr', 'Vi', 'Pg', 'Qg'};
0302       flow_lim_vars = {'Vr', 'Vi'};
0303   else
0304       user_vars = {'Va', 'Vm', 'Pg', 'Qg'};
0305       nodal_balance_vars = {'Va', 'Vm', 'Pg', 'Qg'};
0306       flow_lim_vars = {'Va', 'Vm'};
0307   end
0308   ycon_vars = {'Pg', 'Qg', 'y'};
0309 
0310   %% nonlinear constraint functions
0311   if mpopt.opf.current_balance
0312     mis_cons = {'rImis', 'iImis'};
0313     fcn_mis = @(x)opf_current_balance_fcn(x, mpc, Ybus, mpopt);
0314     hess_mis = @(x, lam)opf_current_balance_hess(x, lam, mpc, Ybus, mpopt);
0315   else
0316     mis_cons = {'Pmis', 'Qmis'};
0317     fcn_mis = @(x)opf_power_balance_fcn(x, mpc, Ybus, mpopt);
0318     hess_mis = @(x, lam)opf_power_balance_hess(x, lam, mpc, Ybus, mpopt);
0319   end
0320   fcn_flow = @(x)opf_branch_flow_fcn(x, mpc, Yf(il, :), Yt(il, :), il, mpopt);
0321   hess_flow = @(x, lam)opf_branch_flow_hess(x, lam, mpc, Yf(il, :), Yt(il, :), il, mpopt);
0322   if vcart
0323     fcn_vref = @(x)opf_vref_fcn(x, mpc, refs, mpopt);
0324     hess_vref = @(x, lam)opf_vref_hess(x, lam, mpc, refs, mpopt);
0325     veq = find(mpc.bus(:, VMIN) == mpc.bus(:, VMAX));
0326     viq = find(mpc.bus(:, VMIN) ~= mpc.bus(:, VMAX));
0327     nveq = length(veq);
0328     nvlims = length(viq);
0329     if nveq
0330       fcn_veq = @(x)opf_veq_fcn(x, mpc, veq, mpopt);
0331       hess_veq = @(x, lam)opf_veq_hess(x, lam, mpc, veq, mpopt);
0332     end
0333     fcn_vlim = @(x)opf_vlim_fcn(x, mpc, viq, mpopt);
0334     hess_vlim = @(x, lam)opf_vlim_hess(x, lam, mpc, viq, mpopt);
0335     fcn_ang = @(x)opf_branch_ang_fcn(x, Aang, lang, uang);
0336     hess_ang = @(x, lam)opf_branch_ang_hess(x, lam, Aang, lang, uang);
0337   end
0338   
0339   %% nonlinear cost functions
0340   if ~isempty(ip3)
0341     cost_Pg = @(x)opf_gen_cost_fcn(x, baseMVA, pcost, ip3, mpopt);
0342   end
0343   if ~isempty(qcost) && ~isempty(iq3)
0344     cost_Qg = @(x)opf_gen_cost_fcn(x, baseMVA, qcost, iq3, mpopt);
0345   end
0346 end
0347 
0348 %% basin constraints for piece-wise linear gen cost variables
0349 if (strcmp(alg, 'PDIPM') && mpopt.pdipm.step_control) || strcmp(alg, 'TRALM')
0350   %% SC-PDIPM or TRALM, no CCV cost vars
0351   ny = 0;
0352   Ay = sparse(0, ng+nq);
0353   by =[];
0354 else
0355   ipwl = find(gencost(:, MODEL) == PW_LINEAR);  %% piece-wise linear costs
0356   ny = size(ipwl, 1);   %% number of piece-wise linear cost vars
0357   [Ay, by] = makeAy(baseMVA, ng, gencost, 1, q1, 1+ng+nq);
0358 end
0359 if any(gencost(:, MODEL) ~= POLYNOMIAL & gencost(:, MODEL) ~= PW_LINEAR)
0360     error('opf_setup: some generator cost rows have invalid MODEL value');
0361 end
0362 
0363 %% more problem dimensions
0364 nx    = nb+nv + ng+nq;  %% number of standard OPF control variables
0365 if nlin
0366   nz = size(mpc.A, 2) - nx; %% number of user z variables
0367   if nz < 0
0368     error('opf_setup: user supplied A matrix must have at least %d columns.', nx);
0369   end
0370 else
0371   nz = 0;               %% number of user z variables
0372   if nw                 %% still need to check number of columns of N
0373     if size(mpc.N, 2) ~= nx;
0374       error('opf_setup: user supplied N matrix must have %d columns.', nx);
0375     end
0376   end
0377 end
0378 
0379 %% construct OPF model object
0380 om = opf_model(mpc);
0381 if ~isempty(pwl1)
0382   om.userdata.pwl1 = pwl1;
0383 end
0384 om.userdata.iang = iang;
0385 if dc
0386   %% user data
0387   om.userdata.Bf = Bf;
0388   om.userdata.Pfinj = Pfinj;
0389 
0390   %% optimization variables
0391   om.add_var('Va', nb, Va, Val, Vau);
0392   om.add_var('Pg', ng, Pg, Pmin, Pmax);
0393 
0394   %% linear constraints
0395   om.add_lin_constraint('Pmis', Amis, bmis, bmis, {'Va', 'Pg'});    %% nb
0396   om.add_lin_constraint('Pf',  Bf(il,:), -upt, upf, {'Va'});        %% nl2
0397   om.add_lin_constraint('ang', Aang, lang, uang, {'Va'});           %% nang
0398 
0399   %% quadratic generator costs
0400   if ~isempty(cpg)
0401     om.add_quad_cost('polPg', Qpg, cpg, kpg, {'Pg'});
0402   end
0403 else
0404   %% user data
0405   om.userdata.Apqdata = Apqdata;
0406 
0407   %% optimization variables
0408   if vcart
0409       Vclim = 1.1 * bus(:, VMAX);
0410       om.add_var('Vr', nb, Vr, -Vclim, Vclim);
0411       om.add_var('Vi', nb, Vi, -Vclim, Vclim);
0412   else
0413       om.add_var('Va', nb, Va, Val, Vau);
0414       om.add_var('Vm', nb, Vm, bus(:, VMIN), bus(:, VMAX));
0415   end
0416   om.add_var('Pg', ng, Pg, Pmin, Pmax);
0417   om.add_var('Qg', ng, Qg, Qmin, Qmax);
0418 
0419   %% nonlinear constraints
0420   om.add_nln_constraint(mis_cons, [nb;nb], 1, fcn_mis, hess_mis, nodal_balance_vars);
0421   if legacy_formulation
0422     om.add_nln_constraint({'Sf', 'St'}, [nl;nl], 0, fcn_flow, hess_flow, flow_lim_vars);
0423   else
0424     om.add_nln_constraint({'Sf', 'St'}, [nl2;nl2], 0, fcn_flow, hess_flow, flow_lim_vars);
0425   end
0426   if vcart
0427     om.userdata.veq = veq;  %% buses with voltage magnitude equality constraints
0428     om.userdata.viq = viq;  %% buses with voltage magnitude limits
0429     om.add_nln_constraint('Vref', length(refs), 1, fcn_vref, hess_vref, {'Vr', 'Vi'});
0430     if nveq
0431       om.add_nln_constraint('Veq', nveq, 1, fcn_veq, hess_veq, {'Vr', 'Vi'});
0432     end
0433     om.add_nln_constraint({'Vmin', 'Vmax'}, [nvlims;nvlims], 0, fcn_vlim, hess_vlim, {'Vr', 'Vi'});
0434     om.add_nln_constraint({'angL', 'angU'}, [nang;nang], 0, fcn_ang, hess_ang, {'Vr', 'Vi'});
0435   end
0436 
0437   %% linear constraints
0438   om.add_lin_constraint('PQh', Apqh, [], ubpqh, {'Pg', 'Qg'});      %% npqh
0439   om.add_lin_constraint('PQl', Apql, [], ubpql, {'Pg', 'Qg'});      %% npql
0440   om.add_lin_constraint('vl',  Avl, lvl, uvl,   {'Pg', 'Qg'});      %% nvl
0441   if ~vcart
0442     om.add_lin_constraint('ang', Aang, lang, uang, {'Va'});         %% nang
0443   end
0444 
0445   %% polynomial generator costs
0446   if ~legacy_formulation
0447     %% quadratic/linear generator costs
0448     if ~isempty(cpg)
0449       om.add_quad_cost('polPg', Qpg, cpg, kpg, {'Pg'});
0450     end
0451     if ~isempty(cqg)
0452       om.add_quad_cost('polQg', Qqg, cqg, kqg, {'Qg'});
0453     end
0454 
0455     %% higher order polynomial generator costs
0456     if ~isempty(ip3)
0457       om.add_nln_cost('polPg', 1, cost_Pg, {'Pg'});
0458     end
0459     if ~isempty(qcost) && ~isempty(iq3)
0460       om.add_nln_cost('polQg', 1, cost_Qg, {'Qg'});
0461     end
0462   end
0463 end
0464 
0465 %% y vars, constraints for piece-wise linear gen costs
0466 if ny > 0
0467   om.add_var('y', ny);
0468   om.add_lin_constraint('ycon', Ay, [], by, ycon_vars);             %% ncony
0469   if dc || ~legacy_formulation
0470     om.add_quad_cost('pwl', [], ones(ny, 1), 0, {'y'});
0471   end
0472 end
0473 
0474 %% add user vars, constraints and costs (as specified via A, ..., N, ...)
0475 if nz > 0
0476   om.add_var('z', nz, z0, zl, zu);
0477   user_vars{end+1} = 'z';
0478 end
0479 if nlin
0480   om.add_lin_constraint('usr', mpc.A, lbu, ubu, user_vars);         %% nlin
0481 end
0482 if nnle
0483   for k = 1:length(mpc.user_constraints.nle)
0484     nlc = mpc.user_constraints.nle{k};
0485     fcn  = eval(['@(x)' nlc{3} '(x, nlc{6}{:})']);
0486     hess = eval(['@(x, lam)' nlc{4} '(x, lam, nlc{6}{:})']);
0487     om.add_nln_constraint(nlc{1:2}, 1, fcn, hess, nlc{5});
0488   end
0489 end
0490 if nnli
0491   for k = 1:length(mpc.user_constraints.nli)
0492     nlc = mpc.user_constraints.nli{k};
0493     fcn  = eval(['@(x)' nlc{3} '(x, nlc{6}{:})']);
0494     hess = eval(['@(x, lam)' nlc{4} '(x, lam, nlc{6}{:})']);
0495     om.add_nln_constraint(nlc{1:2}, 0, fcn, hess, nlc{5});
0496   end
0497 end
0498 if nw
0499   user_cost.N = mpc.N;
0500   user_cost.Cw = Cw;
0501   if ~isempty(fparm)
0502     user_cost.dd = fparm(:, 1);
0503     user_cost.rh = fparm(:, 2);
0504     user_cost.kk = fparm(:, 3);
0505     user_cost.mm = fparm(:, 4);
0506   end
0507   if ~isempty(H)
0508     user_cost.H = H;
0509   end
0510   om.add_legacy_cost('usr', user_cost, user_vars);
0511 end
0512 
0513 %% execute userfcn callbacks for 'formulation' stage
0514 om = run_userfcn(userfcn, 'formulation', om, mpopt);
0515 
0516 %% implement legacy user costs using quadratic or general non-linear costs
0517 cp = om.params_legacy_cost();   %% construct/fetch the parameters
0518 [N, H, Cw, rh, mm] = deal(cp.N, cp.H, cp.Cw, cp.rh, cp.mm);
0519 [nw, nx] = size(N);
0520 if nw
0521     if any(cp.dd ~= 1) || any(cp.kk)    %% not simple quadratic form
0522         if dc                           %% (includes "dead zone" or
0523             if any(cp.dd ~= 1)          %%  quadratic "penalty")
0524                 error('opf_setup: DC OPF can only handle legacy user-defined costs with d = 1');
0525             end
0526             if any(cp.kk)
0527                 error('opf_setup: DC OPF can only handle legacy user-defined costs with no "dead zone", i.e. k = 0');
0528             end
0529         elseif ~legacy_formulation
0530             %% use general nonlinear cost to implement legacy user cost
0531             legacy_cost_fcn = @(x)opf_legacy_user_cost_fcn(x, cp);
0532             om.add_nln_cost('usr', 1, legacy_cost_fcn);
0533         end
0534     else                                %% simple quadratic form
0535         %% use a quadratic cost to implement legacy user cost
0536         if dc || ~legacy_formulation
0537             %% f = 1/2 * w'*H*w + Cw'*w, where w = diag(mm)*(N*x - rh)
0538             %% Let: MN = diag(mm)*N
0539             %%      MR = M * rh
0540             %%      HMR  = H  * MR;
0541             %%      HtMR = H' * MR;
0542             %%  =>   w = MN*x - MR
0543             %% f = 1/2 * (MN*x - MR)'*H*(MN*x - MR) + Cw'*(MN*x - MR)
0544             %%   = 1/2 * x'*MN'*H*MN*x +
0545             %%          (Cw'*MN - 1/2 * MR'*(H+H')*MN)*x +
0546             %%          1/2 * MR'*H*MR - Cw'*MR
0547             %%   = 1/2 * x'*Q*w + c'*x + k
0548             
0549             M    = sparse(1:nw, 1:nw, mm, nw, nw);
0550             MN   = M * N;
0551             MR   = M * rh;
0552             HMR  = H  * MR;
0553             HtMR = H' * MR;
0554             Q = MN' * H * MN;
0555             c = full(MN' * (Cw - 1/2*(HMR+HtMR)));
0556             k = (1/2 * HtMR - Cw)' * MR;
0557             om.add_quad_cost('usr', Q, c, k);
0558         end
0559     end
0560 end
opf_setup

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SOURCE CODE
