Description of cpf

0001 function [max_lambda, predicted_list, corrected_list, combined_list, success, et] = cpf(casedata, loadvarloc, sigmaForLambda, sigmaForVoltage)
0002 %CPF  Run continuation power flow (CPF) solver.
0003 %   [INPUT PARAMETERS]
0004 %   loadvarloc: load variation location(in external bus numbering). Single
0005 %               bus supported so far.
0006 %   sigmaForLambda: stepsize for lambda
0007 %   sigmaForVoltage: stepsize for voltage
0008 %   [OUTPUT PARAMETERS]
0009 %   max_lambda: the lambda in p.u. w.r.t. baseMVA at (or near) the nose
0010 %               point of PV curve
0011 %   NOTE: the first column in return parameters 'predicted_list,
0012 %   corrected_list, combined_list' is bus number; the last row is lambda.
0013 %   created by Rui Bo on 2007/11/12
0014 
0015 %   MATPOWER
0016 %   $Id: cpf.m,v 1.7 2010/04/26 19:45:26 ray Exp $
0017 %   by Rui Bo
0018 %   and Ray Zimmerman, PSERC Cornell
0019 %   Copyright (c) 1996-2010 by Power System Engineering Research Center (PSERC)
0020 %   Copyright (c) 2009-2010 by Rui Bo
0021 %
0022 %   This file is part of MATPOWER.
0023 %   See http://www.pserc.cornell.edu/matpower/ for more info.
0024 %
0025 %   MATPOWER is free software: you can redistribute it and/or modify
0026 %   it under the terms of the GNU General Public License as published
0027 %   by the Free Software Foundation, either version 3 of the License,
0028 %   or (at your option) any later version.
0029 %
0030 %   MATPOWER is distributed in the hope that it will be useful,
0031 %   but WITHOUT ANY WARRANTY; without even the implied warranty of
0032 %   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
0033 %   GNU General Public License for more details.
0034 %
0035 %   You should have received a copy of the GNU General Public License
0036 %   along with MATPOWER. If not, see <http://www.gnu.org/licenses/>.
0037 %
0038 %   Additional permission under GNU GPL version 3 section 7
0039 %
0040 %   If you modify MATPOWER, or any covered work, to interface with
0041 %   other modules (such as MATLAB code and MEX-files) available in a
0042 %   MATLAB(R) or comparable environment containing parts covered
0043 %   under other licensing terms, the licensors of MATPOWER grant
0044 %   you additional permission to convey the resulting work.
0045 
0046 %% define named indices into bus, gen, branch matrices
0047 [PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ...
0048     VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus;
0049 [F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, ...
0050     RATE_C, TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST] = idx_brch;
0051 [GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, ...
0052     GEN_STATUS, PMAX, PMIN, MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN] = idx_gen;
0053 
0054 %% assign default parameters
0055 if nargin < 3
0056     sigmaForLambda = 0.1;       % stepsize for lambda
0057     sigmaForVoltage = 0.025;    % stepsize for voltage
0058 end
0059 
0060 %% options
0061 max_iter = 400;                 % depends on selection of stepsizes
0062 verbose = 1;
0063 %% instead of using condition number as criteria for switching between
0064 %% modes...
0065 % % condNumThresh_Phase1 = 0.2e-2;  % condition number shreshold for voltage prediction-correction (lambda increasing)
0066 % % condNumThresh_Phase2 = 0.2e-2;  % condition number shreshold for lambda prediction-correction
0067 % % condNumThresh_Phase3 = 0.1e-5;  % condition number shreshold for voltage prediction-correction in backward direction (lambda decreasing)
0068 %% ...we use PV curve slopes as the criteria for switching modes
0069 slopeThresh_Phase1 = 0.5;       % PV curve slope shreshold for voltage prediction-correction (with lambda increasing)
0070 slopeThresh_Phase2 = 0.3;       % PV curve slope shreshold for lambda prediction-correction
0071 
0072 %% load the case & convert to internal bus numbering
0073 [baseMVA, bus, gen, branch] = loadcase(casedata);
0074 [i2e, bus, gen, branch] = ext2int(bus, gen, branch);
0075 e2i = sparse(max(i2e), 1);
0076 e2i(i2e) = (1:size(bus, 1))';
0077 loadvarloc_i = e2i(loadvarloc);
0078 
0079 %% get bus index lists of each type of bus
0080 [ref, pv, pq] = bustypes(bus, gen);
0081 
0082 %% generator info
0083 on = find(gen(:, GEN_STATUS) > 0);      %% which generators are on?
0084 gbus = gen(on, GEN_BUS);                %% what buses are they at?
0085 
0086 %% form Ybus matrix
0087 [Ybus, Yf, Yt] = makeYbus(baseMVA, bus, branch);
0088 
0089 %% set up indexing
0090 npv = length(pv);
0091 npq = length(pq);
0092 pv_bus = ~isempty(find(pv == loadvarloc_i));
0093 
0094 %% initialize parameters
0095 % set lambda to be increasing
0096 flag_lambdaIncrease = true;  % flag indicating lambda is increasing or decreasing
0097 if bus(loadvarloc_i, PD) == 0
0098     initQPratio = 0;
0099     fprintf('\t[Warning]:\tLoad real power at bus %d is 0. Q/P ratio will be fixed at 0.\n', loadvarloc);
0100 else
0101     initQPratio = bus(loadvarloc_i, QD)./bus(loadvarloc_i, PD);
0102 end
0103 lambda0 = 0; 
0104 lambda = lambda0;
0105 Vm = ones(size(bus, 1), 1);          %% flat start
0106 Va = bus(ref(1), VA) * Vm;
0107 V  = Vm .* exp(sqrt(-1) * pi/180 * Va);
0108 V(gbus) = gen(on, VG) ./ abs(V(gbus)).* V(gbus);
0109 
0110 pointCnt = 0;
0111 
0112 %% do voltage correction (ie, power flow) to get initial voltage profile
0113 lambda_predicted = lambda;
0114 V_predicted = V;
0115 [V, lambda, success, iterNum] = cpf_correctVoltage(baseMVA, bus, gen, Ybus, V_predicted, lambda_predicted, initQPratio, loadvarloc_i);
0116 %% record data
0117 pointCnt = pointCnt + 1;
0118 corrected_list(:, pointCnt) = [V;lambda];
0119 
0120 %%------------------------------------------------
0121 %% do cpf prediction-correction iterations
0122 %%------------------------------------------------
0123 t0 = clock;
0124 %% --- Start Phase 1: voltage prediction-correction (lambda increasing)
0125 if verbose > 0
0126     fprintf('Start Phase 1: voltage prediction-correction (lambda increasing).\n');
0127 end
0128 i = 0;
0129 while i < max_iter
0130     %% update iteration counter
0131     i = i + 1;
0132     
0133     % save good data
0134     V_saved = V;
0135     lambda_saved = lambda;
0136     
0137     %% do voltage prediction to find predicted point (predicting voltage)
0138     [V_predicted, lambda_predicted, J] = cpf_predict(Ybus, ref, pv, pq, V, lambda, sigmaForLambda, 1, initQPratio, loadvarloc, flag_lambdaIncrease);
0139     
0140     %% do voltage correction to find corrected point
0141     [V, lambda, success, iterNum] = cpf_correctVoltage(baseMVA, bus, gen, Ybus, V_predicted, lambda_predicted, initQPratio, loadvarloc_i);
0142 
0143     %% calculate slope (dP/dLambda) at current point
0144     slope = abs(V(loadvarloc_i) - V_saved(loadvarloc_i))/(lambda - lambda_saved);
0145 
0146     %% instead of using condition number as criteria for switching between
0147     %% modes...
0148     %%    if rcond(J) <= condNumThresh_Phase1 | success == false % Jacobian matrix is ill-conditioned, or correction step fails
0149     %% ...we use PV curve slopes as the criteria for switching modes
0150     if abs(slope) >= slopeThresh_Phase1 | success == false % Approaching nose area of PV curve, or correction step fails
0151         % restore good data
0152         V = V_saved;
0153         lambda = lambda_saved;
0154       
0155         if verbose > 0
0156             fprintf('\t[Info]:\tApproaching nose area of PV curve, or voltage correction fails.\n');
0157         end
0158         break;
0159     else
0160         if verbose > 2
0161             fprintf('\nVm_predicted\tVm_corrected\n');
0162             [[abs(V_predicted);lambda_predicted] [abs(V);lambda]]
0163         end
0164 
0165         %% record data
0166         pointCnt = pointCnt + 1;
0167         predicted_list(:, pointCnt-1) = [V_predicted;lambda_predicted];
0168         corrected_list(:, pointCnt) = [V;lambda];
0169     end
0170 end
0171 pointCnt_Phase1 = pointCnt; % collect number of points obtained at this phase
0172 if verbose > 0
0173     fprintf('\t[Info]:\t%d data points contained in phase 1.\n', pointCnt_Phase1);
0174 end
0175 
0176 %% --- Switch to Phase 2: lambda prediction-correction (voltage decreasing)
0177 if verbose > 0
0178     fprintf('Switch to Phase 2: lambda prediction-correction (voltage decreasing).\n');
0179 end
0180 k = 0;
0181 while k < max_iter
0182     %% update iteration counter
0183     k = k + 1;
0184 
0185     % save good data
0186     V_saved = V;
0187     lambda_saved = lambda;
0188 
0189     %% do lambda prediction to find predicted point (predicting lambda)
0190     [V_predicted, lambda_predicted, J] = cpf_predict(Ybus, ref, pv, pq, V, lambda, sigmaForVoltage, 2, initQPratio, loadvarloc);
0191     %% do lambda correction to find corrected point
0192     Vm_assigned = abs(V_predicted);
0193     [V, lambda, success, iterNum] = cpf_correctLambda(baseMVA, bus, gen, Ybus, Vm_assigned, V_predicted, lambda_predicted, initQPratio, loadvarloc, ref, pv, pq);
0194 
0195     %% calculate slope (dP/dLambda) at current point
0196     slope = abs(V(loadvarloc_i) - V_saved(loadvarloc_i))/(lambda - lambda_saved);
0197 
0198     %% instead of using condition number as criteria for switching between
0199     %% modes...
0200     %%    if rcond(J) >= condNumThresh_Phase2 | success == false % Jacobian matrix is good-conditioned, or correction step fails
0201     %% ...we use PV curve slopes as the criteria for switching modes
0202     if abs(slope) <= slopeThresh_Phase2 | success == false % Leaving nose area of PV curve, or correction step fails
0203         % restore good data
0204         V = V_saved;
0205         lambda = lambda_saved;
0206 
0207         %% ---change to voltage prediction-correction (lambda decreasing)
0208         if verbose > 0
0209             fprintf('\t[Info]:\tLeaving nose area of PV curve, or lambda correction fails.\n');
0210         end
0211         break;
0212     else
0213         if verbose > 2
0214             fprintf('\nVm_predicted\tVm_corrected\n');
0215             [[abs(V_predicted);lambda_predicted] [abs(V);lambda]]
0216         end
0217 
0218         %% record data
0219         pointCnt = pointCnt + 1;
0220         predicted_list(:, pointCnt-1) = [V_predicted;lambda_predicted];
0221         corrected_list(:, pointCnt) = [V;lambda];
0222     end
0223 end
0224 pointCnt_Phase2 = pointCnt - pointCnt_Phase1; % collect number of points obtained at this phase
0225 if verbose > 0
0226     fprintf('\t[Info]:\t%d data points contained in phase 2.\n', pointCnt_Phase2);
0227 end
0228 
0229 %% --- Switch to Phase 3: voltage prediction-correction (lambda decreasing)
0230 if verbose > 0
0231     fprintf('Switch to Phase 3: voltage prediction-correction (lambda decreasing).\n');
0232 end
0233 % set lambda to be decreasing
0234 flag_lambdaIncrease = false; 
0235 i = 0;
0236 while i < max_iter
0237     %% update iteration counter
0238     i = i + 1;
0239     
0240     %% do voltage prediction to find predicted point (predicting voltage)
0241     [V_predicted, lambda_predicted, J] = cpf_predict(Ybus, ref, pv, pq, V, lambda, sigmaForLambda, 1, initQPratio, loadvarloc, flag_lambdaIncrease);
0242     
0243     %% do voltage correction to find corrected point
0244     [V, lambda, success, iterNum] = cpf_correctVoltage(baseMVA, bus, gen, Ybus, V_predicted, lambda_predicted, initQPratio, loadvarloc_i);
0245 
0246     %% calculate slope (dP/dLambda) at current point
0247     slope = abs(V(loadvarloc_i) - V_saved(loadvarloc_i))/(lambda - lambda_saved);
0248 
0249     if lambda < 0 % lambda is less than 0, then stops CPF simulation
0250         if verbose > 0
0251             fprintf('\t[Info]:\tlambda is less than 0.\n\t\t\tCPF finished.\n');
0252         end
0253         break;
0254     end
0255     
0256     %% instead of using condition number as criteria for switching between
0257     %% modes...
0258     %%    if rcond(J) <= condNumThresh_Phase3 | success == false % Jacobian matrix is ill-conditioned, or correction step fails
0259     %% ...we use PV curve slopes as the criteria for switching modes
0260     if success == false % voltage correction step fails.
0261         if verbose > 0
0262             fprintf('\t[Info]:\tVoltage correction step fails..\n');
0263         end
0264         break;
0265     else
0266         if verbose > 2
0267             fprintf('\nVm_predicted\tVm_corrected\n');
0268             [[abs(V_predicted);lambda_predicted] [abs(V);lambda]]
0269         end
0270 
0271         %% record data
0272         pointCnt = pointCnt + 1;
0273         predicted_list(:, pointCnt-1) = [V_predicted;lambda_predicted];
0274         corrected_list(:, pointCnt) = [V;lambda];
0275     end
0276 end
0277 pointCnt_Phase3 = pointCnt - pointCnt_Phase2 - pointCnt_Phase1; % collect number of points obtained at this phase
0278 if verbose > 0
0279     fprintf('\t[Info]:\t%d data points contained in phase 3.\n', pointCnt_Phase3);
0280 end
0281 
0282 et = etime(clock, t0);
0283 
0284 %% combine the prediction and correction data in the sequence of appearance
0285 % NOTE: number of prediction data is one less than that of correction data
0286 predictedCnt = size(predicted_list, 2);
0287 combined_list(:, 1) = corrected_list(:, 1);
0288 for i = 1:predictedCnt
0289     combined_list(:, 2*i)     = predicted_list(:, i);
0290     combined_list(:, 2*i+1)   = corrected_list(:, i+1);
0291 end
0292 
0293 %% convert back to original bus numbering & print results
0294 [bus, gen, branch] = int2ext(i2e, bus, gen, branch);
0295 
0296 %% add bus number as the first column to the prediction, correction, and combined data list
0297 nb          = size(bus, 1);
0298 max_lambda  = max(corrected_list(nb+1, :));
0299 predicted_list = [[bus(:, BUS_I);0] predicted_list];
0300 corrected_list = [[bus(:, BUS_I);0] corrected_list];
0301 combined_list  = [[bus(:, BUS_I);0] combined_list];
0302 
0303 if verbose > 1
0304     Vm_corrected = abs(corrected_list);
0305     Vm_predicted = abs(predicted_list);
0306     Vm_combined  = abs(combined_list);
0307     Vm_corrected
0308     Vm_predicted
0309     Vm_combined
0310     pointCnt_Phase1
0311     pointCnt_Phase2
0312     pointCnt_Phase3
0313     pointCnt
0314 end
cpf

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SOURCE CODE
