# Mathematical Modeling of Transient Responses in a Large Scale Multiconductor System

Authors: Caiwang Sheng; Xiaoqing Zhang
Journal: Journal of Applied Mathematics
 (2014)

---

## Abstract

This paper proposes a reduced-order model for the large scale circuits representing the multiconductor systems. It is based on the block Arnoldi algorithm for calculating congruence transformation matrix. By setting up the state equations in frequency domain for the multiconductor systems, the transfer functions are calculated and so the frequency response curves are obtained. A comparison is made between those curves obtained from the circuits with and without reduced-order treatment and a better agreement appears between them. By using inverse Laplace transform, the lightning transient responses in the multiconductor systems are given in time-domain. A better agreement is shown between calculated and experimental results, which confirm the validity of the proposed model.


---

## Body

## 1. Introduction

An external lightning protection system of a structure is usually constituted by numerous interconnected longitudinal and transverse conducting branches. It takes a form of three-dimensional multiconductor system and can be converted into an equivalent passive RLC network. For an actual tall structure, its external lightning protection system is large in physical dimension and complex in construction, which results in a large scale equivalent network. From the view point of lightning protection design, the need exists for obtaining lightning transient responses of the equivalent network. For the sake of saving storage capacity and computation time, a simplification should be made for the equivalent network; otherwise the circuit calculation is difficult to perform. From the existing research on lightning transient responses [1–6], we have almost not found the relevant procedure that can give a great simplification for the large scale equivalent network. Therefore, the order reduction model is necessary for the transient analysis of lightning protection systems. Recently, the investigation on passive network synthesis has been reported in literature [7–10], which is closely related to the order reduction model. The investigation work is of practical significance in analysis of large scale electrical networks. Based on the achievement of the existing work, a reduced-order algorithm (ROA) is introduced in this paper as a Krylov space method which provides a passive, stable, and accurate macromodel compared with the reduction techniques [11, 12]. In order to reduce the matrices directly in the system state equation, the transformation matrix can be calculated by using the block Arnoldi algorithm [13]. In the ROA, we have developed a robust framework by unveiling the connections between the passive reduction technique and other well known Krylov space based processes. In terms of the reduced-order model of the structure under study, amplitude-frequency and phase-frequency curves can be quantitatively calculated and so we can compare those curves with the frequency response curves of the original circuit system. Furthermore, the time function can be obtained with the help of the inverse Laplace transform and the lightning transient response is given by using the convolution theorem. The comparison is also made between calculated and measured results to confirm the validity of the proposed model.


## 2. State Equations of the Equivalent Circuit

Figure1 shows the configuration of a cage-like multiconductor system formed by vertical and horizontal conducting branches. For setting up the state equations of the equivalent network, the electrical parameters of the conducting branches in the multiconductor system are represented by resistances, inductances, and capacitances. Since the distance between the two conducting branches is much greater than their respective radius, the electromagnetic couplings between them can be ignored. Therefore, the formulas for evaluating the electrical parameters are derived as follows.Figure 1
A cage-like multiconductor system.Consider a horizontal conducting branch in the multiconductor system, as shown in Figure2. Its capacitance can be expressed by [14]
(1)Cj=2πε0ljln⁡2hj/a-D1,
where D1=ln⁡1+1+4hj2/lj2+2hj/lj-1+4hj2/lj2+0.307.Figure 2
A horizontal conducting branch.Similarly, the capacitance of a vertical conducting branch, as shown in Figure3, can be expressed by [14]
(2)Ck=2πε0lkln⁡lk/a-D2,
where the coefficient D2 is
(3)D2=1+hklkln⁡4hklk+1+hklkln⁡1-hklk-1+2hklkln⁡1+2hklkhklk<1,0.307+1+hklkln⁡1+lkhk-1+2hklkln⁡1+lk2hkhklk>1.
According to the electromagnetic analogy [15], the product of capacitance and inductance is equal to that of space permeability μ0 and permittivity ε0. This gives the inductances for horizontal and vertical conducting branches
(4)Lj=μ0ε0Cj,Lk=μ0ε0Ck.
Additionally, the conductance of an arbitrary conducting branch is estimated by
(5)Gξ=1ρ·πa2lξ,
where ρ is the resistivity of conducting branch and the subscript ξ takes as j or k.Figure 3
A vertical conducting branch.In consideration of the propagation phenomenon of lightning current, each conducting branch of the multiconductor system is subdivided into a suitable number of segments. The length of each segment is assumed to be less than or equal to one-tenth of the wavelength corresponding to the maximum frequency likely to affect the system transient [16]. Each segment is represented by a π-type circuit, as shown in Figure 4, where M and N are the segment numbers of horizontal and vertical branches, respectively.Figure 4
Aπ-type circuit of a two-segment subsystem.As a result, the multiconductor system is converted into an equivalent network constituting a large number ofπ-type circuits. On the basis of the electrical parameters calculated from (1)–(5), the capacitance, conductance, and inductance matrices can be formed for the equivalent network:
(6)C=diag⁡C1,C2…,CnC,G=diag⁡G1,G2…,GnG,L=diag⁡L1,L2…,LnL,
where the subscripts nC, nG, and nL denote the total numbers of capacitances, conductances, and inductances of the equivalent network, respectively. The circuit components of the equivalent network are numbered in the sequence of capacitances, inductances, and lightning current source. In accordance with the component number, the incidence matrix of the equivalent network is written as
(7)A=AC,AL,AS.
The node voltages (i.e., the voltages across the capacitances) and the currents through the inductances are selected as the state variables:
(8)x=un1,un2…,unβ,iL1,iL2,…,iLnLT,
where the subscript β denotes the total node number of the equivalent network. Thus, the state equation is set up for the equivalent network:
(9)Px˙+Qx=Bist,uot=BTx,
where
(10)P=C^00L,Q=G^AL-ALT0,B=As0,
where C^=ACCACT and G^=ALGALT∈Rβ×β. Note that C^, L, and G^ are symmetric and nondefinite matrices.Next, the so-called transfer function can be calculated, which describes the input-output behavior of a linear dynamical system(9) in frequency domain [17]. Taking the Laplace transformation of (9), the transfer function is given by
(11)Hs=BTQ+sP-1B.The computation ofH(s) is called simulation in frequency domain and of course not done by calculating Q+sP-1. However, it is still expensive for large system. For this reason, the reduced-order model is necessarily proposed in this paper.


## 3. Reduced-Order Model

For simplifying the lightning transient calculation of the large scale RLC network, the ROA is employed here. The algorithm is based on the block Arnoldi method. Therefore, we introduce our notion of block Krylov subspaces for multiple starting vectors.


### 3.1. Model Order Reduction Techniques

A widely used class of model order reduction techniques is the class of Krylov subspace methods. Before block Krylov subspace methods are used, the two matricesQ and P in (11) have to be reduced, denoted by D in the following. This can be done by rewriting (11) as
(12)Hs=BTI+sD-1W,
where
(13)D=Q-1P,W=Q-1B.
Although Q and P are sparse matrices, in general, the matrix D is a dense matrix. However, block Krylov subspace methods involve D only in the form of matrix-vector products Dv and possibly DTw [18].The proper definition of block Krylov subspaces is necessary to be involved and the use of these subspaces results in much more efficient reduced-order modeling techniques. The block Krylov space generated by matricesD∈Rβ+nL×β+nL and W=w1w2⋯wζ∈Rβ+nL×ζ is defined as
(14)KrD,W,q=colspW,DW,D2W,…,DqW.


### 3.2. The Algorithm of Calculating Congruence Transformation Matrix

The congruence transformation matrixV is calculated based on the above definition of Kr(D,W,q). Here, V=v1v2⋯vq∈Rβ+nL×q is an orthonormal matrix and subjected to VTV=I (I∈Rq×q is an identity matrix) [19]. The classical Arnoldi process [20] generates orthonormal basis vectors for the sequence of Krylov subspaces KrD,W,q, q≥1, induced by D and W. In order to saving storage capacity and computation time, a smaller value of the order q is better. In fact, exact deflation at step q of the Arnoldi-type process occurs if and only if v~q=0. Similarly, inexact deflation occurs if and only if v~q≈0, but v~q≠0. Therefore, the order q will be determined according to v~q⩽dtol, where dtol⩾0 is a suitably chosen deflation tolerance. A simple Arnoldi-type algorithm is presented in the following [18].(0)
Setv~i=wi (i=1,2,…,m).Setmc=m.Forn=1,2,…,k do the following:(1)
computev~q and check if the deflation criterion (v~q⩽dtol) is fulfilled.If yes,v~q is deflated by doing the following.Setmc=mc-1. If mc=0, set q=q-1 and stop.Setv~i=v~i+1 for i=q,q+1,…,q+mc-1.Return to step(1).(2)
Sethq,q-mc=v~q and vq=v~q/hq,q-mc.(3)
Computevq+mc=Dvq.(4)
Fori=1,2,…,q do the following:sethi,q=viTv~q+mc and v~q+mc=v~q+mc-vihi,q.(5)
Fori=q-mc+1,q-mc+2,…,q-1 do the following:sethq,i=vqTv~i+mc and v~i+mc=v~i+mc-vqhq,i.


### 3.3. A Reduced-Order Algorithm

After computing the transformation matrixV, we can easily find the ROA matrices. Applying the change of variable x=Vx~ in (9) and multiplying both sides by VT give
(15)VTPVx~˙+VTQVx~=VTBistuot=BTVx~.
For the macromodel, the reduced-order matrices are
(16)P~=VTPV,Q~=VTQV,B~=VTB.These types of transformations are known as congruence transformations. From(15) and (16), the reduced-order transfer function is given by
(17)H~s=B~TQ~+sP~-1B~.
Define
(18)D~=Q~-1P~,W~=Q~-1B~.
The reduced system can be formulated as
(19)D~x~˙+x~=W~istuot=B~Tx~.By applying any eigendecomposition routine onD~, we can obtain
(20)D~=JΛJ-1,
where Λ is a diagonal matrix. Since D~ is a real matrix, any complex eigenvalue or eigenvector has its conjugate. Therefore, we can use real J and Λ matrices in the following transformations:
(21)Jr=JK-1,Λr=KΛK-1,
where
(22)D~=JrΛrJr-1,
and Jr,Λr∈Rq×q.  K∈Cq×q is the transformation matrix that can be defined as follows:ifΛi,i is real, Ki,i=1;ifΛi,i=Λ-i+1,i+1, Ki,iKi,i+1Ki+1,iKi+1,i+1=11j-j.Λr will contain 2 × 2 blocks on the diagonal for complex eigenvalues, and 1 × 1 blocks for the real eigenvalues. By substituting (22) into (19), we can get
(23)JrΛrJr-1x~˙+x~=W~ist,uot=B~Tx~.
The final transfer function is given by
(24)H~s=B~TJr-1Λr-1s+Iq-1JrW~.
As a result, we can obtain the final reduced-order model in the following form:
(25)H~ijs=∑m=1qK~ijms+p~m+d=bijksq+⋯+bij1s+bij0sq+ak-1sq-1+⋯a1s+a0,
where H~ijs=H~jis, 1≤i, j≤N, p~ and K~ are the dominant pole and corresponding residue, respectively, q is the number of reduce-order, and d represents any direct coupling between ports.


## 3.1. Model Order Reduction Techniques

A widely used class of model order reduction techniques is the class of Krylov subspace methods. Before block Krylov subspace methods are used, the two matricesQ and P in (11) have to be reduced, denoted by D in the following. This can be done by rewriting (11) as
(12)Hs=BTI+sD-1W,
where
(13)D=Q-1P,W=Q-1B.
Although Q and P are sparse matrices, in general, the matrix D is a dense matrix. However, block Krylov subspace methods involve D only in the form of matrix-vector products Dv and possibly DTw [18].The proper definition of block Krylov subspaces is necessary to be involved and the use of these subspaces results in much more efficient reduced-order modeling techniques. The block Krylov space generated by matricesD∈Rβ+nL×β+nL and W=w1w2⋯wζ∈Rβ+nL×ζ is defined as
(14)KrD,W,q=colspW,DW,D2W,…,DqW.


## 3.2. The Algorithm of Calculating Congruence Transformation Matrix

The congruence transformation matrixV is calculated based on the above definition of Kr(D,W,q). Here, V=v1v2⋯vq∈Rβ+nL×q is an orthonormal matrix and subjected to VTV=I (I∈Rq×q is an identity matrix) [19]. The classical Arnoldi process [20] generates orthonormal basis vectors for the sequence of Krylov subspaces KrD,W,q, q≥1, induced by D and W. In order to saving storage capacity and computation time, a smaller value of the order q is better. In fact, exact deflation at step q of the Arnoldi-type process occurs if and only if v~q=0. Similarly, inexact deflation occurs if and only if v~q≈0, but v~q≠0. Therefore, the order q will be determined according to v~q⩽dtol, where dtol⩾0 is a suitably chosen deflation tolerance. A simple Arnoldi-type algorithm is presented in the following [18].(0)
Setv~i=wi (i=1,2,…,m).Setmc=m.Forn=1,2,…,k do the following:(1)
computev~q and check if the deflation criterion (v~q⩽dtol) is fulfilled.If yes,v~q is deflated by doing the following.Setmc=mc-1. If mc=0, set q=q-1 and stop.Setv~i=v~i+1 for i=q,q+1,…,q+mc-1.Return to step(1).(2)
Sethq,q-mc=v~q and vq=v~q/hq,q-mc.(3)
Computevq+mc=Dvq.(4)
Fori=1,2,…,q do the following:sethi,q=viTv~q+mc and v~q+mc=v~q+mc-vihi,q.(5)
Fori=q-mc+1,q-mc+2,…,q-1 do the following:sethq,i=vqTv~i+mc and v~i+mc=v~i+mc-vqhq,i.


## 3.3. A Reduced-Order Algorithm

After computing the transformation matrixV, we can easily find the ROA matrices. Applying the change of variable x=Vx~ in (9) and multiplying both sides by VT give
(15)VTPVx~˙+VTQVx~=VTBistuot=BTVx~.
For the macromodel, the reduced-order matrices are
(16)P~=VTPV,Q~=VTQV,B~=VTB.These types of transformations are known as congruence transformations. From(15) and (16), the reduced-order transfer function is given by
(17)H~s=B~TQ~+sP~-1B~.
Define
(18)D~=Q~-1P~,W~=Q~-1B~.
The reduced system can be formulated as
(19)D~x~˙+x~=W~istuot=B~Tx~.By applying any eigendecomposition routine onD~, we can obtain
(20)D~=JΛJ-1,
where Λ is a diagonal matrix. Since D~ is a real matrix, any complex eigenvalue or eigenvector has its conjugate. Therefore, we can use real J and Λ matrices in the following transformations:
(21)Jr=JK-1,Λr=KΛK-1,
where
(22)D~=JrΛrJr-1,
and Jr,Λr∈Rq×q.  K∈Cq×q is the transformation matrix that can be defined as follows:ifΛi,i is real, Ki,i=1;ifΛi,i=Λ-i+1,i+1, Ki,iKi,i+1Ki+1,iKi+1,i+1=11j-j.Λr will contain 2 × 2 blocks on the diagonal for complex eigenvalues, and 1 × 1 blocks for the real eigenvalues. By substituting (22) into (19), we can get
(23)JrΛrJr-1x~˙+x~=W~ist,uot=B~Tx~.
The final transfer function is given by
(24)H~s=B~TJr-1Λr-1s+Iq-1JrW~.
As a result, we can obtain the final reduced-order model in the following form:
(25)H~ijs=∑m=1qK~ijms+p~m+d=bijksq+⋯+bij1s+bij0sq+ak-1sq-1+⋯a1s+a0,
where H~ijs=H~jis, 1≤i, j≤N, p~ and K~ are the dominant pole and corresponding residue, respectively, q is the number of reduce-order, and d represents any direct coupling between ports.


## 4. Time-Domain Simulation of the Reduced-Order Model

For a complete simulation of circuit response, the linear elements should be simulated along with the reduced-order macromodels. Here, we describe a method to include the ROA macromodel. The method is based on ay-parameter description of the macromodel [21]. In order to calculate the y-parameters of the reduced-order system, the eigendecomposition steps [(22) and (24)] are used. After obtaining the final form of the transfer function H~ijs, the time function h~ijt can be calculated by using inverse Laplace transform of partial fraction expansion.According to(25), we need to transform the rational fraction into proper fraction. So, we can get the following:
(26)H~ijs=d+y~ijsuijs,
where the remainder term y~ijs/uijs is a proper fraction.The denominator polynomials can be factored via partial fraction expansion, and the roots of the equationuijs=0 are solved. The roots of the equations have two kinds of situation, such as simple root and conjugate complex.(i) If the equationus=0 has n simple roots, we can set the n simple roots as p1,p2,…,pn. Therefore, the expansion form of H~s is
(27)H~s=K1s-p1+K2s-p2+⋯+Kns-pn,
where Ki=y~s/u′ss=pi, i=1,2,3,…,n.After determining the coefficients of(27), the time function is
(28)ht=∑i=1ny~su′sepit.(ii) If the equationus=0 has two conjugate complex roots (p1=α+jω, p2=α-jω), we thus have
(29)K1=s-α-jωH~ss=α+jω=y~su′ss=α+jω,K2=s-α+jωH~ss=α-jω=y~su′ss=α-jω.
We set K1=K1ejθ1 and K2=K1e-jθ1; thus
(30)ht=K1eα+jωt+K2eα-jωt=K1ejθ1eα+jωt+K1e-jθ1eα-jωt=2K1eαtcos⁡ωt+θ1.Subsequent to the calculation of the time functionhij(t), convolution is needed for time-domain analysis. Hence, the circuit response in time-domain at the ports can be given by
(31)yijt=∑j=1N∫0tujthijt-τdτ,1≤i≤n,
which requires O(T2) complexity, where T is the number of time points during simulation. For this reason, recursive convolution [22] and time-domain y-parameter macromodels [21] were developed, where the complexity is linear with the number of time points. The algorithms of these approaches are used by MATLAB in detail.


## 5. Numerical Example

A reduced-scale structure of lightning protection system is shown in Figure5. The diameter of each steel bar is 8 mm. According to the above model and algorithm, the programs are written by MATLAB. Table 1 summarizes the reduction results of matrix dimension and the computational time by using ROA.Table 1
The dimension of the matrices and the computational time for the original and reduced model.
System
Order
Computational time (s)
Original
69
5.680
ROA
10
1.779Figure 5
The structure diagram of a reduced-scale model.In order to verify the validity of the reduced-order model proposed above, the amplitude-frequency curves and phase-frequency curves of original and reduced-order model are computed, as shown in Figures6 and 7.Figure 6
The amplitude-frequency curves of the circuit.Figure 7
The phase-frequency curves of the circuit.From the two figures above, we can see that a better agreement appears between original and reduced-order model.For the sake of computing the circuit response of the reduced-scale model, the injection lightning current needs to be known. The waveform of the current sourceis(t) is shown in Figure 8. Therefore, the circuit response of reduced-order model is computed by the program we have written. The waveforms of the lightning transient responses i1(t) and uo(t) can be obtained and they are compared with that measured experimentally [23] as shown in Figures 9 and 10. It can be seen from Figures 8 and 9 that very good agreement is also found between computed and measured results.Figure 8
The waveform of injection current.Figure 9
Measured and computed current response curves for the reduced-scale structure.Figure 10
Measured and computed voltage response curves for the reduced-scale structure.


## 6. Conclusions

A novel reduced-order macromodeling algorithm has been proposed for calculating the lightning transient responses in the lightning protection systems of structures. The reduced-order model given here uses the block Arnoldi algorithm to reduce the matrix dimension. On the basis of matrix operation and discretization, the amplitude-frequency curves and phase-frequency curves can be obtained by means of the MATLAB program. These curves under including and excluding the reduced-order treatment are compared with each other and a better agreement is shown between them. The inverse Laplace transform of partial fraction expansion has also been performed for time-domain analysis. The computed lightning transient responses in a lightning protection system conform reasonably with those from experimental measurement, which confirms the validity of the proposed model.



---

*Source: 484189-2014-08-05.xml*