Home > matpower7.1 > lib > newtonpf_S_cart.m

newtonpf_S_cart

PURPOSE ^

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

SYNOPSIS ^

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

DESCRIPTION ^

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

   Solves for bus voltages using a full Newton-Raphson method, using nodal
   power balance equations and cartesian 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_I_POLAR, NEWTONPF_I_CART.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [V, converged, i] = newtonpf_S_cart(Ybus, Sbus, V0, ref, pv, pq, mpopt)
0002 %NEWTONPF_S_CART  Solves power flow using full Newton's method (power/cartesian)
0003 %   [V, CONVERGED, I] = NEWTONPF_S_CART(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 cartesian 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_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 %   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 Vm = abs(V);
0054 Vmpv = Vm(pv);
0055 
0056 %% set up indexing for updating V
0057 npv = length(pv);
0058 npq = length(pq);
0059 j1 = 1;         j2 = npq;           %% j1:j2 - Vr of pq buses
0060 j3 = j2 + 1;    j4 = j2 + npv;      %% j3:j4 - Vr of pv buses
0061 j5 = j4 + 1;    j6 = j4 + npq;      %% j5:j6 - Vi of pq buses
0062 j7 = j6 + 1;    j8 = j6 + npv;      %% j7:j8 - Vi of pv buses
0063 
0064 %% evaluate F(x0)
0065 mis = V .* conj(Ybus * V) - Sbus(Vm);
0066 F = [   real(mis([pq; pv]));
0067         imag(mis(pq));
0068         V(pv) .* conj(V(pv)) - Vmpv.^2  ];
0069 
0070 %% check tolerance
0071 normF = norm(F, inf);
0072 if mpopt.verbose > 1
0073     fprintf('\n it    max P & Q mismatch (p.u.)');
0074     fprintf('\n----  ---------------------------');
0075     fprintf('\n%3d        %10.3e', i, normF);
0076 end
0077 if normF < tol
0078     converged = 1;
0079     if mpopt.verbose > 1
0080         fprintf('\nConverged!\n');
0081     end
0082 end
0083 
0084 %% attempt to pick fastest linear solver, if not specified
0085 if isempty(lin_solver)
0086     nx = length(F);
0087     if nx <= 10 || have_feature('octave')
0088         lin_solver = '\';       %% default \ operator
0089     else    %% MATLAB and nx > 10 or Octave and nx > 2000
0090         lin_solver = 'LU3';     %% LU decomp with 3 output args, AMD ordering
0091     end
0092 end
0093 
0094 %% do Newton iterations
0095 while (~converged && i < max_it)
0096     %% update iteration counter
0097     i = i + 1;
0098 
0099     %% evaluate Jacobian
0100     [dSbus_dVr, dSbus_dVi] = dSbus_dV(Ybus, V, 1);
0101     dV2_dVr = sparse(1:npv, npq+(1:npv), 2*real(V(pv)), npv, npv+npq);
0102     dV2_dVi = sparse(1:npv, npq+(1:npv), 2*imag(V(pv)), npv, npv+npq);
0103 
0104     %% handling of derivatives for voltage dependent loads
0105     %% (not yet implemented) goes here
0106 
0107     j11 = real(dSbus_dVr([pq; pv], [pq; pv]));
0108     j12 = real(dSbus_dVi([pq; pv], [pq; pv]));
0109     j21 = imag(dSbus_dVr(pq, [pq; pv]));
0110     j22 = imag(dSbus_dVi(pq, [pq; pv]));
0111     j31 = dV2_dVr;
0112     j32 = dV2_dVi;
0113 
0114     J = [   j11 j12;
0115             j21 j22;
0116             j31 j32;    ];
0117 
0118     %% compute update step
0119     dx = mplinsolve(J, -F, lin_solver);
0120 
0121     %% update voltage
0122     if npv
0123         V(pv) = V(pv) + dx(j3:j4) + 1j * dx(j7:j8);
0124     end
0125     if npq
0126         V(pq) = V(pq) + dx(j1:j2) + 1j * dx(j5:j6);
0127     end
0128 
0129     %% evalute F(x)
0130     mis = V .* conj(Ybus * V) - Sbus(Vm);
0131     F = [   real(mis([pq; pv]));
0132             imag(mis(pq));
0133             V(pv) .* conj(V(pv)) - Vmpv.^2  ];
0134 
0135     %% check for convergence
0136     normF = norm(F, inf);
0137     if mpopt.verbose > 1
0138         fprintf('\n%3d        %10.3e', i, normF);
0139     end
0140     if normF < tol
0141         converged = 1;
0142         if mpopt.verbose
0143             fprintf('\nNewton''s method power flow (power balance, cartesian) converged in %d iterations.\n', i);
0144         end
0145     end
0146 end
0147 
0148 if mpopt.verbose
0149     if ~converged
0150         fprintf('\nNewton''s method power flow (power balance, cartesian) did not converge in %d iterations.\n', i);
0151     end
0152 end

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