Home > matpower7.1 > lib > newtonpf_S_hybrid.m

newtonpf_S_hybrid

PURPOSE ^

NEWTONPF_S_HYBRID Solves power flow using full Newton's method (power/hybrid)

SYNOPSIS ^

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

DESCRIPTION ^

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

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

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