Description of fmincopf6

0001 function [results, success, raw] = fmincopf6_solver(om, mpopt)
0002 %------------------------------  deprecated  ------------------------------
0003 %   MATLAB 6.x support to be removed in a future version.
0004 %--------------------------------------------------------------------------
0005 %FMINCOPF6_SOLVER  Solves AC optimal power flow using FMINCON (for MATLAB 6.x).
0006 %
0007 %   [RESULTS, SUCCESS, RAW] = FMINCOPF6_SOLVER(OM, MPOPT)
0008 %
0009 %   Inputs are an OPF model object and a MATPOWER options vector.
0010 %
0011 %   Outputs are a RESULTS struct, SUCCESS flag and RAW output struct.
0012 %
0013 %   RESULTS is a MATPOWER case struct (mpc) with the usual baseMVA, bus
0014 %   branch, gen, gencost fields, along with the following additional
0015 %   fields:
0016 %       .order      see 'help ext2int' for details of this field
0017 %       .x          final value of optimization variables (internal order)
0018 %       .f          final objective function value
0019 %       .mu         shadow prices on ...
0020 %           .var
0021 %               .l  lower bounds on variables
0022 %               .u  upper bounds on variables
0023 %           .nln
0024 %               .l  lower bounds on nonlinear constraints
0025 %               .u  upper bounds on nonlinear constraints
0026 %           .lin
0027 %               .l  lower bounds on linear constraints
0028 %               .u  upper bounds on linear constraints
0029 %
0030 %   SUCCESS     1 if solver converged successfully, 0 otherwise
0031 %
0032 %   RAW         raw output in form returned by MINOS
0033 %       .xr     final value of optimization variables
0034 %       .pimul  constraint multipliers
0035 %       .info   solver specific termination code
0036 %       .output solver specific output information
0037 %
0038 %   See also OPF, FMINCON.
0039 
0040 %   MATPOWER
0041 %   $Id: fmincopf6_solver.m,v 1.26 2010/06/09 14:56:58 ray Exp $
0042 %   by Ray Zimmerman, PSERC Cornell
0043 %   and Carlos E. Murillo-Sanchez, PSERC Cornell & Universidad Autonoma de Manizales
0044 %   Copyright (c) 2000-2010 by Power System Engineering Research Center (PSERC)
0045 %
0046 %   This file is part of MATPOWER.
0047 %   See http://www.pserc.cornell.edu/matpower/ for more info.
0048 %
0049 %   MATPOWER is free software: you can redistribute it and/or modify
0050 %   it under the terms of the GNU General Public License as published
0051 %   by the Free Software Foundation, either version 3 of the License,
0052 %   or (at your option) any later version.
0053 %
0054 %   MATPOWER is distributed in the hope that it will be useful,
0055 %   but WITHOUT ANY WARRANTY; without even the implied warranty of
0056 %   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
0057 %   GNU General Public License for more details.
0058 %
0059 %   You should have received a copy of the GNU General Public License
0060 %   along with MATPOWER. If not, see <http://www.gnu.org/licenses/>.
0061 %
0062 %   Additional permission under GNU GPL version 3 section 7
0063 %
0064 %   If you modify MATPOWER, or any covered work, to interface with
0065 %   other modules (such as MATLAB code and MEX-files) available in a
0066 %   MATLAB(R) or comparable environment containing parts covered
0067 %   under other licensing terms, the licensors of MATPOWER grant
0068 %   you additional permission to convey the resulting work.
0069 
0070 %%----- initialization -----
0071 %% define named indices into data matrices
0072 [PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ...
0073     VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus;
0074 [GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, GEN_STATUS, PMAX, PMIN, ...
0075     MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN, PC1, PC2, QC1MIN, QC1MAX, ...
0076     QC2MIN, QC2MAX, RAMP_AGC, RAMP_10, RAMP_30, RAMP_Q, APF] = idx_gen;
0077 [F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, RATE_C, ...
0078     TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST, ...
0079     ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch;
0080 
0081 %% options
0082 verbose = mpopt(31);    %% VERBOSE
0083 
0084 %% unpack data
0085 mpc = get_mpc(om);
0086 [baseMVA, bus, gen, branch] = ...
0087     deal(mpc.baseMVA, mpc.bus, mpc.gen, mpc.branch);
0088 [vv, ll, nn] = get_idx(om);
0089 
0090 %% problem dimensions
0091 nb = size(bus, 1);          %% number of buses
0092 nl = size(branch, 1);       %% number of branches
0093 ny = getN(om, 'var', 'y');  %% number of piece-wise linear costs
0094 
0095 %% bounds on optimization vars
0096 [x0, LB, UB] = getv(om);
0097 
0098 %% linear constraints
0099 [A, l, u] = linear_constraints(om);
0100 
0101 %% so, can we do anything good about lambda initialization?
0102 if all(bus(:, LAM_P) == 0)
0103   bus(:, LAM_P) = (10)*ones(nb, 1);
0104 end
0105 
0106 %% split l <= A*x <= u into less than, equal to, greater than, and
0107 %% doubly-bounded sets
0108 ieq = find( abs(u-l) <= eps );          %% equality
0109 igt = find( u >=  1e10 & l > -1e10 );   %% greater than, unbounded above
0110 ilt = find( l <= -1e10 & u <  1e10 );   %% less than, unbounded below
0111 ibx = find( (abs(u-l) > eps) & (u < 1e10) & (l > -1e10) );
0112 Af  = [ A(ilt, :); -A(igt, :); A(ibx, :); -A(ibx, :) ];
0113 bf  = [ u(ilt);   -l(igt);     u(ibx);    -l(ibx)];
0114 Afeq = A(ieq, :);
0115 bfeq = u(ieq);
0116 
0117 %% build admittance matrices
0118 [Ybus, Yf, Yt] = makeYbus(baseMVA, bus, branch);
0119 
0120 %% tolerances
0121 if mpopt(19) == 0           %% CONSTR_MAX_IT
0122   mpopt(19) = 150 + 2*nb;
0123 end
0124 
0125 %% find branches with flow limits
0126 il = find(branch(:, RATE_A) ~= 0 & branch(:, RATE_A) < 1e10);
0127 nl2 = length(il);           %% number of constrained lines
0128 
0129 %% basic optimset options needed for fmincon
0130 fmoptions = optimset('GradObj', 'on', 'GradConstr', 'on', ...
0131             'MaxIter', mpopt(19), 'TolCon', mpopt(16), ...
0132             'TolX', mpopt(17), 'TolFun', mpopt(18) );
0133 fmoptions.MaxFunEvals = 4 * fmoptions.MaxIter;
0134 if verbose == 0,
0135   fmoptions.Display = 'off';
0136 elseif verbose == 1
0137   fmoptions.Display = 'iter';
0138 else
0139   fmoptions.Display = 'testing';
0140 end
0141 
0142 %% select algorithm
0143 fmoptions = optimset(fmoptions, 'LargeScale', 'off');
0144 Af = full(Af);
0145 Afeq = full(Afeq);
0146 
0147 %%-----  run opf  -----
0148 f_fcn = @opf_costfcn;
0149 gh_fcn = @opf_consfcn;
0150 mpopt(51) = 0;  %% hijack SPARSE_QP to tell OPF_CONSFCN to use full matrices
0151 [x, f, info, Output, Lambda] = ...
0152   fmincon(f_fcn, x0, Af, bf, Afeq, bfeq, LB, UB, gh_fcn, fmoptions, ...
0153          om, Ybus, Yf(il,:), Yt(il,:), mpopt, il);
0154 success = (info > 0);
0155 
0156 %% update solution data
0157 Va = x(vv.i1.Va:vv.iN.Va);
0158 Vm = x(vv.i1.Vm:vv.iN.Vm);
0159 Pg = x(vv.i1.Pg:vv.iN.Pg);
0160 Qg = x(vv.i1.Qg:vv.iN.Qg);
0161 V = Vm .* exp(1j*Va);
0162 
0163 %%-----  calculate return values  -----
0164 %% update voltages & generator outputs
0165 bus(:, VA) = Va * 180/pi;
0166 bus(:, VM) = Vm;
0167 gen(:, PG) = Pg * baseMVA;
0168 gen(:, QG) = Qg * baseMVA;
0169 gen(:, VG) = Vm(gen(:, GEN_BUS));
0170 
0171 %% compute branch flows
0172 Sf = V(branch(:, F_BUS)) .* conj(Yf * V);  %% cplx pwr at "from" bus, p.u.
0173 St = V(branch(:, T_BUS)) .* conj(Yt * V);  %% cplx pwr at "to" bus, p.u.
0174 branch(:, PF) = real(Sf) * baseMVA;
0175 branch(:, QF) = imag(Sf) * baseMVA;
0176 branch(:, PT) = real(St) * baseMVA;
0177 branch(:, QT) = imag(St) * baseMVA;
0178 
0179 %% line constraint is actually on square of limit
0180 %% so we must fix multipliers
0181 muSf = zeros(nl, 1);
0182 muSt = zeros(nl, 1);
0183 if ~isempty(il)
0184     muSf(il) = 2 * Lambda.ineqnonlin(1:nl2)       .* branch(il, RATE_A) / baseMVA;
0185     muSt(il) = 2 * Lambda.ineqnonlin((1:nl2)+nl2) .* branch(il, RATE_A) / baseMVA;
0186 end
0187 
0188 %% update Lagrange multipliers
0189 bus(:, MU_VMAX)  = Lambda.upper(vv.i1.Vm:vv.iN.Vm);
0190 bus(:, MU_VMIN)  = Lambda.lower(vv.i1.Vm:vv.iN.Vm);
0191 gen(:, MU_PMAX)  = Lambda.upper(vv.i1.Pg:vv.iN.Pg) / baseMVA;
0192 gen(:, MU_PMIN)  = Lambda.lower(vv.i1.Pg:vv.iN.Pg) / baseMVA;
0193 gen(:, MU_QMAX)  = Lambda.upper(vv.i1.Qg:vv.iN.Qg) / baseMVA;
0194 gen(:, MU_QMIN)  = Lambda.lower(vv.i1.Qg:vv.iN.Qg) / baseMVA;
0195 bus(:, LAM_P)    = Lambda.eqnonlin(nn.i1.Pmis:nn.iN.Pmis) / baseMVA;
0196 bus(:, LAM_Q)    = Lambda.eqnonlin(nn.i1.Qmis:nn.iN.Qmis) / baseMVA;
0197 branch(:, MU_SF) = muSf / baseMVA;
0198 branch(:, MU_ST) = muSt / baseMVA;
0199 
0200 %% package up results
0201 nlnN = getN(om, 'nln');
0202 nlt = length(ilt);
0203 ngt = length(igt);
0204 nbx = length(ibx);
0205 
0206 %% extract multipliers for nonlinear constraints
0207 kl = find(Lambda.eqnonlin < 0);
0208 ku = find(Lambda.eqnonlin > 0);
0209 nl_mu_l = zeros(nlnN, 1);
0210 nl_mu_u = [zeros(2*nb, 1); muSf; muSt];
0211 nl_mu_l(kl) = -Lambda.eqnonlin(kl);
0212 nl_mu_u(ku) =  Lambda.eqnonlin(ku);
0213 
0214 %% extract multipliers for linear constraints
0215 kl = find(Lambda.eqlin < 0);
0216 ku = find(Lambda.eqlin > 0);
0217 
0218 mu_l = zeros(size(u));
0219 mu_l(ieq(kl)) = -Lambda.eqlin(kl);
0220 mu_l(igt) = Lambda.ineqlin(nlt+(1:ngt));
0221 mu_l(ibx) = Lambda.ineqlin(nlt+ngt+nbx+(1:nbx));
0222 
0223 mu_u = zeros(size(u));
0224 mu_u(ieq(ku)) = Lambda.eqlin(ku);
0225 mu_u(ilt) = Lambda.ineqlin(1:nlt);
0226 mu_u(ibx) = Lambda.ineqlin(nlt+ngt+(1:nbx));
0227 
0228 mu = struct( ...
0229   'var', struct('l', Lambda.lower, 'u', Lambda.upper), ...
0230   'nln', struct('l', nl_mu_l, 'u', nl_mu_u), ...
0231   'lin', struct('l', mu_l, 'u', mu_u) );
0232 
0233 results = mpc;
0234 [results.bus, results.branch, results.gen, ...
0235     results.om, results.x, results.mu, results.f] = ...
0236         deal(bus, branch, gen, om, x, mu, f);
0237 
0238 pimul = [ ...
0239   results.mu.nln.l - results.mu.nln.u;
0240   results.mu.lin.l - results.mu.lin.u;
0241   -ones(ny>0, 1);
0242   results.mu.var.l - results.mu.var.u;
0243 ];
0244 raw = struct('xr', x, 'pimul', pimul, 'info', info, 'output', Output);
fmincopf6_solver

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SOURCE CODE
