Home > matpower7.1 > lib > newtonpf_I_polar.m

newtonpf_I_polar

PURPOSE ^

NEWTONPF_I_POLAR Solves power flow using full Newton's method (current/cartesian)

SYNOPSIS ^

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

DESCRIPTION ^

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

   Solves for bus voltages using a full Newton-Raphson method, using nodal
   current 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, NEWTONPF_S_CART, NEWTONPF_I_CART.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [V, converged, i] = newtonpf_I_polar(Ybus, Sbus, V0, ref, pv, pq, mpopt)
0002 %NEWTONPF_I_POLAR  Solves power flow using full Newton's method (current/cartesian)
0003 %   [V, CONVERGED, I] = NEWTONPF_I_POLAR(YBUS, SBUS, V0, REF, PV, PQ, MPOPT)
0004 %
0005 %   Solves for bus voltages using a full Newton-Raphson method, using nodal
0006 %   current 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, NEWTONPF_S_CART, NEWTONPF_I_CART.
0029 
0030 %   MATPOWER
0031 %   Copyright (c) 1996-2019, Power Systems Engineering Research Center (PSERC)
0032 %   by Ray Zimmerman, PSERC Cornell
0033 %   and Baljinnyam Sereeter, Delft University of Technology
0034 %
0035 %   This file is part of MATPOWER.
0036 %   Covered by the 3-clause BSD License (see LICENSE file for details).
0037 %   See https://matpower.org for more info.
0038 
0039 %% default arguments
0040 if nargin < 7
0041     mpopt = mpoption;
0042 end
0043 
0044 %% options
0045 tol         = mpopt.pf.tol;
0046 max_it      = mpopt.pf.nr.max_it;
0047 lin_solver  = mpopt.pf.nr.lin_solver;
0048 
0049 %% initialize
0050 converged = 0;
0051 i = 0;
0052 V = V0;
0053 Va = angle(V);
0054 Vm = abs(V);
0055 n = length(V0);
0056 
0057 %% set up indexing for updating V
0058 npv = length(pv);
0059 npq = length(pq);
0060 j1 = 1;         j2 = npv;           %% j1:j2 - V angle of pv buses
0061 j3 = j2 + 1;    j4 = j2 + npq;      %% j3:j4 - V angle of pq buses
0062 j5 = j4 + 1;    j6 = j4 + npv;      %% j5:j6 - Q of pv buses
0063 j7 = j6 + 1;    j8 = j6 + npq;      %% j7:j8 - V mag of pq buses
0064 
0065 %% evaluate F(x0)
0066 Sb = Sbus(Vm);
0067 Sb(pv) = real(Sb(pv)) + 1j * imag(V(pv) .* conj(Ybus(pv, :) * V));
0068 mis = Ybus * V - conj(Sb ./ V);
0069 F = [   real(mis([pv; pq]));
0070         imag(mis([pv; pq]))   ];
0071 
0072 %% check tolerance
0073 normF = norm(F, inf);
0074 if mpopt.verbose > 1
0075     fprintf('\n it   max Ir & Ii mismatch (p.u.)');
0076     fprintf('\n----  ---------------------------');
0077     fprintf('\n%3d        %10.3e', i, normF);
0078 end
0079 if normF < tol
0080     converged = 1;
0081     if mpopt.verbose > 1
0082         fprintf('\nConverged!\n');
0083     end
0084 end
0085 
0086 %% attempt to pick fastest linear solver, if not specified
0087 if isempty(lin_solver)
0088     nx = length(F);
0089     if nx <= 10 || have_feature('octave')
0090         lin_solver = '\';       %% default \ operator
0091     else    %% MATLAB and nx > 10 or Octave and nx > 2000
0092         lin_solver = 'LU3';     %% LU decomp with 3 output args, AMD ordering
0093     end
0094 end
0095 
0096 %% do Newton iterations
0097 while (~converged && i < max_it)
0098     %% update iteration counter
0099     i = i + 1;
0100 
0101     %% evaluate Jacobian
0102     dImis_dQ = sparse(pv, pv, 1j./conj(V(pv)), n, n);
0103     [dImis_dVa, dImis_dVm] = dImis_dV(Sb, Ybus, V);
0104     dImis_dVm(:, pv) = dImis_dQ(:, pv);
0105 
0106     %% handling of derivatives for voltage dependent loads
0107     %% (not yet implemented) goes here
0108 
0109     j11 = real(dImis_dVa([pv; pq], [pv; pq]));
0110     j12 = real(dImis_dVm([pv; pq], [pv; pq]));
0111     j21 = imag(dImis_dVa([pv; pq], [pv; pq]));
0112     j22 = imag(dImis_dVm([pv; pq], [pv; pq]));
0113 
0114     J = [   j11 j12;
0115             j21 j22;    ];
0116 
0117 %     %% evaluate Jacobian
0118 %     dImis_dQ = sparse(pv, pv, 1j./conj(V(pv)), n, n);
0119 %     [dImis_dVa, dImis_dVm] = dImis_dV(Sb, Ybus, V);
0120 %
0121 %     %% handling of derivatives for voltage dependent loads
0122 %     %% (not yet implemented) goes here
0123 %
0124 %     j11 = real(dImis_dVa([pv; pq], [pv; pq]));
0125 %     j12 = real(dImis_dQ([pv; pq], pv));
0126 %     j13 = real(dImis_dVm([pv; pq], pq));
0127 %     j21 = imag(dImis_dVa([pv; pq], [pv; pq]));
0128 %     j22 = imag(dImis_dQ([pv; pq], pv));
0129 %     j23 = imag(dImis_dVm([pv; pq], pq));
0130 %
0131 %     J = [   j11 j12 j13;
0132 %             j21 j22 j23;    ];
0133 
0134     %% compute update step
0135     dx = mplinsolve(J, -F, lin_solver);
0136 
0137     %% update voltage
0138     if npv
0139         Va(pv) = Va(pv) + dx(j1:j2);
0140         Sb(pv) = real(Sb(pv)) + 1j * (imag(Sb(pv)) + dx(j5:j6));
0141     end
0142     if npq
0143         Va(pq) = Va(pq) + dx(j3:j4);
0144         Vm(pq) = Vm(pq) + dx(j7:j8);
0145     end
0146     V = Vm .* exp(1j * Va);
0147     Vm = abs(V);            %% update Vm and Va again in case
0148     Va = angle(V);          %% we wrapped around with a negative Vm
0149 
0150     %% evalute F(x)
0151     mis = Ybus * V - conj(Sb ./ V);
0152     F = [   real(mis([pv; pq]));
0153             imag(mis([pv; pq]))   ];
0154 
0155     %% check for convergence
0156     normF = norm(F, inf);
0157     if mpopt.verbose > 1
0158         fprintf('\n%3d        %10.3e', i, normF);
0159     end
0160     if normF < tol
0161         converged = 1;
0162         if mpopt.verbose
0163             fprintf('\nNewton''s method power flow (current balance, polar) converged in %d iterations.\n', i);
0164         end
0165     end
0166 end
0167 
0168 if mpopt.verbose
0169     if ~converged
0170         fprintf('\nNewton''s method power flow (current balance, polar) did not converge in %d iterations.\n', i);
0171     end
0172 end

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