Home > matpower7.1 > lib > newtonpf_I_hybrid.m

newtonpf_I_hybrid

PURPOSE ^

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

SYNOPSIS ^

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

DESCRIPTION ^

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

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

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

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

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