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