Home > matpower7.1 > lib > newtonpf_I_cart.m

newtonpf_I_cart

PURPOSE ^

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

SYNOPSIS ^

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

DESCRIPTION ^

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

   Solves for bus voltages using a full Newton-Raphson method, using nodal
   current 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_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_cart(Ybus, Sbus, V0, ref, pv, pq, mpopt)
0002 %NEWTONPF_I_CART  Solves power flow using full Newton's method (current/cartesian)
0003 %   [V, CONVERGED, I] = NEWTONPF_I_CART(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 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_S_CART, NEWTONPF_I_POLAR.
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 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 - Q of pv buses
0061 j3 = j2 + 1;    j4 = j2 + npq;      %% j3:j4 - Vr of pq buses
0062 j5 = j4 + 1;    j6 = j4 + npv;      %% j5:j6 - Vr of pv buses
0063 j7 = j6 + 1;    j8 = j6 + npq;      %% j7:j8 - Vi of pq buses
0064 j9 = j8 + 1;    j10= j8 + npv;      %% j9:j10- Vi of pv 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         V(pv) .* conj(V(pv)) - Vmpv.^2  ];
0073 
0074 %% check tolerance
0075 normF = norm(F, inf);
0076 if mpopt.verbose > 1
0077     fprintf('\n it   max Ir & Ii mismatch (p.u.)');
0078     fprintf('\n----  ---------------------------');
0079     fprintf('\n%3d        %10.3e', i, normF);
0080 end
0081 if normF < tol
0082     converged = 1;
0083     if mpopt.verbose > 1
0084         fprintf('\nConverged!\n');
0085     end
0086 end
0087 
0088 %% attempt to pick fastest linear solver, if not specified
0089 if isempty(lin_solver)
0090     nx = length(F);
0091     if nx <= 10 || have_feature('octave')
0092         lin_solver = '\';       %% default \ operator
0093     else    %% MATLAB and nx > 10 or Octave and nx > 2000
0094         lin_solver = 'LU3';     %% LU decomp with 3 output args, AMD ordering
0095     end
0096 end
0097 
0098 %% do Newton iterations
0099 while (~converged && i < max_it)
0100     %% update iteration counter
0101     i = i + 1;
0102 
0103     %% evaluate Jacobian
0104     dImis_dQ = sparse(pv, pv, 1j./conj(V(pv)), n, n);
0105     dV2_dVr = sparse(1:npv, npq+(1:npv), 2*real(V(pv)), npv, npv+npq);
0106     dV2_dVi = sparse(1:npv, npq+(1:npv), 2*imag(V(pv)), npv, npv+npq);
0107     [dImis_dVr, dImis_dVi] = dImis_dV(Sb, Ybus, V, 1);
0108 
0109     %% handling of derivatives for voltage dependent loads
0110     %% (not yet implemented) goes here
0111 
0112     j11 = real(dImis_dQ([pv; pq], pv));
0113     j12 = real(dImis_dVr([pv; pq], [pq; pv]));
0114     j13 = real(dImis_dVi([pv; pq], [pq; pv]));
0115     j21 = imag(dImis_dQ([pv; pq], pv));
0116     j22 = imag(dImis_dVr([pv; pq], [pq; pv]));
0117     j23 = imag(dImis_dVi([pv; pq], [pq; pv]));
0118     j31 = sparse(npv, npv);
0119     j32 = dV2_dVr;
0120     j33 = dV2_dVi;
0121 
0122     J = [   j11 j12 j13;
0123             j21 j22 j23;
0124             j31 j32 j33;    ];
0125 
0126     %% compute update step
0127     dx = mplinsolve(J, -F, lin_solver);
0128 
0129     %% update voltage
0130     if npv
0131         V(pv) = V(pv) + dx(j5:j6) + 1j * dx(j9:j10);
0132         Sb(pv) = real(Sb(pv)) + 1j * (imag(Sb(pv)) + dx(j1:j2));
0133     end
0134     if npq
0135         V(pq) = V(pq) + dx(j3:j4) + 1j * dx(j7:j8);
0136     end
0137 
0138     %% evalute F(x)
0139     mis = Ybus * V - conj(Sb ./ V);
0140     F = [   real(mis([pv; pq]));
0141             imag(mis([pv; pq]));
0142             V(pv) .* conj(V(pv)) - Vmpv.^2  ];
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, cartesian) 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, cartesian) did not converge in %d iterations.\n', i);
0160     end
0161 end

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