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    (deprecated) 2*nb+2*nl - Pmis, Qmis, Sf, St
0021 %               .l  lower bounds on nonlinear constraints
0022 %               .u  upper bounds on nonlinear constraints
0023 %           .nle    nonlinear equality constraints
0024 %           .nli    nonlinear inequality constraints
0025 %           .lin
0026 %               .l  lower bounds on linear constraints
0027 %               .u  upper bounds on linear constraints
0028 %
0029 %   SUCCESS     1 if solver converged successfully, 0 otherwise
0030 %
0031 %   RAW         raw output in form returned by MINOS
0032 %       .xr     final value of optimization variables
0033 %       .pimul  constraint multipliers
0034 %       .info   solver specific termination code
0035 %       .output solver specific output information
0036 %
0037 %   See also OPF, IPOPT.
0038 
0039 %   MATPOWER
0040 %   Copyright (c) 2000-2017, Power Systems Engineering Research Center (PSERC)
0041 %   by Ray Zimmerman, PSERC Cornell
0042 %   and Carlos E. Murillo-Sanchez, PSERC Cornell & Universidad Nacional de Colombia
0043 %
0044 %   This file is part of MATPOWER.
0045 %   Covered by the 3-clause BSD License (see LICENSE file for details).
0046 %   See https://matpower.org 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 = om.get_mpc();
0062 [baseMVA, bus, gen, branch, gencost] = ...
0063     deal(mpc.baseMVA, mpc.bus, mpc.gen, mpc.branch, mpc.gencost);
0064 [vv, ll, nne, nni] = om.get_idx();
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 = om.getN('var', 'y');   %% number of piece-wise linear costs
0071 
0072 %% linear constraints
0073 [A, l, u] = om.params_lin_constraint();
0074 
0075 %% bounds on optimization vars
0076 [x0, xmin, xmax] = om.params_var();
0077 
0078 %% build admittance matrices
0079 [Ybus, Yf, Yt] = makeYbus(baseMVA, bus, branch);
0080 
0081 %% try to select an interior initial point, unless requested not to
0082 if mpopt.opf.start < 2
0083     s = 1;                      %% set init point inside bounds by s
0084     lb = xmin; ub = xmax;
0085     lb(xmin == -Inf) = -1e10;   %% replace Inf with numerical proxies
0086     ub(xmax ==  Inf) =  1e10;
0087     x0 = (lb + ub) / 2;         %% set x0 mid-way between bounds
0088     k = find(xmin == -Inf & xmax < Inf);    %% if only bounded above
0089     x0(k) = xmax(k) - s;                    %% set just below upper bound
0090     k = find(xmin > -Inf & xmax == Inf);    %% if only bounded below
0091     x0(k) = xmin(k) + s;                    %% set just above lower bound
0092     Varefs = bus(bus(:, BUS_TYPE) == REF, VA) * (pi/180);
0093     Vmax = min(bus(:, VMAX), 1.5);
0094     Vmin = max(bus(:, VMIN), 0.5);
0095     Vm = (Vmax + Vmin) / 2;
0096     if mpopt.opf.v_cartesian
0097         V = Vm * exp(1j*Varefs(1));
0098         x0(vv.i1.Vr:vv.iN.Vr) = real(V);
0099         x0(vv.i1.Vi:vv.iN.Vi) = imag(V);
0100     else
0101         x0(vv.i1.Va:vv.iN.Va) = Varefs(1);  %% angles set to first reference angle
0102         x0(vv.i1.Vm:vv.iN.Vm) = Vm;         %% voltage magnitudes
0103         if ny > 0
0104             ipwl = find(gencost(:, MODEL) == PW_LINEAR);
0105             c = gencost(sub2ind(size(gencost), ipwl, NCOST+2*gencost(ipwl, NCOST)));    %% largest y-value in CCV data
0106             x0(vv.i1.y:vv.iN.y) = max(c) + 0.1 * abs(max(c));
0107         end
0108     end
0109 end
0110 
0111 %% find branches with flow limits
0112 il = find(branch(:, RATE_A) ~= 0 & branch(:, RATE_A) < 1e10);
0113 nl2 = length(il);           %% number of constrained lines
0114 nx = length(x0);
0115 
0116 %% replace equality variable bounds with an equality constraint
0117 %% (since IPOPT does not return shadow prices on variables that it eliminates)
0118 kk = find(xmin(nb+1:end) == xmax(nb+1:end));    %% all bounds except ref angles
0119 nk = length(kk);
0120 if nk
0121     kk = kk + nb;               %% adjust index for missing ref angles
0122     A = [ A; sparse((1:nk)', kk, 1, nk, nx) ];
0123     l = [ l; xmin(kk) ];
0124     u = [ u; xmax(kk) ];
0125     xmin(kk) = -Inf;
0126     xmax(kk) = Inf;
0127 end
0128 
0129 %%-----  run opf  -----
0130 %% build Jacobian and Hessian structure
0131 nA = size(A, 1);                %% number of original linear constraints
0132 % f = branch(:, F_BUS);                           %% list of "from" buses
0133 % t = branch(:, T_BUS);                           %% list of "to" buses
0134 % Cf = sparse(1:nl, f, ones(nl, 1), nl, nb);      %% connection matrix for line & from buses
0135 % Ct = sparse(1:nl, t, ones(nl, 1), nl, nb);      %% connection matrix for line & to buses
0136 % Cl = Cf + Ct;
0137 % Cb = Cl' * Cl + speye(nb);
0138 % Cl2 = Cl(il, :);
0139 % Cg = sparse(gen(:, GEN_BUS), (1:ng)', 1, nb, ng);
0140 % nz = nx - 2*(nb+ng);
0141 % nxtra = nx - 2*nb;
0142 % Js = [
0143 %     Cb      Cb      Cg              sparse(nb,ng)   sparse(nb,nz);
0144 %     Cb      Cb      sparse(nb,ng)   Cg              sparse(nb,nz);
0145 %     Cl2     Cl2     sparse(nl2, 2*ng)               sparse(nl2,nz);
0146 %     Cl2     Cl2     sparse(nl2, 2*ng)               sparse(nl2,nz);
0147 %     A;
0148 % ];
0149 % [f, df, d2f] = opf_costfcn(x0, om);
0150 % Hs = tril(d2f + [
0151 %     Cb  Cb  sparse(nb,nxtra);
0152 %     Cb  Cb  sparse(nb,nxtra);
0153 %     sparse(nxtra,nx);
0154 % ]);
0155 randx = rand(size(x0));
0156 [h, g, dh, dg] = opf_consfcn(randx, om, Ybus, Yf(il,:), Yt(il,:), mpopt, il);
0157 Js = [dg'; dh'; A];
0158 lam = struct('eqnonlin', rand(size(dg,2),1), 'ineqnonlin', rand(size(dh,2),1) );
0159 Hs = tril(opf_hessfcn(randx, lam, 1, om, Ybus, Yf(il,:), Yt(il,:), mpopt, il));
0160 neq = length(g);
0161 niq = length(h);
0162 
0163 %% set options struct for IPOPT
0164 options.ipopt = ipopt_options([], mpopt);
0165 
0166 %% extra data to pass to functions
0167 options.auxdata = struct( ...
0168     'om',       om, ...
0169     'Ybus',     Ybus, ...
0170     'Yf',       Yf(il,:), ...
0171     'Yt',       Yt(il,:), ...
0172     'mpopt',    mpopt, ...
0173     'il',       il, ...
0174     'A',        A, ...
0175     'nA',       nA, ...
0176     'neqnln',   neq, ...
0177     'niqnln',   niq, ...
0178     'Js',       Js, ...
0179     'Hs',       Hs    );
0180 
0181 % %% check Jacobian and Hessian structure
0182 % xr                  = rand(size(x0));
0183 % lambda              = rand(neq+niq, 1);
0184 % options.auxdata.Js  = jacobian(xr, options.auxdata);
0185 % options.auxdata.Hs  = tril(hessian(xr, 1, lambda, options.auxdata));
0186 % Js1 = options.auxdata.Js;
0187 % options.auxdata.Js = Js;
0188 % Hs1 = options.auxdata.Hs;
0189 % [i1, j1, s] = find(Js);
0190 % [i2, j2, s] = find(Js1);
0191 % if length(i1) ~= length(i2) || norm(i1-i2) ~= 0 || norm(j1-j2) ~= 0
0192 %     error('something''s wrong with the Jacobian structure');
0193 % end
0194 % [i1, j1, s] = find(Hs);
0195 % [i2, j2, s] = find(Hs1);
0196 % if length(i1) ~= length(i2) || norm(i1-i2) ~= 0 || norm(j1-j2) ~= 0
0197 %     error('something''s wrong with the Hessian structure');
0198 % end
0199 
0200 %% define variable and constraint bounds
0201 options.lb = xmin;
0202 options.ub = xmax;
0203 options.cl = [zeros(neq, 1);  -Inf(niq, 1); l];
0204 options.cu = [zeros(neq, 1); zeros(niq, 1); u];
0205 
0206 %% assign function handles
0207 funcs.objective         = @objective;
0208 funcs.gradient          = @gradient;
0209 funcs.constraints       = @constraints;
0210 funcs.jacobian          = @jacobian;
0211 funcs.hessian           = @hessian;
0212 funcs.jacobianstructure = @(d) Js;
0213 funcs.hessianstructure  = @(d) Hs;
0214 %funcs.jacobianstructure = @jacobianstructure;
0215 %funcs.hessianstructure  = @hessianstructure;
0216 
0217 %% run the optimization
0218 if have_fcn('ipopt_auxdata')
0219     [x, info] = ipopt_auxdata(x0,funcs,options);
0220 else
0221     [x, info] = ipopt(x0,funcs,options);
0222 end
0223 
0224 if info.status == 0 || info.status == 1
0225     success = 1;
0226 else
0227     success = 0;
0228 end
0229 if isfield(info, 'iter')
0230     output.iterations = info.iter;
0231 else
0232     output.iterations = [];
0233 end
0234 f = opf_costfcn(x, om);
0235 
0236 %% update solution data
0237 if mpopt.opf.v_cartesian
0238     Vi = x(vv.i1.Vi:vv.iN.Vi);
0239     Vr = x(vv.i1.Vr:vv.iN.Vr);
0240     V = Vr + 1j*Vi;
0241     Va = angle(V);
0242     Vm = abs(V);
0243 else
0244     Va = x(vv.i1.Va:vv.iN.Va);
0245     Vm = x(vv.i1.Vm:vv.iN.Vm);
0246     V = Vm .* exp(1j*Va);
0247 end
0248 Pg = x(vv.i1.Pg:vv.iN.Pg);
0249 Qg = x(vv.i1.Qg:vv.iN.Qg);
0250 
0251 %%-----  calculate return values  -----
0252 %% update voltages & generator outputs
0253 bus(:, VA) = Va * 180/pi;
0254 bus(:, VM) = Vm;
0255 gen(:, PG) = Pg * baseMVA;
0256 gen(:, QG) = Qg * baseMVA;
0257 gen(:, VG) = Vm(gen(:, GEN_BUS));
0258 
0259 %% compute branch flows
0260 Sf = V(branch(:, F_BUS)) .* conj(Yf * V);  %% cplx pwr at "from" bus, p.u.
0261 St = V(branch(:, T_BUS)) .* conj(Yt * V);  %% cplx pwr at "to" bus, p.u.
0262 branch(:, PF) = real(Sf) * baseMVA;
0263 branch(:, QF) = imag(Sf) * baseMVA;
0264 branch(:, PT) = real(St) * baseMVA;
0265 branch(:, QT) = imag(St) * baseMVA;
0266 
0267 %% line constraint is typically on square of limit
0268 %% so we must fix multipliers
0269 muSf = zeros(nl, 1);
0270 muSt = zeros(nl, 1);
0271 if ~isempty(il)
0272     if upper(mpopt.opf.flow_lim(1)) == 'P'
0273         muSf(il) = info.lambda(om.nle.N+(nni.i1.Sf:nni.iN.Sf));
0274         muSt(il) = info.lambda(om.nle.N+(nni.i1.St:nni.iN.St));
0275     else
0276         muSf(il) = 2 * info.lambda(om.nle.N+(nni.i1.Sf:nni.iN.Sf)) .* branch(il, RATE_A) / baseMVA;
0277         muSt(il) = 2 * info.lambda(om.nle.N+(nni.i1.St:nni.iN.St)) .* branch(il, RATE_A) / baseMVA;
0278     end
0279 end
0280 
0281 %% extract shadow prices for equality var bounds converted to eq constraints
0282 %% (since IPOPT does not return shadow prices on variables that it eliminates)
0283 if nk
0284     offset = om.nle.N + om.nli.N + nA - nk;
0285     lam_tmp = info.lambda(offset+(1:nk));
0286     kl = find(lam_tmp < 0);             %% lower bound binding
0287     ku = find(lam_tmp > 0);             %% upper bound binding
0288     info.zl(kk(kl)) = -lam_tmp(kl);
0289     info.zu(kk(ku)) =  lam_tmp(ku);
0290 
0291     info.lambda(offset+(1:nk)) = [];    %% remove these shadow prices
0292     nA = nA - nk;                       %% reduce dimension accordingly
0293 end
0294 
0295 %% update Lagrange multipliers
0296 if mpopt.opf.v_cartesian
0297     if om.userdata.veq
0298         lam = info.lambda(nne.i1.Veq:nne.iN.Veq);
0299         mu_Vmax = zeros(size(lam));
0300         mu_Vmin = zeros(size(lam));
0301         mu_Vmax(lam > 0) =  lam(lam > 0);
0302         mu_Vmin(lam < 0) = -lam(lam < 0);
0303         bus(om.userdata.veq, MU_VMAX) = mu_Vmax;
0304         bus(om.userdata.veq, MU_VMIN) = mu_Vmin;
0305     end
0306     bus(om.userdata.viq, MU_VMAX)  = info.lambda(om.nle.N+(nni.i1.Vmax:nni.iN.Vmax));
0307     bus(om.userdata.viq, MU_VMIN)  = info.lambda(om.nle.N+(nni.i1.Vmin:nni.iN.Vmin));
0308 else
0309     bus(:, MU_VMAX)  = info.zu(vv.i1.Vm:vv.iN.Vm);
0310     bus(:, MU_VMIN)  = info.zl(vv.i1.Vm:vv.iN.Vm);
0311 end
0312 gen(:, MU_PMAX)  = info.zu(vv.i1.Pg:vv.iN.Pg) / baseMVA;
0313 gen(:, MU_PMIN)  = info.zl(vv.i1.Pg:vv.iN.Pg) / baseMVA;
0314 gen(:, MU_QMAX)  = info.zu(vv.i1.Qg:vv.iN.Qg) / baseMVA;
0315 gen(:, MU_QMIN)  = info.zl(vv.i1.Qg:vv.iN.Qg) / baseMVA;
0316 if mpopt.opf.current_balance
0317     %% convert current balance shadow prices to equivalent lamP and lamQ
0318     %% P + jQ = (Vr + jVi) * (M - jN)
0319     %% M = (Vr P + Vi Q) / (Vr^2 + Vi^2)
0320     %% N = (Vi P - Vr Q) / (Vr^2 + Vi^2)
0321     %% lamP = df/dP = df/dM * dM/dP + df/dN + dN/dP
0322     %% lamQ = df/dQ = df/dM * dM/dQ + df/dN + dN/dQ
0323     VV = V ./ (V .* conj(V));   %% V / Vm^2
0324     VVr = real(VV);
0325     VVi = imag(VV);
0326     lamM = info.lambda(nne.i1.rImis:nne.iN.rImis);
0327     lamN = info.lambda(nne.i1.iImis:nne.iN.iImis);
0328     bus(:, LAM_P) = (VVr.*lamM + VVi.*lamN) / baseMVA;
0329     bus(:, LAM_Q) = (VVi.*lamM - VVr.*lamN) / baseMVA;
0330 else
0331     bus(:, LAM_P) = info.lambda(nne.i1.Pmis:nne.iN.Pmis) / baseMVA;
0332     bus(:, LAM_Q) = info.lambda(nne.i1.Qmis:nne.iN.Qmis) / baseMVA;
0333 end
0334 branch(:, MU_SF) = muSf / baseMVA;
0335 branch(:, MU_ST) = muSt / baseMVA;
0336 
0337 %% package up results
0338 nlnN = 2*nb + 2*nl;     %% because muSf and muSt are nl x 1, not nl2 x 1
0339 
0340 %% extract multipliers for nonlinear constraints
0341 kl = find(info.lambda(1:om.nle.N) < 0);
0342 ku = find(info.lambda(1:om.nle.N) > 0);
0343 nl_mu_l = zeros(nlnN, 1);
0344 nl_mu_u = [zeros(2*nb, 1); muSf; muSt];
0345 nl_mu_l(kl) = -info.lambda(kl);
0346 nl_mu_u(ku) =  info.lambda(ku);
0347 
0348 %% extract multipliers for linear constraints
0349 lam_lin = info.lambda(om.nle.N+om.nli.N+(1:nA));    %% lambda for linear constraints
0350 kl = find(lam_lin < 0);                     %% lower bound binding
0351 ku = find(lam_lin > 0);                     %% upper bound binding
0352 mu_l = zeros(nA, 1);
0353 mu_l(kl) = -lam_lin(kl);
0354 mu_u = zeros(nA, 1);
0355 mu_u(ku) = lam_lin(ku);
0356 
0357 mu = struct( ...
0358   'var', struct('l', info.zl, 'u', info.zu), ...
0359   'nln', struct('l', nl_mu_l, 'u', nl_mu_u), ...
0360   'nle', info.lambda(1:om.nle.N), ...
0361   'nli', info.lambda(om.nle.N + (1:om.nli.N)), ...
0362   'lin', struct('l', mu_l, 'u', mu_u) );
0363 
0364 results = mpc;
0365 [results.bus, results.branch, results.gen, ...
0366     results.om, results.x, results.mu, results.f] = ...
0367         deal(bus, branch, gen, om, x, mu, f);
0368 
0369 pimul = [ ...
0370   results.mu.nln.l - results.mu.nln.u;
0371   results.mu.lin.l - results.mu.lin.u;
0372   -ones(ny>0, 1);
0373   results.mu.var.l - results.mu.var.u;
0374 ];
0375 raw = struct('xr', x, 'pimul', pimul, 'info', info.status, 'output', output);
0376 
0377 
0378 %-----  callback functions  -----
0379 function f = objective(x, d)
0380 f = opf_costfcn(x, d.om);
0381 
0382 function df = gradient(x, d)
0383 [f, df] = opf_costfcn(x, d.om);
0384 
0385 function c = constraints(x, d)
0386 [hn, gn] = opf_consfcn(x, d.om, d.Ybus, d.Yf, d.Yt, d.mpopt, d.il);
0387 if isempty(d.A)
0388     c = [gn; hn];
0389 else
0390     c = [gn; hn; d.A*x];
0391 end
0392 
0393 function J = jacobian(x, d)
0394 [hn, gn, dhn, dgn] = opf_consfcn(x, d.om, d.Ybus, d.Yf, d.Yt, d.mpopt, d.il);
0395 J = [dgn'; dhn'; d.A];
0396 
0397 function H = hessian(x, sigma, lambda, d)
0398 lam.eqnonlin   = lambda(1:d.neqnln);
0399 lam.ineqnonlin = lambda(d.neqnln+(1:d.niqnln));
0400 H = tril(opf_hessfcn(x, lam, sigma, d.om, d.Ybus, d.Yf, d.Yt, d.mpopt, d.il));
0401 
0402 % function Js = jacobianstructure(d)
0403 % Js = d.Js;
0404 %
0405 % function Hs = hessianstructure(d)
0406 % Hs = d.Hs;
ipoptopf_solver

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SUBFUNCTIONS

SOURCE CODE
