Description of ktropf

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 vector.
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,v 1.1 2011/06/17 20:28:44 cvs Exp $
0039 %   by Ray Zimmerman, PSERC Cornell
0040 %   and Carlos E. Murillo-Sanchez, PSERC Cornell & Universidad Autonoma de Manizales
0041 %   Copyright (c) 2000-2011 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;   %% generate a KNITRO options file on the fly
0081 verbose = mpopt(31);    %% VERBOSE
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 ll = xmin; uu = xmax;
0117 ll(xmin == -Inf) = -1e10;   %% replace Inf with numerical proxies
0118 uu(xmax ==  Inf) =  1e10;
0119 x0 = (ll + uu) / 2;
0120 Varefs = bus(bus(:, BUS_TYPE) == REF, VA) * (pi/180);
0121 x0(vv.i1.Va:vv.iN.Va) = Varefs(1);  %% angles set to first reference angle
0122 if ny > 0
0123     ipwl = find(gencost(:, MODEL) == PW_LINEAR);
0124 %     PQ = [gen(:, PMAX); gen(:, QMAX)];
0125 %     c = totcost(gencost(ipwl, :), PQ(ipwl));
0126     c = gencost(sub2ind(size(gencost), ipwl, NCOST+2*gencost(ipwl, NCOST)));    %% largest y-value in CCV data
0127     x0(vv.i1.y:vv.iN.y) = max(c) + 0.1 * abs(max(c));
0128 %     x0(vv.i1.y:vv.iN.y) = c + 0.1 * abs(c);
0129 end
0130 
0131 %% find branches with flow limits
0132 il = find(branch(:, RATE_A) ~= 0 & branch(:, RATE_A) < 1e10);
0133 nl2 = length(il);           %% number of constrained lines
0134 
0135 %% build Jacobian and Hessian structure
0136 nA = size(A, 1);                %% number of original linear constraints
0137 nx = length(x0);
0138 f = branch(:, F_BUS);                           %% list of "from" buses
0139 t = branch(:, T_BUS);                           %% list of "to" buses
0140 Cf = sparse(1:nl, f, ones(nl, 1), nl, nb);      %% connection matrix for line & from buses
0141 Ct = sparse(1:nl, t, ones(nl, 1), nl, nb);      %% connection matrix for line & to buses
0142 Cl = Cf + Ct;
0143 Cb = Cl' * Cl + speye(nb);
0144 Cl2 = Cl(il, :);
0145 Cg = sparse(gen(:, GEN_BUS), (1:ng)', 1, nb, ng);
0146 nz = nx - 2*(nb+ng);
0147 nxtra = nx - 2*nb;
0148 Js = [
0149     Cl2     Cl2     sparse(nl2, 2*ng)               sparse(nl2,nz);
0150     Cl2     Cl2     sparse(nl2, 2*ng)               sparse(nl2,nz);
0151     Cb      Cb      Cg              sparse(nb,ng)   sparse(nb,nz);
0152     Cb      Cb      sparse(nb,ng)   Cg              sparse(nb,nz);
0153 ];
0154 [f, df, d2f] = opf_costfcn(x0, om);
0155 Hs = d2f + [
0156     Cb  Cb  sparse(nb,nxtra);
0157     Cb  Cb  sparse(nb,nxtra);
0158     sparse(nxtra,nx);
0159 ];
0160 
0161 %% basic optimset options needed for ktrlink
0162 hess_fcn = @(x, lambda)opf_hessfcn(x, lambda, 1, om, Ybus, Yf(il,:), Yt(il,:), mpopt, il);
0163 fmoptions = optimset('GradObj', 'on', 'GradConstr', 'on', ...
0164                 'Hessian', 'user-supplied', 'HessFcn', hess_fcn, ...
0165                 'JacobPattern', Js, 'HessPattern', Hs );
0166 if use_ktropts_file
0167     if mpopt(58)
0168         opt_fname = sprintf('knitro_user_options_%d.txt', mpopt(58));
0169     else
0170         %% create ktropts file
0171         ktropts.algorithm           = 1;
0172         ktropts.outlev              = verbose;
0173         ktropts.feastol             = mpopt(16);
0174         ktropts.xtol                = mpopt(17);
0175         ktropts.opttol              = mpopt(18);
0176         if mpopt(19) ~= 0
0177             ktropts.maxit           = mpopt(19);
0178         end
0179         ktropts.bar_directinterval  = 0;
0180         opt_fname = write_ktropts(ktropts);
0181     end
0182 else
0183     fmoptions = optimset(fmoptions, 'Algorithm', 'interior-point', ...
0184         'TolCon', mpopt(16), 'TolX', mpopt(17), 'TolFun', mpopt(18) );
0185     if mpopt(19) ~= 0
0186         fmoptions = optimset(fmoptions, 'MaxIter', mpopt(19), ...
0187                     'MaxFunEvals', 4 * mpopt(19));
0188     end
0189     if verbose == 0,
0190       fmoptions.Display = 'off';
0191     elseif verbose == 1
0192       fmoptions.Display = 'iter';
0193     else
0194       fmoptions.Display = 'testing';
0195     end
0196     opt_fname = [];
0197 end
0198 
0199 %%-----  run opf  -----
0200 f_fcn = @(x)opf_costfcn(x, om);
0201 gh_fcn = @(x)opf_consfcn(x, om, Ybus, Yf(il,:), Yt(il,:), mpopt, il);
0202 [x, f, info, Output, Lambda] = ktrlink(f_fcn, x0, Af, bf, Afeq, bfeq, ...
0203                                     xmin, xmax, gh_fcn, fmoptions, opt_fname);
0204 success = (info == 0);
0205 
0206 %% delete ktropts file
0207 if use_ktropts_file && ~mpopt(58)   %% ... but only if I created it
0208     delete(opt_fname);
0209 end
0210 
0211 %% update solution data
0212 Va = x(vv.i1.Va:vv.iN.Va);
0213 Vm = x(vv.i1.Vm:vv.iN.Vm);
0214 Pg = x(vv.i1.Pg:vv.iN.Pg);
0215 Qg = x(vv.i1.Qg:vv.iN.Qg);
0216 V = Vm .* exp(1j*Va);
0217 
0218 %%-----  calculate return values  -----
0219 %% update voltages & generator outputs
0220 bus(:, VA) = Va * 180/pi;
0221 bus(:, VM) = Vm;
0222 gen(:, PG) = Pg * baseMVA;
0223 gen(:, QG) = Qg * baseMVA;
0224 gen(:, VG) = Vm(gen(:, GEN_BUS));
0225 
0226 %% compute branch flows
0227 Sf = V(branch(:, F_BUS)) .* conj(Yf * V);  %% cplx pwr at "from" bus, p.u.
0228 St = V(branch(:, T_BUS)) .* conj(Yt * V);  %% cplx pwr at "to" bus, p.u.
0229 branch(:, PF) = real(Sf) * baseMVA;
0230 branch(:, QF) = imag(Sf) * baseMVA;
0231 branch(:, PT) = real(St) * baseMVA;
0232 branch(:, QT) = imag(St) * baseMVA;
0233 
0234 %% line constraint is actually on square of limit
0235 %% so we must fix multipliers
0236 muSf = zeros(nl, 1);
0237 muSt = zeros(nl, 1);
0238 if ~isempty(il)
0239     muSf(il) = 2 * Lambda.ineqnonlin(1:nl2)       .* branch(il, RATE_A) / baseMVA;
0240     muSt(il) = 2 * Lambda.ineqnonlin((1:nl2)+nl2) .* branch(il, RATE_A) / baseMVA;
0241 end
0242 
0243 %% update Lagrange multipliers
0244 bus(:, MU_VMAX)  = Lambda.upper(vv.i1.Vm:vv.iN.Vm);
0245 bus(:, MU_VMIN)  = -Lambda.lower(vv.i1.Vm:vv.iN.Vm);
0246 gen(:, MU_PMAX)  = Lambda.upper(vv.i1.Pg:vv.iN.Pg) / baseMVA;
0247 gen(:, MU_PMIN)  = -Lambda.lower(vv.i1.Pg:vv.iN.Pg) / baseMVA;
0248 gen(:, MU_QMAX)  = Lambda.upper(vv.i1.Qg:vv.iN.Qg) / baseMVA;
0249 gen(:, MU_QMIN)  = -Lambda.lower(vv.i1.Qg:vv.iN.Qg) / baseMVA;
0250 bus(:, LAM_P)    = Lambda.eqnonlin(nn.i1.Pmis:nn.iN.Pmis) / baseMVA;
0251 bus(:, LAM_Q)    = Lambda.eqnonlin(nn.i1.Qmis:nn.iN.Qmis) / baseMVA;
0252 branch(:, MU_SF) = muSf / baseMVA;
0253 branch(:, MU_ST) = muSt / baseMVA;
0254 
0255 %% package up results
0256 nlnN = getN(om, 'nln');
0257 nlt = length(ilt);
0258 ngt = length(igt);
0259 nbx = length(ibx);
0260 
0261 %% extract multipliers for nonlinear constraints
0262 kl = find(Lambda.eqnonlin < 0);
0263 ku = find(Lambda.eqnonlin > 0);
0264 nl_mu_l = zeros(nlnN, 1);
0265 nl_mu_u = [zeros(2*nb, 1); muSf; muSt];
0266 nl_mu_l(kl) = -Lambda.eqnonlin(kl);
0267 nl_mu_u(ku) =  Lambda.eqnonlin(ku);
0268 
0269 %% extract multipliers for linear constraints
0270 kl = find(Lambda.eqlin < 0);
0271 ku = find(Lambda.eqlin > 0);
0272 
0273 mu_l = zeros(size(u));
0274 mu_l(ieq(kl)) = -Lambda.eqlin(kl);
0275 mu_l(igt) = Lambda.ineqlin(nlt+(1:ngt));
0276 mu_l(ibx) = Lambda.ineqlin(nlt+ngt+nbx+(1:nbx));
0277 
0278 mu_u = zeros(size(u));
0279 mu_u(ieq(ku)) = Lambda.eqlin(ku);
0280 mu_u(ilt) = Lambda.ineqlin(1:nlt);
0281 mu_u(ibx) = Lambda.ineqlin(nlt+ngt+(1:nbx));
0282 
0283 mu = struct( ...
0284   'var', struct('l', -Lambda.lower, 'u', Lambda.upper), ...
0285   'nln', struct('l', nl_mu_l, 'u', nl_mu_u), ...
0286   'lin', struct('l', mu_l, 'u', mu_u) );
0287 
0288 results = mpc;
0289 [results.bus, results.branch, results.gen, ...
0290     results.om, results.x, results.mu, results.f] = ...
0291         deal(bus, branch, gen, om, x, mu, f);
0292 
0293 pimul = [ ...
0294   results.mu.nln.l - results.mu.nln.u;
0295   results.mu.lin.l - results.mu.lin.u;
0296   -ones(ny>0, 1);
0297   results.mu.var.l - results.mu.var.u;
0298 ];
0299 raw = struct('xr', x, 'pimul', pimul, 'info', info, 'output', Output);
0300 
0301 
0302 %%-----  write_ktropts  -----
0303 function fname = write_ktropts(ktropts)
0304 
0305 %% generate file name
0306 fname = sprintf('ktropts_%06d.txt', fix(1e6*rand));
0307 
0308 %% open file
0309 [fd, msg] = fopen(fname, 'wt');     %% write options file
0310 if fd == -1
0311     error('could not create %d : %s', fname, msg);
0312 end
0313 
0314 %% write options
0315 fields = fieldnames(ktropts);
0316 for k = 1:length(fields)
0317     fprintf(fd, '%s %g\n', fields{k}, ktropts.(fields{k}));
0318 end
0319 
0320 %% close file
0321 if fd ~= 1
0322     fclose(fd);
0323 end
ktropf_solver

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SUBFUNCTIONS

SOURCE CODE
