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 vector.

   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 vector.
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 %   $Id: fmincopf_solver.m,v 1.25 2011/06/16 17:46:37 cvs Exp $
0039 %   by Ray Zimmerman, PSERC Cornell
0040 %   and Carlos E. Murillo-Sanchez, PSERC Cornell & Universidad Autonoma de Manizales
0041 %   Copyright (c) 2000-2010 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 
0078 %% options
0079 verbose = mpopt(31);    %% VERBOSE
0080 
0081 %% unpack data
0082 mpc = get_mpc(om);
0083 [baseMVA, bus, gen, branch] = ...
0084     deal(mpc.baseMVA, mpc.bus, mpc.gen, mpc.branch);
0085 [vv, ll, nn] = get_idx(om);
0086 
0087 %% problem dimensions
0088 nb = size(bus, 1);          %% number of buses
0089 nl = size(branch, 1);       %% number of branches
0090 ny = getN(om, 'var', 'y');  %% number of piece-wise linear costs
0091 
0092 %% bounds on optimization vars
0093 [x0, LB, UB] = getv(om);
0094 
0095 %% linear constraints
0096 [A, l, u] = linear_constraints(om);
0097 
0098 %% split l <= A*x <= u into less than, equal to, greater than, and
0099 %% doubly-bounded sets
0100 ieq = find( abs(u-l) <= eps );          %% equality
0101 igt = find( u >=  1e10 & l > -1e10 );   %% greater than, unbounded above
0102 ilt = find( l <= -1e10 & u <  1e10 );   %% less than, unbounded below
0103 ibx = find( (abs(u-l) > eps) & (u < 1e10) & (l > -1e10) );
0104 Af  = [ A(ilt, :); -A(igt, :); A(ibx, :); -A(ibx, :) ];
0105 bf  = [ u(ilt);   -l(igt);     u(ibx);    -l(ibx)];
0106 Afeq = A(ieq, :);
0107 bfeq = u(ieq);
0108 
0109 %% build admittance matrices
0110 [Ybus, Yf, Yt] = makeYbus(baseMVA, bus, branch);
0111 
0112 %% find branches with flow limits
0113 il = find(branch(:, RATE_A) ~= 0 & branch(:, RATE_A) < 1e10);
0114 nl2 = length(il);           %% number of constrained lines
0115 
0116 %% basic optimset options needed for fmincon
0117 fmoptions = optimset('GradObj', 'on', 'GradConstr', 'on', ...
0118             'TolCon', mpopt(16), 'TolX', mpopt(17), 'TolFun', mpopt(18) );
0119 if mpopt(19) ~= 0
0120     fmoptions = optimset(fmoptions, 'MaxIter', mpopt(19), ...
0121                 'MaxFunEvals', 4 * mpopt(19));
0122 end
0123 
0124 if verbose == 0,
0125   fmoptions.Display = 'off';
0126 elseif verbose == 1
0127   fmoptions.Display = 'iter';
0128 else
0129   fmoptions.Display = 'testing';
0130 end
0131 
0132 %% select algorithm
0133 otver = ver('optim');
0134 if str2double(otver.Version(1)) < 4
0135   fmoptions = optimset(fmoptions, 'LargeScale', 'off');
0136   Af = full(Af);
0137   Afeq = full(Afeq);
0138 else
0139   if mpopt(55) == 1           %% active-set
0140     fmoptions = optimset(fmoptions, 'Algorithm', 'active-set');
0141     Af = full(Af);
0142     Afeq = full(Afeq);
0143   elseif mpopt(55) == 2       %% interior-point, w/ default 'bfgs' Hessian approx
0144     fmoptions = optimset(fmoptions, 'Algorithm', 'interior-point');
0145   elseif mpopt(55) == 3       %% interior-point, w/ 'lbfgs' Hessian approx
0146     fmoptions = optimset(fmoptions, 'Algorithm', 'interior-point', 'Hessian','lbfgs');
0147   elseif mpopt(55) == 4       %% interior-point, w/ exact user-supplied Hessian
0148     fmc_hessian = @(x, lambda)opf_hessfcn(x, lambda, 1, om, Ybus, Yf(il,:), Yt(il,:), mpopt, il);
0149     fmoptions = optimset(fmoptions, 'Algorithm', 'interior-point', ...
0150         'Hessian', 'user-supplied', 'HessFcn', fmc_hessian);
0151   elseif mpopt(55) == 5       %% interior-point, w/ finite-diff Hessian
0152     fmoptions = optimset(fmoptions, 'Algorithm', 'interior-point', 'Hessian','fin-diff-grads', 'SubProblem', 'cg');
0153   else
0154     error('fmincopf_solver: unknown algorithm specified in FMC_ALG option');
0155   end
0156 end
0157 % fmoptions = optimset(fmoptions, 'DerivativeCheck', 'on', 'FinDiffType', 'central', 'FunValCheck', 'on');
0158 % fmoptions = optimset(fmoptions, 'Diagnostics', 'on');
0159 
0160 %% try to select an interior initial point for interior point solver
0161 if str2double(otver.Version(1)) >= 4 && strcmp(optimget(fmoptions, 'Algorithm'), 'interior-point')
0162   x0 = zeros(getN(om, 'var'), 1);
0163   x0(vv.i1.Va:vv.iN.Va) = 0;
0164   x0(vv.i1.Vm:vv.iN.Vm)   = 1;
0165   x0(vv.i1.Pg:vv.iN.Pg) = (gen(:, PMIN) + gen(:, PMAX)) / 2 / baseMVA;
0166   x0(vv.i1.Qg:vv.iN.Qg) = (gen(:, QMIN) + gen(:, QMAX)) / 2 / baseMVA;
0167 end
0168 
0169 %%-----  run opf  -----
0170 f_fcn = @(x)opf_costfcn(x, om);
0171 gh_fcn = @(x)opf_consfcn(x, om, Ybus, Yf(il,:), Yt(il,:), mpopt, il);
0172 [x, f, info, Output, Lambda] = ...
0173   fmincon(f_fcn, x0, Af, bf, Afeq, bfeq, LB, UB, gh_fcn, fmoptions);
0174 success = (info > 0);
0175 
0176 %% update solution data
0177 Va = x(vv.i1.Va:vv.iN.Va);
0178 Vm = x(vv.i1.Vm:vv.iN.Vm);
0179 Pg = x(vv.i1.Pg:vv.iN.Pg);
0180 Qg = x(vv.i1.Qg:vv.iN.Qg);
0181 V = Vm .* exp(1j*Va);
0182 
0183 %%-----  calculate return values  -----
0184 %% update voltages & generator outputs
0185 bus(:, VA) = Va * 180/pi;
0186 bus(:, VM) = Vm;
0187 gen(:, PG) = Pg * baseMVA;
0188 gen(:, QG) = Qg * baseMVA;
0189 gen(:, VG) = Vm(gen(:, GEN_BUS));
0190 
0191 %% compute branch flows
0192 Sf = V(branch(:, F_BUS)) .* conj(Yf * V);  %% cplx pwr at "from" bus, p.u.
0193 St = V(branch(:, T_BUS)) .* conj(Yt * V);  %% cplx pwr at "to" bus, p.u.
0194 branch(:, PF) = real(Sf) * baseMVA;
0195 branch(:, QF) = imag(Sf) * baseMVA;
0196 branch(:, PT) = real(St) * baseMVA;
0197 branch(:, QT) = imag(St) * baseMVA;
0198 
0199 %% line constraint is actually on square of limit
0200 %% so we must fix multipliers
0201 muSf = zeros(nl, 1);
0202 muSt = zeros(nl, 1);
0203 if ~isempty(il)
0204     muSf(il) = 2 * Lambda.ineqnonlin(1:nl2)       .* branch(il, RATE_A) / baseMVA;
0205     muSt(il) = 2 * Lambda.ineqnonlin((1:nl2)+nl2) .* branch(il, RATE_A) / baseMVA;
0206 end
0207 
0208 %% update Lagrange multipliers
0209 bus(:, MU_VMAX)  = Lambda.upper(vv.i1.Vm:vv.iN.Vm);
0210 bus(:, MU_VMIN)  = Lambda.lower(vv.i1.Vm:vv.iN.Vm);
0211 gen(:, MU_PMAX)  = Lambda.upper(vv.i1.Pg:vv.iN.Pg) / baseMVA;
0212 gen(:, MU_PMIN)  = Lambda.lower(vv.i1.Pg:vv.iN.Pg) / baseMVA;
0213 gen(:, MU_QMAX)  = Lambda.upper(vv.i1.Qg:vv.iN.Qg) / baseMVA;
0214 gen(:, MU_QMIN)  = Lambda.lower(vv.i1.Qg:vv.iN.Qg) / baseMVA;
0215 bus(:, LAM_P)    = Lambda.eqnonlin(nn.i1.Pmis:nn.iN.Pmis) / baseMVA;
0216 bus(:, LAM_Q)    = Lambda.eqnonlin(nn.i1.Qmis:nn.iN.Qmis) / baseMVA;
0217 branch(:, MU_SF) = muSf / baseMVA;
0218 branch(:, MU_ST) = muSt / baseMVA;
0219 
0220 %% package up results
0221 nlnN = getN(om, 'nln');
0222 nlt = length(ilt);
0223 ngt = length(igt);
0224 nbx = length(ibx);
0225 
0226 %% extract multipliers for nonlinear constraints
0227 kl = find(Lambda.eqnonlin < 0);
0228 ku = find(Lambda.eqnonlin > 0);
0229 nl_mu_l = zeros(nlnN, 1);
0230 nl_mu_u = [zeros(2*nb, 1); muSf; muSt];
0231 nl_mu_l(kl) = -Lambda.eqnonlin(kl);
0232 nl_mu_u(ku) =  Lambda.eqnonlin(ku);
0233 
0234 %% extract multipliers for linear constraints
0235 kl = find(Lambda.eqlin < 0);
0236 ku = find(Lambda.eqlin > 0);
0237 
0238 mu_l = zeros(size(u));
0239 mu_l(ieq(kl)) = -Lambda.eqlin(kl);
0240 mu_l(igt) = Lambda.ineqlin(nlt+(1:ngt));
0241 mu_l(ibx) = Lambda.ineqlin(nlt+ngt+nbx+(1:nbx));
0242 
0243 mu_u = zeros(size(u));
0244 mu_u(ieq(ku)) = Lambda.eqlin(ku);
0245 mu_u(ilt) = Lambda.ineqlin(1:nlt);
0246 mu_u(ibx) = Lambda.ineqlin(nlt+ngt+(1:nbx));
0247 
0248 mu = struct( ...
0249   'var', struct('l', Lambda.lower, 'u', Lambda.upper), ...
0250   'nln', struct('l', nl_mu_l, 'u', nl_mu_u), ...
0251   'lin', struct('l', mu_l, 'u', mu_u) );
0252 
0253 results = mpc;
0254 [results.bus, results.branch, results.gen, ...
0255     results.om, results.x, results.mu, results.f] = ...
0256         deal(bus, branch, gen, om, x, mu, f);
0257 
0258 pimul = [ ...
0259   results.mu.nln.l - results.mu.nln.u;
0260   results.mu.lin.l - results.mu.lin.u;
0261   -ones(ny>0, 1);
0262   results.mu.var.l - results.mu.var.u;
0263 ];
0264 raw = struct('xr', x, 'pimul', pimul, 'info', info, 'output', Output);