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 sdp = strcmp(alg, 'SDPOPF');
0033
0034
0035 om = build_cost_params(om);
0036
0037
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
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
0053 if mpopt.verbose > 0
0054 fprintf(' -- AC Optimal Power Flow\n');
0055 end
0056
0057
0058 if strcmp(alg, 'DEFAULT')
0059 if have_fcn('pdipmopf')
0060 alg = 'PDIPM';
0061 else
0062 alg = 'MIPS';
0063 end
0064 mpopt = mpoption(mpopt, 'opf.ac.solver', alg);
0065 end
0066
0067
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
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
0136 results.gen(:, VG) = results.bus(results.gen(:, GEN_BUS), VM);
0137
0138
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
0147 if mpopt.opf.return_raw_der
0148
0149 if isfield(results, 'dg')
0150 raw.dg = results.dg;
0151 raw.g = results.g;
0152 end
0153
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'];
0160 end
0161
0162 [f, df, d2f] = opf_costfcn(results.x, om);
0163 raw.df = df;
0164 raw.d2f = d2f;
0165 end
0166 end
0167
0168
0169 if isfield(results, 'dg')
0170 rmfield(results, 'dg');
0171 rmfield(results, 'g');
0172 end
0173
0174
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
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
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
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
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
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
0237
0238
0239 pwl1 = userdata(om, 'pwl1');
0240 if ~isempty(pwl1) && ~strcmp(alg, 'TRALM') && ~(strcmp(alg, 'PDIPM') && mpopt.pdipm.step_control)
0241
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