Home > matpower7.1 > lib > newtonpf.m

newtonpf

PURPOSE ^

NEWTONPF Solves power flow using full Newton's method (power/polar)

SYNOPSIS ^

function [V, converged, i] = newtonpf(Ybus, Sbus, V0, ref, pv, pq, mpopt)

DESCRIPTION ^

NEWTONPF  Solves power flow using full Newton's method (power/polar)
   [V, CONVERGED, I] = NEWTONPF(YBUS, SBUS, V0, REF, PV, PQ, MPOPT)

   Solves for bus voltages using a full Newton-Raphson method, using nodal
   power balance equations and polar coordinate representation of
   voltages, given the following inputs:
       YBUS  - full system admittance matrix (for all buses)
       SBUS  - handle to function that returns the complex bus power
               injection vector (for all buses), given the bus voltage
               magnitude vector (for all buses)
       V0    - initial vector of complex bus voltages
       REF   - bus index of reference bus (voltage ang reference & gen slack)
       PV    - vector of bus indices for PV buses
       PQ    - vector of bus indices for PQ buses
       MPOPT - (optional) MATPOWER option struct, used to set the
               termination tolerance, maximum number of iterations, and
               output options (see MPOPTION for details).

   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.

   Returns the final complex voltages, a flag which indicates whether it
   converged or not, and the number of iterations performed.

   See also RUNPF, NEWTONPF_S_CART, NEWTONPF_I_POLAR, NEWTONPF_I_CART.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [V, converged, i] = newtonpf(Ybus, Sbus, V0, ref, pv, pq, mpopt)
0002 %NEWTONPF  Solves power flow using full Newton's method (power/polar)
0003 %   [V, CONVERGED, I] = NEWTONPF(YBUS, SBUS, V0, REF, PV, PQ, MPOPT)
0004 %
0005 %   Solves for bus voltages using a full Newton-Raphson method, using nodal
0006 %   power balance equations and polar coordinate representation of
0007 %   voltages, given the following inputs:
0008 %       YBUS  - full system admittance matrix (for all buses)
0009 %       SBUS  - handle to function that returns the complex bus power
0010 %               injection vector (for all buses), given the bus voltage
0011 %               magnitude vector (for all buses)
0012 %       V0    - initial vector of complex bus voltages
0013 %       REF   - bus index of reference bus (voltage ang reference & gen slack)
0014 %       PV    - vector of bus indices for PV buses
0015 %       PQ    - vector of bus indices for PQ buses
0016 %       MPOPT - (optional) MATPOWER option struct, used to set the
0017 %               termination tolerance, maximum number of iterations, and
0018 %               output options (see MPOPTION for details).
0019 %
0020 %   The bus voltage vector contains the set point for generator
0021 %   (including ref bus) buses, and the reference angle of the swing
0022 %   bus, as well as an initial guess for remaining magnitudes and
0023 %   angles.
0024 %
0025 %   Returns the final complex voltages, a flag which indicates whether it
0026 %   converged or not, and the number of iterations performed.
0027 %
0028 %   See also RUNPF, NEWTONPF_S_CART, NEWTONPF_I_POLAR, NEWTONPF_I_CART.
0029 
0030 %   MATPOWER
0031 %   Copyright (c) 1996-2019, Power Systems Engineering Research Center (PSERC)
0032 %   by Ray Zimmerman, PSERC Cornell
0033 %
0034 %   This file is part of MATPOWER.
0035 %   Covered by the 3-clause BSD License (see LICENSE file for details).
0036 %   See https://matpower.org for more info.
0037 
0038 %% default arguments
0039 if nargin < 7
0040     mpopt = mpoption;
0041 end
0042 
0043 %% options
0044 tol         = mpopt.pf.tol;
0045 max_it      = mpopt.pf.nr.max_it;
0046 lin_solver  = mpopt.pf.nr.lin_solver;
0047 
0048 %% initialize
0049 converged = 0;
0050 i = 0;
0051 V = V0;
0052 Va = angle(V);
0053 Vm = abs(V);
0054 
0055 %% set up indexing for updating V
0056 npv = length(pv);
0057 npq = length(pq);
0058 j1 = 1;         j2 = npv;           %% j1:j2 - V angle of pv buses
0059 j3 = j2 + 1;    j4 = j2 + npq;      %% j3:j4 - V angle of pq buses
0060 j5 = j4 + 1;    j6 = j4 + npq;      %% j5:j6 - V mag of pq buses
0061 
0062 %% evaluate F(x0)
0063 mis = V .* conj(Ybus * V) - Sbus(Vm);
0064 F = [   real(mis([pv; pq]));
0065         imag(mis(pq))   ];
0066 
0067 %% check tolerance
0068 normF = norm(F, inf);
0069 if mpopt.verbose > 1
0070     fprintf('\n it    max P & Q mismatch (p.u.)');
0071     fprintf('\n----  ---------------------------');
0072     fprintf('\n%3d        %10.3e', i, normF);
0073 end
0074 if normF < tol
0075     converged = 1;
0076     if mpopt.verbose > 1
0077         fprintf('\nConverged!\n');
0078     end
0079 end
0080 
0081 %% attempt to pick fastest linear solver, if not specified
0082 if isempty(lin_solver)
0083     nx = length(F);
0084     if nx <= 10 || have_feature('octave')
0085         lin_solver = '\';       %% default \ operator
0086     else    %% MATLAB and nx > 10 or Octave and nx > 2000
0087         lin_solver = 'LU3';     %% LU decomp with 3 output args, AMD ordering
0088     end
0089 end
0090 
0091 %% do Newton iterations
0092 while (~converged && i < max_it)
0093     %% update iteration counter
0094     i = i + 1;
0095 
0096     %% evaluate Jacobian
0097     [dSbus_dVa, dSbus_dVm] = dSbus_dV(Ybus, V);
0098     [dummy, neg_dSd_dVm] = Sbus(Vm);
0099     dSbus_dVm = dSbus_dVm - neg_dSd_dVm;
0100 
0101     j11 = real(dSbus_dVa([pv; pq], [pv; pq]));
0102     j12 = real(dSbus_dVm([pv; pq], pq));
0103     j21 = imag(dSbus_dVa(pq, [pv; pq]));
0104     j22 = imag(dSbus_dVm(pq, pq));
0105 
0106     J = [   j11 j12;
0107             j21 j22;    ];
0108 
0109     %% compute update step
0110     dx = mplinsolve(J, -F, lin_solver);
0111 
0112     %% update voltage
0113     if npv
0114         Va(pv) = Va(pv) + dx(j1:j2);
0115     end
0116     if npq
0117         Va(pq) = Va(pq) + dx(j3:j4);
0118         Vm(pq) = Vm(pq) + dx(j5:j6);
0119     end
0120     V = Vm .* exp(1j * Va);
0121     Vm = abs(V);            %% update Vm and Va again in case
0122     Va = angle(V);          %% we wrapped around with a negative Vm
0123 
0124     %% evalute F(x)
0125     mis = V .* conj(Ybus * V) - Sbus(Vm);
0126     F = [   real(mis([pv; pq]));
0127             imag(mis(pq))   ];
0128 
0129     %% check for convergence
0130     normF = norm(F, inf);
0131     if mpopt.verbose > 1
0132         fprintf('\n%3d        %10.3e', i, normF);
0133     end
0134     if normF < tol
0135         converged = 1;
0136         if mpopt.verbose
0137             fprintf('\nNewton''s method power flow (power balance, polar) converged in %d iterations.\n', i);
0138         end
0139     end
0140 end
0141 
0142 if mpopt.verbose
0143     if ~converged
0144         fprintf('\nNewton''s method power flow (power balance, polar) did not converge in %d iterations.\n', i);
0145     end
0146 end

Generated on Fri 09-Oct-2020 11:21:31 by m2html © 2005