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