Description of fmincopf

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);
fmincopf_solver

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SOURCE CODE
