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

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