newtonpf.rar_Jacobian_matpower newtonpf_matpower潮流_矩阵潮流计算_雅克比

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 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. % % See also RUNPF. % MATPOWER % $Id: newtonpf.m,v 1.14 2011/12/14 17:05:18 cvs Exp $ % by Ray Zimmerman, PSERC Cornell % Copyright (c) 1996-2011 by Power System Engineering Research Center (PSERC) % % This file is part of MATPOWER. % See http://www.pserc.cornell.edu/matpower/ for more info. % % MATPOWER is free software: you can redistribute it and/or modify % it under the terms of the GNU General Public License as published % by the Free Software Foundation, either version 3 of the License, % or (at your option) any later version. % % MATPOWER is distributed in the hope that it will be useful, % but WITHOUT ANY WARRANTY; without even the implied warranty of % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the % GNU General Public License for more details. % % You should have received a copy of the GNU General Public License % along with MATPOWER. If not, see <http://www.gnu.org/licenses/>. % % Additional permission under GNU GPL version 3 section 7 % % If you modify MATPOWER, or any covered work, to interface with % other modules (such as MATLAB code and MEX-files) available in a % MATLAB(R) or comparable environment containing parts covered % under other licensing terms, the licensors of MATPOWER grant % you additional permission to convey the resulting work. %% default arguments if nargin < 7 mpopt = mpoption; end %% options tol = mpopt(2); max_it = mpopt(3); verbose = mpopt(31); %% initialize 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); 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)); 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(1j * 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('\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