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

Generated on Mon 26-Jan-2015 15:21:31 by m2html © 2005