Efﬁcient Statistical Extraction of the Per-Unit-Length Capacitance and Inductance Matrices of Cables with Random Parameters

Cable bundles often exhibit random parameter variations due to uncertain or uncontrollable physical properties and wire positioning. Efﬁcient tools, based on the so-called polynomial chaos, exist to rapidly assess the impact of such variations on the per-unit-length capacitance and inductance matrices, and on the pertinent cable response. Never-theless, the state-of-the-art method for the statistical extraction of the per-unit-length capacitance and inductance matrices of cables suffers from several inefﬁciencies that hinder its applicability to large problems, in terms of number of random parameters and/or conductors. This paper presents an improved methodology that overcomes the aforementioned


Introduction
Cable bundles are commonly employed in a variety of critical applications, including medical, industrial and aerospace equipment.Their precise position, as well as their physical properties, are hardly controllable and a probabilistic assessment of these interconnections is therefore crucial.Indeed, electromagnetic compatibility (EMC) studies comprise the statistical assessment of the coupling (crosstalk) occurring in wire bundles with random parameters [1]- [10].Most of the above contributions aimed at deriving closed-form expressions for the statistics of the random crosstalk.In fact, the brute-force Monte Carlo (MC) method, which is commonly employed for statistical analyses, is known to require a large number of samples to converge, thus becoming inefficient when the calculation involves many conductors and requires numerical methods.
Among the many recent approaches to account for uncertainty in electromagnetic problems (see, e.g., [11,12]), an alternative strategy was proposed to characterize cable responses from a statistical standpoint [13].The approach is based on the so-called polynomial chaos (PC) frame-work [14] and was also applied to printed circuit board (PCB) lines [15].According to PC, the stochastic electrical quantities, i.e., the line voltages, currents and per-unitlength (p.u.l.) parameters, are expanded in series of orthogonal polynomials.Once the series coefficients are determined, pertinent statistical information is readily obtained from these expansions.
The calculation of the PC-expansion (PCE) coefficients of the p.u.l.capacitance and inductance matrices for random cables was specifically addressed in [16].A new numerical scheme was put forward based on the twofold expansion of the charge distributions along the wire peripheries in terms of Fourier and PC series.A deterministic system of equations for the PCE coefficients is constructed by applying a stochastic Galerkin method (SGM) [17] to the classical numerical schemes for circular conductors [18]- [20], and then solved by standard linear system inversion.The technique was successfully used for the statistical assessment of crosstalk in cable bundles [21].However, the approach in [16] rapidly becomes intractable when a large number of random parameters and/or conductors is considered.
In this paper, the aforementioned limitations are pinpointed and overcome by adopting a new strategy to construct the deterministic system of equations for the extraction of the PCE coefficients of the p.u.l.capacitance and inductance.Namely, the stochastic testing (ST), proposed in [22] for statistical circuit simulation, is used in place of the SGM.This consists of a point matching or collocation [23], rather than Galerkin, approach.ST yields a system of equations that is inherently sparse and decoupled, thus allowing a much faster and less memory-demanding computation.
The rest of the paper is organized as follows: Section 2 briefly recalls the deterministic numerical scheme for the calculation of the p.u.l.capacitance and inductance matrices for circular wires.The principles of the PC-based approach and the state-of-the-art stochastic scheme to calculate the PCE coefficients of the p.u.l.parameters are summarized in Section 3. Section 4 outlines the advocated ST-based technique.Numerical validations are provided in Section 5. Finally, conclusions are drawn in Section 6.
Figure 1: System of dielectric-coated wires.The charge distributions ρ k and ρ ′ k on the kth conductor and dielectric peripheries, respectively, are indicated, together with the potential ϕ i and the normal components of the displacement field ⃗ D they produce on conductor i (k, i = 0, . . ., n − 1).

Deterministic Scheme
The classical deterministic numerical scheme for the calculation of the p.u.l.capacitance and inductance matrices for circular conductors [18]- [20] is presented first.The discussion starts from a system of n wires like in Fig. 1, with the zeroth wire being the reference for line voltages and the return path for line currents.Owing to the circular symmetry, the charge distributions ρ and ρ ′ on the conductor and dielectric peripheries, respectively, are expanded as Fourier series: with k = 0, . . ., n − 1.The number of expansion terms N = 1 + 2A defines the accuracy of the calculation.A larger N is required when the charge distributions are strongly nonuniform, which is caused by wire proximity.Hence, the worst-case scenario is the one of nearly touching wires, for which the scheme is found to converge for A ≥ 10 (i.e., N ≥ 21) [18].Specifically, conductor proximity turns out to be more critical than dielectric proximity, which makes the convergence even faster for coated wires.Each of the charge distributions produces an electric potential ϕ as well as an electric field ⃗ E = −∇ϕ.There exist a closed-form expression for the total potential and electric field in a given point in the space, as a result of the superposition of all the charge distributions [24].To solve for the all 2nN Fourier coefficients, two types of boundary conditions are enforced: i) the potential on each conductor interface must be independent of the angular position θ; ii) the difference between the normal components of the electric displacement field ⃗ D = ε ⃗ E across each dielectric interface must be null.Both boundary conditions are enforced at N discrete points on each metallic and dielectric interface of the n wires, yielding a system of 2nN equations in the form where Φ = [Φ 0 , . . ., Φ n−1 ] T , with Φ i = [ϕ i , . . ., ϕ i ] of length N and ϕ i being the potential on the periphery of the ith conductor.Vectors A and A ′ collect the Fourier coefficients of the charge distributions, whilst the contributions of the latter fill matrices D uv (u, v = 1, 2) [24].The entries of the generalized capacitance matrix C are defined as [18] where r wk and r dk are the conductor and dielectric radii of the kth wire, respectively (see Fig. 1).The coefficients a k0 and a ′ k0 are the the leading terms in the Fourier expansions (1), and are computed via the inversion of (2).
The entries of the standard (Maxwellian) p.u.l.capacitance matrix C, as well as those of the p.u.l.inductance matrix L, are readily derived [19].When the wires are lying above an infinite ground plane, a similar system of equations can be obtained by considering the additional contributions from the wire images [19].Finally, circular shields can be accounted for as additional conductors by expanding the charge distribution on their inner peripheries as in (1a).

Polynomial Chaos-Based Analysis
In this section, the state-of-the-art numerical scheme for the extraction of the capacitance and inductance PCE coefficients for circular conductors is reviewed and its limitations are highlighted.Only the equations that are relevant for the present paper are reported.For detailed information, the reader should refer to [16].

The Polynomial Chaos Expansion
Let the problem become stochastic and dependent on d RVs denoted as ξ = [ξ 1 , . . ., ξ d ].For the sake of simplicity, and in the general absence of information on the correlation among wire parameters, only the case of independent RVs is considered in this paper.It is worth mentioning that dependent parameters can be included in a relatively easy fashion, based on Karhunen-Loève expansion [25], only when they are Gaussian distributed (see also [26]).The modeling of dependent and non-Gaussian parameters still represents an active field of research instead [17].Moreover, it is further assumed that, within the given variability, the overlap between wires cannot occur or is extremely unlikely.
Since the p.u.l.parameters depend on the RVs ξ, they become stochastic themselves.Any stochastic quantity is expressed as a PCE [17].For instance, the p.u.l.capacitance is written as where C γ are the pertinent PCE (matrix) coefficients, whilst {φ γ } P γ=0 forms a complete basis of multivariate polynomial functions that are orthonormal with respect to the joint probability density function (PDF) w(ξ) of ξ.
The multivariate polynomials φ γ (ξ) satisfy the orthonormality condition with δ γβ denoting the Kronecker's delta.They are constructed as the product combination of the univariate polynomials that are orthonormal with respect to the univariate PDF of ξ i [17].For Gaussian and uniform RVs, these univariate basis functions are the Hermite and Legendre polynomials, respectively [14].By selecting polynomials of total degree up to p, the number of PCE terms is and choosing p = 2 can already provide satisfactory accuracy [16,21].
Once the coefficients of the PCE are available, statistical information is readily extracted, either analytically or numerically, from (4).For example, the average p.u.l.capacitance matrix µ C corresponds to while its standard deviation σ C is where the square and square root in (8) are not intended as matrix operations, but rather as operations on each individual entry.The derivation of closed-form expressions for higher-order moments becomes increasingly more cumbersome [26].However, this information is easily retrieved numerically by randomly sampling the PCE (4).The rationale of PC-based approaches is that the determination of the coefficients C γ is usually much faster than collecting statistical information using e.g., MC or similar methods.
It is worth mentioning that the knowledge of the PCE coefficients of the capacitance and inductance matrices not only provides pertinent statistics.It also allows to obtain stochastic information concerning the cable response, like for example the crosstalk between the wires.This is achieved through the simulation of an equivalent, deterministic cable with augmented p.u.l.matrices, which are constructed from the PCE coefficients of the capacitance and inductance matrices [13].

Calculation of the Expansion Coefficients
Given the PCE (4) and the definition of the inner product (5), different strategies are available to solve for the unknown PCE coefficients [17,23,26].Sampling-based strategies, including pseudo-spectral or collocation methods, as well as regression-based techniques [27]- [29], sample the responses at a suitable subset of points in the random space.The PC coefficients are then retrieved via quadrature, interpolation or linear regression of these samples, respectively.In particular, the pseudo-spectral method is a projection-based approach that calculates the PCE coefficients according to the classical projection theorem: i.e., by integration of the p.u.l.matrix in the random space.
The integration is carried out for every matrix entry via Gaussian quadrature rules [30].
The sampling-based approaches are intrinsically nonintrusive, but the required number of samples rapidly grows with the number of RVs, even when sparse grids are adopted.It is important to remark that, except for the case of bare and widely-spaced wires, for which close-form result are available [24], the evaluation of the p.u.l.matrices requires the numerical scheme in Section 2. This hinders the efficient calculation of the many samples required by these techniques.
The SGM is an alternative, intrusive approach that requires a modification of the governing equations, but potentially allows for a more accurate simulation [17].Unfortunately, the SGM requires to solve a larger set of coupled equations, which dramatically reduces the efficiency when the problem dimension is large.In [16], this stochastic numerical scheme was integrated in the classical deterministic procedure for the determination of the p.u.l.matrices of cable bundles described in Section 2, but the effectiveness was limited to structures consisting of a limited number of wires and RVs, as explained in the next section.

State-of-the-Art Galerkin-Based Scheme
The stochastic scheme [16] is based on the additional expansion of the (random) charge distributions (1) into PC series.Specifically, each Fourier coefficient becomes ξ-dependent and is represented with a PCE like (4), e.g., and similarly for b km , a ′ km and b ′ km .The potentials ϕ i and the entries of matrices D uv in (2) are also expanded.To obtain the PCE coefficients for matrices D uv , a pseudospectral integration (refer to Section 3.2) is applied to their analytical entries.
Then, a SGM, consisting in weighing the resulting equations using the basis functions φ γ and integrating them with the inner product (5), is applied, leading to the following augmented and deterministic system In (11), Φ is a vector collecting the PCE coefficients ϕ i,γ of the wire potentials, whilst the coefficients of the charge distributions are collected into Â and Â′ .
The PCE coefficients of the generalized p.u.l.capacitance matrix are obtained, after the inversion of (11), similarly to (3), using the PCE coefficients a k0,γ and a ′ k0,γ of the first Fourier terms.The PCE coefficients of the Maxwellian p.u.l.capacitance matrix and those of the p.u.l.inductance matrix are derived via simple matrix manipulations.
It is important to point out that: 1.The numerical integration of the entries of matrices D uv requires (p + 1) d evaluations of such matrices.
2. The PCE coefficients are K for each variable, thus yielding a system (11) of dimension 2nN K.
Since the system is full and does not have any sparsity pattern, its dimension rapidly becomes intractable when the number of wires n and/or RVs d increases.Indeed, the creation of ( 11) is very memory intensive and its inversion is computationally demanding.These are the major issues in the numerical scheme proposed in [16].

Proposed Improved Implementation
This section provides the proposed improved implementation and compares its performance against the previous methodology described in Section 3.

Stochastic Testing-Based Scheme
The departure point is the system (2), which is here rewritten explicitly by emphasizing the dependence on ξ: The ξ-dependent vectors Φ, A and A ′ are again written as PCEs, leading to ) where Φ γ , A γ and A ′ γ collect the γth PCE coefficients of the wire potentials and charge distributions.
Suppose now that a suitable set of K test points {ξ ν } K ν=1 in the random space R d is available.Point matching (13), i.e., enforcing the system to be satisfied at the points ξ ν , yields deterministic equations: ) with ν = 1, . . ., K, and having defined the deterministic scalar and matrix coefficients w νγ = φ γ (ξ ν ) and D uv | ν = D uv (ξ ν ), respectively, i.e., the PC basis functions and the system matrices evaluated in the test points.A proper set of test points ξ ν is generated with the algorithm proposed in [22].
Collecting the equations ( 14) ∀ν = 1, . . ., K leads to the following system where the definitions of Φ, Â and Â′ are as in (11).The matrix quantities are with u, v = 1, 2, and where ⊗ denotes the Kronecker product, W is a K ×K matrix with entries w νγ , and I nN denotes the identity matrix of size nN × nN .
The system (15) has the same size as the SGM-based system (11), i.e., 2nN K, but with an additional and remarkable property: matrices Γ and Λ uv are all inherently sparse, the former stemming from the Kronecker product of an identity matrix, and the latter being block diagonal.
The inversion of the system ( 15) is where The inverse ] −1 (20) with and being the Schur complements [32] of matrices Λ 11 and Λ 22 , respectively.Moreover, since the matrices Λ uv are block diagonal, the matrices Ψ uv are also block diagonal, and each block is computed by applying operations (21) to the individual blocks, i.e., where, consistently with ( 16), the notation | ν indicates that the corresponding block is evaluated for the test point ξ ν .This means that the operations ( 23) can be carried out iteratively for each test point ξ ν , by looping over ν = 1, . . ., K.
The evaluation of D uv | ν is performed according to the traditional deterministic calculation in Section 2, by considering the values of the cable parameters that correspond to the random sample ξ ν .
According to (18) and the definition ( 20), the Fourier-PC coefficients are computed as where, due to the particular structure of matrix Φ, with n columns filled-in by ones and zeros (see [16]), the calculation of the product Γ Φ amounts to the sums of the first n groups of N columns of Γ.The generalized p.u.l.capacitance matrix, the Maxwellian p.u.l.capacitance matrix and the p.u.l.inductance matrix are then calculated via the same procedure used in [16].To further improve the efficency, the remaining SGM-based operations involved in these computations are possibly replaced by analogous ST procedures.

Comparison with the Galerkin-Based Scheme
The ST-based scheme outlined in the previous section has some remarkable computational advantages with respect to the state-of-the-art Galerkin-based approach: 1.The matrices Λ uv are block diagonal, thus allowing to construct the inverse of Λ, to be used in (24), iteratively using relations (23).
2. The matrices Λ uv , Γ and Γ −1 are sparse, thus requiring much less memory consumption and allowing to speed-up matrix operations, in particular the multiplication in (24).Moreover, the matrix Γ −1 is efficiently computed through the inversion of a smaller matrix, i.e., W in (19).
3. The construction of system (15) does not require the numerical integration of matrices D uv .
In summary, the inversion of ( 15) requires to evaluate K matrices of size 2nN × 2nN , to apply operations (23) to them, and to perform the multiplication (24) of three sparse matrices of size 2nN K ×2nN K.As shown in the next section, this turns out to be much faster and less memory intensive than directly building and inverting the full system (11) of size 2nN K × 2nN K, as required by the Galerkin-based approach.Furthermore, it is worth noting that the terms w νγ only depend on the type of basis functions φ γ and on the set of testing points ξ ν .These are in turn only dependent on the distribution type and on the number of RVs d.Therefore, the matrices W for a wide class of problems can be pre-computed and stored into look-up tables.

Validation and Numerical Results
The advocated ST-based numerical scheme has been implemented in MATLAB .In this section, its performance against the previous Galerkin-based technique and the reference MC method is thoroughly assessed.Two different cable structures are considered: a ribbon cable and a shielded cable.
For the MC analysis, an iterative simulation is implemented starting from a set of 100 samples.The number of samples is then doubled at every iteration, until the standard deviations of the p.u.l.parameters converge within a given threshold.As a stopping criterion, a maximum relative difference below 1% is sought between the new and the old estimation of the standard deviations.It is relevant to point out the exponential growth of the sample size is necessary to ensure that the new samples bring enough additional information to the already available data.The first application considers the 5-wire ribbon cable in Fig. 2, whose geometry and parameters are taken from [24].The wire-to-wire separation is 50 mil and the conductor radius is 7.5 mil.The radius of the dielectric coating is 17.5 mil, and its relative permittivity is 3.5.The wires are numbered sequentially from left to right, the leftmost conductor being the reference.Two different scenarios are analyzed: in the first one, d = 4 RVs are considered, i.e., each wire-to-wire separation, assumed to have a Gaussian distribution with a standard deviation of 2 mil.In the second case, an additional set of five RVs, i.e., the radii of the dielectric coatings, with a standard deviation of 1.5 mil, is included, leading to an overall number of d = 9 RVs.It should be noted that an overlap of the wires might in principle occur, but this is extremely unlikely because of the relatively large difference between the gap and the standard deviations of the separation and coat thickness.For the Fourier expansion of the charge distributions, A = 10 is considered.

Ribbon Cable
Tab. 1 collects the standard deviations of the entries of the p.u.l.inductance and capacitance matrices for the two scenarios.Owing to the matrix symmetry, only the upper triangular parts are indicated.The results obtained by means of the proposed PC technique with p = 3 are compared against the estimation from a MC analysis.It is worth noting that no difference is found for the standard deviation of the inductance matrix in the second scenario, because the additional variations considered concerns the dielectrics and as such they affect the capacitance matrix only.3 shows the convergence of the MC simulation, by indicating the progress of the maximum relative difference on the standard deviations of the inductance and capacitance matrices as the sample size is increased.The target accuracy of 1% is achieved after 51200 runs for both the d = 4 and d = 9 scenarios, leading to a simulation time of 56 min and 16 s.As far as the computational efficiency is concerned, Tab. 2 collects the execution times for MC and the advocated ST-based method with first-, secondand third-order PCEs.The simulation time of the previous SGM-based implementation, when not run out of memory, is also provided within parenthesis.The table additionally reports the accuracy of ST, defined as the maximum difference with respect to MC estimations.In the first scenario, using p = 2 already suffices to reach an accuracy below 1%.Moreover, the new numerical scheme provides impressive speed-ups of over three orders of magnitude against MC.For the second scenario, which considers a relatively large number of RVs, the ST-based technique requires a thirdorder expansion to achieve a similar accuracy, nevertheless still providing a remarkable speed-up of 130×.

Shielded Cable
In the second application example, the shielded cable of Fig. 4 is analyzed.The structure is taken from [20].The conductors, the dielectric coatings and the shield have nor- Figure 4: Cross-section of the shielded cable for the second application test case.The radii of the wires, dielectric insulations and shield are of 1, 2 and 10 units, respectively.Wires #2 to #7 are situated radially at 5 units from wire #1.
malized radii of 1, 2 and 10 units of length, respectively.Wire #1 is located at the center of the shield, whilst the remaining six wires lie at 5 units of length from the center.The relative permittivity of the dielectric coatings is 4. Two different scenarios are again considered.In the first case, the radial positions of the wires #2 ... #7 are considered as independent uniform RVs in the range [4.5, 5.5] (hence, d = 6).In the second case, both the horizontal and vertical coordinates of these wires are considered as uniform RVs, with variations of ±10% around their nominal positions.For this structure, wire overlap never occurs within the given uncertainty bounds.For the charge distributions, A = 10 Fourier coefficients are again considered.Tab. 3 shows the standard deviations of the p.u.l.inductance and capacitance matrix entries for the two scenarios, estimated by means of both MC and PC with p = 3. Owing to the symmetry of the structure, the deviations of some entries coincide and are therefore omitted from the table.The MC analysis, whose convergence is shown in Fig. 5, achieves again the target accuracy of 1% after 51200 runs, which in this case correspond to a simulation time of 2 h and 23 min.Tab. 4 collects the main figures concerning the efficiency and the accuracy of the proposed method.A significant speed-up of 37× is achieved even when up to 12 variables are included.For both scenarios, third-order PCEs are needed to approach the MC accuracy.For this example, the SGM implementation only executes for d = 6 and a first-order PCE, for which it takes 14 s.

Conclusions
This paper presents a considerable improvement of the state-of-the-art technique for the extraction of the PCE coefficients of the p.u.l.capacitance and inductance matrices for cables with random parameters.The novelty of the approach lies in the fact that a point matching technique, i.e., ST, is used in place of the traditional SGM.As such, the new equations are decoupled and sparse, thus allowing for a much faster solution and the handling of larger problems, in terms of number of conductors and/or RVs.The technique is validated via the statistical assessment of the p.u.l.capacitance and inductance matrices of a ribbon cable and a shielded cable with random geometrical parameters.

Figure 2 :
Figure 2: Cross-section of the ribbon cable for the first application test case.The deviations indicated are inteded as standard deviations of a Gaussian distribution.

Figure 3 :
Figure 3: Convergence of the MC simulation for the standard deviation of the p.u.l.parameters of the ribbon cable.

Fig.
Fig.3shows the convergence of the MC simulation, by indicating the progress of the maximum relative difference on the standard deviations of the inductance and capacitance matrices as the sample size is increased.The target accuracy of 1% is achieved after 51200 runs for both the d = 4 and d = 9 scenarios, leading to a simulation time of 56 min and 16 s.As far as the computational efficiency is concerned, Tab. 2 collects the execution times for MC and the advocated ST-based method with first-, secondand third-order PCEs.The simulation time of the previous SGM-based implementation, when not run out of memory, is also provided within parenthesis.The table additionally reports the accuracy of ST, defined as the maximum difference with respect to MC estimations.In the first scenario, using p = 2 already suffices to reach an accuracy below 1%.Moreover, the new numerical scheme provides impressive speed-ups of over three orders of magnitude against MC.For the second scenario, which considers a relatively large number of RVs, the ST-based technique requires a thirdorder expansion to achieve a similar accuracy, nevertheless still providing a remarkable speed-up of 130×.

Figure 5 :
Figure 5: Convergence of the MC simulation for the standard deviation of the p.u.l.parameters of the shielded cable.

Table 1 :
Standard deviations of the p.u.l.inductance and capacitance matrix entries for the ribbon cable of Fig.2.

Table 2 :
Simulation times and accuracy for the ribbon cable analysis.

Table 3 :
Standard deviations of the p.u.l.inductance and capacitance matrix entries for the shielded cable of Fig.4.

Table 4 :
Simulation times and accuracy for the shielded cable analysis.