0001 function [results, success, raw] = ktropf_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 use_ktropts_file = 1;
0062 create_ktropts_file = 0;
0063
0064
0065 mpc = get_mpc(om);
0066 [baseMVA, bus, gen, branch, gencost] = ...
0067 deal(mpc.baseMVA, mpc.bus, mpc.gen, mpc.branch, mpc.gencost);
0068 [vv, ll, nn] = get_idx(om);
0069
0070
0071 nb = size(bus, 1);
0072 ng = size(gen, 1);
0073 nl = size(branch, 1);
0074 ny = getN(om, 'var', 'y');
0075
0076
0077 [x0, xmin, xmax] = getv(om);
0078
0079
0080 [A, l, u] = linear_constraints(om);
0081
0082
0083
0084 ieq = find( abs(u-l) <= eps );
0085 igt = find( u >= 1e10 & l > -1e10 );
0086 ilt = find( l <= -1e10 & u < 1e10 );
0087 ibx = find( (abs(u-l) > eps) & (u < 1e10) & (l > -1e10) );
0088 Af = [ A(ilt, :); -A(igt, :); A(ibx, :); -A(ibx, :) ];
0089 bf = [ u(ilt); -l(igt); u(ibx); -l(ibx)];
0090 Afeq = A(ieq, :);
0091 bfeq = u(ieq);
0092
0093
0094 [Ybus, Yf, Yt] = makeYbus(baseMVA, bus, branch);
0095
0096
0097 if mpopt.opf.init_from_mpc ~= 1
0098 ll = xmin; uu = xmax;
0099 ll(xmin == -Inf) = -1e10;
0100 uu(xmax == Inf) = 1e10;
0101 x0 = (ll + uu) / 2;
0102 k = find(xmin == -Inf & xmax < Inf);
0103 x0(k) = xmax(k) - 1;
0104 k = find(xmin > -Inf & xmax == Inf);
0105 x0(k) = xmin(k) + 1;
0106 Varefs = bus(bus(:, BUS_TYPE) == REF, VA) * (pi/180);
0107 x0(vv.i1.Va:vv.iN.Va) = Varefs(1);
0108 if ny > 0
0109 ipwl = find(gencost(:, MODEL) == PW_LINEAR);
0110
0111
0112 c = gencost(sub2ind(size(gencost), ipwl, NCOST+2*gencost(ipwl, NCOST)));
0113 x0(vv.i1.y:vv.iN.y) = max(c) + 0.1 * abs(max(c));
0114
0115 end
0116 end
0117
0118
0119 il = find(branch(:, RATE_A) ~= 0 & branch(:, RATE_A) < 1e10);
0120 nl2 = length(il);
0121
0122
0123 nA = size(A, 1);
0124 nx = length(x0);
0125 f = branch(:, F_BUS);
0126 t = branch(:, T_BUS);
0127 Cf = sparse(1:nl, f, ones(nl, 1), nl, nb);
0128 Ct = sparse(1:nl, t, ones(nl, 1), nl, nb);
0129 Cl = Cf + Ct;
0130 Cb = Cl' * Cl + speye(nb);
0131 Cl2 = Cl(il, :);
0132 Cg = sparse(gen(:, GEN_BUS), (1:ng)', 1, nb, ng);
0133 nz = nx - 2*(nb+ng);
0134 nxtra = nx - 2*nb;
0135 Js = [
0136 Cl2 Cl2 sparse(nl2, 2*ng) sparse(nl2,nz);
0137 Cl2 Cl2 sparse(nl2, 2*ng) sparse(nl2,nz);
0138 Cb Cb Cg sparse(nb,ng) sparse(nb,nz);
0139 Cb Cb sparse(nb,ng) Cg sparse(nb,nz);
0140 ];
0141 [f, df, d2f] = opf_costfcn(x0, om);
0142 Hs = d2f + [
0143 Cb Cb sparse(nb,nxtra);
0144 Cb Cb sparse(nb,nxtra);
0145 sparse(nxtra,nx);
0146 ];
0147
0148
0149 hess_fcn = @(x, lambda)opf_hessfcn(x, lambda, 1, om, Ybus, Yf(il,:), Yt(il,:), mpopt, il);
0150 fmoptions = optimset('GradObj', 'on', 'GradConstr', 'on', ...
0151 'Hessian', 'user-supplied', 'HessFcn', hess_fcn, ...
0152 'JacobPattern', Js, 'HessPattern', Hs );
0153 if use_ktropts_file
0154 if ~isempty(mpopt.knitro.opt_fname)
0155 opt_fname = mpopt.knitro.opt_fname;
0156 elseif mpopt.knitro.opt
0157 opt_fname = sprintf('knitro_user_options_%d.txt', mpopt.knitro.opt);
0158 else
0159
0160 ktropts.algorithm = 1;
0161 ktropts.outlev = mpopt.verbose;
0162 ktropts.feastol = mpopt.opf.violation;
0163 ktropts.xtol = mpopt.knitro.tol_x;
0164 ktropts.opttol = mpopt.knitro.tol_f;
0165 if mpopt.fmincon.max_it ~= 0
0166 ktropts.maxit = mpopt.fmincon.max_it;
0167 end
0168 ktropts.bar_directinterval = 0;
0169 opt_fname = write_ktropts(ktropts);
0170 create_ktropts_file = 1;
0171 end
0172 else
0173 fmoptions = optimset(fmoptions, 'Algorithm', 'interior-point', ...
0174 'TolCon', mpopt.opf.violation, 'TolX', mpopt.knitro.tol_x, ...
0175 'TolFun', mpopt.knitro.tol_f );
0176 if mpopt.fmincon.max_it ~= 0
0177 fmoptions = optimset(fmoptions, 'MaxIter', mpopt.fmincon.max_it, ...
0178 'MaxFunEvals', 4 * mpopt.fmincon.max_it);
0179 end
0180 if mpopt.verbose == 0,
0181 fmoptions.Display = 'off';
0182 elseif mpopt.verbose == 1
0183 fmoptions.Display = 'iter';
0184 else
0185 fmoptions.Display = 'testing';
0186 end
0187 opt_fname = [];
0188 end
0189
0190
0191 f_fcn = @(x)opf_costfcn(x, om);
0192 gh_fcn = @(x)opf_consfcn(x, om, Ybus, Yf(il,:), Yt(il,:), mpopt, il);
0193 if have_fcn('knitromatlab')
0194 [x, f, info, Output, Lambda] = knitromatlab(f_fcn, x0, Af, bf, Afeq, bfeq, ...
0195 xmin, xmax, gh_fcn, [], fmoptions, opt_fname);
0196 else
0197 [x, f, info, Output, Lambda] = ktrlink(f_fcn, x0, Af, bf, Afeq, bfeq, ...
0198 xmin, xmax, gh_fcn, fmoptions, opt_fname);
0199 end
0200 success = (info == 0);
0201
0202
0203 if create_ktropts_file
0204 delete(opt_fname);
0205 end
0206
0207
0208 Va = x(vv.i1.Va:vv.iN.Va);
0209 Vm = x(vv.i1.Vm:vv.iN.Vm);
0210 Pg = x(vv.i1.Pg:vv.iN.Pg);
0211 Qg = x(vv.i1.Qg:vv.iN.Qg);
0212 V = Vm .* exp(1j*Va);
0213
0214
0215
0216 bus(:, VA) = Va * 180/pi;
0217 bus(:, VM) = Vm;
0218 gen(:, PG) = Pg * baseMVA;
0219 gen(:, QG) = Qg * baseMVA;
0220 gen(:, VG) = Vm(gen(:, GEN_BUS));
0221
0222
0223 Sf = V(branch(:, F_BUS)) .* conj(Yf * V);
0224 St = V(branch(:, T_BUS)) .* conj(Yt * V);
0225 branch(:, PF) = real(Sf) * baseMVA;
0226 branch(:, QF) = imag(Sf) * baseMVA;
0227 branch(:, PT) = real(St) * baseMVA;
0228 branch(:, QT) = imag(St) * baseMVA;
0229
0230
0231
0232 muSf = zeros(nl, 1);
0233 muSt = zeros(nl, 1);
0234 if ~isempty(il)
0235 muSf(il) = 2 * Lambda.ineqnonlin(1:nl2) .* branch(il, RATE_A) / baseMVA;
0236 muSt(il) = 2 * Lambda.ineqnonlin((1:nl2)+nl2) .* branch(il, RATE_A) / baseMVA;
0237 end
0238
0239
0240 bus(:, MU_VMAX) = Lambda.upper(vv.i1.Vm:vv.iN.Vm);
0241 bus(:, MU_VMIN) = -Lambda.lower(vv.i1.Vm:vv.iN.Vm);
0242 gen(:, MU_PMAX) = Lambda.upper(vv.i1.Pg:vv.iN.Pg) / baseMVA;
0243 gen(:, MU_PMIN) = -Lambda.lower(vv.i1.Pg:vv.iN.Pg) / baseMVA;
0244 gen(:, MU_QMAX) = Lambda.upper(vv.i1.Qg:vv.iN.Qg) / baseMVA;
0245 gen(:, MU_QMIN) = -Lambda.lower(vv.i1.Qg:vv.iN.Qg) / baseMVA;
0246 bus(:, LAM_P) = Lambda.eqnonlin(nn.i1.Pmis:nn.iN.Pmis) / baseMVA;
0247 bus(:, LAM_Q) = Lambda.eqnonlin(nn.i1.Qmis:nn.iN.Qmis) / baseMVA;
0248 branch(:, MU_SF) = muSf / baseMVA;
0249 branch(:, MU_ST) = muSt / baseMVA;
0250
0251
0252 nlnN = getN(om, 'nln');
0253 nlt = length(ilt);
0254 ngt = length(igt);
0255 nbx = length(ibx);
0256
0257
0258 kl = find(Lambda.eqnonlin < 0);
0259 ku = find(Lambda.eqnonlin > 0);
0260 nl_mu_l = zeros(nlnN, 1);
0261 nl_mu_u = [zeros(2*nb, 1); muSf; muSt];
0262 nl_mu_l(kl) = -Lambda.eqnonlin(kl);
0263 nl_mu_u(ku) = Lambda.eqnonlin(ku);
0264
0265
0266 kl = find(Lambda.eqlin < 0);
0267 ku = find(Lambda.eqlin > 0);
0268
0269 mu_l = zeros(size(u));
0270 mu_l(ieq(kl)) = -Lambda.eqlin(kl);
0271 mu_l(igt) = Lambda.ineqlin(nlt+(1:ngt));
0272 mu_l(ibx) = Lambda.ineqlin(nlt+ngt+nbx+(1:nbx));
0273
0274 mu_u = zeros(size(u));
0275 mu_u(ieq(ku)) = Lambda.eqlin(ku);
0276 mu_u(ilt) = Lambda.ineqlin(1:nlt);
0277 mu_u(ibx) = Lambda.ineqlin(nlt+ngt+(1:nbx));
0278
0279 mu = struct( ...
0280 'var', struct('l', -Lambda.lower, 'u', Lambda.upper), ...
0281 'nln', struct('l', nl_mu_l, 'u', nl_mu_u), ...
0282 'lin', struct('l', mu_l, 'u', mu_u) );
0283
0284 results = mpc;
0285 [results.bus, results.branch, results.gen, ...
0286 results.om, results.x, results.mu, results.f] = ...
0287 deal(bus, branch, gen, om, x, mu, f);
0288
0289 pimul = [ ...
0290 results.mu.nln.l - results.mu.nln.u;
0291 results.mu.lin.l - results.mu.lin.u;
0292 -ones(ny>0, 1);
0293 results.mu.var.l - results.mu.var.u;
0294 ];
0295 raw = struct('xr', x, 'pimul', pimul, 'info', info, 'output', Output);
0296
0297
0298
0299 function fname = write_ktropts(ktropts)
0300
0301
0302 fname = sprintf('ktropts_%06d.txt', fix(1e6*rand));
0303
0304
0305 [fd, msg] = fopen(fname, 'wt');
0306 if fd == -1
0307 error('could not create %d : %s', fname, msg);
0308 end
0309
0310
0311 fields = fieldnames(ktropts);
0312 for k = 1:length(fields)
0313 fprintf(fd, '%s %g\n', fields{k}, ktropts.(fields{k}));
0314 end
0315
0316
0317 if fd ~= 1
0318 fclose(fd);
0319 end