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