Home > matpower5.1 > ipoptopf_solver.m

ipoptopf_solver

PURPOSE ^

IPOPTOPF_SOLVER Solves AC optimal power flow using MIPS.

SYNOPSIS ^

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

DESCRIPTION ^

IPOPTOPF_SOLVER  Solves AC optimal power flow using MIPS.

   [RESULTS, SUCCESS, RAW] = IPOPTOPF_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, IPOPT.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 function [results, success, raw] = ipoptopf_solver(om, mpopt)
0002 %IPOPTOPF_SOLVER  Solves AC optimal power flow using MIPS.
0003 %
0004 %   [RESULTS, SUCCESS, RAW] = IPOPTOPF_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, IPOPT.
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: ipoptopf_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 %% unpack data
0061 mpc = get_mpc(om);
0062 [baseMVA, bus, gen, branch, gencost] = ...
0063     deal(mpc.baseMVA, mpc.bus, mpc.gen, mpc.branch, mpc.gencost);
0064 [vv, ll, nn] = get_idx(om);
0065 
0066 %% problem dimensions
0067 nb = size(bus, 1);          %% number of buses
0068 ng = size(gen, 1);          %% number of gens
0069 nl = size(branch, 1);       %% number of branches
0070 ny = getN(om, 'var', 'y');  %% number of piece-wise linear costs
0071 
0072 %% linear constraints
0073 [A, l, u] = linear_constraints(om);
0074 
0075 %% bounds on optimization vars
0076 [x0, xmin, xmax] = getv(om);
0077 
0078 %% build admittance matrices
0079 [Ybus, Yf, Yt] = makeYbus(baseMVA, bus, branch);
0080 
0081 %% try to select an interior initial point
0082 if mpopt.opf.init_from_mpc ~= 1
0083     ll = xmin; uu = xmax;
0084     ll(xmin == -Inf) = -1e10;   %% replace Inf with numerical proxies
0085     uu(xmax ==  Inf) =  1e10;
0086     x0 = (ll + uu) / 2;         %% set x0 mid-way between bounds
0087     k = find(xmin == -Inf & xmax < Inf);    %% if only bounded above
0088     x0(k) = xmax(k) - 1;                    %% set just below upper bound
0089     k = find(xmin > -Inf & xmax == Inf);    %% if only bounded below
0090     x0(k) = xmin(k) + 1;                    %% set just above lower bound
0091     Varefs = bus(bus(:, BUS_TYPE) == REF, VA) * (pi/180);
0092     x0(vv.i1.Va:vv.iN.Va) = Varefs(1);  %% angles set to first reference angle
0093     if ny > 0
0094         ipwl = find(gencost(:, MODEL) == PW_LINEAR);
0095     %     PQ = [gen(:, PMAX); gen(:, QMAX)];
0096     %     c = totcost(gencost(ipwl, :), PQ(ipwl));
0097         c = gencost(sub2ind(size(gencost), ipwl, NCOST+2*gencost(ipwl, NCOST)));    %% largest y-value in CCV data
0098         x0(vv.i1.y:vv.iN.y) = max(c) + 0.1 * abs(max(c));
0099     %     x0(vv.i1.y:vv.iN.y) = c + 0.1 * abs(c);
0100     end
0101 end
0102 
0103 %% find branches with flow limits
0104 il = find(branch(:, RATE_A) ~= 0 & branch(:, RATE_A) < 1e10);
0105 nl2 = length(il);           %% number of constrained lines
0106 
0107 %%-----  run opf  -----
0108 %% build Jacobian and Hessian structure
0109 nA = size(A, 1);                %% number of original linear constraints
0110 nx = length(x0);
0111 f = branch(:, F_BUS);                           %% list of "from" buses
0112 t = branch(:, T_BUS);                           %% list of "to" buses
0113 Cf = sparse(1:nl, f, ones(nl, 1), nl, nb);      %% connection matrix for line & from buses
0114 Ct = sparse(1:nl, t, ones(nl, 1), nl, nb);      %% connection matrix for line & to buses
0115 Cl = Cf + Ct;
0116 Cb = Cl' * Cl + speye(nb);
0117 Cl2 = Cl(il, :);
0118 Cg = sparse(gen(:, GEN_BUS), (1:ng)', 1, nb, ng);
0119 nz = nx - 2*(nb+ng);
0120 nxtra = nx - 2*nb;
0121 Js = [
0122     Cb      Cb      Cg              sparse(nb,ng)   sparse(nb,nz);
0123     Cb      Cb      sparse(nb,ng)   Cg              sparse(nb,nz);
0124     Cl2     Cl2     sparse(nl2, 2*ng)               sparse(nl2,nz);
0125     Cl2     Cl2     sparse(nl2, 2*ng)               sparse(nl2,nz);
0126     A;
0127 ];
0128 [f, df, d2f] = opf_costfcn(x0, om);
0129 Hs = tril(d2f + [
0130     Cb  Cb  sparse(nb,nxtra);
0131     Cb  Cb  sparse(nb,nxtra);
0132     sparse(nxtra,nx);
0133 ]);
0134 
0135 %% set options struct for IPOPT
0136 options.ipopt = ipopt_options([], mpopt);
0137 
0138 %% extra data to pass to functions
0139 options.auxdata = struct( ...
0140     'om',       om, ...
0141     'Ybus',     Ybus, ...
0142     'Yf',       Yf(il,:), ...
0143     'Yt',       Yt(il,:), ...
0144     'mpopt',    mpopt, ...
0145     'il',       il, ...
0146     'A',        A, ...
0147     'nA',       nA, ...
0148     'neqnln',   2*nb, ...
0149     'niqnln',   2*nl2, ...
0150     'Js',       Js, ...
0151     'Hs',       Hs    );
0152 
0153 % %% check Jacobian and Hessian structure
0154 % xr                  = rand(size(x0));
0155 % lambda              = rand(2*nb+2*nl2, 1);
0156 % options.auxdata.Js  = jacobian(xr, options.auxdata);
0157 % options.auxdata.Hs  = tril(hessian(xr, 1, lambda, options.auxdata));
0158 % Js1 = options.auxdata.Js;
0159 % options.auxdata.Js = Js;
0160 % Hs1 = options.auxdata.Hs;
0161 % [i1, j1, s] = find(Js);
0162 % [i2, j2, s] = find(Js1);
0163 % if length(i1) ~= length(i2) || norm(i1-i2) ~= 0 || norm(j1-j2) ~= 0
0164 %     error('something''s wrong with the Jacobian structure');
0165 % end
0166 % [i1, j1, s] = find(Hs);
0167 % [i2, j2, s] = find(Hs1);
0168 % if length(i1) ~= length(i2) || norm(i1-i2) ~= 0 || norm(j1-j2) ~= 0
0169 %     error('something''s wrong with the Hessian structure');
0170 % end
0171 
0172 %% define variable and constraint bounds
0173 options.lb = xmin;
0174 options.ub = xmax;
0175 options.cl = [zeros(2*nb, 1);  -Inf(2*nl2, 1); l];
0176 options.cu = [zeros(2*nb, 1); zeros(2*nl2, 1); u];
0177 
0178 %% assign function handles
0179 funcs.objective         = @objective;
0180 funcs.gradient          = @gradient;
0181 funcs.constraints       = @constraints;
0182 funcs.jacobian          = @jacobian;
0183 funcs.hessian           = @hessian;
0184 funcs.jacobianstructure = @(d) Js;
0185 funcs.hessianstructure  = @(d) Hs;
0186 %funcs.jacobianstructure = @jacobianstructure;
0187 %funcs.hessianstructure  = @hessianstructure;
0188 
0189 %% run the optimization
0190 if have_fcn('ipopt_auxdata')
0191     [x, info] = ipopt_auxdata(x0,funcs,options);
0192 else
0193     [x, info] = ipopt(x0,funcs,options);
0194 end
0195 
0196 if info.status == 0 || info.status == 1
0197     success = 1;
0198 else
0199     success = 0;
0200 end
0201 if isfield(info, 'iter')
0202     output.iterations = info.iter;
0203 else
0204     output.iterations = [];
0205 end
0206 f = opf_costfcn(x, om);
0207 
0208 %% update solution data
0209 Va = x(vv.i1.Va:vv.iN.Va);
0210 Vm = x(vv.i1.Vm:vv.iN.Vm);
0211 Pg = x(vv.i1.Pg:vv.iN.Pg);
0212 Qg = x(vv.i1.Qg:vv.iN.Qg);
0213 V = Vm .* exp(1j*Va);
0214 
0215 %%-----  calculate return values  -----
0216 %% update voltages & generator outputs
0217 bus(:, VA) = Va * 180/pi;
0218 bus(:, VM) = Vm;
0219 gen(:, PG) = Pg * baseMVA;
0220 gen(:, QG) = Qg * baseMVA;
0221 gen(:, VG) = Vm(gen(:, GEN_BUS));
0222 
0223 %% compute branch flows
0224 Sf = V(branch(:, F_BUS)) .* conj(Yf * V);  %% cplx pwr at "from" bus, p.u.
0225 St = V(branch(:, T_BUS)) .* conj(Yt * V);  %% cplx pwr at "to" bus, p.u.
0226 branch(:, PF) = real(Sf) * baseMVA;
0227 branch(:, QF) = imag(Sf) * baseMVA;
0228 branch(:, PT) = real(St) * baseMVA;
0229 branch(:, QT) = imag(St) * baseMVA;
0230 
0231 %% line constraint is actually on square of limit
0232 %% so we must fix multipliers
0233 muSf = zeros(nl, 1);
0234 muSt = zeros(nl, 1);
0235 if ~isempty(il)
0236     muSf(il) = 2 * info.lambda(2*nb+    (1:nl2)) .* branch(il, RATE_A) / baseMVA;
0237     muSt(il) = 2 * info.lambda(2*nb+nl2+(1:nl2)) .* branch(il, RATE_A) / baseMVA;
0238 end
0239 
0240 %% update Lagrange multipliers
0241 bus(:, MU_VMAX)  = info.zu(vv.i1.Vm:vv.iN.Vm);
0242 bus(:, MU_VMIN)  = info.zl(vv.i1.Vm:vv.iN.Vm);
0243 gen(:, MU_PMAX)  = info.zu(vv.i1.Pg:vv.iN.Pg) / baseMVA;
0244 gen(:, MU_PMIN)  = info.zl(vv.i1.Pg:vv.iN.Pg) / baseMVA;
0245 gen(:, MU_QMAX)  = info.zu(vv.i1.Qg:vv.iN.Qg) / baseMVA;
0246 gen(:, MU_QMIN)  = info.zl(vv.i1.Qg:vv.iN.Qg) / baseMVA;
0247 bus(:, LAM_P)    = info.lambda(nn.i1.Pmis:nn.iN.Pmis) / baseMVA;
0248 bus(:, LAM_Q)    = info.lambda(nn.i1.Qmis:nn.iN.Qmis) / baseMVA;
0249 branch(:, MU_SF) = muSf / baseMVA;
0250 branch(:, MU_ST) = muSt / baseMVA;
0251 
0252 %% package up results
0253 nlnN = getN(om, 'nln');
0254 
0255 %% extract multipliers for nonlinear constraints
0256 kl = find(info.lambda(1:2*nb) < 0);
0257 ku = find(info.lambda(1:2*nb) > 0);
0258 nl_mu_l = zeros(nlnN, 1);
0259 nl_mu_u = [zeros(2*nb, 1); muSf; muSt];
0260 nl_mu_l(kl) = -info.lambda(kl);
0261 nl_mu_u(ku) =  info.lambda(ku);
0262 
0263 %% extract multipliers for linear constraints
0264 lam_lin = info.lambda(2*nb+2*nl2+(1:nA));   %% lambda for linear constraints
0265 kl = find(lam_lin < 0);                     %% lower bound binding
0266 ku = find(lam_lin > 0);                     %% upper bound binding
0267 mu_l = zeros(nA, 1);
0268 mu_l(kl) = -lam_lin(kl);
0269 mu_u = zeros(nA, 1);
0270 mu_u(ku) = lam_lin(ku);
0271 
0272 mu = struct( ...
0273   'var', struct('l', info.zl, 'u', info.zu), ...
0274   'nln', struct('l', nl_mu_l, 'u', nl_mu_u), ...
0275   'lin', struct('l', mu_l, 'u', mu_u) );
0276 
0277 results = mpc;
0278 [results.bus, results.branch, results.gen, ...
0279     results.om, results.x, results.mu, results.f] = ...
0280         deal(bus, branch, gen, om, x, mu, f);
0281 
0282 pimul = [ ...
0283   results.mu.nln.l - results.mu.nln.u;
0284   results.mu.lin.l - results.mu.lin.u;
0285   -ones(ny>0, 1);
0286   results.mu.var.l - results.mu.var.u;
0287 ];
0288 raw = struct('xr', x, 'pimul', pimul, 'info', info.status, 'output', output);
0289 
0290 
0291 %-----  callback functions  -----
0292 function f = objective(x, d)
0293 f = opf_costfcn(x, d.om);
0294 
0295 function df = gradient(x, d)
0296 [f, df] = opf_costfcn(x, d.om);
0297 
0298 function c = constraints(x, d)
0299 [hn, gn] = opf_consfcn(x, d.om, d.Ybus, d.Yf, d.Yt, d.mpopt, d.il);
0300 if isempty(d.A)
0301     c = [gn; hn];
0302 else
0303     c = [gn; hn; d.A*x];
0304 end
0305 
0306 function J = jacobian(x, d)
0307 [hn, gn, dhn, dgn] = opf_consfcn(x, d.om, d.Ybus, d.Yf, d.Yt, d.mpopt, d.il);
0308 J = [dgn'; dhn'; d.A];
0309 
0310 function H = hessian(x, sigma, lambda, d)
0311 lam.eqnonlin   = lambda(1:d.neqnln);
0312 lam.ineqnonlin = lambda(d.neqnln+(1:d.niqnln));
0313 H = tril(opf_hessfcn(x, lam, sigma, d.om, d.Ybus, d.Yf, d.Yt, d.mpopt, d.il));
0314 
0315 % function Js = jacobianstructure(d)
0316 % Js = d.Js;
0317 %
0318 % function Hs = hessianstructure(d)
0319 % Hs = d.Hs;

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