Home > matpower7.0 > lib > 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 https://matpower.org 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 vcart = ~dc && mpopt.opf.v_cartesian;
0034 
0035 %% get indexing
0036 [vv, ll, nne, nni] = om.get_idx();
0037 
0038 if mpopt.verbose > 0
0039     v = mpver('all');
0040     fprintf('\nMATPOWER Version %s, %s', v.Version, v.Date);
0041 end
0042 
0043 %%-----  run DC OPF solver  -----
0044 if dc
0045   if mpopt.verbose > 0
0046     fprintf(' -- DC Optimal Power Flow\n');
0047   end
0048   [results, success, raw] = dcopf_solver(om, mpopt);
0049 else
0050   %%-----  run AC OPF solver  -----
0051   if mpopt.verbose > 0
0052     fprintf(' -- AC Optimal Power Flow\n  AC OPF formulation: ');
0053     if sdp
0054         fprintf('SDP relaxation\n');
0055     else
0056         if vcart
0057             v = 'cartesian';
0058         else
0059             v = 'polar';
0060         end
0061         if mpopt.opf.current_balance
0062             v2 = 'current';
0063         else
0064             v2 = 'power';
0065         end
0066         fprintf('%s voltages, %s balance eqns\n', v, v2);
0067     end
0068   end
0069 
0070   %% ZIP loads?
0071   if (~isempty(mpopt.exp.sys_wide_zip_loads.pw) && ...
0072           any(mpopt.exp.sys_wide_zip_loads.pw(2:3))) || ...
0073           (~isempty(mpopt.exp.sys_wide_zip_loads.qw) && ...
0074           any(mpopt.exp.sys_wide_zip_loads.qw(2:3)))
0075     switch alg
0076     case {'PDIPM', 'TRALM', 'MINOPF', 'SDPOPF'}
0077       warning('opf_execute: ''%s'' solver does not support ZIP load model. Converting to constant power loads.', alg)
0078       mpopt = mpoption(mpopt, 'exp.sys_wide_zip_loads', ...
0079                         struct('pw', [], 'qw', []));
0080     end
0081   end
0082 
0083   %% run specific AC OPF solver
0084   switch alg
0085     case 'MIPS'
0086       [results, success, raw] = mipsopf_solver(om, mpopt);
0087     case 'IPOPT'
0088       if ~have_fcn('ipopt')
0089         error('opf_execute: MPOPT.opf.ac.solver = ''%s'' requires IPOPT (see https://github.com/coin-or/Ipopt)', alg);
0090       end
0091       [results, success, raw] = ipoptopf_solver(om, mpopt);
0092     case 'PDIPM'
0093       if mpopt.pdipm.step_control
0094         if ~have_fcn('scpdipmopf')
0095           error('opf_execute: MPOPT.opf.ac.solver = ''%s'' requires SCPDIPMOPF (see http://www.pserc.cornell.edu/tspopf/)', alg);
0096         end
0097       else
0098         if ~have_fcn('pdipmopf')
0099           error('opf_execute: MPOPT.opf.ac.solver = ''%s'' requires PDIPMOPF (see http://www.pserc.cornell.edu/tspopf/)', alg);
0100         end
0101       end
0102       [results, success, raw] = tspopf_solver(om, mpopt);
0103     case 'TRALM'
0104       if ~have_fcn('tralmopf')
0105         error('opf_execute: MPOPT.opf.ac.solver = ''%s'' requires TRALM (see http://www.pserc.cornell.edu/tspopf/)', alg);
0106       end
0107       [results, success, raw] = tspopf_solver(om, mpopt);
0108     case 'MINOPF'
0109       if ~have_fcn('minopf')
0110         error('opf_execute: MPOPT.opf.ac.solver = ''%s'' requires MINOPF (see http://www.pserc.cornell.edu/minopf/)', alg);
0111       end
0112       [results, success, raw] = mopf_solver(om, mpopt);
0113     case 'FMINCON'
0114       if ~have_fcn('fmincon')
0115         error('opf_execute: MPOPT.opf.ac.solver = ''%s'' requires FMINCON (Optimization Toolbox 2.x or later)', alg);
0116       end
0117       [results, success, raw] = fmincopf_solver(om, mpopt);
0118     case 'KNITRO'
0119       if ~have_fcn('knitro')
0120         error('opf_execute: MPOPT.opf.ac.solver = ''%s'' requires Artelys Knitro (see https://www.artelys.com/solvers/knitro/)', alg);
0121       end
0122       [results, success, raw] = ktropf_solver(om, mpopt);
0123     case 'SDPOPF'
0124       if ~have_fcn('yalmip')
0125         error('opf_execute: MPOPT.opf.ac.solver = ''%s'' requires YALMIP (see https://yalmip.github.io)', alg);
0126       end
0127       [results, success, raw] = sdpopf_solver(om, mpopt);
0128     otherwise
0129       error('opf_execute: MPOPT.opf.ac.solver = ''%s'' is not a valid AC OPF solver selection', alg);
0130   end
0131 end
0132 if ~isfield(raw, 'output') || ~isfield(raw.output, 'alg') || isempty(raw.output.alg)
0133     raw.output.alg = alg;
0134 end
0135 
0136 if success
0137   if ~dc && ~sdp
0138     %% copy bus voltages back to gen matrix
0139     results.gen(:, VG) = results.bus(results.gen(:, GEN_BUS), VM);
0140 
0141     %% gen PQ capability curve multipliers
0142     if ll.N.PQh > 0 || ll.N.PQl > 0
0143       mu_PQh = results.mu.lin.l(ll.i1.PQh:ll.iN.PQh) - results.mu.lin.u(ll.i1.PQh:ll.iN.PQh);
0144       mu_PQl = results.mu.lin.l(ll.i1.PQl:ll.iN.PQl) - results.mu.lin.u(ll.i1.PQl:ll.iN.PQl);
0145       Apqdata = om.get_userdata('Apqdata');
0146       results.gen = update_mupq(results.baseMVA, results.gen, mu_PQh, mu_PQl, Apqdata);
0147     end
0148 
0149     %% compute g, dg, f, df, d2f if requested by opf.return_raw_der = 1
0150     if mpopt.opf.return_raw_der
0151       %% move from results to raw if using v4.0 of MINOPF or TSPOPF
0152       if isfield(results, 'dg')
0153         raw.dg = results.dg;
0154         raw.g = results.g;
0155       end
0156       %% compute g, dg, unless already done by post-v4.0 MINOPF or TSPOPF
0157       if ~isfield(raw, 'dg')
0158         mpc = om.get_mpc();
0159         [Ybus, Yf, Yt] = makeYbus(mpc.baseMVA, mpc.bus, mpc.branch);
0160         [g, geq, dg, dgeq] = opf_consfcn(results.x, om, Ybus, Yf, Yt, mpopt);
0161         raw.g = [ geq; g];
0162         raw.dg = [ dgeq'; dg'];   %% true Jacobian organization
0163       end
0164       %% compute df, d2f
0165       [f, df, d2f] = opf_costfcn(results.x, om);
0166       raw.df = df;
0167       raw.d2f = d2f;
0168     end
0169   end
0170 
0171   %% delete g and dg fieldsfrom results if using v4.0 of MINOPF or TSPOPF
0172   if isfield(results, 'dg')
0173     rmfield(results, 'dg');
0174     rmfield(results, 'g');
0175   end
0176 
0177   %% angle limit constraint multipliers
0178   iang = om.get_userdata('iang');
0179   if ~sdp && length(iang)
0180     if vcart
0181       iang = om.get_userdata('iang');
0182       results.branch(iang, MU_ANGMIN) = results.mu.nli(nni.i1.angL:nni.iN.angL) * pi/180;
0183       results.branch(iang, MU_ANGMAX) = results.mu.nli(nni.i1.angU:nni.iN.angU) * pi/180;
0184     else
0185       iang = om.get_userdata('iang');
0186       results.branch(iang, MU_ANGMIN) = results.mu.lin.l(ll.i1.ang:ll.iN.ang) * pi/180;
0187       results.branch(iang, MU_ANGMAX) = results.mu.lin.u(ll.i1.ang:ll.iN.ang) * pi/180;
0188     end
0189   end
0190 else
0191   %% assign empty g, dg, f, df, d2f if requested by opf.return_raw_der = 1
0192   if ~dc && mpopt.opf.return_raw_der
0193     raw.dg = [];
0194     raw.g = [];
0195     raw.df = [];
0196     raw.d2f = [];
0197   end
0198 end
0199 
0200 if ~sdp
0201   %% assign values and limit shadow prices for variables
0202   om_var_order = om.get('var', 'order');
0203   for k = 1:length(om_var_order)
0204     name = om_var_order(k).name;
0205     if om.getN('var', name)
0206       idx = vv.i1.(name):vv.iN.(name);
0207       results.var.val.(name) = results.x(idx);
0208       results.var.mu.l.(name) = results.mu.var.l(idx);
0209       results.var.mu.u.(name) = results.mu.var.u(idx);
0210     end
0211   end
0212 
0213   %% assign shadow prices for linear constraints
0214   om_lin_order = om.get('lin', 'order');
0215   for k = 1:length(om_lin_order)
0216     name = om_lin_order(k).name;
0217     if om.getN('lin', name)
0218       idx = ll.i1.(name):ll.iN.(name);
0219       results.lin.mu.l.(name) = results.mu.lin.l(idx);
0220       results.lin.mu.u.(name) = results.mu.lin.u(idx);
0221     end
0222   end
0223 
0224   %% assign shadow prices for nonlinear constraints
0225   if ~dc
0226     om_nle_order = om.get('nle', 'order');
0227     for k = 1:length(om_nle_order)
0228       name = om_nle_order(k).name;
0229       if om.getN('nle', name)
0230         results.nle.lambda.(name) = results.mu.nle(nne.i1.(name):nne.iN.(name));
0231       end
0232     end
0233 
0234     om_nli_order = om.get('nli', 'order');
0235     for k = 1:length(om_nli_order)
0236       name = om_nli_order(k).name;
0237       if om.getN('nli', name)
0238         results.nli.mu.(name) = results.mu.nli(nni.i1.(name):nni.iN.(name));
0239       end
0240     end
0241   end
0242 
0243   %% assign values for components of quadratic cost
0244   om_qdc_order = om.get('qdc', 'order');
0245   for k = 1:length(om_qdc_order)
0246     name = om_qdc_order(k).name;
0247     if om.getN('qdc', name)
0248       results.qdc.(name) = om.eval_quad_cost(results.x, name);
0249     end
0250   end
0251 
0252   %% assign values for components of general nonlinear cost
0253   om_nlc_order = om.get('nlc', 'order');
0254   for k = 1:length(om_nlc_order)
0255     name = om_nlc_order(k).name;
0256     if om.getN('nlc', name)
0257       results.nlc.(name) = om.eval_nln_cost(results.x, name);
0258     end
0259   end
0260 
0261   %% assign values for components of legacy user cost
0262   om_cost_order = om.get('cost', 'order');
0263   for k = 1:length(om_cost_order)
0264     name = om_cost_order(k).name;
0265     if om.getN('cost', name)
0266       results.cost.(name) = om.eval_legacy_cost(results.x, name);
0267     end
0268   end
0269 
0270   %% if single-block PWL costs were converted to POLY, insert dummy y into x
0271   %% Note: The "y" portion of x will be nonsense, but everything should at
0272   %%       least be in the expected locations.
0273   pwl1 = om.get_userdata('pwl1');
0274   if ~isempty(pwl1) && ~strcmp(alg, 'TRALM') && ~(strcmp(alg, 'PDIPM') && mpopt.pdipm.step_control)
0275     %% get indexing
0276     vv = om.get_idx();
0277     if dc
0278       nx = vv.iN.Pg;
0279     else
0280       nx = vv.iN.Qg;
0281     end
0282     y = zeros(length(pwl1), 1);
0283     raw.xr = [ raw.xr(1:nx); y; raw.xr(nx+1:end)];
0284     results.x = [ results.x(1:nx); y; results.x(nx+1:end)];
0285   end
0286 end
Generated on Mon 24-Jun-2019 15:58:45 by m2html © 2005