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