Home > matpower5.1 > ktropf_solver.m

ktropf_solver

PURPOSE ^

KTROPF_SOLVER Solves AC optimal power flow using KNITRO.

SYNOPSIS ^

function [results, success, raw] = ktropf_solver(om, mpopt)

DESCRIPTION ^

KTROPF_SOLVER  Solves AC optimal power flow using KNITRO.

   [RESULTS, SUCCESS, RAW] = KTROPF_SOLVER(OM, MPOPT)

   Inputs are an OPF model object and a MATPOWER options struct.

   Outputs are a RESULTS struct, SUCCESS flag and RAW output struct.

   RESULTS is a MATPOWER case struct (mpc) with the usual baseMVA, bus
   branch, gen, gencost fields, along with the following additional
   fields:
       .order      see 'help ext2int' for details of this field
       .x          final value of optimization variables (internal order)
       .f          final objective function value
       .mu         shadow prices on ...
           .var
               .l  lower bounds on variables
               .u  upper bounds on variables
           .nln
               .l  lower bounds on nonlinear constraints
               .u  upper bounds on nonlinear constraints
           .lin
               .l  lower bounds on linear constraints
               .u  upper bounds on linear constraints

   SUCCESS     1 if solver converged successfully, 0 otherwise

   RAW         raw output in form returned by MINOS
       .xr     final value of optimization variables
       .pimul  constraint multipliers
       .info   solver specific termination code
       .output solver specific output information

   See also OPF, KTRLINK.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 function [results, success, raw] = ktropf_solver(om, mpopt)
0002 %KTROPF_SOLVER  Solves AC optimal power flow using 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
0021 %               .l  lower bounds on nonlinear constraints
0022 %               .u  upper bounds on nonlinear constraints
0023 %           .lin
0024 %               .l  lower bounds on linear constraints
0025 %               .u  upper bounds on linear constraints
0026 %
0027 %   SUCCESS     1 if solver converged successfully, 0 otherwise
0028 %
0029 %   RAW         raw output in form returned by MINOS
0030 %       .xr     final value of optimization variables
0031 %       .pimul  constraint multipliers
0032 %       .info   solver specific termination code
0033 %       .output solver specific output information
0034 %
0035 %   See also OPF, KTRLINK.
0036 
0037 %   MATPOWER
0038 %   Copyright (c) 2000-2015 by Power System Engineering Research Center (PSERC)
0039 %   by Ray Zimmerman, PSERC Cornell
0040 %   and Carlos E. Murillo-Sanchez, PSERC Cornell & Universidad Autonoma de Manizales
0041 %
0042 %   $Id: ktropf_solver.m 2644 2015-03-11 19:34:22Z ray $
0043 %
0044 %   This file is part of MATPOWER.
0045 %   Covered by the 3-clause BSD License (see LICENSE file for details).
0046 %   See http://www.pserc.cornell.edu/matpower/ 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 = 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 %% 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 = getN(om, 'var', 'y');  %% number of piece-wise linear costs
0075 
0076 %% bounds on optimization vars
0077 [x0, xmin, xmax] = getv(om);
0078 
0079 %% linear constraints
0080 [A, l, u] = linear_constraints(om);
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
0097 if mpopt.opf.init_from_mpc ~= 1
0098     ll = xmin; uu = xmax;
0099     ll(xmin == -Inf) = -1e10;   %% replace Inf with numerical proxies
0100     uu(xmax ==  Inf) =  1e10;
0101     x0 = (ll + uu) / 2;         %% set x0 mid-way between bounds
0102     k = find(xmin == -Inf & xmax < Inf);    %% if only bounded above
0103     x0(k) = xmax(k) - 1;                    %% set just below upper bound
0104     k = find(xmin > -Inf & xmax == Inf);    %% if only bounded below
0105     x0(k) = xmin(k) + 1;                    %% set just above lower bound
0106     Varefs = bus(bus(:, BUS_TYPE) == REF, VA) * (pi/180);
0107     x0(vv.i1.Va:vv.iN.Va) = Varefs(1);  %% angles set to first reference angle
0108     if ny > 0
0109         ipwl = find(gencost(:, MODEL) == PW_LINEAR);
0110     %     PQ = [gen(:, PMAX); gen(:, QMAX)];
0111     %     c = totcost(gencost(ipwl, :), PQ(ipwl));
0112         c = gencost(sub2ind(size(gencost), ipwl, NCOST+2*gencost(ipwl, NCOST)));    %% largest y-value in CCV data
0113         x0(vv.i1.y:vv.iN.y) = max(c) + 0.1 * abs(max(c));
0114     %     x0(vv.i1.y:vv.iN.y) = c + 0.1 * abs(c);
0115     end
0116 end
0117 
0118 %% find branches with flow limits
0119 il = find(branch(:, RATE_A) ~= 0 & branch(:, RATE_A) < 1e10);
0120 nl2 = length(il);           %% number of constrained lines
0121 
0122 %% build Jacobian and Hessian structure
0123 nA = size(A, 1);                %% number of original linear constraints
0124 nx = length(x0);
0125 f = branch(:, F_BUS);                           %% list of "from" buses
0126 t = branch(:, T_BUS);                           %% list of "to" buses
0127 Cf = sparse(1:nl, f, ones(nl, 1), nl, nb);      %% connection matrix for line & from buses
0128 Ct = sparse(1:nl, t, ones(nl, 1), nl, nb);      %% connection matrix for line & to buses
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 %% basic optimset options needed for ktrlink
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         %% create ktropts file
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;    %% make a note that I created it
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 %%-----  run opf  -----
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 %% delete ktropts file
0203 if create_ktropts_file  %% ... but only if I created it
0204     delete(opt_fname);
0205 end
0206 
0207 %% update solution data
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 %%-----  calculate return values  -----
0215 %% update voltages & generator outputs
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 %% compute branch flows
0223 Sf = V(branch(:, F_BUS)) .* conj(Yf * V);  %% cplx pwr at "from" bus, p.u.
0224 St = V(branch(:, T_BUS)) .* conj(Yt * V);  %% cplx pwr at "to" bus, p.u.
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 %% line constraint is actually on square of limit
0231 %% so we must fix multipliers
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 %% update Lagrange multipliers
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 %% package up results
0252 nlnN = getN(om, 'nln');
0253 nlt = length(ilt);
0254 ngt = length(igt);
0255 nbx = length(ibx);
0256 
0257 %% extract multipliers for nonlinear constraints
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 %% extract multipliers for linear constraints
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 %%-----  write_ktropts  -----
0299 function fname = write_ktropts(ktropts)
0300 
0301 %% generate file name
0302 fname = sprintf('ktropts_%06d.txt', fix(1e6*rand));
0303 
0304 %% open file
0305 [fd, msg] = fopen(fname, 'wt');     %% write options file
0306 if fd == -1
0307     error('could not create %d : %s', fname, msg);
0308 end
0309 
0310 %% write options
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 %% close file
0317 if fd ~= 1
0318     fclose(fd);
0319 end

Generated on Fri 20-Mar-2015 18:23:34 by m2html © 2005