Description of ktropf

0001 function [results, success, raw] = ktropf_solver(om, mpopt)
0002 %KTROPF_SOLVER  Solves AC optimal power flow using Artelys Knitro.
0003 %
0004 %   [RESULTS, SUCCESS, RAW] = KTROPF_SOLVER(OM, MPOPT)
0005 %
0006 %   Inputs are an OPF model object and a MATPOWER options struct.
0007 %
0008 %   Outputs are a RESULTS struct, SUCCESS flag and RAW output struct.
0009 %
0010 %   RESULTS is a MATPOWER case struct (mpc) with the usual baseMVA, bus
0011 %   branch, gen, gencost fields, along with the following additional
0012 %   fields:
0013 %       .order      see 'help ext2int' for details of this field
0014 %       .x          final value of optimization variables (internal order)
0015 %       .f          final objective function value
0016 %       .mu         shadow prices on ...
0017 %           .var
0018 %               .l  lower bounds on variables
0019 %               .u  upper bounds on variables
0020 %           .nln    (deprecated) 2*nb+2*nl - Pmis, Qmis, Sf, St
0021 %               .l  lower bounds on nonlinear constraints
0022 %               .u  upper bounds on nonlinear constraints
0023 %           .nle    nonlinear equality constraints
0024 %           .nli    nonlinear inequality constraints
0025 %           .lin
0026 %               .l  lower bounds on linear constraints
0027 %               .u  upper bounds on linear constraints
0028 %
0029 %   SUCCESS     1 if solver converged successfully, 0 otherwise
0030 %
0031 %   RAW         raw output in form returned by MINOS
0032 %       .xr     final value of optimization variables
0033 %       .pimul  constraint multipliers
0034 %       .info   solver specific termination code
0035 %       .output solver specific output information
0036 %
0037 %   See also OPF, KTRLINK.
0038 
0039 %   MATPOWER
0040 %   Copyright (c) 2000-2017, Power Systems Engineering Research Center (PSERC)
0041 %   by Ray Zimmerman, PSERC Cornell
0042 %   and Carlos E. Murillo-Sanchez, PSERC Cornell & Universidad Nacional de Colombia
0043 %
0044 %   This file is part of MATPOWER.
0045 %   Covered by the 3-clause BSD License (see LICENSE file for details).
0046 %   See https://matpower.org for more info.
0047 
0048 %%----- initialization -----
0049 %% define named indices into data matrices
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 %% options
0061 use_ktropts_file = 1;       %% use a Knitro options file to pass options
0062 create_ktropts_file = 0;    %% generate a Knitro options file on the fly
0063 
0064 %% unpack data
0065 mpc = om.get_mpc();
0066 [baseMVA, bus, gen, branch, gencost] = ...
0067     deal(mpc.baseMVA, mpc.bus, mpc.gen, mpc.branch, mpc.gencost);
0068 [vv, ll, nne, nni] = om.get_idx();
0069 
0070 %% problem dimensions
0071 nb = size(bus, 1);          %% number of buses
0072 ng = size(gen, 1);          %% number of gens
0073 nl = size(branch, 1);       %% number of branches
0074 ny = om.getN('var', 'y');   %% number of piece-wise linear costs
0075 
0076 %% linear constraints
0077 [A, l, u] = om.params_lin_constraint();
0078 
0079 %% bounds on optimization vars
0080 [x0, xmin, xmax] = om.params_var();
0081 
0082 %% split l <= A*x <= u into less than, equal to, greater than, and
0083 %% doubly-bounded sets
0084 ieq = find( abs(u-l) <= eps );          %% equality
0085 igt = find( u >=  1e10 & l > -1e10 );   %% greater than, unbounded above
0086 ilt = find( l <= -1e10 & u <  1e10 );   %% less than, unbounded below
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 %% build admittance matrices
0094 [Ybus, Yf, Yt] = makeYbus(baseMVA, bus, branch);
0095 
0096 %% try to select an interior initial point, unless requested not to
0097 if mpopt.opf.start < 2
0098     s = 1;                      %% set init point inside bounds by s
0099     lb = xmin; ub = xmax;
0100     lb(xmin == -Inf) = -1e10;   %% replace Inf with numerical proxies
0101     ub(xmax ==  Inf) =  1e10;
0102     x0 = (lb + ub) / 2;         %% set x0 mid-way between bounds
0103     k = find(xmin == -Inf & xmax < Inf);    %% if only bounded above
0104     x0(k) = xmax(k) - s;                    %% set just below upper bound
0105     k = find(xmin > -Inf & xmax == Inf);    %% if only bounded below
0106     x0(k) = xmin(k) + s;                    %% set just above lower bound
0107     Varefs = bus(bus(:, BUS_TYPE) == REF, VA) * (pi/180);
0108     Vmax = min(bus(:, VMAX), 1.5);
0109     Vmin = max(bus(:, VMIN), 0.5);
0110     Vm = (Vmax + Vmin) / 2;
0111     if mpopt.opf.v_cartesian
0112         V = Vm * exp(1j*Varefs(1));
0113         x0(vv.i1.Vr:vv.iN.Vr) = real(V);
0114         x0(vv.i1.Vi:vv.iN.Vi) = imag(V);
0115     else
0116         x0(vv.i1.Va:vv.iN.Va) = Varefs(1);  %% angles set to first reference angle
0117         x0(vv.i1.Vm:vv.iN.Vm) = Vm;         %% voltage magnitudes
0118         if ny > 0
0119             ipwl = find(gencost(:, MODEL) == PW_LINEAR);
0120             c = gencost(sub2ind(size(gencost), ipwl, NCOST+2*gencost(ipwl, NCOST)));    %% largest y-value in CCV data
0121             x0(vv.i1.y:vv.iN.y) = max(c) + 0.1 * abs(max(c));
0122         end
0123     end
0124 end
0125 
0126 %% find branches with flow limits
0127 il = find(branch(:, RATE_A) ~= 0 & branch(:, RATE_A) < 1e10);
0128 nl2 = length(il);           %% number of constrained lines
0129 
0130 %% build Jacobian and Hessian structure
0131 % nA = size(A, 1);                %% number of original linear constraints
0132 % nx = length(x0);
0133 % f = branch(:, F_BUS);                           %% list of "from" buses
0134 % t = branch(:, T_BUS);                           %% list of "to" buses
0135 % Cf = sparse(1:nl, f, ones(nl, 1), nl, nb);      %% connection matrix for line & from buses
0136 % Ct = sparse(1:nl, t, ones(nl, 1), nl, nb);      %% connection matrix for line & to buses
0137 % Cl = Cf + Ct;
0138 % Cb = Cl' * Cl + speye(nb);
0139 % Cl2 = Cl(il, :);
0140 % Cg = sparse(gen(:, GEN_BUS), (1:ng)', 1, nb, ng);
0141 % nz = nx - 2*(nb+ng);
0142 % nxtra = nx - 2*nb;
0143 % [h, g, dh, dg] = opf_consfcn(x0, om, Ybus, Yf(il,:), Yt(il,:), mpopt, il);
0144 % Js = [dh'; dg'];
0145 % Js = [
0146 %     Cl2     Cl2     sparse(nl2, 2*ng)               sparse(nl2,nz);
0147 %     Cl2     Cl2     sparse(nl2, 2*ng)               sparse(nl2,nz);
0148 %     Cb      Cb      Cg              sparse(nb,ng)   sparse(nb,nz);
0149 %     Cb      Cb      sparse(nb,ng)   Cg              sparse(nb,nz);
0150 % ];
0151 % lam = struct('eqnonlin', ones(size(dg,2),1), 'ineqnonlin', ones(size(dh,2),1) );
0152 % Hs = opf_hessfcn(x0, lam, 1, om, Ybus, Yf(il,:), Yt(il,:), mpopt, il);
0153 % [f, df, d2f] = opf_costfcn(x0, om);
0154 % Hs = d2f + [
0155 %     Cb  Cb  sparse(nb,nxtra);
0156 %     Cb  Cb  sparse(nb,nxtra);
0157 %     sparse(nxtra,nx);
0158 % ];
0159 randx = rand(size(x0));
0160 [h, g, dh, dg] = opf_consfcn(randx, om, Ybus, Yf(il,:), Yt(il,:), mpopt, il);
0161 Js = [dh'; dg'];
0162 lam = struct('eqnonlin', rand(size(dg,2),1), 'ineqnonlin', rand(size(dh,2),1) );
0163 Hs = opf_hessfcn(randx, lam, 1, om, Ybus, Yf(il,:), Yt(il,:), mpopt, il);
0164 neq = length(g);
0165 niq = length(h);
0166 
0167 %% basic optimset options needed for ktrlink
0168 hess_fcn = @(x, lambda)opf_hessfcn(x, lambda, 1, om, Ybus, Yf(il,:), Yt(il,:), mpopt, il);
0169 fmoptions = optimset('GradObj', 'on', 'GradConstr', 'on', ...
0170                 'Hessian', 'user-supplied', 'HessFcn', hess_fcn, ...
0171                 'JacobPattern', Js, 'HessPattern', Hs );
0172 if use_ktropts_file
0173     if ~isempty(mpopt.knitro.opt_fname)
0174         opt_fname = mpopt.knitro.opt_fname;
0175     elseif mpopt.knitro.opt
0176         opt_fname = sprintf('knitro_user_options_%d.txt', mpopt.knitro.opt);
0177     else
0178         %% create ktropts file
0179         ktropts.algorithm           = 1;
0180         ktropts.outlev              = mpopt.verbose;
0181         ktropts.feastol             = mpopt.opf.violation;
0182         ktropts.xtol                = mpopt.knitro.tol_x;
0183         ktropts.opttol              = mpopt.knitro.tol_f;
0184         if mpopt.knitro.maxit ~= 0
0185             ktropts.maxit           = mpopt.knitro.maxit;
0186         end
0187         ktropts.bar_directinterval  = 0;
0188         opt_fname = write_ktropts(ktropts);
0189         create_ktropts_file = 1;    %% make a note that I created it
0190     end
0191 else
0192     fmoptions = optimset(fmoptions, 'Algorithm', 'interior-point', ...
0193         'TolCon', mpopt.opf.violation, 'TolX', mpopt.knitro.tol_x, ...
0194         'TolFun', mpopt.knitro.tol_f );
0195     if mpopt.fmincon.max_it ~= 0
0196         fmoptions = optimset(fmoptions, 'MaxIter', mpopt.fmincon.max_it, ...
0197                     'MaxFunEvals', 4 * mpopt.fmincon.max_it);
0198     end
0199     if mpopt.verbose == 0,
0200       fmoptions.Display = 'off';
0201     elseif mpopt.verbose == 1
0202       fmoptions.Display = 'iter';
0203     else
0204       fmoptions.Display = 'testing';
0205     end
0206     opt_fname = [];
0207 end
0208 
0209 %%-----  run opf  -----
0210 f_fcn = @(x)opf_costfcn(x, om);
0211 gh_fcn = @(x)opf_consfcn(x, om, Ybus, Yf(il,:), Yt(il,:), mpopt, il);
0212 if have_fcn('knitromatlab')
0213     [x, f, info, Output, Lambda] = knitromatlab(f_fcn, x0, Af, bf, Afeq, bfeq, ...
0214                                     xmin, xmax, gh_fcn, [], fmoptions, opt_fname);
0215 else
0216     [x, f, info, Output, Lambda] = ktrlink(f_fcn, x0, Af, bf, Afeq, bfeq, ...
0217                                     xmin, xmax, gh_fcn, fmoptions, opt_fname);
0218 end
0219 success = (info == 0);
0220 
0221 %% delete ktropts file
0222 if create_ktropts_file  %% ... but only if I created it
0223     delete(opt_fname);
0224 end
0225 
0226 %% update solution data
0227 if mpopt.opf.v_cartesian
0228     Vi = x(vv.i1.Vi:vv.iN.Vi);
0229     Vr = x(vv.i1.Vr:vv.iN.Vr);
0230     V = Vr + 1j*Vi;
0231     Va = angle(V);
0232     Vm = abs(V);
0233 else
0234     Va = x(vv.i1.Va:vv.iN.Va);
0235     Vm = x(vv.i1.Vm:vv.iN.Vm);
0236     V = Vm .* exp(1j*Va);
0237 end
0238 Pg = x(vv.i1.Pg:vv.iN.Pg);
0239 Qg = x(vv.i1.Qg:vv.iN.Qg);
0240 
0241 %%-----  calculate return values  -----
0242 %% update voltages & generator outputs
0243 bus(:, VA) = Va * 180/pi;
0244 bus(:, VM) = Vm;
0245 gen(:, PG) = Pg * baseMVA;
0246 gen(:, QG) = Qg * baseMVA;
0247 gen(:, VG) = Vm(gen(:, GEN_BUS));
0248 
0249 %% compute branch flows
0250 Sf = V(branch(:, F_BUS)) .* conj(Yf * V);  %% cplx pwr at "from" bus, p.u.
0251 St = V(branch(:, T_BUS)) .* conj(Yt * V);  %% cplx pwr at "to" bus, p.u.
0252 branch(:, PF) = real(Sf) * baseMVA;
0253 branch(:, QF) = imag(Sf) * baseMVA;
0254 branch(:, PT) = real(St) * baseMVA;
0255 branch(:, QT) = imag(St) * baseMVA;
0256 
0257 %% line constraint is typically on square of limit
0258 %% so we must fix multipliers
0259 muSf = zeros(nl, 1);
0260 muSt = zeros(nl, 1);
0261 if ~isempty(il)
0262     if upper(mpopt.opf.flow_lim(1)) == 'P'
0263         muSf(il) = Lambda.ineqnonlin(nni.i1.Sf:nni.iN.Sf);
0264         muSt(il) = Lambda.ineqnonlin(nni.i1.St:nni.iN.St);
0265     else
0266         muSf(il) = 2 * Lambda.ineqnonlin(nni.i1.Sf:nni.iN.Sf) .* branch(il, RATE_A) / baseMVA;
0267         muSt(il) = 2 * Lambda.ineqnonlin(nni.i1.St:nni.iN.St) .* branch(il, RATE_A) / baseMVA;
0268     end
0269 end
0270 
0271 %% fix Lambdas
0272 %% (shadow prices on equality variable bounds can come back on wrong limit)
0273 kl = find(Lambda.lower > 0 & Lambda.upper == 0);
0274 Lambda.upper(kl) = Lambda.lower(kl);
0275 Lambda.lower(kl) = 0;
0276 ku = find(Lambda.upper < 0 & Lambda.lower == 0);
0277 Lambda.lower(ku) = Lambda.upper(ku);
0278 Lambda.upper(ku) = 0;
0279 
0280 %% update Lagrange multipliers
0281 if mpopt.opf.v_cartesian
0282     if om.userdata.veq
0283         lam = Lambda.eqnonlin(nne.i1.Veq:nne.iN.Veq);
0284         mu_Vmax = zeros(size(lam));
0285         mu_Vmin = zeros(size(lam));
0286         mu_Vmax(lam > 0) =  lam(lam > 0);
0287         mu_Vmin(lam < 0) = -lam(lam < 0);
0288         bus(om.userdata.veq, MU_VMAX) = mu_Vmax;
0289         bus(om.userdata.veq, MU_VMIN) = mu_Vmin;
0290     end
0291     bus(om.userdata.viq, MU_VMAX) = Lambda.ineqnonlin(nni.i1.Vmax:nni.iN.Vmax);
0292     bus(om.userdata.viq, MU_VMIN) = Lambda.ineqnonlin(nni.i1.Vmin:nni.iN.Vmin);
0293 else
0294     bus(:, MU_VMAX)  = Lambda.upper(vv.i1.Vm:vv.iN.Vm);
0295     bus(:, MU_VMIN)  = -Lambda.lower(vv.i1.Vm:vv.iN.Vm);
0296 end
0297 gen(:, MU_PMAX)  = Lambda.upper(vv.i1.Pg:vv.iN.Pg) / baseMVA;
0298 gen(:, MU_PMIN)  = -Lambda.lower(vv.i1.Pg:vv.iN.Pg) / baseMVA;
0299 gen(:, MU_QMAX)  = Lambda.upper(vv.i1.Qg:vv.iN.Qg) / baseMVA;
0300 gen(:, MU_QMIN)  = -Lambda.lower(vv.i1.Qg:vv.iN.Qg) / baseMVA;
0301 if mpopt.opf.current_balance
0302     %% convert current balance shadow prices to equivalent lamP and lamQ
0303     %% P + jQ = (Vr + jVi) * (M - jN)
0304     %% M = (Vr P + Vi Q) / (Vr^2 + Vi^2)
0305     %% N = (Vi P - Vr Q) / (Vr^2 + Vi^2)
0306     %% lamP = df/dP = df/dM * dM/dP + df/dN + dN/dP
0307     %% lamQ = df/dQ = df/dM * dM/dQ + df/dN + dN/dQ
0308     VV = V ./ (V .* conj(V));   %% V / Vm^2
0309     VVr = real(VV);
0310     VVi = imag(VV);
0311     lamM = Lambda.eqnonlin(nne.i1.rImis:nne.iN.rImis);
0312     lamN = Lambda.eqnonlin(nne.i1.iImis:nne.iN.iImis);
0313     bus(:, LAM_P) = (VVr.*lamM + VVi.*lamN) / baseMVA;
0314     bus(:, LAM_Q) = (VVi.*lamM - VVr.*lamN) / baseMVA;
0315 else
0316     bus(:, LAM_P) = Lambda.eqnonlin(nne.i1.Pmis:nne.iN.Pmis) / baseMVA;
0317     bus(:, LAM_Q) = Lambda.eqnonlin(nne.i1.Qmis:nne.iN.Qmis) / baseMVA;
0318 end
0319 branch(:, MU_SF) = muSf / baseMVA;
0320 branch(:, MU_ST) = muSt / baseMVA;
0321 
0322 %% package up results
0323 nlnN = 2*nb + 2*nl;     %% because muSf and muSt are nl x 1, not nl2 x 1
0324 nlt = length(ilt);
0325 ngt = length(igt);
0326 nbx = length(ibx);
0327 
0328 %% extract multipliers for nonlinear constraints
0329 kl = find(Lambda.eqnonlin < 0);
0330 ku = find(Lambda.eqnonlin > 0);
0331 nl_mu_l = zeros(nlnN, 1);
0332 nl_mu_u = [zeros(2*nb, 1); muSf; muSt];
0333 nl_mu_l(kl) = -Lambda.eqnonlin(kl);
0334 nl_mu_u(ku) =  Lambda.eqnonlin(ku);
0335 
0336 %% extract multipliers for linear constraints
0337 kl = find(Lambda.eqlin < 0);
0338 ku = find(Lambda.eqlin > 0);
0339 
0340 mu_l = zeros(size(u));
0341 mu_l(ieq(kl)) = -Lambda.eqlin(kl);
0342 mu_l(igt) = Lambda.ineqlin(nlt+(1:ngt));
0343 mu_l(ibx) = Lambda.ineqlin(nlt+ngt+nbx+(1:nbx));
0344 
0345 mu_u = zeros(size(u));
0346 mu_u(ieq(ku)) = Lambda.eqlin(ku);
0347 mu_u(ilt) = Lambda.ineqlin(1:nlt);
0348 mu_u(ibx) = Lambda.ineqlin(nlt+ngt+(1:nbx));
0349 
0350 mu = struct( ...
0351   'var', struct('l', -Lambda.lower, 'u', Lambda.upper), ...
0352   'nln', struct('l', nl_mu_l, 'u', nl_mu_u), ...
0353   'nle', Lambda.eqnonlin, ...
0354   'nli', Lambda.ineqnonlin, ...
0355   'lin', struct('l', mu_l, 'u', mu_u) );
0356 
0357 results = mpc;
0358 [results.bus, results.branch, results.gen, ...
0359     results.om, results.x, results.mu, results.f] = ...
0360         deal(bus, branch, gen, om, x, mu, f);
0361 
0362 pimul = [ ...
0363   results.mu.nln.l - results.mu.nln.u;
0364   results.mu.lin.l - results.mu.lin.u;
0365   -ones(ny>0, 1);
0366   results.mu.var.l - results.mu.var.u;
0367 ];
0368 raw = struct('xr', x, 'pimul', pimul, 'info', info, 'output', Output);
0369 
0370 
0371 %%-----  write_ktropts  -----
0372 function fname = write_ktropts(ktropts)
0373 
0374 %% generate file name
0375 fname = sprintf('ktropts_%06d.txt', fix(1e6*rand));
0376 
0377 %% open file
0378 [fd, msg] = fopen(fname, 'wt');     %% write options file
0379 if fd == -1
0380     error('could not create %d : %s', fname, msg);
0381 end
0382 
0383 %% write options
0384 fields = fieldnames(ktropts);
0385 for k = 1:length(fields)
0386     fprintf(fd, '%s %g\n', fields{k}, ktropts.(fields{k}));
0387 end
0388 
0389 %% close file
0390 if fd ~= 1
0391     fclose(fd);
0392 end
ktropf_solver

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SUBFUNCTIONS

SOURCE CODE
