Home > matpower4.1 > extras > cpf > cpf_correctLambda.m

cpf_correctLambda

PURPOSE ^

CPF_CORRECTLAMBDA Correct lambda in correction step near load point.

SYNOPSIS ^

function [V, lambda, converged, iterNum] = cpf_correctLambda(baseMVA, bus, gen, Ybus, Vm_assigned, V_predicted, lambda_predicted, initQPratio, loadvarloc, ref, pv, pq)

DESCRIPTION ^

CPF_CORRECTLAMBDA  Correct lambda in correction step near load point.
   function: correct lambda(ie, real power of load) in cpf correction step
   near the nose point. Use NR's method to solve the nonlinear equations
   [INPUT PARAMETERS]
   loadvarloc: (in internal bus numbering)
   created by Rui Bo on 2007/11/12

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [V, lambda, converged, iterNum] = cpf_correctLambda(baseMVA, bus, gen, Ybus, Vm_assigned, V_predicted, lambda_predicted, initQPratio, loadvarloc, ref, pv, pq)
0002 %CPF_CORRECTLAMBDA  Correct lambda in correction step near load point.
0003 %   function: correct lambda(ie, real power of load) in cpf correction step
0004 %   near the nose point. Use NR's method to solve the nonlinear equations
0005 %   [INPUT PARAMETERS]
0006 %   loadvarloc: (in internal bus numbering)
0007 %   created by Rui Bo on 2007/11/12
0008 
0009 %   MATPOWER
0010 %   $Id: cpf_correctLambda.m,v 1.4 2010/04/26 19:45:26 ray Exp $
0011 %   by Rui Bo
0012 %   and Ray Zimmerman, PSERC Cornell
0013 %   Copyright (c) 1996-2010 by Power System Engineering Research Center (PSERC)
0014 %   Copyright (c) 2009-2010 by Rui Bo
0015 %
0016 %   This file is part of MATPOWER.
0017 %   See http://www.pserc.cornell.edu/matpower/ for more info.
0018 %
0019 %   MATPOWER is free software: you can redistribute it and/or modify
0020 %   it under the terms of the GNU General Public License as published
0021 %   by the Free Software Foundation, either version 3 of the License,
0022 %   or (at your option) any later version.
0023 %
0024 %   MATPOWER is distributed in the hope that it will be useful,
0025 %   but WITHOUT ANY WARRANTY; without even the implied warranty of
0026 %   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
0027 %   GNU General Public License for more details.
0028 %
0029 %   You should have received a copy of the GNU General Public License
0030 %   along with MATPOWER. If not, see <http://www.gnu.org/licenses/>.
0031 %
0032 %   Additional permission under GNU GPL version 3 section 7
0033 %
0034 %   If you modify MATPOWER, or any covered work, to interface with
0035 %   other modules (such as MATLAB code and MEX-files) available in a
0036 %   MATLAB(R) or comparable environment containing parts covered
0037 %   under other licensing terms, the licensors of MATPOWER grant
0038 %   you additional permission to convey the resulting work.
0039 
0040 %% define named indices into bus, gen, branch matrices
0041 [PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ...
0042     VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus;
0043 
0044 %% options
0045 tol     = 1e-5; % mpopt(2);
0046 max_it  = 100;  % mpopt(3);
0047 verbose = 0;    % mpopt(31);
0048 
0049 %% initialize
0050 j = sqrt(-1);
0051 converged = 0;
0052 i = 0;
0053 V = V_predicted;
0054 lambda = lambda_predicted;
0055 Va = angle(V);
0056 Vm = abs(V);
0057 
0058 %% set up indexing for updating V
0059 npv = length(pv);
0060 npq = length(pq);
0061 j1 = 1;         j2 = npv;           %% j1:j2 - V angle of pv buses
0062 j3 = j2 + 1;    j4 = j2 + npq;      %% j3:j4 - V angle of pq buses
0063 j5 = j4 + 1;    j6 = j4 + npq;      %% j5:j6 - V mag of pq buses
0064 j7 = j6 + 1;                        %% j7 - lambda
0065 
0066 pv_bus = ~isempty(find(pv == loadvarloc));
0067 
0068 %% set load as lambda indicates
0069 bus(loadvarloc, PD) = lambda*baseMVA;
0070 bus(loadvarloc, QD) = lambda*baseMVA*initQPratio;
0071 
0072 %% compute complex bus power injections (generation - load)
0073 SbusInj = makeSbus(baseMVA, bus, gen);
0074 
0075 %% evalute F(x0)
0076 mis = V .* conj(Ybus * V) - SbusInj;
0077 mis = -mis; % NOTE: use reverse mismatch and correspondingly use '(-)Jacobian" obtained from dSbus_dV
0078 F = [   real(mis(pv));
0079         real(mis(pq));
0080         imag(mis(pq));
0081         abs(V(loadvarloc)) - Vm_assigned(loadvarloc);   ];
0082     
0083 %% do Newton iterations
0084 while (~converged & i < max_it)
0085     %% update iteration counter
0086     i = i + 1;
0087 
0088     %% evaluate Jacobian
0089     [dSbus_dVm, dSbus_dVa] = dSbus_dV(Ybus, V);
0090 
0091     j11 = real(dSbus_dVa([pv; pq], [pv; pq]));
0092     j12 = real(dSbus_dVm([pv; pq], pq));
0093     j21 = imag(dSbus_dVa(pq, [pv; pq]));
0094     j22 = imag(dSbus_dVm(pq, pq));
0095 
0096     J = [   j11 j12;
0097             j21 j22;    ];
0098 
0099     %% evaluate dDeltaP/dLambda, dDeltaQ/dLambda, dDeltaVm/dLambda,
0100     %% dDeltaVm/dVa, dDeltaVm/dVm
0101     dDeltaP_dLambda = zeros(npv+npq, 1);
0102     dDeltaQ_dLambda = zeros(npq, 1);
0103     if pv_bus % pv bus
0104         dDeltaP_dLambda(find(pv == loadvarloc)) = -1;                         % corresponding to deltaP
0105     else % pq bus
0106         dDeltaP_dLambda(npv + find(pq == loadvarloc)) = -1;                   % corresponding to deltaP
0107         dDeltaQ_dLambda(find(pq == loadvarloc)) = -initQPratio;   % corresponding to deltaQ
0108     end
0109     dDeltaVm_dLambda = zeros(1, 1);
0110     dDeltaVm_dVa = zeros(1, npv+npq);
0111     dDeltaVm_dVm = zeros(1, npq);
0112     dDeltaVm_dVm(1, find(pq == loadvarloc)) = -1;
0113 
0114     %% form augmented Jacobian
0115     J12 = [dDeltaP_dLambda; 
0116            dDeltaQ_dLambda];
0117     J21 = [dDeltaVm_dVa   dDeltaVm_dVm];
0118     J22 = dDeltaVm_dLambda;
0119     augJ = [    -J   J12;
0120                 J21 J22;    ];
0121 
0122     %% compute update step
0123     dx = -(augJ \ F);
0124 
0125     %% update voltage.
0126     % NOTE: voltage magnitude of pv buses, voltage magnitude
0127     % and angle of reference bus are not updated, so they keep as constants
0128     % (ie, the value as in the initial guess)
0129     if npv
0130         Va(pv) = Va(pv) + dx(j1:j2);
0131     end
0132     if npq
0133         Va(pq) = Va(pq) + dx(j3:j4);
0134         Vm(pq) = Vm(pq) + dx(j5:j6);
0135     end
0136     lambda = lambda + dx(j7);
0137 
0138     V = Vm .* exp(j * Va); % NOTE: angle is in radians in pf solver, but in degree in case data
0139     Vm = abs(V);            %% update Vm and Va again in case
0140     Va = angle(V);          %% we wrapped around with a negative Vm
0141 
0142     %% set load as lambda indicates
0143     bus(loadvarloc, PD) = lambda*baseMVA;
0144     bus(loadvarloc, QD) = lambda*baseMVA*initQPratio;
0145 
0146     %% compute complex bus power injections (generation - load)
0147     SbusInj = makeSbus(baseMVA, bus, gen);
0148     
0149     %% evalute F(x)
0150     mis = V .* conj(Ybus * V) - SbusInj;
0151     mis = -mis;
0152     F = [   real(mis(pv));
0153             real(mis(pq));
0154             imag(mis(pq));
0155             abs(V(loadvarloc)) - Vm_assigned(loadvarloc);   ];
0156 
0157     %% check for convergence
0158     normF = norm(F, inf);
0159     if verbose > 1
0160         fprintf('\niteration [%3d]\t\tnorm of mismatch: %10.3e', i, normF);
0161     end
0162     if normF < tol
0163         converged = 1;
0164     end
0165 end
0166 
0167 iterNum = i;

Generated on Mon 26-Jan-2015 15:00:13 by m2html © 2005