0001 function [results, success, raw] = opf_execute(om, mpopt)
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
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
0029
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
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
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
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
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
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
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
0136 end
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
0144 results.gen(:, VG) = results.bus(results.gen(:, GEN_BUS), VM);
0145
0146
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
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
0161 if mpopt.opf.return_raw_der
0162
0163 if isfield(results, 'dg')
0164 raw.dg = results.dg;
0165 raw.g = results.g;
0166 end
0167
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'];
0173 end
0174
0175 [f, df, d2f] = nlp_costfcn(om, results.x);
0176 raw.df = df;
0177 raw.d2f = d2f;
0178 end
0179 end
0180
0181
0182 if isfield(results, 'dg')
0183 rmfield(results, 'dg');
0184 rmfield(results, 'g');
0185 end
0186
0187
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
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
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
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
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
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
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
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
0281
0282
0283 pwl1 = om.get_userdata('pwl1');
0284 if ~isempty(pwl1) && ~strcmp(alg, 'TRALM') && ~(strcmp(alg, 'PDIPM') && mpopt.pdipm.step_control)
0285
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