gusucode.com > matpower工具箱源码程序 > matpower工具箱源码程序/MP2_0/newtonpf.m
function [V, converged, i] = newtonpf(Ybus, Sbus, V0, ref, pv, pq, mpopt) %NEWTONPF Solves the power flow using a full Newton's method. % [V, converged, i] = newtonpf(Ybus, Sbus, V0, ref, pv, pq, mpopt) % solves for bus voltages given the full system admittance matrix (for % all buses), the complex bus power injection vector (for all buses), % the initial vector of complex bus voltages, and column vectors with % the lists of bus indices for the swing bus, PV buses, and PQ buses, % respectively. The bus voltage vector contains the set point for % generator (including ref bus) buses, and the reference angle of the % swing bus, as well as an initial guess for remaining magnitudes and % angles. mpopt is a MATPOWER options vector which can be used to % set the termination tolerance, maximum number of iterations, and % output options (see 'help mpoption' for details). Uses default % options if this parameter is not given. Returns the final complex % voltages, a flag which indicates whether it converged or not, and % the number of iterations performed. % MATPOWER Version 2.0 % by Ray Zimmerman, PSERC Cornell 12/24/97 % Copyright (c) 1996, 1997 by Power System Engineering Research Center (PSERC) % See http://www.pserc.cornell.edu/ for more info. %% default arguments if nargin < 7 mpopt = mpoption; end %% options tol = mpopt(2); max_it = mpopt(3); verbose = mpopt(31); %% initialize j = sqrt(-1); converged = 0; i = 0; V = V0; Va = angle(V); Vm = abs(V); %% set up indexing for updating V npv = length(pv); npq = length(pq); j1 = 1; j2 = npv; %% j1:j2 - V angle of pv buses j3 = j2 + 1; j4 = j2 + npq; %% j3:j4 - V angle of pq buses j5 = j4 + 1; j6 = j4 + npq; %% j5:j6 - V mag of pq buses %% evaluate F(x0) mis = V .* conj(Ybus * V) - Sbus; F = [ real(mis([pv; pq])); imag(mis(pq)) ]; %% check tolerance normF = norm(F, inf); if verbose > 1 fprintf('\n it max P & Q mismatch (p.u.)'); fprintf('\n---- ---------------------------'); fprintf('\n%3d %10.3e', i, normF); end if normF < tol converged = 1; if verbose > 1 fprintf('\nConverged!\n'); end end %% do Newton iterations while (~converged & i < max_it) %% update iteration counter i = i + 1; %% evaluate Jacobian [dSbus_dVm, dSbus_dVa] = dSbus_dV(Ybus, V); %% selecting a subset of rows of a large sparse matrix is very slow %% in Matlab 5 (but not Matlab 4 ... go figure), but selecting a %% subset of the columns is fast, and so is transposing, so instead %% of doing this ... % j11 = real(dSbus_dVa([pv; pq], [pv; pq])); % j12 = real(dSbus_dVm([pv; pq], pq)); % j21 = imag(dSbus_dVa(pq, [pv; pq])); % j22 = imag(dSbus_dVm(pq, pq)); %% ... we do the equivalent thing using %% a temporary matrix and transposing temp = real(dSbus_dVa(:, [pv; pq]))'; j11 = temp(:, [pv; pq])'; temp = real(dSbus_dVm(:, pq))'; j12 = temp(:, [pv; pq])'; temp = imag(dSbus_dVa(:, [pv; pq]))'; j21 = temp(:, pq)'; temp = imag(dSbus_dVm(:, pq))'; j22 = temp(:, pq)'; J = [ j11 j12; j21 j22; ]; %% compute update step dx = -(J \ F); %% update voltage Va(pv) = Va(pv) + dx(j1:j2); Va(pq) = Va(pq) + dx(j3:j4); Vm(pq) = Vm(pq) + dx(j5:j6); V = Vm .* exp(j * Va); %% evalute F(x) mis = V .* conj(Ybus * V) - Sbus; F = [ real(mis(pv)); real(mis(pq)); imag(mis(pq)) ]; %% check for convergence normF = norm(F, inf); if verbose > 1 fprintf('\n%3d %10.3e', i, normF); end if normF < tol converged = 1; if verbose fprintf('\nNewton''s method power flow converged in %d iterations.\n', i); end end end if verbose if ~converged fprintf('\nNewton''s method power did not converge in %d iterations.\n', i); end end