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