0001 function [results, success, raw] = ipoptopf_solver(om, mpopt)
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027
0028
0029
0030
0031
0032
0033
0034
0035
0036
0037
0038
0039
0040
0041
0042
0043
0044
0045
0046
0047
0048 [PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ...
0049 VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus;
0050 [GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, GEN_STATUS, PMAX, PMIN, ...
0051 MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN, PC1, PC2, QC1MIN, QC1MAX, ...
0052 QC2MIN, QC2MAX, RAMP_AGC, RAMP_10, RAMP_30, RAMP_Q, APF] = idx_gen;
0053 [F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, RATE_C, ...
0054 TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST, ...
0055 ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch;
0056 [PW_LINEAR, POLYNOMIAL, MODEL, STARTUP, SHUTDOWN, NCOST, COST] = idx_cost;
0057
0058
0059 mpc = get_mpc(om);
0060 [baseMVA, bus, gen, branch, gencost] = ...
0061 deal(mpc.baseMVA, mpc.bus, mpc.gen, mpc.branch, mpc.gencost);
0062 [vv, ll, nn] = get_idx(om);
0063
0064
0065 nb = size(bus, 1);
0066 ng = size(gen, 1);
0067 nl = size(branch, 1);
0068 ny = getN(om, 'var', 'y');
0069
0070
0071 [A, l, u] = linear_constraints(om);
0072
0073
0074 [x0, xmin, xmax] = getv(om);
0075
0076
0077 [Ybus, Yf, Yt] = makeYbus(baseMVA, bus, branch);
0078
0079
0080 if mpopt.opf.init_from_mpc ~= 1
0081 ll = xmin; uu = xmax;
0082 ll(xmin == -Inf) = -1e10;
0083 uu(xmax == Inf) = 1e10;
0084 x0 = (ll + uu) / 2;
0085 k = find(xmin == -Inf & xmax < Inf);
0086 x0(k) = xmax(k) - 1;
0087 k = find(xmin > -Inf & xmax == Inf);
0088 x0(k) = xmin(k) + 1;
0089 Varefs = bus(bus(:, BUS_TYPE) == REF, VA) * (pi/180);
0090 x0(vv.i1.Va:vv.iN.Va) = Varefs(1);
0091 if ny > 0
0092 ipwl = find(gencost(:, MODEL) == PW_LINEAR);
0093
0094
0095 c = gencost(sub2ind(size(gencost), ipwl, NCOST+2*gencost(ipwl, NCOST)));
0096 x0(vv.i1.y:vv.iN.y) = max(c) + 0.1 * abs(max(c));
0097
0098 end
0099 end
0100
0101
0102 il = find(branch(:, RATE_A) ~= 0 & branch(:, RATE_A) < 1e10);
0103 nl2 = length(il);
0104 nx = length(x0);
0105
0106
0107
0108 kk = find(xmin(nb+1:end) == xmax(nb+1:end));
0109 nk = length(kk);
0110 if nk
0111 kk = kk + nb;
0112 A = [ A; sparse((1:nk)', kk, 1, nk, nx) ];
0113 l = [ l; xmin(kk) ];
0114 u = [ u; xmax(kk) ];
0115 xmin(kk) = -Inf;
0116 xmax(kk) = Inf;
0117 end
0118
0119
0120
0121 nA = size(A, 1);
0122 f = branch(:, F_BUS);
0123 t = branch(:, T_BUS);
0124 Cf = sparse(1:nl, f, ones(nl, 1), nl, nb);
0125 Ct = sparse(1:nl, t, ones(nl, 1), nl, nb);
0126 Cl = Cf + Ct;
0127 Cb = Cl' * Cl + speye(nb);
0128 Cl2 = Cl(il, :);
0129 Cg = sparse(gen(:, GEN_BUS), (1:ng)', 1, nb, ng);
0130 nz = nx - 2*(nb+ng);
0131 nxtra = nx - 2*nb;
0132 Js = [
0133 Cb Cb Cg sparse(nb,ng) sparse(nb,nz);
0134 Cb Cb sparse(nb,ng) Cg sparse(nb,nz);
0135 Cl2 Cl2 sparse(nl2, 2*ng) sparse(nl2,nz);
0136 Cl2 Cl2 sparse(nl2, 2*ng) sparse(nl2,nz);
0137 A;
0138 ];
0139 [f, df, d2f] = opf_costfcn(x0, om);
0140 Hs = tril(d2f + [
0141 Cb Cb sparse(nb,nxtra);
0142 Cb Cb sparse(nb,nxtra);
0143 sparse(nxtra,nx);
0144 ]);
0145
0146
0147 options.ipopt = ipopt_options([], mpopt);
0148
0149
0150 options.auxdata = struct( ...
0151 'om', om, ...
0152 'Ybus', Ybus, ...
0153 'Yf', Yf(il,:), ...
0154 'Yt', Yt(il,:), ...
0155 'mpopt', mpopt, ...
0156 'il', il, ...
0157 'A', A, ...
0158 'nA', nA, ...
0159 'neqnln', 2*nb, ...
0160 'niqnln', 2*nl2, ...
0161 'Js', Js, ...
0162 'Hs', Hs );
0163
0164
0165
0166
0167
0168
0169
0170
0171
0172
0173
0174
0175
0176
0177
0178
0179
0180
0181
0182
0183
0184 options.lb = xmin;
0185 options.ub = xmax;
0186 options.cl = [zeros(2*nb, 1); -Inf(2*nl2, 1); l];
0187 options.cu = [zeros(2*nb, 1); zeros(2*nl2, 1); u];
0188
0189
0190 funcs.objective = @objective;
0191 funcs.gradient = @gradient;
0192 funcs.constraints = @constraints;
0193 funcs.jacobian = @jacobian;
0194 funcs.hessian = @hessian;
0195 funcs.jacobianstructure = @(d) Js;
0196 funcs.hessianstructure = @(d) Hs;
0197
0198
0199
0200
0201 if have_fcn('ipopt_auxdata')
0202 [x, info] = ipopt_auxdata(x0,funcs,options);
0203 else
0204 [x, info] = ipopt(x0,funcs,options);
0205 end
0206
0207 if info.status == 0 || info.status == 1
0208 success = 1;
0209 else
0210 success = 0;
0211 end
0212 if isfield(info, 'iter')
0213 output.iterations = info.iter;
0214 else
0215 output.iterations = [];
0216 end
0217 f = opf_costfcn(x, om);
0218
0219
0220 Va = x(vv.i1.Va:vv.iN.Va);
0221 Vm = x(vv.i1.Vm:vv.iN.Vm);
0222 Pg = x(vv.i1.Pg:vv.iN.Pg);
0223 Qg = x(vv.i1.Qg:vv.iN.Qg);
0224 V = Vm .* exp(1j*Va);
0225
0226
0227
0228 bus(:, VA) = Va * 180/pi;
0229 bus(:, VM) = Vm;
0230 gen(:, PG) = Pg * baseMVA;
0231 gen(:, QG) = Qg * baseMVA;
0232 gen(:, VG) = Vm(gen(:, GEN_BUS));
0233
0234
0235 Sf = V(branch(:, F_BUS)) .* conj(Yf * V);
0236 St = V(branch(:, T_BUS)) .* conj(Yt * V);
0237 branch(:, PF) = real(Sf) * baseMVA;
0238 branch(:, QF) = imag(Sf) * baseMVA;
0239 branch(:, PT) = real(St) * baseMVA;
0240 branch(:, QT) = imag(St) * baseMVA;
0241
0242
0243
0244 muSf = zeros(nl, 1);
0245 muSt = zeros(nl, 1);
0246 if ~isempty(il)
0247 muSf(il) = 2 * info.lambda(2*nb+ (1:nl2)) .* branch(il, RATE_A) / baseMVA;
0248 muSt(il) = 2 * info.lambda(2*nb+nl2+(1:nl2)) .* branch(il, RATE_A) / baseMVA;
0249 end
0250
0251
0252
0253 if nk
0254 lam_tmp = info.lambda(2*nb+2*nl2+nA-nk+(1:nk));
0255 kl = find(lam_tmp < 0);
0256 ku = find(lam_tmp > 0);
0257 info.zl(kk(kl)) = -lam_tmp(kl);
0258 info.zu(kk(ku)) = lam_tmp(ku);
0259
0260 info.lambda(2*nb+2*nl2+nA-nk+(1:nk)) = [];
0261 nA = nA - nk;
0262 end
0263
0264
0265
0266 bus(:, MU_VMAX) = info.zu(vv.i1.Vm:vv.iN.Vm);
0267 bus(:, MU_VMIN) = info.zl(vv.i1.Vm:vv.iN.Vm);
0268 gen(:, MU_PMAX) = info.zu(vv.i1.Pg:vv.iN.Pg) / baseMVA;
0269 gen(:, MU_PMIN) = info.zl(vv.i1.Pg:vv.iN.Pg) / baseMVA;
0270 gen(:, MU_QMAX) = info.zu(vv.i1.Qg:vv.iN.Qg) / baseMVA;
0271 gen(:, MU_QMIN) = info.zl(vv.i1.Qg:vv.iN.Qg) / baseMVA;
0272 bus(:, LAM_P) = info.lambda(nn.i1.Pmis:nn.iN.Pmis) / baseMVA;
0273 bus(:, LAM_Q) = info.lambda(nn.i1.Qmis:nn.iN.Qmis) / baseMVA;
0274 branch(:, MU_SF) = muSf / baseMVA;
0275 branch(:, MU_ST) = muSt / baseMVA;
0276
0277
0278 nlnN = getN(om, 'nln');
0279
0280
0281 kl = find(info.lambda(1:2*nb) < 0);
0282 ku = find(info.lambda(1:2*nb) > 0);
0283 nl_mu_l = zeros(nlnN, 1);
0284 nl_mu_u = [zeros(2*nb, 1); muSf; muSt];
0285 nl_mu_l(kl) = -info.lambda(kl);
0286 nl_mu_u(ku) = info.lambda(ku);
0287
0288
0289 lam_lin = info.lambda(2*nb+2*nl2+(1:nA));
0290 kl = find(lam_lin < 0);
0291 ku = find(lam_lin > 0);
0292 mu_l = zeros(nA, 1);
0293 mu_l(kl) = -lam_lin(kl);
0294 mu_u = zeros(nA, 1);
0295 mu_u(ku) = lam_lin(ku);
0296
0297 mu = struct( ...
0298 'var', struct('l', info.zl, 'u', info.zu), ...
0299 'nln', struct('l', nl_mu_l, 'u', nl_mu_u), ...
0300 'lin', struct('l', mu_l, 'u', mu_u) );
0301
0302 results = mpc;
0303 [results.bus, results.branch, results.gen, ...
0304 results.om, results.x, results.mu, results.f] = ...
0305 deal(bus, branch, gen, om, x, mu, f);
0306
0307 pimul = [ ...
0308 results.mu.nln.l - results.mu.nln.u;
0309 results.mu.lin.l - results.mu.lin.u;
0310 -ones(ny>0, 1);
0311 results.mu.var.l - results.mu.var.u;
0312 ];
0313 raw = struct('xr', x, 'pimul', pimul, 'info', info.status, 'output', output);
0314
0315
0316
0317 function f = objective(x, d)
0318 f = opf_costfcn(x, d.om);
0319
0320 function df = gradient(x, d)
0321 [f, df] = opf_costfcn(x, d.om);
0322
0323 function c = constraints(x, d)
0324 [hn, gn] = opf_consfcn(x, d.om, d.Ybus, d.Yf, d.Yt, d.mpopt, d.il);
0325 if isempty(d.A)
0326 c = [gn; hn];
0327 else
0328 c = [gn; hn; d.A*x];
0329 end
0330
0331 function J = jacobian(x, d)
0332 [hn, gn, dhn, dgn] = opf_consfcn(x, d.om, d.Ybus, d.Yf, d.Yt, d.mpopt, d.il);
0333 J = [dgn'; dhn'; d.A];
0334
0335 function H = hessian(x, sigma, lambda, d)
0336 lam.eqnonlin = lambda(1:d.neqnln);
0337 lam.ineqnonlin = lambda(d.neqnln+(1:d.niqnln));
0338 H = tril(opf_hessfcn(x, lam, sigma, d.om, d.Ybus, d.Yf, d.Yt, d.mpopt, d.il));
0339
0340
0341
0342
0343
0344