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