Description of opf

0001 function [f, df, d2f] = opf_costfcn(x, om, varargin)
0002 %OPF_COSTFCN  Evaluates objective function, gradient and Hessian for OPF.
0003 %   [F, DF, D2F] = OPF_COSTFCN(X, OM)
0004 %
0005 %   Objective function evaluation routine for AC optimal power flow,
0006 %   suitable for use with MIPS or FMINCON. Computes objective function value,
0007 %   gradient and Hessian.
0008 %
0009 %   Inputs:
0010 %     X : optimization vector
0011 %     OM : OPF model object
0012 %
0013 %   Outputs:
0014 %     F   : value of objective function
0015 %     DF  : (optional) gradient of objective function (column vector)
0016 %     D2F : (optional) Hessian of objective function (sparse matrix)
0017 %
0018 %   Examples:
0019 %       f = opf_costfcn(x, om);
0020 %       [f, df] = opf_costfcn(x, om);
0021 %       [f, df, d2f] = opf_costfcn(x, om);
0022 %
0023 %   See also OPF_CONSFCN, OPF_HESSFCN.
0024 
0025 %   MATPOWER
0026 %   $Id: opf_costfcn.m 1635 2010-04-26 19:45:26Z ray $
0027 %   by Carlos E. Murillo-Sanchez, PSERC Cornell & Universidad Autonoma de Manizales
0028 %   and Ray Zimmerman, PSERC Cornell
0029 %   Copyright (c) 1996-2010 by Power System Engineering Research Center (PSERC)
0030 %
0031 %   This file is part of MATPOWER.
0032 %   See http://www.pserc.cornell.edu/matpower/ for more info.
0033 %
0034 %   MATPOWER is free software: you can redistribute it and/or modify
0035 %   it under the terms of the GNU General Public License as published
0036 %   by the Free Software Foundation, either version 3 of the License,
0037 %   or (at your option) any later version.
0038 %
0039 %   MATPOWER is distributed in the hope that it will be useful,
0040 %   but WITHOUT ANY WARRANTY; without even the implied warranty of
0041 %   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
0042 %   GNU General Public License for more details.
0043 %
0044 %   You should have received a copy of the GNU General Public License
0045 %   along with MATPOWER. If not, see <http://www.gnu.org/licenses/>.
0046 %
0047 %   Additional permission under GNU GPL version 3 section 7
0048 %
0049 %   If you modify MATPOWER, or any covered work, to interface with
0050 %   other modules (such as MATLAB code and MEX-files) available in a
0051 %   MATLAB(R) or comparable environment containing parts covered
0052 %   under other licensing terms, the licensors of MATPOWER grant
0053 %   you additional permission to convey the resulting work.
0054 
0055 %%----- initialize -----
0056 %% define named indices into data matrices
0057 [GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, GEN_STATUS, PMAX, PMIN, ...
0058     MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN, PC1, PC2, QC1MIN, QC1MAX, ...
0059     QC2MIN, QC2MAX, RAMP_AGC, RAMP_10, RAMP_30, RAMP_Q, APF] = idx_gen;
0060 [F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, RATE_C, ...
0061     TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST, ...
0062     ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch;
0063 [PW_LINEAR, POLYNOMIAL, MODEL, STARTUP, SHUTDOWN, NCOST, COST] = idx_cost;
0064 
0065 %% unpack data
0066 mpc = get_mpc(om);
0067 [baseMVA, gen, gencost] = deal(mpc.baseMVA, mpc.gen, mpc.gencost);
0068 cp = get_cost_params(om);
0069 [N, Cw, H, dd, rh, kk, mm] = deal(cp.N, cp.Cw, cp.H, cp.dd, ...
0070                                     cp.rh, cp.kk, cp.mm);
0071 vv = get_idx(om);
0072 
0073 %% problem dimensions
0074 ng = size(gen, 1);          %% number of dispatchable injections
0075 ny = getN(om, 'var', 'y');  %% number of piece-wise linear costs
0076 nxyz = length(x);           %% total number of control vars of all types
0077 
0078 %% grab Pg & Qg
0079 Pg = x(vv.i1.Pg:vv.iN.Pg);  %% active generation in p.u.
0080 Qg = x(vv.i1.Qg:vv.iN.Qg);  %% reactive generation in p.u.
0081 
0082 %%----- evaluate objective function -----
0083 %% polynomial cost of P and Q
0084 % use totcost only on polynomial cost; in the minimization problem
0085 % formulation, pwl cost is the sum of the y variables.
0086 ipol = find(gencost(:, MODEL) == POLYNOMIAL);   %% poly MW and MVAr costs
0087 xx = [ Pg; Qg ] * baseMVA;
0088 if ~isempty(ipol)
0089   f = sum( totcost(gencost(ipol, :), xx(ipol)) );  %% cost of poly P or Q
0090 else
0091   f = 0;
0092 end
0093 
0094 %% piecewise linear cost of P and Q
0095 if ny > 0
0096   ccost = full(sparse(ones(1,ny), vv.i1.y:vv.iN.y, ones(1,ny), 1, nxyz));
0097   f = f + ccost * x;
0098 else
0099   ccost = zeros(1, nxyz);
0100 end
0101 
0102 %% generalized cost term
0103 if ~isempty(N)
0104     nw = size(N, 1);
0105     r = N * x - rh;                 %% Nx - rhat
0106     iLT = find(r < -kk);            %% below dead zone
0107     iEQ = find(r == 0 & kk == 0);   %% dead zone doesn't exist
0108     iGT = find(r > kk);             %% above dead zone
0109     iND = [iLT; iEQ; iGT];          %% rows that are Not in the Dead region
0110     iL = find(dd == 1);             %% rows using linear function
0111     iQ = find(dd == 2);             %% rows using quadratic function
0112     LL = sparse(iL, iL, 1, nw, nw);
0113     QQ = sparse(iQ, iQ, 1, nw, nw);
0114     kbar = sparse(iND, iND, [   ones(length(iLT), 1);
0115                                 zeros(length(iEQ), 1);
0116                                 -ones(length(iGT), 1)], nw, nw) * kk;
0117     rr = r + kbar;                  %% apply non-dead zone shift
0118     M = sparse(iND, iND, mm(iND), nw, nw);  %% dead zone or scale
0119     diagrr = sparse(1:nw, 1:nw, rr, nw, nw);
0120     
0121     %% linear rows multiplied by rr(i), quadratic rows by rr(i)^2
0122     w = M * (LL + QQ * diagrr) * rr;
0123 
0124     f = f + (w' * H * w) / 2 + Cw' * w;
0125 end
0126 
0127 %%----- evaluate cost gradient -----
0128 if nargout > 1
0129   %% index ranges
0130   iPg = vv.i1.Pg:vv.iN.Pg;
0131   iQg = vv.i1.Qg:vv.iN.Qg;
0132 
0133   %% polynomial cost of P and Q
0134   df_dPgQg = zeros(2*ng, 1);        %% w.r.t p.u. Pg and Qg
0135   df_dPgQg(ipol) = baseMVA * polycost(gencost(ipol, :), xx(ipol), 1);
0136   df = zeros(nxyz, 1);
0137   df(iPg) = df_dPgQg(1:ng);
0138   df(iQg) = df_dPgQg((1:ng) + ng);
0139 
0140   %% piecewise linear cost of P and Q
0141   df = df + ccost';  % As in MINOS, the linear cost row is additive wrt
0142                      % any nonlinear cost.
0143 
0144   %% generalized cost term
0145   if ~isempty(N)
0146     HwC = H * w + Cw;
0147     AA = N' * M * (LL + 2 * QQ * diagrr);
0148     df = df + AA * HwC;
0149     
0150     %% numerical check
0151     if 0    %% 1 to check, 0 to skip check
0152       ddff = zeros(size(df));
0153       step = 1e-7;
0154       tol  = 1e-3;
0155       for k = 1:length(x)
0156         xx = x;
0157         xx(k) = xx(k) + step;
0158         ddff(k) = (opf_costfcn(xx, om) - f) / step;
0159       end
0160       if max(abs(ddff - df)) > tol
0161         idx = find(abs(ddff - df) == max(abs(ddff - df)));
0162         fprintf('\nMismatch in gradient\n');
0163         fprintf('idx             df(num)         df              diff\n');
0164         fprintf('%4d%16g%16g%16g\n', [ 1:length(df); ddff'; df'; abs(ddff - df)' ]);
0165         fprintf('MAX\n');
0166         fprintf('%4d%16g%16g%16g\n', [ idx'; ddff(idx)'; df(idx)'; abs(ddff(idx) - df(idx))' ]);
0167         fprintf('\n');
0168       end
0169     end     %% numeric check
0170   end
0171 
0172   %% ---- evaluate cost Hessian -----
0173   if nargout > 2
0174     pcost = gencost(1:ng, :);
0175     if size(gencost, 1) > ng
0176         qcost = gencost(ng+1:2*ng, :);
0177     else
0178         qcost = [];
0179     end
0180     
0181     %% polynomial generator costs
0182     d2f_dPg2 = sparse(ng, 1);               %% w.r.t. p.u. Pg
0183     d2f_dQg2 = sparse(ng, 1);               %% w.r.t. p.u. Qg
0184     ipolp = find(pcost(:, MODEL) == POLYNOMIAL);
0185     d2f_dPg2(ipolp) = baseMVA^2 * polycost(pcost(ipolp, :), Pg(ipolp)*baseMVA, 2);
0186     if ~isempty(qcost)          %% Qg is not free
0187         ipolq = find(qcost(:, MODEL) == POLYNOMIAL);
0188         d2f_dQg2(ipolq) = baseMVA^2 * polycost(qcost(ipolq, :), Qg(ipolq)*baseMVA, 2);
0189     end
0190     i = [iPg iQg]';
0191     d2f = sparse(i, i, [d2f_dPg2; d2f_dQg2], nxyz, nxyz);
0192 
0193     %% generalized cost
0194     if ~isempty(N)
0195         d2f = d2f + AA * H * AA' + 2 * N' * M * QQ * sparse(1:nw, 1:nw, HwC, nw, nw) * N;
0196     end
0197   end
0198 end
opf_costfcn

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SOURCE CODE
