fmincopf_solver

FMINCOPF_SOLVER  Solves AC optimal power flow using FMINCON.

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

0001 function [results, success, raw] = fmincopf_solver(om, mpopt)
0002 %FMINCOPF_SOLVER  Solves AC optimal power flow using FMINCON.
0003 %
0004 %   [RESULTS, SUCCESS, RAW] = FMINCOPF_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, FMINCON.
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 
0057 %% unpack data
0058 mpc = get_mpc(om);
0059 [baseMVA, bus, gen, branch] = ...
0060     deal(mpc.baseMVA, mpc.bus, mpc.gen, mpc.branch);
0061 [vv, ll, nn] = get_idx(om);
0062 
0063 %% problem dimensions
0064 nb = size(bus, 1);          %% number of buses
0065 nl = size(branch, 1);       %% number of branches
0066 ny = getN(om, 'var', 'y');  %% number of piece-wise linear costs
0067 
0068 %% bounds on optimization vars
0069 [x0, LB, UB] = getv(om);
0070 
0071 %% linear constraints
0072 [A, l, u] = linear_constraints(om);
0073 
0074 %% split l <= A*x <= u into less than, equal to, greater than, and
0075 %% doubly-bounded sets
0076 ieq = find( abs(u-l) <= eps );          %% equality
0077 igt = find( u >=  1e10 & l > -1e10 );   %% greater than, unbounded above
0078 ilt = find( l <= -1e10 & u <  1e10 );   %% less than, unbounded below
0079 ibx = find( (abs(u-l) > eps) & (u < 1e10) & (l > -1e10) );
0080 Af  = [ A(ilt, :); -A(igt, :); A(ibx, :); -A(ibx, :) ];
0081 bf  = [ u(ilt);   -l(igt);     u(ibx);    -l(ibx)];
0082 Afeq = A(ieq, :);
0083 bfeq = u(ieq);
0084 
0085 %% build admittance matrices
0086 [Ybus, Yf, Yt] = makeYbus(baseMVA, bus, branch);
0087 
0088 %% find branches with flow limits
0089 il = find(branch(:, RATE_A) ~= 0 & branch(:, RATE_A) < 1e10);
0090 nl2 = length(il);           %% number of constrained lines
0091 
0092 %% basic optimset options needed for fmincon
0093 fmoptions = optimset('GradObj', 'on', 'GradConstr', 'on', ...
0094             'TolCon', mpopt.opf.violation, 'TolX', mpopt.fmincon.tol_x, ...
0095             'TolFun', mpopt.fmincon.tol_f );
0096 if mpopt.fmincon.max_it ~= 0
0097     fmoptions = optimset(fmoptions, 'MaxIter', mpopt.fmincon.max_it, ...
0098             'MaxFunEvals', 4 * mpopt.fmincon.max_it);
0099 end
0100 
0101 if mpopt.verbose == 0,
0102   fmoptions.Display = 'off';
0103 elseif mpopt.verbose == 1
0104   fmoptions.Display = 'iter';
0105 else
0106   fmoptions.Display = 'testing';
0107 end
0108 
0109 %% select algorithm
0110 if have_fcn('fmincon_ipm')
0111   switch mpopt.fmincon.alg
0112     case 1              %% active-set (does not use sparse matrices, not suitable for large problems)
0113       fmoptions = optimset(fmoptions, 'Algorithm', 'active-set');
0114       Af = full(Af);
0115       Afeq = full(Afeq);
0116     case 2              %% interior-point, w/ default 'bfgs' Hessian approx
0117       fmoptions = optimset(fmoptions, 'Algorithm', 'interior-point');
0118     case 3              %% interior-point, w/ 'lbfgs' Hessian approx
0119       fmoptions = optimset(fmoptions, 'Algorithm', 'interior-point', 'Hessian','lbfgs');
0120     case 4              %% interior-point, w/ exact user-supplied Hessian
0121       fmc_hessian = @(x, lambda)opf_hessfcn(x, lambda, 1, om, Ybus, Yf(il,:), Yt(il,:), mpopt, il);
0122       fmoptions = optimset(fmoptions, 'Algorithm', 'interior-point', ...
0123           'Hessian', 'user-supplied', 'HessFcn', fmc_hessian);
0124     case 5              %% interior-point, w/ finite-diff Hessian
0125       fmoptions = optimset(fmoptions, 'Algorithm', 'interior-point', 'Hessian','fin-diff-grads', 'SubProblem', 'cg');
0126     case 6              %% sqp (does not use sparse matrices, not suitable for large problems)
0127       fmoptions = optimset(fmoptions, 'Algorithm', 'sqp');
0128       Af = full(Af);
0129       Afeq = full(Afeq);
0130     otherwise
0131       error('fmincopf_solver: unknown algorithm specified in ''fmincon.alg'' option');
0132   end
0133 else
0134   fmoptions = optimset(fmoptions, 'LargeScale', 'off');
0135   Af = full(Af);
0136   Afeq = full(Afeq);
0137 end
0138 % fmoptions = optimset(fmoptions, 'DerivativeCheck', 'on', 'FinDiffType', 'central', 'FunValCheck', 'on');
0139 % fmoptions = optimset(fmoptions, 'Diagnostics', 'on');
0140 
0141 %%-----  run opf  -----
0142 f_fcn = @(x)opf_costfcn(x, om);
0143 gh_fcn = @(x)opf_consfcn(x, om, Ybus, Yf(il,:), Yt(il,:), mpopt, il);
0144 [x, f, info, Output, Lambda] = ...
0145   fmincon(f_fcn, x0, Af, bf, Afeq, bfeq, LB, UB, gh_fcn, fmoptions);
0146 success = (info > 0);
0147 
0148 %% update solution data
0149 Va = x(vv.i1.Va:vv.iN.Va);
0150 Vm = x(vv.i1.Vm:vv.iN.Vm);
0151 Pg = x(vv.i1.Pg:vv.iN.Pg);
0152 Qg = x(vv.i1.Qg:vv.iN.Qg);
0153 V = Vm .* exp(1j*Va);
0154 
0155 %%-----  calculate return values  -----
0156 %% update voltages & generator outputs
0157 bus(:, VA) = Va * 180/pi;
0158 bus(:, VM) = Vm;
0159 gen(:, PG) = Pg * baseMVA;
0160 gen(:, QG) = Qg * baseMVA;
0161 gen(:, VG) = Vm(gen(:, GEN_BUS));
0162 
0163 %% compute branch flows
0164 Sf = V(branch(:, F_BUS)) .* conj(Yf * V);  %% cplx pwr at "from" bus, p.u.
0165 St = V(branch(:, T_BUS)) .* conj(Yt * V);  %% cplx pwr at "to" bus, p.u.
0166 branch(:, PF) = real(Sf) * baseMVA;
0167 branch(:, QF) = imag(Sf) * baseMVA;
0168 branch(:, PT) = real(St) * baseMVA;
0169 branch(:, QT) = imag(St) * baseMVA;
0170 
0171 %% line constraint is actually on square of limit
0172 %% so we must fix multipliers
0173 muSf = zeros(nl, 1);
0174 muSt = zeros(nl, 1);
0175 if ~isempty(il)
0176     muSf(il) = 2 * Lambda.ineqnonlin(1:nl2)       .* branch(il, RATE_A) / baseMVA;
0177     muSt(il) = 2 * Lambda.ineqnonlin((1:nl2)+nl2) .* branch(il, RATE_A) / baseMVA;
0178 end
0179 
0180 %% fix Lambdas
0181 %% (shadow prices on equality variable bounds can come back on wrong limit)
0182 kl = find(Lambda.lower < 0 & Lambda.upper == 0);
0183 Lambda.upper(kl) = -Lambda.lower(kl);
0184 Lambda.lower(kl) = 0;
0185 ku = find(Lambda.upper < 0 & Lambda.lower == 0);
0186 Lambda.lower(ku) = -Lambda.upper(ku);
0187 Lambda.upper(ku) = 0;
0188 
0189 %% update Lagrange multipliers
0190 bus(:, MU_VMAX)  = Lambda.upper(vv.i1.Vm:vv.iN.Vm);
0191 bus(:, MU_VMIN)  = Lambda.lower(vv.i1.Vm:vv.iN.Vm);
0192 gen(:, MU_PMAX)  = Lambda.upper(vv.i1.Pg:vv.iN.Pg) / baseMVA;
0193 gen(:, MU_PMIN)  = Lambda.lower(vv.i1.Pg:vv.iN.Pg) / baseMVA;
0194 gen(:, MU_QMAX)  = Lambda.upper(vv.i1.Qg:vv.iN.Qg) / baseMVA;
0195 gen(:, MU_QMIN)  = Lambda.lower(vv.i1.Qg:vv.iN.Qg) / baseMVA;
0196 bus(:, LAM_P)    = Lambda.eqnonlin(nn.i1.Pmis:nn.iN.Pmis) / baseMVA;
0197 bus(:, LAM_Q)    = Lambda.eqnonlin(nn.i1.Qmis:nn.iN.Qmis) / baseMVA;
0198 branch(:, MU_SF) = muSf / baseMVA;
0199 branch(:, MU_ST) = muSt / baseMVA;
0200 
0201 %% package up results
0202 nlnN = getN(om, 'nln');
0203 nlt = length(ilt);
0204 ngt = length(igt);
0205 nbx = length(ibx);
0206 
0207 %% extract multipliers for nonlinear constraints
0208 kl = find(Lambda.eqnonlin < 0);
0209 ku = find(Lambda.eqnonlin > 0);
0210 nl_mu_l = zeros(nlnN, 1);
0211 nl_mu_u = [zeros(2*nb, 1); muSf; muSt];
0212 nl_mu_l(kl) = -Lambda.eqnonlin(kl);
0213 nl_mu_u(ku) =  Lambda.eqnonlin(ku);
0214 
0215 %% extract multipliers for linear constraints
0216 kl = find(Lambda.eqlin < 0);
0217 ku = find(Lambda.eqlin > 0);
0218 
0219 mu_l = zeros(size(u));
0220 mu_l(ieq(kl)) = -Lambda.eqlin(kl);
0221 mu_l(igt) = Lambda.ineqlin(nlt+(1:ngt));
0222 mu_l(ibx) = Lambda.ineqlin(nlt+ngt+nbx+(1:nbx));
0223 
0224 mu_u = zeros(size(u));
0225 mu_u(ieq(ku)) = Lambda.eqlin(ku);
0226 mu_u(ilt) = Lambda.ineqlin(1:nlt);
0227 mu_u(ibx) = Lambda.ineqlin(nlt+ngt+(1:nbx));
0228 
0229 mu = struct( ...
0230   'var', struct('l', Lambda.lower, 'u', Lambda.upper), ...
0231   'nln', struct('l', nl_mu_l, 'u', nl_mu_u), ...
0232   'lin', struct('l', mu_l, 'u', mu_u) );
0233 
0234 results = mpc;
0235 [results.bus, results.branch, results.gen, ...
0236     results.om, results.x, results.mu, results.f] = ...
0237         deal(bus, branch, gen, om, x, mu, f);
0238 
0239 pimul = [ ...
0240   results.mu.nln.l - results.mu.nln.u;
0241   results.mu.lin.l - results.mu.lin.u;
0242   -ones(ny>0, 1);
0243   results.mu.var.l - results.mu.var.u;
0244 ];
0245 raw = struct('xr', x, 'pimul', pimul, 'info', info, 'output', Output);