Industrial PLC Network Modeling and Parameter Identification Using Sensitivity Analysis and Mean Field Variational Inference

A multiple input multiple output (MIMO) power line communication (PLC) model for industrial facilities was developed that uses the physics of a bottom-up model but can be calibrated like top-down models. The PLC model considers 4-conductor cables (three-phase conductors and a ground conductor) and has several load types, including motor loads. The model is calibrated to data using mean field variational inference with a sensitivity analysis to reduce the parameter space. The results show that the inference method can accurately identify many of the model parameters, and the model is accurate even when the network is modified.


Introduction
Power line communications (PLC) systems are becoming increasingly common and are a key component of smart grids, as they allow grid usage to be monitored without adding extensive infrastructure. Currently, PLC systems are commonly used for automatic meter reading and communication between electrical substations, but there is a wide range of potential applications, including dynamic load management, appliance and lighting control, device-specific billing, and power outage detection [1][2][3][4][5][6]. There are three classifications of PLC systems based on bandwidth: ultra narrowband (UNB) systems have frequencies ranging from 0.3-3 kHz, narrowband (NB) systems have a frequency range of 3-500 kHz, and broadband (BB) systems have frequencies of 1.8-250 MHz [7]. Lower transmission frequencies have slower data transmission rates but are less affected by losses and are, therefore, more effective over longer distances. In contrast, broadband systems are limited to applications within the home or industrial facilities but have much faster data transmission rates.
Power line networks are far from ideal for signal transmission. The loads, wires, junctions, and network configuration all affect the signal transmission, resulting in systems with high variability and noise. PLC systems must, therefore, be customized for individual systems to effectively transmit signals, which requires developing models that characterize transmission pathways for individual networks.
Two types of PLC models have been developed: top-down and bottom-up models [8,9]. The top-down approach generates a parametric model based on the transmission line (TL) theory [10] but does not incorporate information about the network topology and loads. The model parameters are found by fitting the model output to the data. The topdown modeling approach is best used when the full network information is unavailable, but measurements of the network are available. Because the network topology is not incorporated into the model, the model can no longer be used if the network components or the topology are changed.

PLC Network Model
The PLC network model used in this paper is a modified version of the networks used in Refs. [12,29], which are models of home networks. There are a number of differences between this work and Ref. [29]: this work models 4-conductor cables to account for threephase wiring typical of industrial facilities instead of 2-conductor cables, the load models are different and account for motors, and the networks contain a larger number of separate circuits and loads.

Network Topology
The PLC network consists of a number of rooms, each on a separate circuit connected to the service panel, as shown in Figure 1. For this work, we use six rooms, each containing four outlets. Each circuit is connected to a junction box, from which cables leading to the outlets are connected. The outlets are organized in a bus connection, where only two cables are connected directly to the junction box. The service panel and junction boxes are assumed to be ideal connections with zero impedance. Two outlets in the facility are connected to the transmitter and receiver; the rest are connected to loads. The loads in the rooms that are not at wire terminations are connected in parallel to the wiring.

Cable Model
The cables are modeled using TL theory [10,30]. It is assumed that all four conductors in the cable have identical diameters and materials, and the conductors have a square layout, shown in Figure 2. The resistance, inductance, capacitance, and conductance are designated R, L, C, and G with values: where r w is the conductor radius, µ 0 = 4π · 10 −7 [H/m] is the vacuum permeability, σ = 5.8 · 10 7 [S/m] is the conductivity of copper, = 3.19 · 10 −11 [F/m] is dielectric constant of PVC, f is the frequency. The loss tangent tan(δ L ) is assumed to be zero, so G = 0. The distances between conductors are d 0, where d c is given by the relation: This relation is from Ref. [29] and was found through a linear fit between conductor radius r w and conductor-to-conductor distance for Thermoplastic High Heat-resistant Nylon (THHN) wires between gauge sizes 6 and 16. The skin depth δ is given by

Load Models
Each load consists of six component loads connecting each conductor in the threephase cable with each other conductor, as shown in Figure 3. The components of each load connected to the ground cable are designated Z G1 , Z G2 , and Z G3 , and the three phase-tophase components are designated Z 12 , Z 13 , and Z 23 .
There are three types of loads: constant impedance, double RLC loads, and motor loads. Constant impedance loads are parameterized by a resistance R const and leakage capacitance C c,leak . The three phase-to-phase components are equal to R const , and the three ground-to-phase components are the impedance of the leakage capacitance. The load equations for a constant impedance load are where ω = 2π f . The constant impedance load is a model for simple loads such as incandescent lights and battery chargers. A double RLC load has a series RLC circuit connected to a parallel RLC circuit, as shown in Figure 4. Unlike the constant impedance load, there is an asymmetry in the phase-to-phase components. The double RLC load has 9 parameters, which are R s , ω 0s , ζ s , R p , ω 0p , ζ p , ∆ 1 , ∆ 2 , and C d,leak . The load equations are where the parameters ∆ 1 and ∆ 2 provide asymmetry in the load. ∆ 2 Note that each RLC circuit is converted from R, L, and C to natural frequency ω 0 and damping factor ζ using the formula: The double RLC load can be considered a model for many loads with more complexity than constant impedance loads. For example, Refs. [31][32][33] took measurements of a variety of appliances and showed that many of them had resonant frequencies. Ref. [34] modeled a number of household appliances using combinations of parallel and series RLC loads, and the double RLC model is a generalization for many of the appliance models. Industrial networks often contain large numbers of motors, so in this work, a separate model is used for motors [35]. This model is shown in Figure 5 and has been used in a variety of works, including Refs. [36][37][38][39][40] for three-phase induction motors. The motor model has 5 parameters: C m , L m , R m1 , R m2 , and C m,leak . However, the model is insensitive to R m2 except at low frequencies, so its value is fixed at R m2 = 5 Ω. The motor load formula can be reduced to where Z = jωL m + R m2 .

Network Generation
The cable lengths l w and conductor radius r w are generated with the parameter ranges in Table 1. The conductor radii for wires connected to the service panel are for wires with gauges between 12 and 6, while all other wires have ranges for gauges 14 to 10. Each load in the network is generated with a random load type, where each load type has equal probability. The load parameters are generated over the ranges given in Table 1 and using either a uniform distribution or a log-uniform distribution.
When performing the inference, all of the parameters are allowed to vary over the ranges defined in Table 1 except cable length l w and conductor radius r w . The conductor radius is fixed to the original value from when the parameters were assigned, and the cable length can vary ±1 m from its original value. This was done because it was found that performing inference when all the parameters were allowed to vary over large ranges gave poor results. The wire parameters are more likely to be known in a facility than the specific load parameters, so we chose to reduce the uncertainties associated with the wire parameters while keeping the load parameter ranges large.

Transmitter and Receiver Models
The transmitter and receiver models are adapted from Ref. [20] for three-phase systems. The transmitter model is shown in Figure 6 and consists of one voltage source connected to each of the three-phase conductors. There is a single resistor between each pair of conductors as well as one after each voltage source. The receiver model is shown in Figure 7 and consists of a single resistor on each line.  The network is solved by reducing the entire network except for the transmitter and receiver to a single transmission parameter matrix, then computing the transfer function. The equations modeling the transfer function through the system are derived in the Appendix A.

Model Solver
The PLC network model is solved in Python using 4-port transmission parameter matrices. Each transmission line and load can be formulated as a 6 × 6 matrix, shown in Figure 8 and Equation (19).
Each load can be put in the form of an admittance matrix Y L using the equation: The admittance matrix of a load connected in parallel to transmission lines can be converted to a transmission parameter matrix using Equation (21), and the transmission parameter matrix for loads connected in series is given by Equation (22).

Mean Field Variational Inference
To estimate model parameters, we use a variational Bayes approach. We seek a posterior distribution that summarizes what we know about model parameters based on data and prior knowledge. In general, the posterior is calculated according to the Bayes rule.
where d ∈ R N d are the observed data, θ are the model parameters, and p( ) designates a probability. p(θ) is called the prior probability. The prior probability is based upon established knowledge in the area of a particular application. p(d|θ), the likelihood, measures agreement between model output and the simulated data for given parameter values. Hence, parameters that align with prior beliefs and result in a close match with the data have the highest posterior probability. For this work, the prior probability is modeled with a uniform distribution where the bounds are set in Table 1. We define our likelihood using a Student's t-distribution. Using a Student's t-distribution resulted in faster convergence than a Gaussian function.
where y ∈ R N d represents our model output and T (µ, σ, ν) denotes a Student's t-distribution with mean µ, variance σ, and degrees of freedom ν. For our purposes, we set σ and ν to be 0.4 and 1, respectively. To apply Bayesian inference, p(θ|d) must be approximated as it can rarely be computed directly. Markov chain Monte Carlo (MCMC) methods are generally used to draw samples from the distribution of p(θ|d). In MCMC, samples are drawn using a Markov chain which is constructed to have stationary distribution p(θ|d). Unfortunately, when N θ is large, MCMC methods often take a prohibitive amount of time, and convergence can be difficult to determine.
Variational inference methods provide an alternate approach to approximate otherwise intractable posterior distributions and are faster in most cases [41,42]. Assuming our posterior distribution comes from a specific family of distributions, we can cast the inference problem as an optimization problem. Variational inference methods most commonly use Kullback-Leibler divergence to measure the difference between distributions. We call a realization of our family of distributions q φ (θ).
Prior to performing inference, all network parameters are converted from the ranges shown in Table 1 to the range [0, 1] to ensure that they are all assigned equal weights. Accordingly, the variational family is chosen to be a Gaussian distribution truncated to be between 0 and 1. Here, we start with a standard Gaussian distribution and use rejection sampling to truncate the distribution according to our bounds. Parameter independence is the mean-field assumption, equivalent to setting the covariance matrix to the identity matrix. In general, the goal is to choose a variational family that is simple enough for efficient optimization but can capture a density close to p(θ, d).
Minimizing the KL divergence between our variational family and the posterior is equivalent to maximizing the Evidence Lower Bound (ELBO) with respect to q φ (θ).
For details on the derivation of the ELBO, see Ref. [41].
To maximize our objective function, it is desirable to have gradients of Equation (25). Finding an exact gradient is typically difficult as this would require gradients of both the model and variational distribution. We use a stochastic gradient estimator as we cannot compute gradients of arbitrarily complicated PLC models and variational distributions. To optimize the ELBO with stochastic optimization, we need an estimator of the gradient, which can be computed using samples from the variational distribution.
This is called the score function gradient estimator, and its derivation can be found in Ref. [43]. As the derivative is taken with respect to the parameters of q, this approach does not require the calculation of the derivative of p(d, θ). Due to the high variance of the score function estimator, using this alone often leads to poor results. It is common to use this estimator combined with various variance reduction techniques. Using a higher number of samples will lead to improved estimates of the gradient and faster convergence, although it also increases computational costs. In our work, we use a control variate strategy, where terms are added to the estimator that has zero expectation but is specifically chosen to reduce variance. Essentially, terms of the form log q φ (θ)g(θ) are replaced with log q φ (θ)(g(θ) − b) where b is a constant and g is not to be differentiated with respect to φ. The identity E q φ (θ) [∇ φ (log q φ (θ) × b)] ensures this does not change the mean of our gradient estimator. A carefully chosen b can reduce the variance of our estimator. The value b is referred to as a baseline. In this work, our baseline is constructed using a decaying running average of samples from the objective function. For further discussion of baselines, see Ref. [44].
To perform inference on our model, we use the probabilistic software package called pyro [45,46]. The score function gradient estimator is computed using 14 samples, which are computed in parallel. An AdaGrad optimizer with an initial learning rate of 0.6 converges samples to the posterior distribution [47].

Results
Model calibration relies upon data simulated using the PLC model with randomly generated parameters to compute a transfer function for subsequent inference of network parameters. This also allows a direct comparison between the randomly generated parameters and the inferred parameters. 1.0 Decibel of noise is added to the transfer function before applying mean field variational inference.
Each network model can have over 100 parameters, and a parameter space this highdimensional can make inferring parameters challenging. Furthermore, many of the network parameters are from cables or loads far from the main path between transmitter and receiver and have very little effect on signals between them. Therefore, this study uses sensitivity analysis to determine which parameters contribute least to the variance of the transfer function, similar to Ref. [29]. Computing a full sensitivity analysis as was done in Ref. [29] is computationally expensive, so this study instead uses an approximate sensitivity analysis.
To dramatically reduce the number of model evaluations that would otherwise be required by a global sensitivity analysis, a one-at-a-time (OAT) approach is employed. Varying each parameter ten times at values distributed throughout their domain while fixing all other parameters at their nominal values, the variation, defined as the twonorm of the difference with the nominal transfer function, is calculated with respect to each parameter. The total variation is defined as the sum of all of the variations due to each parameter. An estimate of the relative sensitivities of the parameters is obtained by normalizing each parameter's variation by the total variation. Then, based on an adaptively selected threshold, parameters are kept or discarded based on their contribution. It was found through a series of tests that a one-half percent threshold often reduced the parameter space by more than half and, correspondingly, the runtime by a significant margin while achieving virtually equivalent results to that of the full parameter model.

Demonstration on Single Network
As an initial demonstration, a network is generated with random load types, but all the network parameters are assigned a "true" value of 0.25, so the reader can easily identify mispredicted parameters. This demonstration is to quantify the accuracy of the inference methodology so the loads are fixed to their true types.
The mean and variance of each sensitive parameter are shown in each iteration step in Figure 9. The network has a total of 133 parameters, and 75 of them are sensitive parameters. All but 6 of them converged to mean values between 0.08 and 0.4, showing that mean field variational inference can accurately infer most of the sensitive parameters. Table 4 shows the mean parameter values at the end of the run ordered by decreasing sensitivity. The most sensitive parameters are inferred very accurately, while the accuracy noticeably decreases for less sensitive parameters.   Figure 10 shows the resulting transfer function component H 1,1 ( f ) along with the 95% confidence interval computed from the push-forward posterior. The push-forward posterior is the distribution of the transfer function computed from samples drawn from the posterior. The figure shows that the inference methodology accurately fits the model to the true transfer function, and the uncertainties are highest at frequencies where the transfer function has a peak. The focus of this work is on the calibration methodology, so for brevity, we refrain from a detailed study of the full transfer function matrix. However, the transfer function was found to be full rank across the entire frequency range, justifying its usage for MIMO applications.
One of the main advantages of a bottom-up model is that the model can still be used if the network changes. Here, a new network is generated with random parameter values over the range [0, 1] and then calibrated with mean field variational inference. The network is then modified by powering off all of the the loads in the two central rooms of the network. This is done by assigning them as constant loads with the resistance of R const = 10 6 Ω and parasitic capacitance C c,leak = 0.1 nF. Figure 11 shows the transfer function of the model with the loads turned off compared with the true model. While there are slight discrepancies, the calibrated model is very accurate despite not having been trained using data from the modified network.

Demonstration on Multiple Networks
The results in Section 4.1 showed that the inference method could calibrate the model for a single realization of the network. Next, the method is repeated 10 times for different network realizations to demonstrate that the method is consistently able to calibrate the model. However, instead of setting the true values of the parameters to 0.25, all network parameters have randomly generated parameters over the range [0, 1]. Figure 12 shows the mean transfer function and true transfer function for each realization. Slight inaccuracies can be observed, particularly at frequency peaks, but generally, there is good agreement between the data and calibrated model.  Table 5 shows the magnitude of the average error for parameters in order of decreasing sensitivity. Each error in the table is averaged across the indicated parameters and the 10 realizations. The parameters with the highest sensitivity to the 10th highest sensitivity have an average error magnitude of 0.029, but the errors increase as the sensitivity decreases. Note that for parameters chosen randomly in [0, 1], randomly guessing the parameters would give an average error magnitude of 0.333. Since all the error magnitudes in Table 5 are significantly lower than 0.333, the inference method can consistently identify parameters, particularly for parameters with high sensitivity. Table 5. Error magnitude of inferred parameters averaged across 10 realizations and the parameters in the given range. Parameters are ordered with decreasing sensitivity.

Parameters
Average Error Magnitude

Inferring Load Types
A natural extension of this analysis is to consider the case where there are unknown load types that may be desirable to infer. In a home kitchen, for example, knowledge of a fridge and microwave may be known, while additional appliances may not be as clear, although a set of likely appliances may be assumed, given the setting. It follows that one may be interested in resolving the identity of loads; however, it must first be assessed whether this is a feasible task provided the selection of potential load types, the model to which they are applied, and the method of inferencing.
In the previous sections, the load types (constant, double RLC, or motors) during inference were identical to the true load types. In reality, it may not always be possible to know all of the loads in the networks, especially as loads may change during operation. It may, therefore, be important to infer load types in the network from data.
We performed an analysis to identify how much the transfer function changed based on load types. For brevity, we summarize the method and findings. Individual loads in a network were changed to a different load type, then inference was performed to see if a set of parameters could be found that resulted in an accurate transfer function even though one of the load types was wrong. It was found that if a motor load was replaced with a different load type or vice versa, it was generally impossible for the calibrated model to have an accurate transfer function. In contrast, constant loads and double RLC loads could generally be interchanged, and the calibrated model output would still be accurate. This suggests that constant and double RLC loads are too similar to identify from each other but that motor loads can be distinguished from the other load types.
We treat the constant and double RLC load types as a single group represented by the constant impedance load, as it has fewer parameters, and a slightly modified version of the algorithm offered in [29] is employed to predict load types. The algorithm initializes by assigning random load types to each load. This assignment is evaluated by performing a full parameter inference and recording the likelihood function. Then, three loads are randomly changed in type, and the change is accepted or rejected based on whether or not the subsequent full parameter inference yields better agreement than that of the previous assignment. This process continues for a predetermined number of steps, although it may be readily apparent to the user that changes are no longer being accepted, indicating potential convergence. The results are illustrated in Table 6.
The inference successfully identified all twenty-two of the unknown load groupings. While succeeding in this case, it should be noted that it may require a large number of iterations for the random triplets to happen upon the exact load identities. For this problem, with a limited number of iterations, it is more often the case that all but a random few will be inferred correctly. It should also be mentioned that positionally symmetric loads may interfere with the accuracy of load type inference. Positionally symmetric loads, in this context, are loads that are connected independently and identically to a network. As an example, two loads branching off of the same node by wires of equivalent length and matching properties are positionally symmetric loads. In particular, if the loads were to be swapped, the transfer function would be unaffected, as the network has effectively not changed. Consequently, in the presence of positionally symmetric loads, this method of load type inference would achieve equivalent objective function values for predicting the correct load types or the swapped, correct load types, which may be incorrect. Table 6. True and predicted load types. The load names are room number and outlet number. The * indicates that the constant impedance load was used as the representative member during the computations.

Load
True Predicted The success in our case is also attributable to the load types under consideration. When considering motor versus constant impedance loads, the impact of the two load types is substantially different, allowing for accurate resolution via tangible differences in the transfer function, even after inferring optimal load parameters. Additionally, accounting for only two load types simplifies the problem-more robust algorithms would likely be required if operating on a larger set of load types.

Conclusions
This work developed a PLC model for three-phase industrial networks and used mean field variational inference with sensitivity analysis to calibrate the model. The calibrated model has the advantages of both bottom-up and top-down models as it incorporates full knowledge of the network topology but can also be tuned to match data. The results showed that many of the parameters could be inferred accurately. When averaged across different network realizations, the 10 most sensitive parameters had an average error of 2.9% after calibration, and the 40 most sensitive parameters had an average error of 7.1%. This work also found good agreement between the transfer function produced by the calibrated model and a network modified by powering offloads, suggesting that the model maintains accuracy when the network undergoes changes such as loads being powered off.
This study used mean field variational inference for calibration, whereas our previous work (Ref. [29]) used transitional Markov chain Monte Carlo (TMCMC). The key difference between the methods is that mean field variational inference assumes the distribution of each parameter is independent of all other parameters and follows an assumed distribution profile, which in this case is Gaussian. While these assumptions allow mean field variational inference to converge to a solution faster than TMCMC and work better in high dimensional parameter spaces, care must be taken that it is only applied to problems where the assumptions are valid. Our previous work found that very few network parameters were strongly correlated, and generally, the distributions had a single mode. That, together with the accurate results obtained here, justifies using mean field variational inference.
Our analysis of inferring load types found that motor loads could be distinguished from constant or double RLC loads. In a 6-room network, all 22 motor loads were inferred correctly, but constant loads could not be distinguished from double RLC loads. This raises the question of what kinds of loads are needed for a PLC model. It is possible that a reasonable PLC model could contain one of those two types of loads rather than both.
The accuracy of the network structure and load models should be considered when analyzing PLC models. In this and many other previous works, both the network structure and load models were generated using intuition and generalizations rather than from statistics in actual applications. While this is done due to the lack of available statistics and does not discredit such works, it must be taken into account that the models may not reflect PLC networks in actual applications.
There are several directions future works could explore. One area where more research is needed is accurate load modeling. Most load models in PLC network models are unrealistic due to the lack of data from loads common in households and industrial networks. More work is needed to infer load types or cable parameters from data. Here and in Ref. [29], a greedy search algorithm is used, but that method is expensive and can be inaccurate. Future work could further explore the potential for MIMO applications in industrial facilities. Finally, research is needed to measure facilities and develop calibrated models custom to each facility. Data Availability Statement: Data will be made available by the corresponding author upon reasonable request.

Acknowledgments:
We are grateful for advice and suggestions by Monica Jaramillo and Ricky Tang. Affiliation: Sandia National Laboratories, New Mexico, Albuquerque, NM 87185.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript:

Appendix A. Transfer Function Computation
A schematic of the network with the transmitter and receiver is shown in Figure A1, where the entire network except for the transmitter and receiver is reduced to a single transmission parameter matrix Φ nw . Three locations a, b, and c are designated by the blue dashed lines. The transfer function is defined as We define transfer functions through a part of the network as Using the chain rule, H = H nw H trans (A4) Figure A1. Schematic of the network reduced to a transmission parameter matrix.
We first compute H nw by expanding the network parameter matrix Φ nw into submatrices and then substituting in the receiver load: V c I c = φ nw,11 φ nw,12 φ nw,21 φ nw,22 We can also define submatrices of the inverse matrix The receiver load matrix is defined as Z rec where, We assume that the power grid has a fixed voltage and, therefore, does not carry signals from the transmitter. For the particular receiver in Figure A1, the transmitter load matrix is computed as follows: V c,1 = Z R1 I c,1 + (I c,1 + I c,2 + I c,3 )Z RG V c,2 = Z R2 I c,1 + (I c,1 + I c,2 + I c,3 )Z RG V c,3 = Z R3 I c,1 + (I c,1 + I c,2 + I c,3 )Z RG (A8) Substituting Equation (A7) into the first row of Equation (A6) yields: Reorganizing gives: H nw = (φ −1 nw,11 + φ −1 nw,12 Z −1 rec ) −1 (A11) It is also useful to compute an effective admittance Y nw of the network and receiver. Equation (A7) can be substituted into Equation (A5) to obtain: where: Y nw = −(φ nw,12 − Z rec φ nw,22 ) −1 (φ nw,11 − Z rec φ nw21 ) (A13)

Appendix A.2. Transmitter Computation
Three equations can be formed by applying Krichoff's current law to the three conductors in the transmitter: The impedance law across resistors Z T0 is Equations (A15) and (A12) are substituted into (A14) to eliminate the variables I a and I b .