Home > matpower6.0 > 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-2016, Power Systems Engineering Research Center (PSERC)
0039 %   by Ray Zimmerman, PSERC Cornell
0040 %   and Carlos E. Murillo-Sanchez, PSERC Cornell & Universidad Nacional de Colombia
0041 %
0042 %   This file is part of MATPOWER.
0043 %   Covered by the 3-clause BSD License (see LICENSE file for details).
0044 %   See http://www.pserc.cornell.edu/matpower/ for more info.
0045 
0046 %%----- initialization -----
0047 %% define named indices into data matrices
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 %% options
0059 use_ktropts_file = 1;       %% use a KNITRO options file to pass options
0060 create_ktropts_file = 0;    %% generate a KNITRO options file on the fly
0061 
0062 %% unpack data
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 %% problem dimensions
0069 nb = size(bus, 1);          %% number of buses
0070 ng = size(gen, 1);          %% number of gens
0071 nl = size(branch, 1);       %% number of branches
0072 ny = getN(om, 'var', 'y');  %% number of piece-wise linear costs
0073 
0074 %% bounds on optimization vars
0075 [x0, xmin, xmax] = getv(om);
0076 
0077 %% linear constraints
0078 [A, l, u] = linear_constraints(om);
0079 
0080 %% split l <= A*x <= u into less than, equal to, greater than, and
0081 %% doubly-bounded sets
0082 ieq = find( abs(u-l) <= eps );          %% equality
0083 igt = find( u >=  1e10 & l > -1e10 );   %% greater than, unbounded above
0084 ilt = find( l <= -1e10 & u <  1e10 );   %% less than, unbounded below
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 %% build admittance matrices
0092 [Ybus, Yf, Yt] = makeYbus(baseMVA, bus, branch);
0093 
0094 %% try to select an interior initial point
0095 if mpopt.opf.init_from_mpc ~= 1
0096     ll = xmin; uu = xmax;
0097     ll(xmin == -Inf) = -1e10;   %% replace Inf with numerical proxies
0098     uu(xmax ==  Inf) =  1e10;
0099     x0 = (ll + uu) / 2;         %% set x0 mid-way between bounds
0100     k = find(xmin == -Inf & xmax < Inf);    %% if only bounded above
0101     x0(k) = xmax(k) - 1;                    %% set just below upper bound
0102     k = find(xmin > -Inf & xmax == Inf);    %% if only bounded below
0103     x0(k) = xmin(k) + 1;                    %% set just above lower bound
0104     Varefs = bus(bus(:, BUS_TYPE) == REF, VA) * (pi/180);
0105     x0(vv.i1.Va:vv.iN.Va) = Varefs(1);  %% angles set to first reference angle
0106     if ny > 0
0107         ipwl = find(gencost(:, MODEL) == PW_LINEAR);
0108     %     PQ = [gen(:, PMAX); gen(:, QMAX)];
0109     %     c = totcost(gencost(ipwl, :), PQ(ipwl));
0110         c = gencost(sub2ind(size(gencost), ipwl, NCOST+2*gencost(ipwl, NCOST)));    %% largest y-value in CCV data
0111         x0(vv.i1.y:vv.iN.y) = max(c) + 0.1 * abs(max(c));
0112     %     x0(vv.i1.y:vv.iN.y) = c + 0.1 * abs(c);
0113     end
0114 end
0115 
0116 %% find branches with flow limits
0117 il = find(branch(:, RATE_A) ~= 0 & branch(:, RATE_A) < 1e10);
0118 nl2 = length(il);           %% number of constrained lines
0119 
0120 %% build Jacobian and Hessian structure
0121 nA = size(A, 1);                %% number of original linear constraints
0122 nx = length(x0);
0123 f = branch(:, F_BUS);                           %% list of "from" buses
0124 t = branch(:, T_BUS);                           %% list of "to" buses
0125 Cf = sparse(1:nl, f, ones(nl, 1), nl, nb);      %% connection matrix for line & from buses
0126 Ct = sparse(1:nl, t, ones(nl, 1), nl, nb);      %% connection matrix for line & to buses
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 %% basic optimset options needed for ktrlink
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         %% create ktropts file
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;    %% make a note that I created it
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 %%-----  run opf  -----
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 %% delete ktropts file
0201 if create_ktropts_file  %% ... but only if I created it
0202     delete(opt_fname);
0203 end
0204 
0205 %% update solution data
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 %%-----  calculate return values  -----
0213 %% update voltages & generator outputs
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 %% compute branch flows
0221 Sf = V(branch(:, F_BUS)) .* conj(Yf * V);  %% cplx pwr at "from" bus, p.u.
0222 St = V(branch(:, T_BUS)) .* conj(Yt * V);  %% cplx pwr at "to" bus, p.u.
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 %% line constraint is actually on square of limit
0229 %% so we must fix multipliers
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 %% fix Lambdas
0238 %% (shadow prices on equality variable bounds can come back on wrong limit)
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 %% update Lagrange multipliers
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 %% package up results
0259 nlnN = getN(om, 'nln');
0260 nlt = length(ilt);
0261 ngt = length(igt);
0262 nbx = length(ibx);
0263 
0264 %% extract multipliers for nonlinear constraints
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 %% extract multipliers for linear constraints
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 %%-----  write_ktropts  -----
0306 function fname = write_ktropts(ktropts)
0307 
0308 %% generate file name
0309 fname = sprintf('ktropts_%06d.txt', fix(1e6*rand));
0310 
0311 %% open file
0312 [fd, msg] = fopen(fname, 'wt');     %% write options file
0313 if fd == -1
0314     error('could not create %d : %s', fname, msg);
0315 end
0316 
0317 %% write options
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 %% close file
0324 if fd ~= 1
0325     fclose(fd);
0326 end

Generated on Fri 16-Dec-2016 12:45:37 by m2html © 2005