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%   $Id: newtonpf.m,v 1.6 2005/01/14 17:22:23 ray Exp $%   by Ray Zimmerman, PSERC Cornell%   Copyright (c) 1996-2005 by Power System Engineering Research Center (PSERC)%   See http://www.pserc.cornell.edu/matpower/ for more info.%% default argumentsif nargin < 7    mpopt = mpoption;end%% optionstol     = mpopt(2);max_it  = mpopt(3);verbose = mpopt(31);%% initializej = sqrt(-1);converged = 0;i = 0;V = V0;Va = angle(V);Vm = abs(V);%% set up indexing for updating Vnpv = length(pv);npq = length(pq);j1 = 1;         j2 = npv;           %% j1:j2 - V angle of pv busesj3 = j2 + 1;    j4 = j2 + npq;      %% j3:j4 - V angle of pq busesj5 = 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 tolerancenormF = norm(F, inf);if verbose > 1    fprintf('\n it    max P & Q mismatch (p.u.)');    fprintf('\n----  ---------------------------');    fprintf('\n%3d        %10.3e', i, normF);endif normF < tol    converged = 1;    if verbose > 1        fprintf('\nConverged!\n');    endend%% do Newton iterationswhile (~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    if npv        Va(pv) = Va(pv) + dx(j1:j2);    end    if npq        Va(pq) = Va(pq) + dx(j3:j4);        Vm(pq) = Vm(pq) + dx(j5:j6);    end    V = Vm .* exp(j * Va);    Vm = abs(V);            %% update Vm and Va again in case    Va = angle(V);          %% we wrapped around with a negative Vm    %% 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('\n牛頓拉夫遜法潮流計算經過 %d 次迭代達到收斂.\n', i);        end    endendif verbose    if ~converged        fprintf('\n牛頓拉夫遜法潮流計算經過 %d 次迭代達不到收斂.\n', i);    endend