Home > matpower5.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 %   $Id: ktropf_solver.m 2438 2014-12-05 14:32:42Z ray $
0039 %   by Ray Zimmerman, PSERC Cornell
0040 %   and Carlos E. Murillo-Sanchez, PSERC Cornell & Universidad Autonoma de Manizales
0041 %   Copyright (c) 2000-2013 by Power System Engineering Research Center (PSERC)
0042 %
0043 %   This file is part of MATPOWER.
0044 %   See http://www.pserc.cornell.edu/matpower/ for more info.
0045 %
0046 %   MATPOWER is free software: you can redistribute it and/or modify
0047 %   it under the terms of the GNU General Public License as published
0048 %   by the Free Software Foundation, either version 3 of the License,
0049 %   or (at your option) any later version.
0050 %
0051 %   MATPOWER is distributed in the hope that it will be useful,
0052 %   but WITHOUT ANY WARRANTY; without even the implied warranty of
0053 %   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
0054 %   GNU General Public License for more details.
0055 %
0056 %   You should have received a copy of the GNU General Public License
0057 %   along with MATPOWER. If not, see <http://www.gnu.org/licenses/>.
0058 %
0059 %   Additional permission under GNU GPL version 3 section 7
0060 %
0061 %   If you modify MATPOWER, or any covered work, to interface with
0062 %   other modules (such as MATLAB code and MEX-files) available in a
0063 %   MATLAB(R) or comparable environment containing parts covered
0064 %   under other licensing terms, the licensors of MATPOWER grant
0065 %   you additional permission to convey the resulting work.
0066 
0067 %%----- initialization -----
0068 %% define named indices into data matrices
0069 [PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ...
0070     VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus;
0071 [GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, GEN_STATUS, PMAX, PMIN, ...
0072     MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN, PC1, PC2, QC1MIN, QC1MAX, ...
0073     QC2MIN, QC2MAX, RAMP_AGC, RAMP_10, RAMP_30, RAMP_Q, APF] = idx_gen;
0074 [F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, RATE_C, ...
0075     TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST, ...
0076     ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch;
0077 [PW_LINEAR, POLYNOMIAL, MODEL, STARTUP, SHUTDOWN, NCOST, COST] = idx_cost;
0078 
0079 %% options
0080 use_ktropts_file = 1;       %% use a KNITRO options file to pass options
0081 create_ktropts_file = 0;    %% generate a KNITRO options file on the fly
0082 
0083 %% unpack data
0084 mpc = get_mpc(om);
0085 [baseMVA, bus, gen, branch, gencost] = ...
0086     deal(mpc.baseMVA, mpc.bus, mpc.gen, mpc.branch, mpc.gencost);
0087 [vv, ll, nn] = get_idx(om);
0088 
0089 %% problem dimensions
0090 nb = size(bus, 1);          %% number of buses
0091 ng = size(gen, 1);          %% number of gens
0092 nl = size(branch, 1);       %% number of branches
0093 ny = getN(om, 'var', 'y');  %% number of piece-wise linear costs
0094 
0095 %% bounds on optimization vars
0096 [x0, xmin, xmax] = getv(om);
0097 
0098 %% linear constraints
0099 [A, l, u] = linear_constraints(om);
0100 
0101 %% split l <= A*x <= u into less than, equal to, greater than, and
0102 %% doubly-bounded sets
0103 ieq = find( abs(u-l) <= eps );          %% equality
0104 igt = find( u >=  1e10 & l > -1e10 );   %% greater than, unbounded above
0105 ilt = find( l <= -1e10 & u <  1e10 );   %% less than, unbounded below
0106 ibx = find( (abs(u-l) > eps) & (u < 1e10) & (l > -1e10) );
0107 Af  = [ A(ilt, :); -A(igt, :); A(ibx, :); -A(ibx, :) ];
0108 bf  = [ u(ilt);   -l(igt);     u(ibx);    -l(ibx)];
0109 Afeq = A(ieq, :);
0110 bfeq = u(ieq);
0111 
0112 %% build admittance matrices
0113 [Ybus, Yf, Yt] = makeYbus(baseMVA, bus, branch);
0114 
0115 %% try to select an interior initial point
0116 if mpopt.opf.init_from_mpc ~= 1
0117     ll = xmin; uu = xmax;
0118     ll(xmin == -Inf) = -1e10;   %% replace Inf with numerical proxies
0119     uu(xmax ==  Inf) =  1e10;
0120     x0 = (ll + uu) / 2;         %% set x0 mid-way between bounds
0121     k = find(xmin == -Inf & xmax < Inf);    %% if only bounded above
0122     x0(k) = xmax(k) - 1;                    %% set just below upper bound
0123     k = find(xmin > -Inf & xmax == Inf);    %% if only bounded below
0124     x0(k) = xmin(k) + 1;                    %% set just above lower bound
0125     Varefs = bus(bus(:, BUS_TYPE) == REF, VA) * (pi/180);
0126     x0(vv.i1.Va:vv.iN.Va) = Varefs(1);  %% angles set to first reference angle
0127     if ny > 0
0128         ipwl = find(gencost(:, MODEL) == PW_LINEAR);
0129     %     PQ = [gen(:, PMAX); gen(:, QMAX)];
0130     %     c = totcost(gencost(ipwl, :), PQ(ipwl));
0131         c = gencost(sub2ind(size(gencost), ipwl, NCOST+2*gencost(ipwl, NCOST)));    %% largest y-value in CCV data
0132         x0(vv.i1.y:vv.iN.y) = max(c) + 0.1 * abs(max(c));
0133     %     x0(vv.i1.y:vv.iN.y) = c + 0.1 * abs(c);
0134     end
0135 end
0136 
0137 %% find branches with flow limits
0138 il = find(branch(:, RATE_A) ~= 0 & branch(:, RATE_A) < 1e10);
0139 nl2 = length(il);           %% number of constrained lines
0140 
0141 %% build Jacobian and Hessian structure
0142 nA = size(A, 1);                %% number of original linear constraints
0143 nx = length(x0);
0144 f = branch(:, F_BUS);                           %% list of "from" buses
0145 t = branch(:, T_BUS);                           %% list of "to" buses
0146 Cf = sparse(1:nl, f, ones(nl, 1), nl, nb);      %% connection matrix for line & from buses
0147 Ct = sparse(1:nl, t, ones(nl, 1), nl, nb);      %% connection matrix for line & to buses
0148 Cl = Cf + Ct;
0149 Cb = Cl' * Cl + speye(nb);
0150 Cl2 = Cl(il, :);
0151 Cg = sparse(gen(:, GEN_BUS), (1:ng)', 1, nb, ng);
0152 nz = nx - 2*(nb+ng);
0153 nxtra = nx - 2*nb;
0154 Js = [
0155     Cl2     Cl2     sparse(nl2, 2*ng)               sparse(nl2,nz);
0156     Cl2     Cl2     sparse(nl2, 2*ng)               sparse(nl2,nz);
0157     Cb      Cb      Cg              sparse(nb,ng)   sparse(nb,nz);
0158     Cb      Cb      sparse(nb,ng)   Cg              sparse(nb,nz);
0159 ];
0160 [f, df, d2f] = opf_costfcn(x0, om);
0161 Hs = d2f + [
0162     Cb  Cb  sparse(nb,nxtra);
0163     Cb  Cb  sparse(nb,nxtra);
0164     sparse(nxtra,nx);
0165 ];
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.fmincon.max_it ~= 0
0185             ktropts.maxit           = mpopt.fmincon.max_it;
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 Va = x(vv.i1.Va:vv.iN.Va);
0228 Vm = x(vv.i1.Vm:vv.iN.Vm);
0229 Pg = x(vv.i1.Pg:vv.iN.Pg);
0230 Qg = x(vv.i1.Qg:vv.iN.Qg);
0231 V = Vm .* exp(1j*Va);
0232 
0233 %%-----  calculate return values  -----
0234 %% update voltages & generator outputs
0235 bus(:, VA) = Va * 180/pi;
0236 bus(:, VM) = Vm;
0237 gen(:, PG) = Pg * baseMVA;
0238 gen(:, QG) = Qg * baseMVA;
0239 gen(:, VG) = Vm(gen(:, GEN_BUS));
0240 
0241 %% compute branch flows
0242 Sf = V(branch(:, F_BUS)) .* conj(Yf * V);  %% cplx pwr at "from" bus, p.u.
0243 St = V(branch(:, T_BUS)) .* conj(Yt * V);  %% cplx pwr at "to" bus, p.u.
0244 branch(:, PF) = real(Sf) * baseMVA;
0245 branch(:, QF) = imag(Sf) * baseMVA;
0246 branch(:, PT) = real(St) * baseMVA;
0247 branch(:, QT) = imag(St) * baseMVA;
0248 
0249 %% line constraint is actually on square of limit
0250 %% so we must fix multipliers
0251 muSf = zeros(nl, 1);
0252 muSt = zeros(nl, 1);
0253 if ~isempty(il)
0254     muSf(il) = 2 * Lambda.ineqnonlin(1:nl2)       .* branch(il, RATE_A) / baseMVA;
0255     muSt(il) = 2 * Lambda.ineqnonlin((1:nl2)+nl2) .* branch(il, RATE_A) / baseMVA;
0256 end
0257 
0258 %% update Lagrange multipliers
0259 bus(:, MU_VMAX)  = Lambda.upper(vv.i1.Vm:vv.iN.Vm);
0260 bus(:, MU_VMIN)  = -Lambda.lower(vv.i1.Vm:vv.iN.Vm);
0261 gen(:, MU_PMAX)  = Lambda.upper(vv.i1.Pg:vv.iN.Pg) / baseMVA;
0262 gen(:, MU_PMIN)  = -Lambda.lower(vv.i1.Pg:vv.iN.Pg) / baseMVA;
0263 gen(:, MU_QMAX)  = Lambda.upper(vv.i1.Qg:vv.iN.Qg) / baseMVA;
0264 gen(:, MU_QMIN)  = -Lambda.lower(vv.i1.Qg:vv.iN.Qg) / baseMVA;
0265 bus(:, LAM_P)    = Lambda.eqnonlin(nn.i1.Pmis:nn.iN.Pmis) / baseMVA;
0266 bus(:, LAM_Q)    = Lambda.eqnonlin(nn.i1.Qmis:nn.iN.Qmis) / baseMVA;
0267 branch(:, MU_SF) = muSf / baseMVA;
0268 branch(:, MU_ST) = muSt / baseMVA;
0269 
0270 %% package up results
0271 nlnN = getN(om, 'nln');
0272 nlt = length(ilt);
0273 ngt = length(igt);
0274 nbx = length(ibx);
0275 
0276 %% extract multipliers for nonlinear constraints
0277 kl = find(Lambda.eqnonlin < 0);
0278 ku = find(Lambda.eqnonlin > 0);
0279 nl_mu_l = zeros(nlnN, 1);
0280 nl_mu_u = [zeros(2*nb, 1); muSf; muSt];
0281 nl_mu_l(kl) = -Lambda.eqnonlin(kl);
0282 nl_mu_u(ku) =  Lambda.eqnonlin(ku);
0283 
0284 %% extract multipliers for linear constraints
0285 kl = find(Lambda.eqlin < 0);
0286 ku = find(Lambda.eqlin > 0);
0287 
0288 mu_l = zeros(size(u));
0289 mu_l(ieq(kl)) = -Lambda.eqlin(kl);
0290 mu_l(igt) = Lambda.ineqlin(nlt+(1:ngt));
0291 mu_l(ibx) = Lambda.ineqlin(nlt+ngt+nbx+(1:nbx));
0292 
0293 mu_u = zeros(size(u));
0294 mu_u(ieq(ku)) = Lambda.eqlin(ku);
0295 mu_u(ilt) = Lambda.ineqlin(1:nlt);
0296 mu_u(ibx) = Lambda.ineqlin(nlt+ngt+(1:nbx));
0297 
0298 mu = struct( ...
0299   'var', struct('l', -Lambda.lower, 'u', Lambda.upper), ...
0300   'nln', struct('l', nl_mu_l, 'u', nl_mu_u), ...
0301   'lin', struct('l', mu_l, 'u', mu_u) );
0302 
0303 results = mpc;
0304 [results.bus, results.branch, results.gen, ...
0305     results.om, results.x, results.mu, results.f] = ...
0306         deal(bus, branch, gen, om, x, mu, f);
0307 
0308 pimul = [ ...
0309   results.mu.nln.l - results.mu.nln.u;
0310   results.mu.lin.l - results.mu.lin.u;
0311   -ones(ny>0, 1);
0312   results.mu.var.l - results.mu.var.u;
0313 ];
0314 raw = struct('xr', x, 'pimul', pimul, 'info', info, 'output', Output);
0315 
0316 
0317 %%-----  write_ktropts  -----
0318 function fname = write_ktropts(ktropts)
0319 
0320 %% generate file name
0321 fname = sprintf('ktropts_%06d.txt', fix(1e6*rand));
0322 
0323 %% open file
0324 [fd, msg] = fopen(fname, 'wt');     %% write options file
0325 if fd == -1
0326     error('could not create %d : %s', fname, msg);
0327 end
0328 
0329 %% write options
0330 fields = fieldnames(ktropts);
0331 for k = 1:length(fields)
0332     fprintf(fd, '%s %g\n', fields{k}, ktropts.(fields{k}));
0333 end
0334 
0335 %% close file
0336 if fd ~= 1
0337     fclose(fd);
0338 end

Generated on Mon 26-Jan-2015 15:21:31 by m2html © 2005