Home > matpower5.0 > fmincopf_solver.m

fmincopf_solver

PURPOSE ^

FMINCOPF_SOLVER Solves AC optimal power flow using FMINCON.

SYNOPSIS ^

function [results, success, raw] = fmincopf_solver(om, mpopt)

DESCRIPTION ^

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.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

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