Open Collections

UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Time-domain based steady-state initialization for the EMTP Perkins, Brian K. 1993

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata

Download

Media
831-ubc_1993_spring_perkins_brian.pdf [ 2.29MB ]
Metadata
JSON: 831-1.0064915.json
JSON-LD: 831-1.0064915-ld.json
RDF/XML (Pretty): 831-1.0064915-rdf.xml
RDF/JSON: 831-1.0064915-rdf.json
Turtle: 831-1.0064915-turtle.txt
N-Triples: 831-1.0064915-rdf-ntriples.txt
Original Record: 831-1.0064915-source.json
Full Text
831-1.0064915-fulltext.txt
Citation
831-1.0064915.ris

Full Text

TIME-DOMAIN BASED STEADY-STATE INITIALIZATION FOR THE EMTP by BRIAN KENNETH PERKINS B.A.Sc., The Univeristy of Toronto, 1991  A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF APPLIED SCIENCE  in THE FACULTY OF GRADUATE STUDIES DEPARTMENT OF ELECTRICAL ENGINEERING  We accept this thesis as conforming to the required standard  THE UNIVERSITY OF BRITISH COLUMBIA April 1993 © Brian Kenneth Perkins, 1993  In presenting this thesis in partial fulfilment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission.  (Signature)  Department of  L7<e 71rieal 647;49er/ill  The University of British Columbia Vancouver, Canada  Date  DE-6 (2/88)  ..i/r-11 2/, /993  Abstract The existing steady-state initialization algorithm in the electromagnetic transients program EMTP is based on a frequency-domain admittance method. This yields a simple, efficient solution for linear, lumped element networks. However, the basic algorithm cannot be easily extended to handle nonlinear and variable structure networks. A time-domain based steady-state initialization algorithm is formulated for the three classes of lumped element networks, namely, linear, nonlinear and variable structure. The primary contribution of this work is that each network element does not require a dedicated steadystate model and a separate steady-state solution algorithm. Furthermore, as the proposed initialization process is consistent with the time-domain character of the transient simulations themselves, the proposed algorithm can be incorporated into existing electromagnetic transient programs with minimal programming effort.  ii  Table of Contents ii  Abstract ^ List of Tables ^  vi  List of Figures ^  vii  List of Algorithms ^  ix  Chapter 1^Introduction ^  1  1^Background ^  1  2^Existing methods ^  1  3^Scope of the work ^  2  4^Thesis organization ^  2  Chapter 2^Steady-state initialization: State-space formulation ^ 4 1^Introduction ^  4  2^Motivating Example ^  4  3^Formulating the Fixed-Point Problem ^ 7 4^State-space steady-state initialization algorithm ^ 9 5^Evaluating the Jacobian Matrix ^  10  6^Application to the EMTP ^  10  Chapter 3^Steady-state initialization: Linear Networks ^ 11 1^Introduction ^  11  2^Formulation of the network equations ^ 11 3^Triangularization ^  12  4^Solution method ^  12  5^Numerical example: Transient analysis ^ 13 6^Steady-state initialization ^  15  7^Computing the Jacobian Matrix ^  16  8^Steady-state initialization algorithm ^ 17 iii  9^Numerical example: Steady-state initialization ^ 18 10^Application: Capacitive and Inductive Coupling into fence ^ 20 Chapter 4^Steady-state initialization: Nonlinear Networks ^ 23 1^Introduction ^  23  2^Representation of nonlinear elements ^ 23 3^Formulation of the network equations ^ 24 4^Network equivalents ^  25  5^Solution method ^  26  6^Numerical example: Transient analysis ^ 29 7^Steady-state initialization ^  31  8^Computing the Jacobian Matrix ^  32  9^Steady-state initialization algorithm ^ 33 10^Numerical example: Steady-state initialization ^ 34 11^Application: Harmonics due to Transformer Saturation 36 Chapter 5^Steady-state initialization: Variable-structure networks 38 1^Introduction ^  38  2^Representation of ideal-switch devices ^ 38 3^Formulation of the network equations ^ 39 4^Network equivalents ^  39  5^Solution method ^  40  6^Numerical example: Transient analysis ^ 41 7^Steady-state initialization algorithm ^ 44 8^Numerical example: Steady-state initialization ^ 45 9^Application: Initialization of single and three-phase bridge networks ^ iv  46  Chapter 6^Conclusions ^  49  1^Summary ^  49  References ^  50  List of Tables Table 2.1^Progresssion of the fixed-point iteration ^ 9 Table 3.2^Phasor representation of steady-state solution for continuous-time and discrete-time network ^ 20 Table 3.3^Computed steady-state values of the induced voltage on the fence. ^  22  Table 4.4^Phasor representation of voltage and current harmonics in nonlinear inductor (At = 312.5ps) ^  vi  37  List of Figures Figure 2.1^Simple RLC circuit. ^  4  Figure 2.2^RLC circuit response from zero initial conditions ^ 6 Figure 2.3^RLC circuit response from steady-state initial conditions . . ^ 7 Figure 3.4^Simple linear, lumped element network (a) continuous-time network (b) discrete-time network ^  14  Figure 3.5^Solution of sample network from zero initial conditions . . ^ 15 Figure 3.6^Solution of sample network from steady-state initial conditions ^  19  Figure 3.7^Electrostatic and electromagnetic coupling from power line to fence ^  21  Figure 4.8^Convergence and divergence behaviour of fixed-point iteration ^  28  Figure 4.9^Solution of sample network from zero initial conditions . . ^ 31 Figure 4.10^Progression of the steady-state initialization algorithm from zero initial conditions ^  35  Figure 4.11^Progression of the steady-state initialization algorithm from linearized network initial conditions ^ 35 Figure 4.12^Steady-state waveform for the sample network (At = 312.5ps) ^  36  Figure 5.13^State transition diagram for an ideal diode ^ 38 Figure 5.14^Half-bridge rectifier with freewheeling diode (a) continuous-time network (b) discrete-time network ^ 42 Figure 5.15^Solution of sample network from zero initial conditions . . ^ 43 Figure 5.16^Solution of sample network from steady-state initial conditions ^ Figure 5.17^Single-phase diode bridge rectifier ^ vii  46 47  Figure 5.18^Steady-state waveforms for single-phase diode bridge rectifier ^ Figure 5.19^Three-phase diode bridge rectifier ^  47 48  Figure 5.20^Steady-state waveforms for three-phase diode bridge rectifier ^  viii  48  ^  List of Algorithms ^3.1^Transient solve code for linear networks ^ 13 3.2^Time-domain steady-state initialization algorithm for linear networks ^  18  4.3^Transient solve code for nonlinear networks ^ 29 4.4^Time-domain steady-state initialization algorithm for nonlinear networks ^  34  5.5^Transient solve code for variable structure networks ^ 41 5.6^Time-domain steady-state initialization algorithm for variable structure networks ^  ix  45  Chapter 1 Introduction 1.1 Background The EMTP is a digital computer-based program for the solution of transients in large scale networks typical of power systems [1]. The ordinary differential equations (ODEs) describing the network are formulated in nodal form. Solution methods are based on the trapezoidal rule of integration, an implicit numerical method which insures stable solutions of networks with widely seperated eigenvalues. 1 EMTP development can be classified into two major categories: network component modelling and solution methods. Much work has gone into the developement of computerbased models of power system components and their implementation within the basic EMTP framework. Examples are transmission line models, synchronous machine models, and HVdc converter models, to name just a few. The basic EMTP solution methods, however, have not undergone any significant changes thanks in large part to the inherent robustness of the original formulation. There have, nonetheless, been significant contributions in this area. Solution methods for networks with nonlinear or time-varying elements [2] and the critical damping adjustment scheme (CDA) [3] are especially noteworthy. TACS is a contribution that merits membership in both of the above categories. 1.2 Existing methods The periodic steady-state solution of a network is of particular interest to the power systems engineer. Typically, transient studies are executed from the periodic steadystate. Alternatively, the steady-state solution may be of interest in itself for power flow calculations and harmonic analysis. The existing steady-state initialization algorithm is based on a frequency-domain admittance matrix approach [1]. It is limited to the solution of linear, lumped element networks. Harmonic balance methods have been applied to networks containing nonlinear elements Such systems are said to be stiff.  [4] [5]. However, these methods apply to only a small class of networks and have yet to find widespread implementation. The main drawback of frequency-domain methods for steady-state initialization is that they make little or no use of the time-domain solve code common to electromagnetic transient programs; they are essentially standalone programs that interface with the EMTP. 1.3 Scope of the work Various time-domain approaches for the periodic steady-state solution of systems of ordinary differential equations have been proposed [6]. These approaches are invariably based on a state-space formulation of the ODEs. Among them, the autonomous shooting method is by far the simplest and most fundamental. Reference [7] details the application of such a method to the periodic steady-state solution of nonlinear networks with periodic inputs. The proposed steady-state initialization method based on the seminal work of Aprille and Trick [7] offers two significant advantages: 1. the basic formulation applies globally to linear, nonlinear and variable structure network s 2 2. the method is consistent with the time-domain character of the transient simulations themselves 1.4 Thesis organization Chapter 1 serves as a general introduction. Chapter 2 introduces the state-space formulation for a continuous-time-domain based steady-state initialization algorithm. Chapters 3, 4 and 5 detail the development of the time-domain based steady-state initialization algorithm within the framework of the EMTP solve code for three classes of lumped element networks: 1. multiphase linear, lumped element networks (Chapter 3) 2. nonlinear, lumped element networks (Chapter 4) 3. variable structure networks (Chapter 5). 2  ^  Networks with ideal-switch devices (e.g. ideal diodes).  2  Chapter 6 concludes the thesis.  3  Chapter 2 Steady-state initialization: State-space formulation 2.1 Introduction In the analysis of networks driven by periodic inputs, one may be interested in determining the steady-state periodic response. Even for transient simulations, one typically begins at steady-state initial conditions. Unless the network is linear, time-invariant (LTI) for which simple frequency domain analysis yields the steady-state response, one generally has no alternative but to simulate the system over a sufficiently long interval of time for the transient waveform to die out. This procedure is satisfactory if the transient decays rapidly. However, for lightly damped networks typical of power systems the transient will decay very slowly and it becomes prohibitively expensive for the computer to simulate over such a long transient regime. Before presenting an algorithm for overcoming this difficulty, it is instructive to examine a simple LTI example to motivate the problem. 2.2 Motivating Example Consider the simple RLC circuit shown in Fig. 2.1.  R  G C 0S 4.1  Figure 2.1 Simple RLC circuit.  4  The state equations are given by 1  Ti t) c(t) =^L(t) 1  (2.1)  L(t) = (E cos wt — Ri L (t) — vc(t))  dt The complete response is given by  vc (i) = (k1 cosh fit + k2 sinh fit) exp at + Vc cos Pt + 0)  iL(t) = (k3 cosh fit + k4 sinh fit) exp at — wCl/c sin (wt + 0) where  (2.2)  E  VC =  1 — w 2 L C ) 2 + (wC R) 2  9 = — arctan  wC R 1 1_1 — w 2 LCj  (2.3)  R\2 1 a =^=^ 2L^) LC The four arbitrary constants are determined by the initial states vc(0) and iL(0): k 1 = vc(0) — Vc cos 0 k 2 = [3 -1 [C -1 i L(0) + wi7c sin 9 — avc(0) + aVc cos 0]^(2.4)  k 3 = C[ak i + /3k2], k4 = C[ak2 i3k1] It follows from (2.4) that k1 = k2 = k3 = k4 = 0 if and only if v c (0) = Vc cos 0 iL(0) = —u CVc sin 0  (2.5)  ,  With this choice of initial states the complete solution is given by vc (t) = Vc cos (cot + 0) ,  iL(t) = —wCVc sin (wt + 0)  (2.6)  both of which are periodic waveforms of period T = . Hence it is possible to suppress the transient component if the initial states are chosen as above. Numerical example: Consider a lightly damped case of the simple RLC circuit with  parameters E = 1 w = 1007r R = 0.04 ohms C = 200 uF L 1 mH  = —20 .^= j2236 5  (2.7)  Fig. 2.2 depicts the circuit response from zero initial conditions. The transient component takes over ten cycles to decay. Fig. 2.3 depicts the circuit response with steady-state initial conditions specified as in (2.5) corresponding to the steady-state solution vc (t) = 1.0201 cos (10Ort — 0.0026) ,  (2.8)  L(t) = —0.0641 sin (10Ort — 0.0026) §  2 v(t) 1.5  I  1  0.5  i i (t)  ll  0  AP"  I  -0.5  II  -1.5  0  ^  1  ■  0.02^0.04^0.06^0.08^0.1^0.12^0.14^0.16^0.18^0 2 t (s)  Figure 2.2 RLC circuit response from zero initial conditions  6  Figure 2.3 RLC circuit response from steady-state initial conditions  The objective is to generalize the concept of supressing the transient component and devise an algorithm for determining the appropriate initial state x o = x(0) such that the complete solution to the state equations dx dt^f (x ' t)  (2.9)  has no transient component under the assumption that all inputs are periodic in period T f(x,t)^f(x,t T)^  (2.10)  where f(.) is possibly nonlinear. 2.3 Formulating the Fixed-Point Problem Let x(t) be the solution to (2.9) with initial state xo. Integrating both sides of (2.9) from time 0 to time t we obtain x(t)  = f f(x(t),t)dt xo = x(t, xo) 7  (2.11)  The right hand side of (2.11) is denoted by x(t, ro) to emphasize that the solution at any time depends on the initial state. The goal is to find x o such that at t T, x(T, x o ) x o . Since this relationship is seldom satisfied for an arbitrary choice of x o , define the function F(xo) =  x(T, xo)  = f f (x(t),t)dt xo  (2.12)  0  Since we are fixing the time t^T, F(xo) is a function of xo and is independent of t. In the scalar case the function can be represented in the F versus x o plane. It follows from (2.12) that the initial state x o = xo which gives rise to a periodic solution 4(0 must satisfy the equation xo  = F(x 0 )^  ( 2.13)  which is identical to the standard form for a fixed-point iteration [8]. Hence the solution can be found by applying the fixed point iteration algorithm PH-1 ) F (x(0 ) 0  (2.14)  Invoking the assumption f (x ,t) f (x ,t T) it follows that the (i +1)th iterate is given by (i+1)T  (i) x o( i+1) = Fx 0 )— =^(x(t),t)dt x 0( °)  (2.15)  0  This iteration must be repeated until F (x o( i) ) — x 0( i)  <E  (2.16)  where c is some prespecified quantity. The fixed-point iteration algorithm is equivalent to integrating (2.9) until all transient components have decayed. Numerical example: Consider the simple RLC circuit of Section 2.2 with the states defined  as x(t) = [ v e (t) iL(t)] T . The progression of the fixed-point iteration can be calculated from the response (2.2). Table 2.1 summarizes the progression of the fixed-point iteration. 8  Note that in the LTI case, the speed of convergence is dependent on the eigenvalues of the system.§ i  11 F ( 4) ) — x (i) 11 o  0  5.5786e-1 7.9206e-2 1.3779e-3 1.4568e-7 1.7792e-16  5 10 20 50  Table 2.1 Progresssion of the fixed-point iteration  2.4 State-space steady-state initialization algorithm A more efficient algorithm is obtained by formulating (2.13) as a two-point boundary value problem [7]. Rewrite (2.13) in terms of an objective function defined H(x0) = F(xo) — xo^  (2.17)  for which the periodic steady-state initial conditions yield H^= 0. Convergence to x'01` is accelerated by a Newton-Raphson iteration on (2.17) 4+1) x02)  [DF (431 — IT] [F  (x0 i)  ) — x :( ) ]^(2.18)  where D F x o( z ) is the Jacobian matrix of F(xo) evaluated at x o( i) . Thus we obtain quadratic convergence to the periodic steady-state as opposed to the possibly very slow convergence of the fixed-point iteration algorithm. In fact, for a LTI system convergence is realized in a single Newton-Raphson iteration. Numerical example: For the simple RLC network of Section 2.2 the above algorithm  yields 41) = [ 1.0201 —0.06411 T^(2.19)  on the first iteration. The above is equivalent to the initial conditions obtained by (2.5).§ 9  2.5 Evaluating the Jacobian Matrix In general, one cannot obtain a closed form expression for the Jacobian matrix as in the LT' case. For a system of dimension n the Jacobian matrix is given by axi(T,s0)^axi (T,x0)  axi(o)^axn(0) DF(xo) =^.^.^(2.20) • ax„(T,x0)^ax„(,x0) ax (o)^axn(o)  i  Evaluating the Jacobian matrix by finite differences is inefficient and prone to numerical error as one must solve the possibly nonlinear differential equations (n + 1) times. A more accurate and efficient method involves the solution of related sensitivity systems which describe the linearization of the nonlinear system along the trajectory of the solution. This method involves the solution of the nonlinear system and its n associated linear sensitivity systems. The method is formulated as follows. Given that the possibly nonlinear system dx f(x,t)^x(0) = x o^(2.21) dt^ has solution x(t, x0), then the n associated sensitivity systems are defined by  dx ^a^a o f(x,t) ax ( 0) _ [0 • •^1 • • • OF {j = 1..n}^(2.22) ^axj(0) dt^axj(0) whose solutions at time T represents the j = 1..n columns of the Jacobian matrix.  ^a  Rearranging (2.22) the sensitivity systems become  where  ax  d ax^a f, ^ax (2.23 x(t =^ dt ax 3 (0)^ax^'^axj(0)^  )  )  ^the first variation (linearization) of f along the solution trajectory.  In the application of implicit integration methods for the solution of (2.9), -g8 . is readily available.  2.6 Application to the EMTP The state-space formulation presented in this chapter is fine for pedagogical purposes but impractical for large scale networks. The work details a formulation of the steady-state initialization algorithm within the framework of fixed time-step electromagnetic transients programs (EMTP). The algorithm is formulated for the three classes of lumped element networks, namely linear, nonlinear, and variable structure. 10  Chapter 3 Steady-state initialization: Linear Networks 3.1 Introduction A time-domain based steady-state initialization algorithm for linear, lumped element networks within the framework of electromagnetic transient programs is presented in this chapter. The formulation and solution of such networks is also presented to facilitate the derivation of the steady-state initialization algorithm.  3.2 Formulation of the network equations The linear equations of the discrete-time network can be formulated as follows [1] ^[G][V (t)] = [J (t)] = [A h ][Jh (t)] + [A i ][Ji(t)]^(3.24)  where [G] is the symmetric nodal conductance matrix, [V(t)] is the vector of node voltages, [At)] is a vector containing contributions from Norton equivalent history sources [Jh (0]  and injected current sources [J3 (t)]. [Ah] and [A 3 ] are incidence matrices mapping branch current contributions onto the nodal current contributions. The vector of node voltages is partitioned into unknown [VA (t)] and known [VB (t)] voltages. Subdividing the matrices accordingly, the unknown node voltages at each time-step are given by [GAA][VA(t)]^[JA (t)]^[GAB][VB(t)]  where [JA (t ] = [AhA] [Jh ( t )] + [A3 ] [J, ( t )]  ^  ^  )  (3.25)  (3.26)  The Norton equivalent history sources are updated recursively as follows 4 ^[Jh(t + At)] = [LI h]([J h(t)] + 2[G h][A,][V (t)]) ^(3.27) 4  ^  This formulation is applicable to both monophase and multiphase networks.  11  where [uhi diag ( +[U], [Jh] associated with [L] — [U], [Jh ] associated with [C] )  (3.28)  f [L]^[A] associated with [L]) Ph ] = dia g o 2 -67[C], [Jh] associated with [C] )  (3.29)  [Ay] = — [Ah] T^(3.30)  3.3 Triangularization The network equations are simplified prior to the first time step by upper triangularization of [GAA]. The triangularization proceeds by downward operations as follows [GAA GAB AhA AjA] ^[-1:][GAA GAB AhA AjA ]  ^  (3.31)  where [P] represents the permutation matrix that upper triangularizes [GAA]. 3.4 Solution method For the above formulation a solution method akin to that of the EMTP is presented in Algorithm 3.1. Note that the algorithm presented differs slightly from that of the EMTP in that the history sources are updated at the end of the solve code loop as opposed to the beginning. This permits the initial conditions to be specified in terms of Norton equivalent history terms as opposed to branch voltages and currents.  solve code( t +— to, [A]  4—  [Jho]  while t < tf, build rhs( —  Algorithm 3.1 Transient solve code for linear networks (Continued ... ) 12  ^  compute([VB], [rhs] +  —  [AhA] [A]  + [A  ^ -  [  GAB ][ vB ]  solve ( [VA] 4— backsub([GAA], [ rhs])  update_history_terms( [Jd <— [Uh] • GA] + 2[Gh] [Ad [V]  t^t + At  end while  Algorithm 3.1 Transient solve code for linear networks  3.5 Numerical example: Transient analysis Consider the simple network shown in Fig. 3.4(a). Fig. 3.4(b) depicts the corresponding  discrete-time network with inductances and capacitances represented as equivalent conductances shunted by Norton equivalent history sources. Formulation of the network equations: With nodes and history sources defined as in Fig.  3.4(b) the discrete-time network equations become gR s [  —  vi(t)]  -11  [0 —1  gL,^—gLs gL s gLs^gC gL ni i[v2(t)  0 —1^1  [  with history sources updated recursively as follows -^. ^0 ^1^0 phi (t + At) lhi(t) ^0 ih2 (t + At) ih 2 (t) + 2[Gh ] ^1^0 [Uh] At) 1^—1 0 (t h3 113(0  +  _  ,  ih,(t) 1'2(0 3 11 3 (t)  [ 6][v3(t)] —^ 0R  .  (3.32)  v 1 (t )l ) [v2t  v3 (t)  (3.33) [Uh] = diag( [1^—1  1 ])  = diag([gLm^gc gLs 1) 13  ^E(  ^/ se)? ( /0011'11) '  O- Oo 4to 8 _a.  ^4^ E(4) 0./.3 c --= /. 6 777F Zfrt= 3/5.3 m  E (79  =L  c2C pn.s  yks =^024/5;0 980 zis?'^  az,  1c  ac ^ = .2-56  8077  96.3C e-3  (b) Figure 3.4 Simple linear, lumped element network (a) continuous-time network (b) discrete-time network  Triangularization: For the example at hand with At = 1.25 ms downward operations  transform the augmented matrix [ GAA GAB AhA Am] to 249.9057 —4.8077 —245.0980 ^0^0^—1 1 0^7.2772^—4.7152^—1 —1 0.9808  (3.34)  Solution method: Algorithm 3.1 with parameters to = 0 and tf = 200 ms is applied to the  sample network. Fig. 3.5 represents the progression of the voltage and currrent waveforms 14  v2(t) and (t) from zero initial conditions. Note that even after ten cycles the solution has  not yet settled to the periodic steady-state.  vit)  Figure 3.5 Solution of sample network from zero initial conditions  3.6 Steady-state initialization The Norton equivalent history sources represent the state of a discrete-time network. Define [F(Jho)] = [Jh(T; Jho)]^ (3.35)  as the state at time-step T on the solution trajectory uniquely defined by [Jho ] and computed by simulating the network over one period. A linear network with sources of period T has a steady-state solution of period T. The discrete-time network has reached the periodic steady-state when [F(Jho)]^[Jho]^ 15  (3.36)  In general, (3.36) is not satisfied by an arbitrary choice of [Jh o ]. We thus define the objective function [H Jho)] = [F(Jho)] — [ A d  ^  (  (3.37)  for which the steady-state initial condition [4 0 ] yields the solution (zero of the objective function). The two-point boundary value problem can be solved by Newton-Raphson on (3.37). For a linear network only a single iteration is required for convergence to the steady-state initial conditions. We thus have [40] = [Jho] — [DF(Jho) — u] i F(Jho) — Jho] —  ^  [  (3.38)  Since [Jho] is some arbitrary initial state, choose [Jho] = [0]. The steady-state initialization problem thus reduces to solving the linear system [DF(o)—^= —[ F (0)]  ^  (3.39)  3.7 Computing the Jacobian Matrix The proposed approach requires the efficient and accurate computation of [DF(Jh o )]. The brute force approach involves perturbing each of the  T1  states (Norton equivalent history  sources) in turn and computing the Jacobian matrix by finite-differences. This approach is computationally inefficient and prone to numerical innaccuracy. The sensitivity networks approach computes the Jacobian matrix by solving n linear sensitivity networks that are intrinsically related to the discrete-time network. The sensitivity networks approach is motivated as follows. [F(Jh o )] is computed by solving the discrete-time network (3.25) from t = 0 to t = T with history sources updated recursively as in (3.27). Analogously, [DF(Jho)] = [ 8.jhTi ;j" ) 1 is computed by solving , ,[31/A(t)]  [Li-AAJ  = [AhA] [ajh(t)1  aJho^aJho 16  (3.40)  ^  from t = 0 to t = T with history sources updated recursively as follows  rah(0)]  rh(t + At)]  ^[U  =Rib^t = 0 aJho^ [avA(t)-i h l ([aJh(t)] 2[Gh][A,A]  aJho^DA°^aJho  (3.41) t > At  Note that the solution of the n sensitivity networks is simply the matrix-valued solution of the discrete-time network with the time-dependent sources zeroed. The sensitivity networks  approach is a computationally efficient and accurate method for computing the Jacobian matrix. 3.8 Steady-state initialization algorithm Algorithm 3.2 depicts the proposed steady-state initialization algorithm in pseudocode form. Note that the algorithm is implemented within the framework of the solve code presented in Section 3.4. (The additional code over and above that in Algorithm 3.1 is marked with a double-dagger symbol (t)).  steady-state_init ( t  4-  0 tf4 ,  —  T, [Jh]^[0], [D^[U]  while t < t f, build rhs( compute([VB], [J3 ]) [rhs]^[AhA] [Jhi^[AAA] [J3] [Drhs]  4—  —  [G AB ] [VB]  [AhA][DJh] t  solve ( [VA] 4- backsub([GAA], [r hs]) [DVA ] 4— backsub^AA],[Drhs])  Algorithm 3.2 Time-domain steady-state initialization algorithm for linear networks (Continued ^)  17  ^  update_history_terms ( [Jh] <— [Uh] • Ph] + 2 [Gh][Av][V])  [DJh] 4— [Uh] • ([DJh] 2 [Gh][AvA][DVAD  ^t  ^t At  end while newtonraphson_iteration( [Ad^[DA  —  U] -1 • [— A]  t  AlgSrithm 3.2 Time-domain steady-state initialization algorithm for linear networks  3.9 Numerical example: Steady-state initialization The  algorithm can be used by itself to initialize the discrete-time network for transient  studies. A time-domain representation of the steady-state waveforms requires transient simulation over one period from the steady-state initial conditions. For the example of Section 3.5, [Jh(0)] is computed by both the existing frequency domain and the proposed time-domain steady-state initialization algorithm. [Jh (T)] represents the state of the network after an additional period of simulation. Fig. 3.6 represents the solution from the steadystate initial conditions computed by the time-domain method.  18  ^  Figure 3.6 Solution of sample network from steady-state initial conditions  With the existing frequency-domain initialization ^[-0.011231^—0.011242 [Jh(0)] =^0.55838^[A(T)] =^0.56912 0.55327^ 0.58210 1Jh(T) — [Jh(0)]11 2 = 3.0772e — 2  (3.42)  Note that the computed initial condition does not satisfy the periodicity constraint of the discrete-time network. For the proposed time-domain method [-0.011083 -^—0.011083 [A(0)] =^0.57793^[A(T)] =^0.57793 0.55327^ 0.55327  (3.43)  11 ,4(T) — Jh( 0 )112 = 0  the periodicity constraint is in fact satisfied. It is convenient, particularly for power flow calculations, to express the steady-state solution in phasor form. For example, given the sequence [v] = [v(0) v(At) v(2At) • • • v((N — 1)0t)]^(3.44) 19  where N  = ot , the phasor representation (in rectangular coordinates) is N-1  27r V =^v(nAt)exp (-jn Tv-)^  (3.45)  n=0  for which the polar form is expressed asVLq, V = (4)mag(V), = arg(V). Table 3.2 itemizes the phasor representation of the continuous-time solution and the discretetime network solution of the example at hand. The continuous-time solution is the exact Continuous-time network  Discrete-time network (At = 1.25 ms)  V 1 (V rms) V2 (V rms) V3 (V rms) Ii (A rms)  0.7778 L 0.7938 L  -  -  90.1175°  0.7778 L  90.1175°  0.7956 L - 90.1245°  -  90.1245°  90.0°  0.7778 L - 90.0°  0.7778 L -  0.0079 L179.8825°  0.0076 L179.8755°  Table 3.2 Phasor representation of steady-state solution for continuous-time and discrete-time network  solution of the network. The discrete-time network solution shows some discrepancy from the exact solution as the discrete-time network only approximates the continuous-time network.  3.10 Application: Capacitive and Inductive Coupling into fences A fence runs parallel to a three-phase power line as depicted in Fig. 3.7. By simply treating the fence as an additional conductor the following [C], [R] and [L] symmetrical matrices for a nominal multiphase ir -circuit are obtained: 3.78545 0.81330 0.81520 0.08440  1 - [C] =  2  [1=]  5  ^  =  0.8108 0.1148 0.1148 0.1162  3.65440 0.41745 0.13790  0.8108 0.1148 0.1162  This case is provided by Dr. H.W. Dommel  20  3.64995 0.0.5945  0.8108 0.1162  nF  (3.46)  3.48635  Il 3.7214  (3.47)  [L] =  • 5.2304•^• 2.2627 5.2304^•^• 2.2627 1.9852 5.2304^• 1.6807 1.7495 1.6149 5.2802  mH  (3.48)  Power line conductors: R = 0.348 ohms/km Xa = 0.386 ohms/km (60 Hz) (reactance at 1 foot spacing) diameter = 12.7 mm frequency = 60 Hz Fence: R = 1.802 ohms/km solid conductor (nonmagnetic) diameter = 4.064 mm length = 2 km  Figure 3.7 Electrostatic and electromagnetic coupling from power line to fence  Capacitive coupling: Under the assumption that the fence is insulated and is nowhere  grounded, the steady-state voltage on the fence due to capacitive coupling is calculated for two cases of power line energization: 1. symmetrical operation at 345 kV rms line-to-line 2. phase A grounded because of a nearby single line-to-ground fault, phases B and C at normal voltage. Inductive coupling: Under the assumption that the 2 km long fence is grounded at one end  and open at the other end (insulated in between), the steady-state voltage at the open end of the fence is calculated for two cases: 1. symmetrical operation with 'ph ase = 1 kA rms 2. IA = 10 kA, IB = lc =0 because of a nearby single line-to-ground fault. 21  The preceeding cases were solved by the existing EMTP frequency-domain based steadystate initialization algorithm as well as the proposed time-domain based steady-state initialization algorithm (for two vales of At). Table 3.3 summarizes the computed values for the induced voltage on the fence.  CASE  EMTP frequency-domain method  Capacitive coupling (case 1) Capacitive coupling (case 2) Inductive coupling (case 1) Inductive coupling (case 2)  3.9664 L — 101.87° kV  Proposed Proposed time-domain time-domain method method (At = 2.0833 ms) (At = 0.52083 ms) 3.9663 L — 101.87° 3.9663 L — 101.87° kV kV 6.8444 L — 145.45° 6.8443 L — 145.45° kV kV 42.907 L — 179.85° 45.112 L — 179.85° V V  6.8446 L — 145.45° kV 42.782 L — 179.87° V 6.4416 L — 100.39° kV  6.4618 L kV  —  100.36  °  6.7834 L — 99.86° kV  Table 3.3 Computed steady-state values of the induced voltage on the fence.  22  Chapter 4 Steady-state initialization: Nonlinear Networks 4.1 Introduction A general method for the formulation, solution, and steady-state initialization of nonlinear lumped element networks is presented in this chapter. The proposed formulation is limited to current-controlled nonlinear elements 6 although typically voltage-controlled nonlinear elements 7 can be expressed as current-controlled nonlinear elements.  4.2 Representation of nonlinear elements A continuously differentiable scalar function of current is expressed in linearized form as follows f (i(t))^m(t) i(t)^b(t)^  where m(t)  = af  )^ b(t)^f (i(t)) — m(t) • i(t) 1  (4.49)  (4.50)  Note that this permits a convenient representation of piecewise-linear functions. Currentcontrolled nonlinearities are thus represented as follows: Nonlinear resistance: A nonlinear resistance with branch characteristic v(t) = f(i(t)) is  expressed as v(t) = r(t) • i(t)^ek(t)^  where r(t) = m(t) represents a time-varying resistance and  e k (t)  (4.51)  represents a time-varying  'knee' voltage. 6^Current-controlled  nonlinear elements are those in which there is a one-to-one correspondence between the current and voltage (nonlinear resistance) or between the current and the state (nonlinear inductance). 7^Nonlinear  conductances and capacitances.  23  Nonlinear inductance: Upon the application of the trapezoidal rule of integration, a  nonlinear inductance with branch characteristic v(t) = as^  dd  (t)  ,  A(t)^f (i(t)) is expressed  v(t) = r(t) • i(t)^ek(t) + eh (t)  (4.52)  eh(t) = —v(t — At) — r(t — At) • i(t — At) — ek(t — At)  where r(t) im(t) represents a time-varying resistance and ek(t) = b(t) represents a time-varing 'knee' voltage. Note that the discretized nonlinear inductance characteristic is simply a scaled and shifted version of the nonlinear resistance characteristic. A network may contain any number of nonlinear resistances and inductances. We thus introduce the notation [V(/(t))] = [R(t)][/(t)] ^[Ek (t)]^[Ehoet)1^(4.53)  where [V(/(t))] and [I (t)] represent nonlinear branch voltages and currents respectively, and^ [R(t)]^diag([ • • il mi(t) • •^rni(t) • • • D  [Ek (t)]  = [• •^  -ib3 (t)^•^bi (t)^• • • i T^(4.54)  [Eh(t)]^[ • • •^eh,(t)^... ]T represent a diagonal matrix of time-varying resistances, a vector of time-varying 'knee'  voltages, and a vector of Thevenin equivalent history sources respectively. Thevenin equivalent history sources are updated as follows [Eh(t + At)] = [ U 0] GAni[V (t)] + [R(t)][i (t)] + Ek(t)) —  ^  (4.55)  where [ U 0] indicates that the computation proceeds only for branches with nonlinear —  inductances. [A n ] = [C] is an incidence matrix mapping the node voltages onto the nonlinear branch voltages, where [C] is defined in the following section. 4.3 Formulation of the network equations Current-controlled nonlinearities are incorporated into the nodal formulation of Section 3.2 by a compensation method as follows [9]  I  G B V (t)] _^J (t) [ C^0^[ I (t)^1_1/ (I (t)) 24  1  (4.56)  where  + 1, positive incidence  [C] = [c( v )] ,^c(u)^— 1, negative incidence  (4.57)  0, no incidence  maps the linear node voltages onto the nonlinear branch voltages and [B] = [C] T^(4.58)  maps the compensation currents onto the linear part of the network. As in Section 3.2 the vector of node voltages is partitioned into unknown and known node voltages. Partitioning the matrices accordingly, the unknown node voltages and the compensation currents are given by [GAA^B Al[VA(t)] CA^0 I (t)  [^JA(t) V(I(t))  [GAB] r vB (t)]^ (4.59) CB I -  The Norton equivalent history sources associated with linear capacitances and inductances are updated recursively as in Section 3.2. Thevenin equivalent history sources associated with nonlinear inductances are updated as in (4.55).  4.4 Network equivalents A straightforward application of an iterative method for solving the nonlinear system of equations at each time-step is prohibitively computationally intensive for large scale networks. If, however, [VA (t)] is 'decoupled' from [I (t)] by forming a network equivalent of the linear part of the network, an iterative solution method can be applied to the subset of nonlinear equations. The linear part of the network can then be solved by back substitution. A multidimensional Thevenin equivalent with respect to the nonlinear branches is synthesized prior to the first time-step as follows  [ P —K  GAA GAB  BA  AhA  0 CB 0 ][GAA  -T  A hA  U  CA  GAB CB  BA 0  A Ai i  AhA^AjAi 0^0  (4.60)  where [P] represents the permutation matrix that upper triangulanzes [GAA] and [K] = [C A] [GAA] -1 . In practice, the LHS of (4.60) is synthesized by downward operations that 25  leaves [GAA] upper triangularized and [CA] = [0]. We now have network equations of the form  [  ^ G AA BA][VA(t)^[^JA(t)^[GAB rT (4.61) Bkt)i ^CB1Lv —T^I (t)]^V  0^  (I (t)) + 4(01  where [Jit' (t)] =[ 11' hA] [Jh(t)]^  [AAA]  [Ji(t)]. Note that the nonlinear elements are now  coupled to a Thevenin equivalent of the linear part of the network, namely [V  (0]  —  [ T][I (0] = [ V (I (0)1^(4.62)  where [V,;(t)] = CB][VB(t)] — PA' (t)] represents the open circuit voltage. [  4.5 Solution method The proposed solution method is based on a piecewise-linear fixed-point iteration that yields Newton-Raphson-like convergence. Expressing [17(/(t))] as in (4.53) yields ^[vo(t)] — [T][I(t)] = [R(t)][I(t)] + [-Ek(t )]  ^  (4.63)  where [V0 (t)] = [V,;(t)] — [ Eho(t) ] . Casting (4.63) in the form of a fixed-point iteration [/(0] = [P(/(t))], we define ^[F'(I (t))] = [T R(t)] -1 [Vo (t) — E k(t)]^(4.64)  where the iteration  [i(z+ 1) (01^[F'([(a )(0)]  ^  (4.65)  converges to a unique fixed-point [P(t)] such that [P(t)] = [F 1 (I* (t))]. [I* (0] represents the intersection of the linear network Thevenin equivalent and the nonlinear branch characteristic as depicted in Fig 4.8(a) for the scalar case. The fixed-point iteration is 'seeded' with the solution from the previous time-step, namely [/( 3 )(t)] = [P(t — At)]. The nonlinear problem has been reduced to finding the segment of the piecewise-linear branch characteristics that uniquely solves the network at each time-step. We now discuss the convergence and divergence properties of the fixed-point iteration. For [/( ° )(t)] 26  sufficiently close to [/*(t)] the iterates will converge quadratically upon the solution. Figures 4.8(b) and (c) depict typical convergence behaviour in the scalar case. Divergence is characterized by a cyclic path among the piecewise-linear segments [8]. Fig. 4.8(d) depicts a typical divergence mechanism in the scalar case. Note the cyclic path between segments 1 and 3 — the iterates will never converge to the solution which lies on segment 2. An anologous divergence mechanism occurs in the multidimensional case among piecewise-linear hyperplanes (as opposed to segments). Cyclic paths must be detected and corrective action taken. The proposed fixed-point iteration algorithm for the solution of the currents in the nonlinear branches is based on the following paradigms: 1. Convergence to the solution is achieved if two successive iterates visit the same piecewise-linear segment (hyperplane). 2. Divergence is detected if two non-successive iterates visit the same piecewise-linear segment (hyperplane). Corrective action consists of 're-seeding' the iteration on a piecewise-linear segment (hyperplane) that lies between the two iterated segments (hyperplanes).  27  Figure 4.8 Convergence and divergence behaviour of fixed-point iteration  The unknown node voltages associated with the linear part of the network are found by  back-substitution. Such a scheme can be readily incorporated into the algorithm presented in Section 3.3 as presented in Algorithm 4.3.  nonlinear_^_ solve code (  t 4— tO, [A]^[AO], [Eh] <— [Eho]  while t < t f , build_rhs( compute([17/3], [J.7 ]) [rhs]^[AhA][Jh] + [A3 A]r- Ji  GAB] [VB]  Algorithm 4.3 Transient solve code for nonlinear networks (Continued ... )  28  solve_ compensation _ currents( [JA]  4—  {A /hA][Jh] + [A3 A][JI]  [V0] = [GB][VB] [JA] — { Eoh [I, R, Ek] <—fixed_point_iteration(/ , T, V0 ) [rhs]^[rhs] — [B A ] [I] )  solve (  [ VA ] 4— backsub([GAA], [rhs]) update_history_terms ( [Jh] 4— [Uh] • ([Jh] + 2[Gh][A„][V]) [Eh] +-- [ — U 0 ] ([A71][ 11^[R][I]^[Ek]) t^t + At  end while  Algorithm 4.3 Transient solve code for nonlinear networks  4.6 Numerical example: Transient analysis Consider the simple network of Section 3.5 with the linear inductance L n, replaced by a piecewise-linear inductance L 7, defined by {318.3 mH, i < — isat I, = 1.591 mH, — i sat < l < isat  318.3 mH, i^ -I > i-sat where i sa t  =  0.01 A.  1  (4.66)  Formulation of the network equations: With nodes defined as before and the Norton  equivalent representation of L m with history source jb, replaced by the Thevenin equivalent 29  representation of L i, with history source e hi , the discrete-time network equations become  [ g R, + gL3^ —  —  g L,^gL,  gLs  + gc  0^1  0^vi  (01^  [-YR, = JA(t) 0 l[v3(t)] 1^v2(t) lv(ii(t)) 0^i1(t)^ 0  (4.67)  where  1  [JA(t)] - [^ (t ) (4.68) • WO^ ex (t)^ehi(t)^ v(i1(t)) = r(t) )  ehl(t) = —(v2(t — At) + r(t — At) • i i (t — At) + ek(t —  At))  and Norton equivalent history sources updated recursively. Network equivalents For the example at hand with At = 1.25 ms downward operations  transform the augmented matrix (4.60) to [ 249.9057 0  —4.8077 7.2752  0  0  —245.0980 —4.7152 0.6481  0 1 —0.1375  0 —1 0.1375  —1 0.9808 —0.1348  (4.69)  Solution method: Fig. 4.9 represents the progression of the voltage and currrent waveforms  v 2 (t) and i1 (t) from zero initial conditions. Note that even after ten cycles the solution has not yet settled to the periodic steady-state.  30  2 v 2 (t) 1 0 -1 -2  i  0.02  ^  0.04  ^  0.06  ^  0.08  ^  0.1  ^  0.12  ^  0.14  ^  0.16  ^  0.18^0 2 t (s)  2.5 2 1.5  0.5 0  0.02^0.04^0.06^0.08^0.1^0.12^0.14^0.16^0.18^0 2 t (s) Figure 4.9 Solution of sample network from zero initial conditions  4.7 Steady-state initialization  The Norton and Thevenin equivalent history sources represent the state of a discrete-time network with current-controlled nonlinearities defined by[h(t)] = A(t)  Analogous to  [^  Eh(t)]  .  the development of Section 3.6, define [F(ho)] = [Jh(T; ho)]^  (4.70)  as the state at time-step T on the solution trajectory, which is uniquely defined by [ho] and computed by simulating the network over one period. A nonlinear network with sources of period T may exhibit periodic, subharmonic, quasiperiodic, or even chaotic solutions [6]. The periodic steady state is, however, of particular interest. The discrete-time network has reached the periodic steady state when [F(hN)] = [ kJ]  ^  is solved uniquely by [4], the steady-state initial conditions. 31  (4.71)  The two-point boundary value problem is solved by the Newton-Raphson iteration [h,(ji+1)]^[4 j) ]  — [DF (h o( i) ) — U1 -1 [F (h o( i) ) — Ul^(4.72)  This approach yields quadratic convergence to^at the additional expense of calculating the Jacobian matrix [DF(h o( i) )1 at each iteration.  4.8 Computing the Jacobian Matrix The proposed approach requires the efficient and accurate computation of the Jacobian matrix aF(ho)  (4.73)  [DF(h0)] [ a 81(4 ) a E h0  Along the same lines as Section 3.7, the Jacobian matrix is computed by the solution of n linear, time-varying sensitivity networks related to the variational network8 [G AA^BA  r ava,o(t)  AhAa  aE(t)^A' h Jh (t) 0^—T — R(t) . a I (t) ^ah o^aho^'hA ah o  (4.74)  —  from t = 0 to t = T with history sources updated as follows [  a 4(o) ahoUi 0 a E h( 0) _ 0 U2  t = 0^(4.75)  ah o  ravAm i) Iwo +^aaJh(01 +2[Gh][AvAi ^ _ [Uh] ah o aho J) fOEh(t Ain I_^Oho^j  0] ([AnA]  [OV aAh(t)]  ^[R(t)]  [0a/i(t0)1)  (4.76)  (4.77)  Thus the solution of the sensitivity networks is simply the matrix-valued solution of the variational network with independent sources and 'knee' voltage sources zeroed. 8  ^  The variational networks follows by substituting (4.53) in (4.61).  32  4.9 Steady-state initialization algorithm Algorithm 4.4 depicts the proposed steady-state initialization algorithm for networks with current-controlled nonlinearities in pseudocode form.  steady-state_init( while [A] _ [Jho > €, E rilzi 2 t 4- 0 tf^T, [Jh]^[Jho], [Eh]^[Eho] ,  ] [DJh] 4- [ Ul 0 ], [ DEh^ [ 0 U2 while t < tf, build rhs( compute([VB], [JI ]) [rhs]^[AhA][Jh] + [ A7A][J2] — [GAB][VB] [Drhs]^[AhA][DA] solve compensation currents( [JA] <— [A lhA] [Jh] + [ A7A][J7] [Vid = [CB][VB] [JA] — [ oh ]  oh] [DVo ] 4— —[A'hA ][DJh]—[ D E [I, R, Ek] <—fixed_pointiteration(I ,T , V0 ) [DI] 4-- [T R] -1 • [DVo ] [rhs]^[rhs] — [B A ][I] [Drhs]^[Drhs] — [B A ] [DI]  solve ( [VA] 4— backsubUGAAL [rhs]) [DVA ]^backsub^AA ], [Drhs]) Algorithm 4.4 Time-domain steady-state initialization algorithm for nonlinear networks (Continued^)  33  update_history_terms (  [A]^Rid ([A] + 2 [Gh][Av][ 17]) [DJh]  Rid ([DJh] + 2 1GOAvAHDVAD  [Eh] <-- [ — U 0] ([An] [V] + [R][ 1 ] + [Ek]) [DEh] < [ —  —  U 0] GAnit][DVA] + [R][DI])  t^t At  end_while newton_raphson_iteration( -1 ([D [AO] [AO D Eh ho^ Eho  0u [U1 E 0^U2])  (1A^ 1^1 ,jol) uEhi LEho_u  end_while  Algorithm 4.4 Time-domain steady-state initialization algorithm for nonlinear networks  4.10 Numerical example: Steady-state initialization  For the example of Section 4.6, Fig. 4.10 and Fig. 4.11 show the progression of the steadystate initialization algorithm from zero initial conditions and from the steady-state initial conditions of the network linearized about [I] = [0]. Note that only 5 and 2 iterations, respectively, were required for convergence to the periodic steady-state. 34  0.01^0.02^0.03^0.04^0.05^0.06^0.07  0.08^0.09^0.1 t (s)  Figure 4.10 Progression of the steady-state initialization algorithm from zero initial conditions  0.4 i 1 (t) 0.2  -0.2 -0.4 ^ 0^0.005^0.01^0.015^0.02^0.025^0.03^0.035  0.04 t (s)  Figure 4.11 Progression of the steady-state initialization algorithm from linearized network initial conditions 35  ^  4.11 Application: Harmonics due to Transformer Saturation The example of Section 4.6 is a simplified representation of a nonlinear transformer terminating an open-circuited line with the nonlinear magnetizing characteristic defined piecewise-linear [10]. Harmonics due to transformer saturation are generally difficult to predict as they are highly dependent on the magnitude of the applied voltage. Figure 4.12 shows the nonlinear inductor voltage and current waveforms for At =  312.5ps  obtained by the steady-state initialization algorithm. Table 4.4 depicts the corresponding phasor representation of the harmonics of interest. Note that the choice of At imposes a bandwidth constraint on the harmonic content of the steady-state waveform. 2  v  lo  1  0 -1 -2  0.002 0.004 0.006 0.008^0.01  0.012 0.014 0.016 0.018^0.02 t (s)  0.4  i (r) 1  0.2  -0.2 -0.4 ^^ 0.002 0.004 0.006 0.008^0.01^0.012 0.014 0.016 0.018^0.02 0  t (s) Figure 4.12 Steady-state waveform for the sample network (At = 312.5p.^)  Harmonic number 1  Voltage Harmonic (V rms)  Current Harmonic (A rms)  7.9129e — 1L — 90.10°  6.8418e — 2L — 179.82°  3  7.8842e — 3L88.36°  5.1181e — 2L — 179.31°  Table 4.4 Phasor representation of voltage and current harmonics in nonlinear inductor (At = 312.5its)^(Continued ... )  36  3.7372e — 2L — 178.68°  5  1.6733e  7  6.9259e — 2L — 79.00°  2.1392e — 2L — 177.57°  9  3.3725e — 3L — 83.86°  7.7646e — 3L — 174.51°  11  2.0951e — 4L55.97°  8.6580e — 4L — 34.26°  13  6.3113e — 4L86.79°  3.7703e — 3L — 3.3157°  2L88.91°  Table 4.4 Phasor representation of voltage and current harmonics in nonlinear inductor (At = 312.5p,^)  37  ^  Chapter 5 Steady-state initialization: Variable-structure networks 5.1 Introduction A general method for the formulation, solution, and steady-state initialization of variable structure networks 9 is presented in this chapter. 5.2 Representation of ideal-switch devices An ideal-switch device is a two state switch governed by a set of Boolean equations of the terminal characteristics. The present discussion is restricted to an ideal diode. Define { 0, switch off  f  s^  1, switch on  (5.78)  as the state of an ideal diode. Furthermore, define the boolean variables sv =  { 1, v > 0  s  0, otherwise  { 1, i > 0 0, otherwise  (5.79)  where v represents the anode-to-cathode voltage and i represents the anode-to-cathode current. The switch transitions, governed by changes in branch current and voltage are succinctly expressed in Fig. 5.13.  Figure 5.13 State transition diagram for an ideal diode 9^Variable structure networks are those containing ideal-switch devices, such as ideal diodes.  38  5.3 Formulation of the network equations Ideal-switch devices are incorporated into the nodal formulation of the network equations by a compensation method as follows  G Bi 1 V (01 _ rAtri Ls.c^I(t) —^0  (5.80)  where [I(t)] is a vector of switch branch currents and  [S] = diag([ • • •^33^• • •]) [S] =  ^  diag([ • • • Sj • • •])  (5.81)  are diagonal matrices of the switch state (and its conjugate) introduced in the previous section. Again partitioning the vector of node voltages into unknown and known voltages yields [CAA B1[17,4(1 [JAM] [GAB j CB [VB (0]^(5.82) S • CA^ 37^ I(t)^S •  For variable structure networks, the absolutely stable backward Euler rule of integration is employed with the time-step halved [3]. The Norton equivalent history sources are thus updated as follows  [Jh(t + At)] = [Uh][Gh][A,][V(t)] + 2[U + Uh][Jh(t)]^(5.83) where [Uh], [Gh] and [71„] are previously defined.  5.4 Network equivalents Akin to Section 4.4, [S • CA ] is eliminated as follows  [  GAA^GAB^BA^AhA 0^CB^—T*^AhA  1°  P^0 ][GAA^GAB^BA —K*^U^S • CA^S • CB^—S  [  to^In practice the network equivalent is synthesized by downward operations.  39  AjAl  A;A  (5.84)  AhA  0  0 AjAl  where [P] represents the permutation matrix that upper triangularizes [GAA] and [K*] = [S. CA][GAA]. Note that the network equivalent is dependent on the switch state N. Thus  the network equivalent must be re-synthesized after each switch transition n . We now have network equations of the form A B A ][VA(t)1^[^(t)^[G A*131 rT7B(t)] (t)^—AM^CB Lv —T* 1^J [GA  0  (5.85)  where [JAM] = [AI A ][Jh(t)] [A; A ][Ji (t)]. Note that the ideal-switch branches are now coupled to a pseudo-Thevenin equivalent of the linear part of the network [Vo* (t)] — [T*][I(t)] = [0]^ (5.86)  where [V,,*(t)] = [C13,1[VB(t)] [JAM] represents the open circuit voltage. Eq. (5.85) is dubbed a pseudo-Thevenin equivalent since the solution [/(t)] is a vector containing short circuit currents and zero currents depending on the switch states. 5.5 Solution method  A straightforward non-iterative solution method simply solves the network with the switch states computed at the previous time-step. This is consistent with the existing solution method in the EMTP. A new network equivalent is synthesized only if a switch transition occurs. variable structure solve code( -  t  4-  to, [Jh] 4- [AO]  while t < tf, build rhs( compute([/B1, [Ji ]) [rhs]^[AhA] [Jh] + [A^— [G A 8][VB] Algorithm 5.5 Transient solve code for variable structure networks (Continued ... ) it^  [G AA ], however, need only be triangularized once prior to the first time-step.  40  solve_compensation_currents( [J:4] <— [A;;A][Jh + [A;A][`Iil [Vol^[ CI] [VB] + [JA] [I]^[T 1 -1 • [ Vol [rhs] 4-- [rhs] — [B A] [I] )  solve( [VA] 4— backsubaGAA1, [rhs]) )  switch transition( [S, transition]^switchflA n • [ 17], [I])  if transition, [CB, T* , AhA A; A ] 4— net_equiv(GAA , GAB, BA, AhA, AjA, 8) ,  end )  up date_history_terms(  [A]^Rid [Gh] [Av] [V] + 0.5[U + Uh] ) t^t  end_while  Algorithm 5.5 Transient solve code for variable structure networks  5.6 Numerical example: Transient analysis Consider the single-phase half-bridge rectifier with freewheeling diode shown in Fig. 5.14(a). Fig. 5.14(b) depicts the corresponding discrete-time network with the diodes represented as ideal switches with states si and s2. 41  E (7Z)=10 . iao 0020rr rr 1) R— --1  50)  Z. =  -  zr  02  ms le 9k= Q 5R, /00 yz..=^6041c2 e— zic A /01/7  (6  )  Figure 5.14 Half-bridge rectifier with freewheeling diode (a) continuous-time network (b) discrete-time network  Formulation of the network equations: With nodes and history sources defined as in Fig.  5.14(b) the discrete-time network equations become gL  + gR  —g L  0 0  —  gL  gL  0 —81 —8 2  0  0  v i (t )  0 YR .  —1 0  —1 1  0  sl  0  82  0  S2  v 2 (t) v3(t) ii(t) i2(t)  0  ih  =  (  0 0 0  i  )  0 0 0 81 0  [v4(t)]  (5.87)  where j h (t) is updated as in (5.83). Triangularization and network equivalents: The initial switch states s i = s 2 = 0 are  assumed thus only [CAA] need be upper triangularized prior to the first time-step. For the example at hand with At = 1.25ms and switches initially open downward operations 42  transform the augmented matrix (5.84) to 0.2003^—0.0003^0 0^0.0003^0 0^0^100.0 0^0^0 0^0^0  0 0 0 0 0  0 —1 0 1 0  0 —1 1 0 1  1 —0.9987 0 0 0  (5.88)  Solution method: Algorithm 5.5 with parameters to = 0 and tf = 200ms is applied to  the sample network. Fig. 5.15 represents the progression of the load current and node 2 voltage waveforms iL(t) and v2(t) from zero initial conditions. Note that even after twelve cycles the solution has not yet settled to the periodic steady-state.  it)  5 4  3 2  0^0.02^0.04^0.06^0.08^0.1^0.12^0.14^0.16^0.18^0 2 t (s) 200 v 2 (t) 100  -100  0^0.02^0.04^0.06^0.08^0.1^0.12^0.14  ^  Figure 5.15 Solution of sample network from zero initial conditions  43  0.16  ^  0.18^02 t (s)  5.7 Steady-state initialization algorithm Based on the development introduced in Chapter 3, a steady-state initialization algorithm is implemented within the variable structure network solve code as in Algorithm 5.6. steady_state_init( t^0 ,t f 4  —  T, [A]^[0], [D,Ih] 4— [U]  while t < t f , build rhs( compute([VB], [JA) [rhs]^[AhA] [hil + [ A  ] [ Ji] —  [GAB] [17.13]  [Drhs]^[Aha][-DA] solve compensation currents(  [JA]^[AIA][Jhl + [A;A][Ji] [vol^ri3l[vB] + [I] 4— [ 71 -1 ' [ Vol [DI] 4-- [T*] -1 [AI A ][DJh] [rhs]^[rhs] — [B A ][I] [Drhs]^[Drhs] — [B A ][DI]  solve( [VA] 4—  backsub([GAA], [rhs])  [DVA] 4— backsub([GAA], [Drhs])  switch transition( [S,transition] 4— switchGA,2 ] • [V], [I])  if transition, [C73,T*, A hA, A; A ] +— net_equiv(G A A , GAB, BA, AhA, AjA, S) Algorithm 5.6 Time-domain steady-state initialization algorithm for variable structure networks (Continued 44  )  end update_history_terms( [Ai^[Uh][Gh][Av][V] + 0.5[U + Uh]{A} [DA] 4-- [Uh][Gh][AvADVA] + 0.5[U + Uh}[DJh} )  t 4- t  end while newton_raphson_iteration( [Ai  4-  [DA - U] -1 • [ Al -  Algorithm 5.6 Time-domain steady-state initialization algorithm for variable structure networks  5.8 Numerical example: Steady-state initialization  Fig. 5.16 represents the solution from the steady-state initial conditions computed by the time-domain method. Note the scaling of the axis in the first plot. The waveform represents a dc current with ripple. 45  10.7 ^  i L(t) 10.6 10.5 10.4  0  0.002^0.004^0.006^0.008^0.01^0.012^0.014^0.016^0.018  t (s)  200 ^ v2 (t) 150 100 50 0 -50 ^ 0  0.002^0.004^0.006^0.008^0.01^0.012^0.014^0.016^0.018  t (s)  Figure 5.16 Solution of sample network from steady-state initial conditions  5.9 Application: Initialization of single and three-phase bridge networks The following are sample cases from [111. Figure 5.17 depicts a single-phase diode bridge  rectifier. Figure 5.18 shows the steady-state waveforms obtained from Algorithm 5.6 for At = 130.2083ps. Figure 5.20 depicts the steady-state waveforms for the three-phase diode bridge rectifier (Fig. 5.19). Note that there is no frequency-domain method for the solution of such networks. 46  a 0-n.  Figure 5.17 Single-phase diode bridge rectifier  200  -50 -100 -150 -200 ^ 0  0.002^0.004^0.006^0.008^0.01^0.012^0.014  0.016^0.018  t (s)  Figure 5.18 Steady-state waveforms for single-phase diode bridge rectifier  47  Figure 5.19 Three-phase diode bridge rectifier  Figure 5.20 Steady-state waveforms for three-phase diode bridge rectifier  Chapter 6 Conclusions 6.1 Summary A time-domain based steady-state initialization algorithm implementable within the solve code typical of electromagnetic transient programs (EMTP) has been proposed. The algorithm is developed for the three classes of lumped element networks: linear, nonlinear, and variable structure. The primary contribution of this work is that each network element does not require a dedicated steady-state model and a separate steady-state solution algorithm. The proposed algorithm can be implemented in existing electromagnetic transients programs with minimal programming effort as the steady-state initialization process is consistent with the timedomain character of the transient simulations themselves. Future efforts might extend the time-domain method to handle the initilization of distributed parameter elements (line models), TACS models, and machine models.  49  References [1] H. Dommel, "Digital computer solution of electromagnetic transients in single and multiphase networks," IEEE Transactions on Power Apparatus and Systems, vol. PAS88, no. 4, pp. 388-399, April 1969. [2] H. Dommel, "Nonlinear and time-varying elements in digital simulation of electromagnetic transients," IEEE Transactions on Power Apparatus and Systems, vol. PAS-90, no. 6, pp. 2561-2567, Nov./Dec. 1971. [3] J. R. Marti and J. Lin, "Supression of numerical oscilations in the emtp," IEEE Transactions on Power Systems, vol. PAS-4, no. 2, pp. 739-747, May 1989.  [4] H. D. A. Yan and S. Wei, "Harmonics from transformer saturation," IEEE Transactions on Power Delivery, vol. PWRD-1, no. 2, pp. 209-214, April 1986.  [5] A. S. E. Acha and J. Arrilaga, "Newton-type algorithms for the harmonic analysis of nonlinear power circuits in periodical steady state with special reference to magnetic nonlinearities," IEEE Transactions on Power Delivery, vol. PWRD-3, no. 3, pp. 10901098, July 1988. [6] T. Parker and L. Chua, Practical Numerical Algorithms for Chaotic Systems. New York: Springer-Verlag, 1989. [7] T. Aprille and T. Trick, "Steady-state analysis of nonlinear circuits with periodic inputs," IEEE Proceedings, vol. 60, no. 1, pp. 108-114, January 1972.  [8] L. Chua and P.-M. Lin, Computer-aided Analysis of Electronic Circuits. Englewood Cliffs, New Jersey: Prentice Hall, 1975. [9] A. R. Chung-Wen Ho and P. Brennen, "The modified nodal approach to network analysis," IEEE Transactions on Circuits and Systems, vol. CAS-22, no. 6, pp. 504— 509, June 1975. [10]J. Usaola and J. G. Mayordomo, "Fast steady-state technique for harmonic analysis," IEEE Transactions on Power Delivery, vol. PWD-6, no. 4, October 1991. 50  [1 1}N. Mohan, Computer Exercises for Power Electronics Education. Minneapolis, MN: University of Minnesota, 1990.  51  

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

Usage Statistics

Share

Embed

Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                        
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            src="{[{embed.src}]}"
                            data-item="{[{embed.item}]}"
                            data-collection="{[{embed.collection}]}"
                            data-metadata="{[{embed.showMetadata}]}"
                            data-width="{[{embed.width}]}"
                            async >
                            </script>
                            </div>
                        
                    
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0064915/manifest

Comment

Related Items