Home > matpower6.0 > opf_execute.m

opf_execute

PURPOSE ^

OPF_EXECUTE Executes the OPF specified by an OPF model object.

SYNOPSIS ^

function [results, success, raw] = opf_execute(om, mpopt)

DESCRIPTION ^

OPF_EXECUTE  Executes the OPF specified by an OPF model object.
   [RESULTS, SUCCESS, RAW] = OPF_EXECUTE(OM, MPOPT)

   RESULTS are returned with internal indexing, all equipment
   in-service, etc.

   See also OPF, OPF_SETUP.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [results, success, raw] = opf_execute(om, mpopt)
0002 %OPF_EXECUTE  Executes the OPF specified by an OPF model object.
0003 %   [RESULTS, SUCCESS, RAW] = OPF_EXECUTE(OM, MPOPT)
0004 %
0005 %   RESULTS are returned with internal indexing, all equipment
0006 %   in-service, etc.
0007 %
0008 %   See also OPF, OPF_SETUP.
0009 
0010 %   MATPOWER
0011 %   Copyright (c) 2009-2016, Power Systems Engineering Research Center (PSERC)
0012 %   by Ray Zimmerman, PSERC Cornell
0013 %
0014 %   This file is part of MATPOWER.
0015 %   Covered by the 3-clause BSD License (see LICENSE file for details).
0016 %   See http://www.pserc.cornell.edu/matpower/ for more info.
0017 
0018 %% define named indices into data matrices
0019 [PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ...
0020     VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus;
0021 [GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, GEN_STATUS, PMAX, PMIN, ...
0022     MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN, PC1, PC2, QC1MIN, QC1MAX, ...
0023     QC2MIN, QC2MAX, RAMP_AGC, RAMP_10, RAMP_30, RAMP_Q, APF] = idx_gen;
0024 [F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, RATE_C, ...
0025     TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST, ...
0026     ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch;
0027 
0028 %%-----  setup  -----
0029 %% options
0030 dc  = strcmp(upper(mpopt.model), 'DC');
0031 alg = upper(mpopt.opf.ac.solver);
0032 sdp = strcmp(alg, 'SDPOPF');
0033 
0034 %% build user-defined costs
0035 om = build_cost_params(om);
0036 
0037 %% get indexing
0038 [vv, ll, nn] = get_idx(om);
0039 
0040 if mpopt.verbose > 0
0041     v = mpver('all');
0042     fprintf('\nMATPOWER Version %s, %s', v.Version, v.Date);
0043 end
0044 
0045 %%-----  run DC OPF solver  -----
0046 if dc
0047   if mpopt.verbose > 0
0048     fprintf(' -- DC Optimal Power Flow\n');
0049   end
0050   [results, success, raw] = dcopf_solver(om, mpopt);
0051 else
0052   %%-----  run AC OPF solver  -----
0053   if mpopt.verbose > 0
0054     fprintf(' -- AC Optimal Power Flow\n');
0055   end
0056 
0057   %% if opf.ac.solver not set, choose best available option
0058   if strcmp(alg, 'DEFAULT')
0059     if have_fcn('pdipmopf')
0060       alg = 'PDIPM';            %% PDIPM
0061     else
0062       alg = 'MIPS';             %% MIPS
0063     end
0064     mpopt = mpoption(mpopt, 'opf.ac.solver', alg);
0065   end
0066 
0067   %% ZIP loads?
0068   if (~isempty(mpopt.exp.sys_wide_zip_loads.pw) && ...
0069           any(mpopt.exp.sys_wide_zip_loads.pw(2:3))) || ...
0070           (~isempty(mpopt.exp.sys_wide_zip_loads.qw) && ...
0071           any(mpopt.exp.sys_wide_zip_loads.qw(2:3)))
0072     switch alg
0073     case {'PDIPM', 'TRALM', 'MINOPF', 'SDPOPF'}
0074       warning('opf_execute: ''%s'' solver does not support ZIP load model. Converting to constant power loads.', alg)
0075       mpopt = mpoption(mpopt, 'exp.sys_wide_zip_loads', ...
0076                         struct('pw', [], 'qw', []));
0077     end
0078   end
0079 
0080   %% run specific AC OPF solver
0081   switch alg
0082     case 'MIPS'
0083       [results, success, raw] = mipsopf_solver(om, mpopt);
0084     case 'IPOPT'
0085       if ~have_fcn('ipopt')
0086         error('opf_execute: MPOPT.opf.ac.solver = ''%s'' requires IPOPT (see http://www.coin-or.org/projects/Ipopt.xml)', alg);
0087       end
0088       [results, success, raw] = ipoptopf_solver(om, mpopt);
0089     case 'PDIPM'
0090       if mpopt.pdipm.step_control
0091         if ~have_fcn('scpdipmopf')
0092           error('opf_execute: MPOPT.opf.ac.solver = ''%s'' requires SCPDIPMOPF (see http://www.pserc.cornell.edu/tspopf/)', alg);
0093         end
0094       else
0095         if ~have_fcn('pdipmopf')
0096           error('opf_execute: MPOPT.opf.ac.solver = ''%s'' requires PDIPMOPF (see http://www.pserc.cornell.edu/tspopf/)', alg);
0097         end
0098       end
0099       [results, success, raw] = tspopf_solver(om, mpopt);
0100     case 'TRALM'
0101       if ~have_fcn('tralmopf')
0102         error('opf_execute: MPOPT.opf.ac.solver = ''%s'' requires TRALM (see http://www.pserc.cornell.edu/tspopf/)', alg);
0103       end
0104       [results, success, raw] = tspopf_solver(om, mpopt);
0105     case 'MINOPF'
0106       if ~have_fcn('minopf')
0107         error('opf_execute: MPOPT.opf.ac.solver = ''%s'' requires MINOPF (see http://www.pserc.cornell.edu/minopf/)', alg);
0108       end
0109       [results, success, raw] = mopf_solver(om, mpopt);
0110     case 'FMINCON'
0111       if ~have_fcn('fmincon')
0112         error('opf_execute: MPOPT.opf.ac.solver = ''%s'' requires FMINCON (Optimization Toolbox 2.x or later)', alg);
0113       end
0114       [results, success, raw] = fmincopf_solver(om, mpopt);
0115     case 'KNITRO'
0116       if ~have_fcn('knitro')
0117         error('opf_execute: MPOPT.opf.ac.solver = ''%s'' requires KNITRO (see http://www.ziena.com/)', alg);
0118       end
0119       [results, success, raw] = ktropf_solver(om, mpopt);
0120     case 'SDPOPF'
0121       if ~have_fcn('yalmip')
0122         error('opf_execute: MPOPT.opf.ac.solver = ''%s'' requires YALMIP (see http://users.isy.liu.se/johanl/yalmip/)', alg);
0123       end
0124       [results, success, raw] = sdpopf_solver(om, mpopt);
0125     otherwise
0126       error('opf_execute: MPOPT.opf.ac.solver = ''%s'' is not a valid AC OPF solver selection', alg);
0127   end
0128 end
0129 if ~isfield(raw, 'output') || ~isfield(raw.output, 'alg') || isempty(raw.output.alg)
0130     raw.output.alg = alg;
0131 end
0132 
0133 if success
0134   if ~dc && ~sdp
0135     %% copy bus voltages back to gen matrix
0136     results.gen(:, VG) = results.bus(results.gen(:, GEN_BUS), VM);
0137 
0138     %% gen PQ capability curve multipliers
0139     if ll.N.PQh > 0 || ll.N.PQl > 0
0140       mu_PQh = results.mu.lin.l(ll.i1.PQh:ll.iN.PQh) - results.mu.lin.u(ll.i1.PQh:ll.iN.PQh);
0141       mu_PQl = results.mu.lin.l(ll.i1.PQl:ll.iN.PQl) - results.mu.lin.u(ll.i1.PQl:ll.iN.PQl);
0142       Apqdata = userdata(om, 'Apqdata');
0143       results.gen = update_mupq(results.baseMVA, results.gen, mu_PQh, mu_PQl, Apqdata);
0144     end
0145 
0146     %% compute g, dg, f, df, d2f if requested by opf.return_raw_der = 1
0147     if mpopt.opf.return_raw_der
0148       %% move from results to raw if using v4.0 of MINOPF or TSPOPF
0149       if isfield(results, 'dg')
0150         raw.dg = results.dg;
0151         raw.g = results.g;
0152       end
0153       %% compute g, dg, unless already done by post-v4.0 MINOPF or TSPOPF
0154       if ~isfield(raw, 'dg')
0155         mpc = get_mpc(om);
0156         [Ybus, Yf, Yt] = makeYbus(mpc.baseMVA, mpc.bus, mpc.branch);
0157         [g, geq, dg, dgeq] = opf_consfcn(results.x, om, Ybus, Yf, Yt, mpopt);
0158         raw.g = [ geq; g];
0159         raw.dg = [ dgeq'; dg'];   %% true Jacobian organization
0160       end
0161       %% compute df, d2f
0162       [f, df, d2f] = opf_costfcn(results.x, om);
0163       raw.df = df;
0164       raw.d2f = d2f;
0165     end
0166   end
0167 
0168   %% delete g and dg fieldsfrom results if using v4.0 of MINOPF or TSPOPF
0169   if isfield(results, 'dg')
0170     rmfield(results, 'dg');
0171     rmfield(results, 'g');
0172   end
0173 
0174   %% angle limit constraint multipliers
0175   if ~sdp && ll.N.ang > 0
0176     iang = userdata(om, 'iang');
0177     results.branch(iang, MU_ANGMIN) = results.mu.lin.l(ll.i1.ang:ll.iN.ang) * pi/180;
0178     results.branch(iang, MU_ANGMAX) = results.mu.lin.u(ll.i1.ang:ll.iN.ang) * pi/180;
0179   end
0180 else
0181   %% assign empty g, dg, f, df, d2f if requested by opf.return_raw_der = 1
0182   if ~dc && mpopt.opf.return_raw_der
0183     raw.dg = [];
0184     raw.g = [];
0185     raw.df = [];
0186     raw.d2f = [];
0187   end
0188 end
0189 
0190 if ~sdp
0191   %% assign values and limit shadow prices for variables
0192   om_var_order = get(om, 'var', 'order');
0193   for k = 1:length(om_var_order)
0194     name = om_var_order(k).name;
0195     if getN(om, 'var', name)
0196       idx = vv.i1.(name):vv.iN.(name);
0197       results.var.val.(name) = results.x(idx);
0198       results.var.mu.l.(name) = results.mu.var.l(idx);
0199       results.var.mu.u.(name) = results.mu.var.u(idx);
0200     end
0201   end
0202 
0203   %% assign shadow prices for linear constraints
0204   om_lin_order = get(om, 'lin', 'order');
0205   for k = 1:length(om_lin_order)
0206     name = om_lin_order(k).name;
0207     if getN(om, 'lin', name)
0208       idx = ll.i1.(name):ll.iN.(name);
0209       results.lin.mu.l.(name) = results.mu.lin.l(idx);
0210       results.lin.mu.u.(name) = results.mu.lin.u(idx);
0211     end
0212   end
0213 
0214   %% assign shadow prices for nonlinear constraints
0215   if ~dc
0216     om_nln_order = get(om, 'nln', 'order');
0217     for k = 1:length(om_nln_order)
0218       name = om_nln_order(k).name;
0219       if getN(om, 'nln', name)
0220         idx = nn.i1.(name):nn.iN.(name);
0221         results.nln.mu.l.(name) = results.mu.nln.l(idx);
0222         results.nln.mu.u.(name) = results.mu.nln.u(idx);
0223       end
0224     end
0225   end
0226 
0227   %% assign values for components of user cost
0228   om_cost_order = get(om, 'cost', 'order');
0229   for k = 1:length(om_cost_order)
0230     name = om_cost_order(k).name;
0231     if getN(om, 'cost', name)
0232       results.cost.(name) = compute_cost(om, results.x, name);
0233     end
0234   end
0235 
0236   %% if single-block PWL costs were converted to POLY, insert dummy y into x
0237   %% Note: The "y" portion of x will be nonsense, but everything should at
0238   %%       least be in the expected locations.
0239   pwl1 = userdata(om, 'pwl1');
0240   if ~isempty(pwl1) && ~strcmp(alg, 'TRALM') && ~(strcmp(alg, 'PDIPM') && mpopt.pdipm.step_control)
0241     %% get indexing
0242     vv = get_idx(om);
0243     if dc
0244       nx = vv.iN.Pg;
0245     else
0246       nx = vv.iN.Qg;
0247     end
0248     y = zeros(length(pwl1), 1);
0249     raw.xr = [ raw.xr(1:nx); y; raw.xr(nx+1:end)];
0250     results.x = [ results.x(1:nx); y; results.x(nx+1:end)];
0251   end
0252 end
Generated on Fri 16-Dec-2016 12:45:37 by m2html © 2005