Open AccessArticle

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

Raelynn Wonnacott

David S. Ching

John Chilleri

Cosmin Safta

Lee Rashkin

and

Thomas A. Reichardt

Sandia National Laboratories, 7011 East Avenue, Livermore, CA 94550, USA

Author to whom correspondence should be addressed.

Sensors 2023, 23(5), 2416; https://doi.org/10.3390/s23052416

Submission received: 12 January 2023 / Revised: 8 February 2023 / Accepted: 9 February 2023 / Published: 22 February 2023

(This article belongs to the Special Issue Power Line Communication Technologies for Smart Grids)

Download

Browse Figures

Figure 1
Network topology. Blue circles represent loads, and black squares are junction boxes. "> Figure 2
Schematic of the cable cross-section. "> Figure 3
Schematic of load connections to wiring. "> Figure 4
Double RLC phase-to-phase load circuit. "> Figure 5
Motor model. (a) Phase-to-phase circuit. (b) Phase-to-ground circuit. "> Figure 6
Transmitter model. "> Figure 7
Receiver model. "> Figure 8
Transmission matrix definition. "> Figure 9
Mean value of parameters vs iteration (left) and variance vs iteration (right). Each line corresponds to a different parameter. "> Figure 10
True transfer function (without added noise) and mean of transfer function computed from posterior with 95% confidence bounds. "> Figure 11
Transfer function from the original network and network modified by powering offloads in two rooms. "> Figure 12
True transfer function component <math display="inline"><semantics> <mrow> <msub> <mi>H</mi> <mrow> <mn>1</mn> <mo>,</mo> <mn>1</mn> </mrow> </msub> <mrow> <mo>(</mo> <mi>f</mi> <mo>)</mo> </mrow> </mrow> </semantics></math> and inferred transfer function for 10 different realizations of the network. Each color corresponds to a different realization. The dots are the true transfer function and lines are the inferred transfer function. "> Figure A1
Schematic of the network reduced to a transmission parameter matrix. ">

Versions Notes

Abstract

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.

Keywords:

power line communications; mean field variational inference; sensitivity analysis; MIMO

1. 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 top-down 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.

The bottom-up model uses TL theory in combination with the network topology and conductive pathways to computing signal throughput [11]. Unlike top-down models, they can be easily modified if network components are changed. There are a number of previous works constructing bottom-up models for different applications [8,9,12,13,14,15]. The disadvantages of bottom-up models are that they require full knowledge of the network topology and pathways and that they cannot easily be modified to match measurements.

Most PLC models have been developed for home networks and incorporate either 2-conductor cables or 3-conductor cables [8,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26], although several models have been extended to 4-conductor cables [27,28]. In this work, the target networks are industrial networks, which generally have 3-phase power throughout the facilities, so we model the network using 4-conductor cables (3-phase conductors and a ground conductor). Industrial facilities typically have a large number of motors present, unlike home networks, so we incorporate motor models in our network.

Our previous work (Ref. [29]) combined top-down and bottom-up models to gain the advantages of both. It used a Bayesian inference and sampling method called Transitional Markov Chain Monte Carlo (TMCMC) to calibrate model parameters to data. TMCMC effectively calibrates up to approximately 20 parameters, but it became computationally intractable to calibrate more, and a large PLC network model can contain many more unknown parameters. This work instead uses a Bayesian inference method called mean-field variational inference to calibrate model parameters and uses a more complex industrial network model instead of the home network model in Ref. [29]. Our efforts begin with the analysis of a single network, and then expand to multiple network realizations, then address the challenge of inferring load types. This paper is organized as follows: Section 2 describes the PLC network model for an industrial network, Section 3 describes mean-field variational inference, Section 4 contains results, and Section 5 contains conclusions and suggestions for future work.

2. 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 three-phase 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.

2.1. 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.

2.2. 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:

\begin{matrix} R & = [\begin{matrix} 2 r & r & r \\ r & 2 r & r \\ r & r & 2 r \end{matrix}] w i t h r = \{\begin{matrix} \frac{1}{σ π r_{w}^{2}}, & if r_{w} \leq 2 δ \\ \frac{1}{2 r_{w}} \sqrt{\frac{μ_{0} f}{π σ}}, & if r_{w} > 2 δ \end{matrix} \end{matrix}

(1)

\begin{matrix} L & = [\begin{matrix} \frac{μ_{0}}{π} l o g (\frac{d_{1, 0}}{r_{w}}) & \frac{μ_{0}}{2 π} l o g (\frac{d_{1, 0} d_{2, 0}}{d_{1, 2} r_{w}}) & \frac{μ_{0}}{2 π} l o g (\frac{d_{1, 0} d_{3, 0}}{d_{1, 3} r_{w}}) \\ \frac{μ_{0}}{2 π} l o g (\frac{d_{1, 0} d_{2, 0}}{d_{1, 2} r_{w}}) & \frac{μ_{0}}{π} l o g (\frac{d_{2, 0}}{r_{w}}) & \frac{μ_{0}}{2 π} l o g (\frac{d_{2, 0} d_{3, 0}}{d_{2, 3} r_{w}}) \\ \frac{μ_{0}}{2 π} l o g (\frac{d_{1, 0} d_{3, 0}}{d_{1, 3} r_{w}}) & \frac{μ_{0}}{2 π} l o g (\frac{d_{2, 0} d_{3, 0}}{d_{2, 3} r_{w}}) & \frac{μ_{0}}{π} l o g (\frac{d_{3, 0}}{r_{w}}) \end{matrix}] \end{matrix}

(2)

\begin{matrix} C & = μ_{0} ϵ L^{- 1} \end{matrix}

(3)

\begin{matrix} G & = 2 π f t a n (δ_{L}) C \end{matrix}

(4)

where

r_{w}

is the conductor radius,

μ_{0} = 4 π \cdot 10^{- 7}

[H/m] is the vacuum permeability,

σ = 5.8 \cdot 10^{7}

[S/m] is the conductivity of copper,

ϵ = 3.19 \cdot 10^{- 11}

[F/m] is dielectric constant of PVC, f is the frequency. The loss tangent

t a n (δ_{L})

is assumed to be zero, so

G = 0

. The distances between conductors are

d_{0, 1} = d_{0, 3} = d_{1, 2} = d_{2, 3} = d_{c}

and

d_{0, 2} = d_{1, 3} = \sqrt{2} d_{c}

, where

d_{c}

is given by the relation:

d_{c} = 4 \cdot 10^{- 4} + 3.02 r_{w}

(5)

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

δ = \frac{1}{\sqrt{π μ_{0} σ f}}

(6)

Figure 2. Schematic of the cable cross-section.

2.3. Load Models

Each load consists of six component loads connecting each conductor in the three-phase cable with each other conductor, as shown in Figure 3. The components of each load connected to the ground cable are designated

Z_{G 1}

Z_{G 2}

, and

Z_{G 3}

, and the three phase-to-phase 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_{c o n s t}

and leakage capacitance

C_{c, l e a k}

. The three phase-to-phase components are equal to

R_{c o n s t}

, and the three ground-to-phase components are the impedance of the leakage capacitance. The load equations for a constant impedance load are

\begin{matrix} Z_{12} & = Z_{13} = Z_{23} = R_{c o n s t} \end{matrix}

(7)

\begin{matrix} Z_{G 1} & = Z_{G 2} = Z_{G 3} = \frac{1}{j ω C_{c, l e a k}} \end{matrix}

(8)

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}

ω_{0 s}

ζ_{s}

R_{p}

ω_{0 p}

ζ_{p}

Δ_{1}

Δ_{2}

, and

C_{d, l e a k}

. The load equations are

\begin{matrix} Z_{12} & = \frac{2 j \frac{ω}{ω_{0 p}} R_{p} ζ_{p}}{1 + 2 j \frac{ω}{ω_{0 p}} R_{p} ζ_{p} - \frac{ω^{2}}{ω_{0 p}^{2}}} + R_{s} + 2 j R_{s} ζ_{s} (\frac{ω}{ω_{0 s}} - \frac{ω_{0 s}}{ω}) \end{matrix}

(9)

\begin{matrix} Z_{13} & = Z_{12} (1 + Δ_{1}) \end{matrix}

(10)

\begin{matrix} Z_{23} & = Z_{12} (1 + Δ_{2}) \end{matrix}

(11)

\begin{matrix} Z_{G 1} & = Z_{G 2} = Z_{G 3} = \frac{1}{j ω C_{d, l e a k}} \end{matrix}

(12)

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:

\begin{matrix} ω_{0} & = \frac{1}{\sqrt{L C}} \end{matrix}

(13)

\begin{matrix} ζ & = \frac{R}{2} \sqrt{\frac{C}{L}} \end{matrix}

(14)

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_{m 1}

R_{m 2}

, and

C_{m, l e a k}

. However, the model is insensitive to

R_{m 2}

except at low frequencies, so its value is fixed at

R_{m 2} = 5 Ω

. The motor load formula can be reduced to

\begin{matrix} Z_{12} & = \frac{1}{3} \frac{1}{j ω (C_{m} + \frac{C_{m, l e a k}}{2}) + \frac{1}{R_{m 1}} + \frac{1}{Z^{'}}} \end{matrix}

(15)

\begin{matrix} Z_{13} & = Z_{23} = Z_{12} \end{matrix}

(16)

\begin{matrix} Z_{G 1} & = \frac{1}{3} \frac{1}{j ω C_{m, l e a k} + \frac{1}{j ω C_{m, l e a k} + \frac{R_{m 1} Z^{'}}{j ω C_{m} Z^{'} R_{m 1} + Z^{'} + R_{m 1}}}} \end{matrix}

(17)

\begin{matrix} Z_{G 2} & = Z_{G 3} = Z_{G 1} \end{matrix}

(18)

where

Z^{'} = j ω L_{m} + R_{m 2}

2.4. 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

\pm 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.

2.5. 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.

2.6. 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 \times 6

matrix, shown in Figure 8 and Equation (19).

[\begin{matrix} V_{1, i n} \\ V_{2, i n} \\ V_{3, i n} \\ I_{1, i n} \\ I_{2, i n} \\ I_{3, i n} \end{matrix}] = Φ [\begin{matrix} V_{1, o u t} \\ V_{2, o u t} \\ V_{3, o u t} \\ I_{1, o u t} \\ I_{2, o u t} \\ I_{3, o u t} \end{matrix}]

(19)

Each load can be put in the form of an admittance matrix

Y_{L}

using the equation:

Y_{L} = [\begin{matrix} \frac{1}{Z_{G 1}} + \frac{1}{Z_{12}} + \frac{1}{Z_{13}} & - \frac{1}{Z_{12}} & - \frac{1}{Z_{13}} \\ - \frac{1}{Z_{12}} & \frac{1}{Z_{G 2}} + \frac{1}{Z_{12}} + \frac{1}{Z_{23}} & - \frac{1}{Z_{23}} \\ - \frac{1}{Z_{13}} & - \frac{1}{Z_{23}} & \frac{1}{Z_{G 3}} + \frac{1}{Z_{13}} + \frac{1}{Z_{23}} \end{matrix}]

(20)

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).

Φ_{p a r a l l e l L o a d} = [\begin{matrix} 1 & 0 \\ - Y_{L} & 1 \end{matrix}]

(21)

Φ_{s e r i e s L o a d} = [\begin{matrix} 1 & - Z_{L} \\ 0 & 1 \end{matrix}]

(22)

where

Z_{L} = Y_{L}^{- 1}

Figure 8. Transmission matrix definition.

The network transfer function is computed over the frequency range [150 kHz, 30 MHz] for 200 discrete frequencies with logarithmic spacing. To obtain more simulated data, the resistances in the transmitter and receiver are varied, and the transfer function is computed for each combination of resistances, similarly to Ref. [29]. Table 2 and Table 3 show the resistance combinations used to compute transfer functions.

3. 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.

p (θ | d) \propto p (d | θ) p (θ)

(23)

where

d \in 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.

p (d | θ) = \prod_{i = 0}^{N_{d}} T (d_{i} - y_{i}, σ, ν)

(24)

where

y \in 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_{ϕ} (θ)

E L B O = E_{q_{ϕ} (θ)} [log p (d, θ) - log q_{ϕ} (θ)]

(25)

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.

E_{q_{ϕ} (θ)} [\nabla_{ϕ} q_{ϕ} (θ) (log p (d, θ) - log q_{ϕ} (θ))]

(26)

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_{ϕ} (θ)} [\nabla_{ϕ} (log q_{ϕ} (θ) \times 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].

4. 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 high-dimensional 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 two-norm 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.

4.1. 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_{c o n s t} = 10^{6} Ω

and parasitic capacitance

C_{c, l e a k} = 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.

4.2. 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.

4.3. 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.

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.

5. 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.

Author Contributions

Conceptualization, T.A.R.; methodology, D.S.C., R.W., J.C., C.S. and L.R.; software, D.S.C., R.W. and J.C.; validation, D.S.C.; formal analysis, D.S.C., R.W., and J.C.; writing—original draft preparation, D.S.C., R.W. and J.C.; writing—review and editing, D.S.C., R.W. and J.C.; supervision, D.S.C., C.S. and T.A.R.; project administration, T.A.R.; funding acquisition, T.A.R. All authors have read and agreed to the published version of the manuscript.

Funding

This work was supported by Defense Nuclear Nonproliferation Research and Development. Sandia National Laboratories is a multimission laboratory managed and operated by National Technology and Engineering Solutions of Sandia, LLC., a wholly owned subsidiary of Honeywell International, Inc., for the U.S. Department of Energy’s National Nuclear Security Administration under contract DE-NA-0003525. This paper describes objective technical results and analysis. Any subjective views or opinions that might be expressed in the paper do not necessarily represent the views of the U.S. Department of Energy or the United States Government.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

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:

BB	Broadband
ELBO	Evidence lower bound
MIMO	Multiple input multiple output
NB	Narrowband
PLC	Power line communications
RLC	Resistor-inductor-capacitor
TMCMC	Transisional Markov chain Monte Carlo
TL	Transmission line
UNB	Ultra narrowband

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

{\vec{V}}_{c} = H {\vec{V}}_{a}

(A1)

We define transfer functions through a part of the network as

\begin{matrix} {\vec{V}}_{b} & = H_{trans} {\vec{V}}_{a} \end{matrix}

(A2)

\begin{matrix} {\vec{V}}_{c} & = H_{nw} {\vec{V}}_{b} \end{matrix}

(A3)

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

Φ_{n w}

into submatrices and then substituting in the receiver load:

[\begin{matrix} \vec{V_{c}} \\ \vec{I_{c}} \end{matrix}] = [\begin{matrix} ϕ_{n w, 11} & ϕ_{n w, 12} \\ ϕ_{n w, 21} & ϕ_{n w, 22} \end{matrix}] [\begin{matrix} \vec{V_{b}} \\ \vec{I_{b}} \end{matrix}]

(A5)

We can also define submatrices of the inverse matrix

[\begin{matrix} \vec{V_{b}} \\ \vec{I_{b}} \end{matrix}] = [\begin{matrix} ϕ_{n w, 11}^{- 1} & ϕ_{n w, 12}^{- 1} \\ ϕ_{n w, 21}^{- 1} & ϕ_{n w, 22}^{- 1} \end{matrix}] [\begin{matrix} \vec{V_{c}} \\ \vec{I_{c}} \end{matrix}]

(A6)

Appendix A.1. Receiver Load

The receiver load matrix is defined as

Z_{rec}

where,

\vec{V_{c}} = Z_{rec} \vec{I_{c}}

(A7)

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:

\begin{matrix} \begin{matrix} V_{c, 1} = Z_{R 1} I_{c, 1} + (I_{c, 1} + I_{c, 2} + I_{c, 3}) Z_{R G} \\ V_{c, 2} = Z_{R 2} I_{c, 1} + (I_{c, 1} + I_{c, 2} + I_{c, 3}) Z_{R G} \\ V_{c, 3} = Z_{R 3} I_{c, 1} + (I_{c, 1} + I_{c, 2} + I_{c, 3}) Z_{R G} \end{matrix} \end{matrix}

(A8)

Z_{rec} = [\begin{matrix} Z_{R G} + Z_{R 1} & Z_{R G} & Z_{R G} \\ Z_{R G} & Z_{R G} + Z_{R 2} & Z_{R G} \\ Z_{R G} & Z_{R G} & Z_{R G} + Z_{R 3} \end{matrix}]

(A9)

Substituting Equation (A7) into the first row of Equation (A6) yields:

\vec{V_{b}} = ϕ_{n w, 11}^{- 1} \vec{V_{c}} + ϕ_{n w, 12}^{- 1} Z_{rec}^{- 1} \vec{V_{c}}

(A10)

Reorganizing gives:

H_{nw} = {(ϕ_{n w, 11}^{- 1} + ϕ_{n w, 12}^{- 1} Z_{rec}^{- 1})}^{- 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:

Y_{nw} \vec{V_{b}} = \vec{I_{b}}

(A12)

where:

Y_{nw} = - {(ϕ_{n w, 12} - Z_{rec} ϕ_{n w, 22})}^{- 1} (ϕ_{n w, 11} - Z_{rec} ϕ_{n w 21})

(A13)

Appendix A.2. Transmitter Computation

Three equations can be formed by applying Krichoff’s current law to the three conductors in the transmitter:

\begin{matrix} \begin{matrix} I_{b 1} = & I_{a 1} - \frac{V_{b 1}}{Z_{T G 1}} - \frac{V_{b 1} - V_{b 2}}{Z_{T 12}} - \frac{V_{b 1} - V_{b 3}}{Z_{T 13}} \\ = & I_{a 1} - \frac{Z_{T 12} Z_{T 13} + Z_{T G 1} Z_{T 13} + Z_{T G 1} Z_{T 12}}{Z_{T G 1} Z_{T 12} Z_{T 13}} V_{b 1} + \frac{1}{Z_{T 12}} V_{b 2} + \frac{1}{Z_{T 13}} V_{b 3} \\ I_{b 2} = & I_{a 2} - \frac{V_{b 2}}{Z_{T G 2}} - \frac{V_{a 2} - V_{a 1}}{Z_{T 12}} - \frac{V_{b 2} - V_{b 3}}{Z_{T 23}} \\ = & I_{a 2} + \frac{1}{Z_{T 12}} V_{b 1} - \frac{Z_{T G 2} Z_{T 12} + Z_{T G 2} Z_{T 23} + Z_{T 12} Z_{T 23}}{Z_{T G 2} Z_{T 12} Z_{T 23}} V_{b 2} + \frac{1}{Z_{T 23}} V_{b 3} \\ I_{b 3} = & I_{a 3} - \frac{V_{b 3}}{Z_{T G 3}} - \frac{V_{b 3} - V_{b 1}}{Z_{T 13}} - \frac{V_{b 3} - V_{b 2}}{Z_{T 23}} \\ = & I_{a 3} + \frac{1}{Z_{T 13}} V_{b 1} + \frac{1}{Z_{T 23}} V_{b 2} - \frac{Z_{T G 3} Z_{T 13} + Z_{T G 3} Z_{T 23} + Z_{T 13} Z_{T 23}}{Z_{T G 3} Z_{T 13} Z_{T 23}} V_{b 3} \end{matrix} \end{matrix}

(A14)

The impedance law across resistors

Z_{T 0}

\frac{\vec{V_{a}} - \vec{V_{b}}}{Z_{T 0}} = \vec{I_{a}}

(A15)

Equations (A15) and (A12) are substituted into (A14) to eliminate the variables

\vec{I_{a}}

and

\vec{I_{b}}

\begin{matrix} \begin{matrix} Y_{n w, 11} Y_{b 1} + Y_{n w, 12} Y_{b 2} + Y_{n w, 13} Y_{b 3} = & \frac{V_{a 1} - V_{b 1}}{Z_{T 0}} - Y_{T α} V_{b 1} + \frac{1}{Z_{T 12}} V_{b 2} + \frac{1}{Z_{T 13}} V_{b 3} \\ Y_{n w, 21} Y_{b 1} + Y_{n w, 22} Y_{b 2} + Y_{n w, 23} Y_{b 3} = & \frac{V_{a 2} - V_{b 2}}{Z_{T 0}} - \frac{1}{Z_{T 12}} V_{b 1} - Y_{T β} V_{b 2} + \frac{1}{Z_{T 23}} V_{b 3} \\ Y_{n w, 31} Y_{b 1} + Y_{n w, 32} Y_{b 2} + Y_{n w, 33} Y_{b 3} = & \frac{V_{a 3} - V_{b 3}}{Z_{T 0}} + \frac{1}{Z_{T 13}} V_{b 1} + \frac{1}{Z_{T 23}} V_{b 2} - Y_{T γ} V_{b 3} \end{matrix} \end{matrix}

(A16)

where

\begin{matrix} \begin{matrix} Y_{T α} & = \frac{Z_{T 12} Z_{T 13} + Z_{T G 1} Z_{T 13} + Z_{T G 1} Z_{T 12}}{Z_{T G 1} Z_{T 12} Z_{T 13}} \\ Y_{T β} & = \frac{Z_{T G 2} Z_{T 12} + Z_{T G 2} Z_{T 23} + Z_{T 12} Z_{T 23}}{Z_{T G 2} Z_{T 12} Z_{T 23}} \\ Y_{T γ} & = \frac{Z_{T G 3} Z_{T 13} + Z_{T G 3} Z_{T 23} + Z_{T 13} Z_{T 23}}{Z_{T G 3} Z_{T 13} Z_{T 23}} \end{matrix} \end{matrix}

(A17)

This can be reordered to

[\begin{matrix} 1 + Z_{T 0} Y_{n w, 11} + Z_{T 0} Y_{T α} & Z_{T 0} Y_{n w, 12} - \frac{Z_{T 0}}{Z_{T 12}} & Z_{T 0} Y_{n w, 13} - \frac{Z_{T 0}}{Z_{T 13}} \\ Z_{T 0} Y_{n w, 21} - \frac{Z_{T 0}}{Z_{T 21}} & 1 + Z_{T 0} Y_{n w, 22} + Z_{T 0} Y_{T β} & Z_{T 0} Y_{n w, 23} - \frac{Z_{T 0}}{Z_{T 23}} \\ Z_{T 0} Y_{n w, 31} - \frac{Z_{T 0}}{Z_{T 31}} & Z_{T 0} Y_{n w, 32} - \frac{Z_{T 0}}{Z_{T 32}} & 1 + Z_{T 0} Y_{n w, 33} + Z_{T 0} Y_{T γ} \end{matrix}] \vec{V_{b}} = \vec{V_{a}}

(A18)

By comparing to Equation (A2), it is evident that the matrix in Equation (A18) is the inverse of

H_{trans}

. The transfer function is then computed using the chain rule in Equation (A4) with Equation (A11) and Equation (A18).

References

Galli, S.; Scaglione, A.; Wang, Z. For the grid and through the grid: The role of power line communications in the smart grid. Proc. IEEE 2011, 99, 998–1027. [Google Scholar] [CrossRef] [Green Version]
Nassar, M.; Lin, J.; Mortazavi, Y.; Dabak, A.; Kim, I.H.; Evans, B.L. Local utility power line communications in the 3–500 kHz band: Channel impairments, noise, and standards. IEEE Signal Process. Mag. 2012, 29, 116–127. [Google Scholar] [CrossRef]
Sonmez, M.A.; Zehir, M.A.; Bagriyanik, M.; Nak, O. Impulsive noise survey on power line communication networks up to 125 kHz for smart metering infrastructure in systems with solar inverters in Turkey. In Proceedings of the 2013 International Conference on Renewable Energy Research and Applications (ICRERA), Madrid, Spain, 20–23 October 2013; pp. 705–710. [Google Scholar]
Coelho, P.; Gomes, M.; Moreira, C. Smart metering technology. In Microgrids Design and Implementation; Springer: Berlin/Heidelberg, Germany, 2019; pp. 97–137. [Google Scholar]
Van de Kaa, G.; Fens, T.; Rezaei, J.; Kaynak, D.; Hatun, Z.; Tsilimeni-Archangelidi, A. Realizing smart meter connectivity: Analyzing the competing technologies power line communication, mobile telephony, and radio frequency using the best worst method. Renew. Sustain. Energy Rev. 2019, 103, 320–327. [Google Scholar] [CrossRef]
Mlỳnek, P.; Rusz, M.; Benešl, L.; Sláčik, J.; Musil, P. Possibilities of broadband power line communications for smart home and smart building applications. Sensors 2021, 21, 240. [Google Scholar] [CrossRef]
Masood, B.; Khan, M.A.; Baig, S.; Song, G.; Rehman, A.U.; Rehman, S.U.; Asif, R.M.; Rasheed, M.B. Investigation of deterministic, statistical and parametric NB-PLC channel modeling techniques for advanced metering infrastructure. Energies 2020, 13, 3098. [Google Scholar] [CrossRef]
Canete, F.J.; Cortes, J.A.; Diez, L.; Entrambasaguas, J.T. A channel model proposal for indoor power line communications. IEEE Commun. Mag. 2011, 49, 166–174. [Google Scholar] [CrossRef]
Tonello, A.M.; Versolatto, F. Bottom-up statistical PLC channel modeling—Part I: Random topology model and efficient transfer function computation. IEEE Trans. Power Deliv. 2011, 26, 891–898. [Google Scholar] [CrossRef]
Paul, C.R. Analysis of Multiconductor Transmission Lines; John Wiley & Sons: Hoboken, NJ, USA, 2007. [Google Scholar]
Zhu, W.; Zhu, X.; Lim, E.; Huang, Y. State-of-art power line communications channel modelling. Procedia Comput. Sci. 2013, 17, 563–570. [Google Scholar] [CrossRef] [Green Version]
Marrocco, G.; Statovci, D.; Trautmann, S. A PLC broadband channel simulator for indoor communications. In Proceedings of the 2013 IEEE 17th International Symposium on Power Line Communications and Its Applications, Johannesburg, South Africa, 24–27 March 2013; pp. 321–326. [Google Scholar]
Fang, X.; Wang, N.; Gulliver, T.A. A PLC channel model for home area networks. Energies 2018, 11, 3344. [Google Scholar] [CrossRef] [Green Version]
Aalamifar, F.; Schlögl, A.; Harris, D.; Lampe, L. Modelling power line communication using network simulator-3. In Proceedings of the 2013 IEEE Global Communications Conference (GLOBECOM), Atlanta, GA, USA, 9–13 December 2013; pp. 2969–2974. [Google Scholar]
Passerini, F.; Righini, D.; Tonello, A.M. A bottom-up PLC channel model that includes radiation effects. In Proceedings of the 2018 IEEE International Symposium on Power Line Communications and its Applications (ISPLC), Manchester, UK, 8–11 April 2018; pp. 1–6. [Google Scholar]
Tonello, A.M.; Versolatto, F.; Pittolo, A. In-home power line communication channel: Statistical characterization. IEEE Trans. Commun. 2014, 62, 2096–2106. [Google Scholar] [CrossRef]
Versolatto, F.; Tonello, A.M. A MIMO PLC random channel generator and capacity analysis. In Proceedings of the 2011 IEEE International Symposium on Power Line Communications and Its Applications, Udine, Italy, 3–6 April 2011; pp. 66–71. [Google Scholar]
Banwell, T.; Galli, S. A novel approach to the modeling of the indoor power line channel part I: Circuit analysis and companion model. IEEE Trans. Power Deliv. 2005, 20, 655–663. [Google Scholar] [CrossRef]
Galli, S.; Banwell, T. A novel approach to the modeling of the indoor power line channel-Part II: Transfer function and its properties. IEEE Trans. Power Deliv. 2005, 20, 1869–1878. [Google Scholar] [CrossRef]
Corchado, J.A.; Cortés, J.A.; Canete, F.J.; Díez, L. An MTL-based channel model for indoor broadband MIMO power line communications. IEEE J. Sel. Areas Commun. 2016, 34, 2045–2055. [Google Scholar] [CrossRef]
Anatory, J.; Theethayi, N.; Thottappillil, R. Power-line communication channel model for interconnected networks—Part I: Two-conductor system. IEEE Trans. Power Deliv. 2008, 24, 118–123. [Google Scholar] [CrossRef]
Anatory, J.; Theethayi, N.; Thottappillil, R. Power-line communication channel model for interconnected networks—Part II: Multiconductor system. IEEE Trans. Power Deliv. 2008, 24, 124–128. [Google Scholar] [CrossRef]
Franek, L.; Fiedler, P. A multiconductor model of power line communication in medium-voltage lines. Energies 2017, 10, 816. [Google Scholar] [CrossRef] [Green Version]
Versolatto, F.; Tonello, A.M. An MTL theory approach for the simulation of MIMO power-line communication channels. IEEE Trans. Power Deliv. 2011, 26, 1710–1717. [Google Scholar] [CrossRef]
Zhai, M.Y. Transmission characteristics of low-voltage distribution networks in China under the smart grids environment. IEEE Trans. Power Deliv. 2010, 26, 173–180. [Google Scholar] [CrossRef]
Tomasoni, A.; Riva, R.; Bellini, S. Spatial correlation analysis and model for in-home MIMO power line channels. In Proceedings of the 2012 IEEE International Symposium on Power Line Communications and Its Applications, Beijing, China, 27–30 March 2012; pp. 286–291. [Google Scholar]
Berger, L.T.; Schwager, A.; Pagani, P.; Schneider, D.M. MIMO power line communications. IEEE Commun. Surv. Tutor. 2014, 17, 106–124. [Google Scholar] [CrossRef]
Hao, L.; Guo, J. A MIMO-OFDM scheme over coupled multi-conductor power-line communication channel. In Proceedings of the 2007 IEEE International Symposium on Power Line Communications and Its Applications, Pisa, Italy, 26–28 March 2007; pp. 198–203. [Google Scholar]
Ching, D.S.; Safta, C.; Reichardt, T.A. Sensitivity-informed bayesian inference for home PLC network models with unknown parameters. Energies 2021, 14, 2402. [Google Scholar] [CrossRef]
Ferreira, H.C.; Lampe, L.; Newbury, J.; Swart, T.G. Power Line Communications: Theory and Applications for Narrowband and Broadband Communications Over Power Lines; John Wiley & Sons: Hoboken, NJ, USA, 2010. [Google Scholar]
Antoniali, M.; Tonello, A.M. Measurement and characterization of load impedances in home power line grids. IEEE Trans. Instrum. Meas. 2013, 63, 548–556. [Google Scholar] [CrossRef]
Kharraz, M.A.O.; Lavenu, C.; Jensen, P.; Picard, D.; Serhir, M. Characterization of the input impedance of household appliances in the FCC frequency band. In Proceedings of the 2017 IEEE International Symposium on Power Line Communications and Its Applications (ISPLC), Madrid, Spain, 3–5 April 2017; pp. 1–6. [Google Scholar]
Pasdar, A.M.; Cavdar, I.H.; Sozer, Y. Power-line impedance estimation at FCC band based on intelligent home appliances status detection algorithm through their individual energy and impedance signatures. IEEE Trans. Power Deliv. 2014, 29, 1407–1416. [Google Scholar] [CrossRef]
Sudiarto, B.; Widyanto, A.N.; Hirsch, H. The mutual influence of appliances on the disturbance in the frequency range of 9–150 kHz produced by household appliances in the low voltage network. In Proceedings of the 2017 International Symposium on Electromagnetic Compatibility-EMC EUROPE, Angers, France, 4–7 September 2017; pp. 1–6. [Google Scholar]
Moreno, Y.; Almandoz, G.; Egea, A.; Arribas, B.; Urdangarin, A. Analysis of Permanent Magnet Motors in High Frequency—A Review. Appl. Sci. 2021, 11, 6334. [Google Scholar] [CrossRef]
Boglietti, A.; Carpaneto, E. Induction motor high frequency model. In Proceedings of the Conference Record of the 1999 IEEE Industry Applications Conference. Thirty-Forth IAS Annual Meeting (Cat. No. 99CH36370), Phoenix, AZ, USA, 3–7 October 1999; Volume 3, pp. 1551–1558. [Google Scholar]
Boglietti, A.; Cavagnino, A.; Lazzari, M. Experimental high-frequency parameter identification of AC electrical motors. IEEE Trans. Ind. Appl. 2007, 43, 23–29. [Google Scholar] [CrossRef] [Green Version]
Miloudi, H.; Bendaoud, A.; Miloudi, M.; Dickmann, S.; Schenke, S. Common mode and differential mode characteristics of AC motor for EMC analysis. In Proceedings of the 2016 International Symposium on Electromagnetic Compatibility-EMC EUROPE, Wroclaw, Poland, 25–29 July 2016; pp. 765–769. [Google Scholar]
Miloudi, H.; Bendaoud, A.; Miloudi, M. A method for modeling a common-mode impedance for the AC motor. Elektrotehniski Vestn. 2017, 84, 241–246. [Google Scholar]
Miloudi, H.; Miloudi, M.; Gourbi, A.; Bermaki, M.; Bendaoud, A.; Zeghoudi, A. A high-frequency modeling of AC motor in a frequency range from 40 Hz to 110 MHz. Electr. Eng. Electromech. 2022, 6, 3–7. [Google Scholar] [CrossRef]
Blei, D.M.; Kucukelbir, A.; McAuliffe, J.D. Variational Inference: A Review for Statisticians. J. Am. Stat. Assoc. 2017, 112, 859–877. [Google Scholar] [CrossRef] [Green Version]
Salimans, T.; Kingma, D.P.; Welling, M. Markov Chain Monte Carlo and Variational Inference: Bridging the Gap. In Proceedings of the International Conference on Machine Learning, Lille, France, 6–11 July 2015. [Google Scholar]
Ranganath, R.; Gerrish, S.; Blei, D.M. Black Box Variational Inference. In Proceedings of the 17th International Conference on Artificial Intelligence and Statistics, Reykjavik, Iceland, 22–25 April 2014. [Google Scholar]
Mnih, A.; Gregor, K. Neural Variational Inference and Learning in Belief Networks. In Proceedings of the International Conference on Machine Learning, PMLR, Beijing, China, 21–26 June 2014. [Google Scholar]
Bingham, E.; Chen, J.P.; Jankowiak, M.; Obermeyer, F.; Pradhan, N.; Karaletsos, T.; Singh, R.; Szerlip, P.; Horsfall, P.; Goodman, N.D. Pyro: Deep Universal Probabilistic Programming. J. Mach. Learn. Res. 2018, 20, 973–978. [Google Scholar]
Phan, D.; Pradhan, N.; Jankowiak, M. Composable Effects for Flexible and Accelerated Probabilistic Programming in NumPyro. arXiv 2019, arXiv:1912.11554. [Google Scholar]
Duchi, J.; Hazan, E.; Singer, Y. Adaptive subgradient methods for online learning and stochastic optimization. J. Mach. Learn. Res. 2011, 12, 2121–2159. [Google Scholar]

Figure 1. Network topology. Blue circles represent loads, and black squares are junction boxes.

Figure 3. Schematic of load connections to wiring.

Figure 4. Double RLC phase-to-phase load circuit.

Figure 5. Motor model. (a) Phase-to-phase circuit. (b) Phase-to-ground circuit.

Figure 6. Transmitter model.

Figure 7. Receiver model.

Figure 9. Mean value of parameters vs iteration (left) and variance vs iteration (right). Each line corresponds to a different parameter.

Figure 10. True transfer function (without added noise) and mean of transfer function computed from posterior with 95% confidence bounds.

Figure 11. Transfer function from the original network and network modified by powering offloads in two rooms.

Figure 12. True transfer function component

H_{1, 1} (f)

and inferred transfer function for 10 different realizations of the network. Each color corresponds to a different realization. The dots are the true transfer function and lines are the inferred transfer function.

Figure 12. True transfer function component

H_{1, 1} (f)

Table 1. Parameter ranges and distributions.

Parameter	Minimum	Maximum	Distribution Type
$r_{w}$ [mm] connected to service panel	1.03	2.06	uniform
$r_{w}$ [mm] in rooms	0.81	1.29	uniform
$l_{w}$ [m]	2	20	uniform
$R_{c o n s t}$ [ $Ω$ ]	10	200	log-uniform
$C_{c, l e a k}$ [nF]	0.1	2.0	uniform
$R_{s}$ [ $Ω$ ]	10	3000	log-uniform
$ω_{0 s}$ [Mrad/s]	0	30	uniform
$ζ_{s}$	0.1	2	uniform
$R_{p}$ [ $Ω$ ]	10	3000	log-uniform
$ω_{0 p}$ [Mrad/s]	0	30	uniform
$ζ_{p}$	0.1	2	uniform
$Δ_{1}$	−0.1	0.1	uniform
$Δ_{2}$	−0.1	0.1	uniform
$C_{d, l e a k}$ [F]	0.1	2.0	uniform
$C_{m}$ [nF]	0.1	1.0	uniform
$L_{m}$ [mH]	5	20	uniform
$R_{m 1}$ [ $Ω$ ]	2000	15,000	uniform
$C_{m, l e a k}$ [nF]	0.2	5.0	uniform

Table 2. Transmitter resistance combinations.

Parameter	Transmitter 1	Transmitter 2
$Z_{T 0}$	$50 Ω$	$100 Ω$
$Z_{T G 1}$	$50 Ω$	$100 Ω$
$Z_{T G 2}$	$50 Ω$	$100 Ω$
$Z_{T G 3}$	$50 Ω$	$50 Ω$
$Z_{T 12}$	$100 Ω$	$50 Ω$
$Z_{T 13}$	$100 Ω$	$100 Ω$
$Z_{T 23}$	$100 Ω$	$100 Ω$

Table 3. Receiver resistance combinations.

Parameter	Receiver 1	Receiver 2
$Z_{R G}$	$50 Ω$	$100 Ω$
$Z_{R 1}$	$50 Ω$	$100 Ω$
$Z_{R 2}$	$50 Ω$	$100 Ω$
$Z_{R 3}$	$50 Ω$	$50 Ω$

Table 4. Mean parameter values from posterior in order of decreasing sensitivity, ordered left to right, then down. The true values are 0.25.

0.259	0.251	0.250	0.245	0.256	0.202
0.253	0.252	0.256	0.342	0.174	0.249
0.236	0.261	0.249	0.240	0.284	0.233
0.264	0.239	0.199	0.242	0.243	0.308
0.247	0.267	0.244	0.259	0.257	0.243
0.253	0.245	0.247	0.251	0.344	0.257
0.251	0.237	0.252	0.259	0.241	0.113
0.260	0.233	0.253	0.209	0.272	0.249
0.237	0.295	0.163	0.161	0.245	0.119
0.414	0.403	0.199	0.744	0.455	0.367
0.254	0.183	0.207	0.154	0.352	0.098
0.596	0.507	0.136	0.117	0.246	0.140

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
1–10	0.029
10–20	0.080
20–40	0.088
40–60	0.176

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
R1-O2	Constant	{Constant *, Double RLC}
R1-O3	Motor	Motor
R1-O4	Constant	{Constant *, Double RLC}
R2-O1	Motor	Motor
R2-O2	Motor	Motor
R2-O3	Motor	Motor
R2-O4	Motor	Motor
R3-O1	Constant	{Constant *, Double RLC}
R3-O2	Constant	{Constant *, Double RLC}
R3-O3	Constant	{Constant *, Double RLC}
R3-O4	Motor	Motor
R4-O1	Double RLC	{Constant *, Double RLC}
R4-O2	Motor	Motor
R4-O3	Motor	Motor
R4-O4	Double RLC	{Constant *, Double RLC}
R5-O1	Motor	Motor
R5-O2	Constant	{Constant *, Double RLC}
R5-O3	Double RLC	{Constant *, Double RLC}
R5-O4	Motor	Motor
R6-O1	Constant	{Constant *, Double RLC}
R6-O2	Motor	Motor
R6-O3	Motor	Motor

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

© 2023 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

Share and Cite

MDPI and ACS Style

Wonnacott, R.; Ching, D.S.; Chilleri, J.; Safta, C.; Rashkin, L.; Reichardt, T.A. Industrial PLC Network Modeling and Parameter Identification Using Sensitivity Analysis and Mean Field Variational Inference. Sensors 2023, 23, 2416. https://doi.org/10.3390/s23052416

AMA Style

Wonnacott R, Ching DS, Chilleri J, Safta C, Rashkin L, Reichardt TA. Industrial PLC Network Modeling and Parameter Identification Using Sensitivity Analysis and Mean Field Variational Inference. Sensors. 2023; 23(5):2416. https://doi.org/10.3390/s23052416

Chicago/Turabian Style

Wonnacott, Raelynn, David S. Ching, John Chilleri, Cosmin Safta, Lee Rashkin, and Thomas A. Reichardt. 2023. "Industrial PLC Network Modeling and Parameter Identification Using Sensitivity Analysis and Mean Field Variational Inference" Sensors 23, no. 5: 2416. https://doi.org/10.3390/s23052416

APA Style

Wonnacott, R., Ching, D. S., Chilleri, J., Safta, C., Rashkin, L., & Reichardt, T. A. (2023). Industrial PLC Network Modeling and Parameter Identification Using Sensitivity Analysis and Mean Field Variational Inference. Sensors, 23(5), 2416. https://doi.org/10.3390/s23052416

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Menu

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

Abstract

1. Introduction

2. PLC Network Model

2.1. Network Topology

2.2. Cable Model

2.3. Load Models

2.4. Network Generation

2.5. Transmitter and Receiver Models

2.6. Model Solver

3. Mean Field Variational Inference

4. Results

4.1. Demonstration on Single Network

4.2. Demonstration on Multiple Networks

4.3. Inferring Load Types

5. Conclusions

Author Contributions

Funding

Institutional Review Board Statement

Informed Consent Statement

Data Availability Statement

Acknowledgments

Conflicts of Interest

Abbreviations

Appendix A. Transfer Function Computation

Appendix A.1. Receiver Load

Appendix A.2. Transmitter Computation

References

Share and Cite

Article Metrics

Article Access Statistics

Further Information

Guidelines

MDPI Initiatives

Follow MDPI