Description of ipoptopf

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;
ipoptopf_solver

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SUBFUNCTIONS

SOURCE CODE
