Conversion of descriptor representations to state-space form: an extension of the shuffle algorithm

ABSTRACT This paper proposes a systematic procedure for the determination of state-space models from an available descriptor representation of a linear dynamic system. The goal is to determine a state equation, a set of algebraic equations and an output equation in terms of the state and input variables. It is shown that standard methods may fail to convert the descriptor representation to state-space form, even for simple electrical circuit models obtained from Kirchoff's laws and constitutive element equations. A novel procedure to address this problem is then proposed as an extension of the classic shuffle algorithm combined with a singular value decomposition approach. In addition to an illustrative example involving a simple electrical circuit, the proposed method is employed in a case study involving the modelling of three-dimensional RLC networks with a large number of components.


Introduction
Descriptor representations naturally arise in the modelling of dynamic systems with algebraic constraints from first principles (Luenberger, 1977;Müller, 2000). Examples can be found in several application areas such as robotics, power systems, microelectromechanical devices and many others (Wong & Chu, 2008). In chemical processes, for example, differential equations are used to describe the dynamic balances of mass and energy, whereas algebraic equations are used to express thermodynamic equilibrium and steady-state assumptions. In the case of electrical circuits, descriptor representations follow directly from the algebraic constraints imposed by Kirchoff 's laws and the differential relations between voltage and current in reactive circuit elements. For this reason, much research has been conducted on the analysis and design methods for descriptor models (Benner, Sima, & Voigt, 2012;Cobb, 1981;Kazantzidou & Ntogramatzidis, 2016).
In this context, an important problem consists of obtaining a state-space model in standard form from a given descriptor representation. Such a procedure may be useful to facilitate the simulation and analysis of dynamic properties of the system under consideration, as well as to employ existing state-space control design techniques and software.
An indirect conversion to state-space form can be carried out by first computing the transfer matrix CONTACT Sillas Hadjiloucas s.hadjiloucas@reading.ac.uk associated to the descriptor model (Varga, 1989;Misra, 1989) and then using a canonical state-space realisation (Chen, 1984). However, the relation between the state vector and the original descriptor variables may not be straightforward, which can be an inconvenience for the interpretation of results obtained with such model. A more direct approach consists of using elementary matrix operations to determine a linear transformation that relates the descriptor vector with a suitable state vector (Luenberger, 1977). A classic method for this purpose is the so-called 'shuffle algorithm' proposed in Luenberger (1978). Alternatively, singular value decomposition (SVD) can be used to obtain modified descriptor equations from which purely algebraic dependencies can be extracted (Bender & Laub, 1987;Geromel & Palhares, 2011;Safonov, Chiang, & Limebeer, 1987). This allows rewriting a subset of the equations in the form of a statespace representation.
A key condition required in both the shuffle algorithm and the SVD-based method is the invertibility of one of the transformed matrices. The present paper shows that this invertibility condition is not always satisfied for electrical circuit models obtained from the systematic use of Kirchoff's laws, even in the case of very simple networks. The reason for this shortcoming is discussed, and a novel method is proposed as an extension of the shuffle algorithm combined with the SVD-based approach.
In addition to an illustrative example involving a simple circuit, the proposed method is employed in a case study concerning the modelling of three-dimensional (3-D) RLC networks with a large number of components. Within this scope, the present contribution also bridges a gap with respect to previous work concerning the modelling of 3-D RC networks (Galvao, Hadjiloucas, Kienitz, Paiva, & Afonso, 2013;Galvão et al., 2013;Jacyntho et al., 2016). Indeed, the results obtained herein indicate that the standard shuffle and SVD-based methods for conversion to state-space form may fail in the case of 3-D RLC networks.
The remainder of this paper is organised as follows. Section 2 states the problem of conversion from descriptor representation to state-space form and describes the standard shuffle and SVD-based methods that can be used for this purpose. Section 3 presents two examples involving simple RC and RL electrical circuits in order to illustrate a case in which the standard methods can be applied, as well as a case in which they cannot be applied. Section 4 introduces the proposed method, which is then illustrated with the example in which the standard methods could not be employed. The 3-D RLC network case study is presented in Section 5, along with a brief discussion of possible applications associated with this type of network. Finally, concluding remarks are given in Section 6 .

Background
Consider a linear dynamic system represented by a descriptor equation of the form and an output equation where x(t ) ∈ N , u(t ) ∈ m and y(t ) ∈ r are the descriptor, input and output vectors, respectively, and A, The problem under consideration consists of determining a state equation, a set of algebraic equations and an output equation in terms of the state and input variables. The state equation and set of algebraic equations shall be equivalent to the available descriptor representation in the sense that the descriptor vector can be readily calculated from the input vector, the solution to the state equation and the solution to the set of algebraic equations.
This problem can be addressed by using the shuffle algorithm, which was originally proposed in Luenberger (1978) for the discrete-time case, but can also be applied to continuous-time models (Gerding, 2004). Initially, elementary row operations (i.e. linear combinations of rows) are carried out in order to rewrite Equation (1) in the form whereÊ 1 is an (n × N) matrix of rank n andÂ 1 ,Â 2 ,B 1 andB 2 are matrices of dimensions (n × N), which can be differentiated with respect to time to obtain A 2ẋ (t ) = −B 2u (t ). This expression can be inserted in (3) in place of the last N − n rows to obtain This is not a state equation in standard form, in view of the dependence on the input derivative termu(t ). However, a state-space model can be obtained by choosing a state vector z(t ) ∈ n as Indeed, from (4) and (7), one arrives at and thus where Q and R are matrices of dimensions (N × n) and (N × (N − n)) obtained as In addition, from (3) and (7), it follows thaṫ From (9) and (11), a state equation for z(t) can be written asż where F =Â 1 Q, G =B 1 −Â 1 RB 2 . Finally, from (2) and (9), the output equation becomes with H = CQ, L = D − CRB 2 . The procedure described above demands the inver- T . An alternative method, which can be more numerically reliable for large N, involves the transformation of the descriptor representation into an SVD coordinate system (Bender & Laub, 1987;Geromel & Palhares, 2011;Ishihara & Terra, 2002;Safonov et al., 1987). For this purpose, E is decomposed as where 1 is an (n × n) diagonal matrix of non-zero singular values, p = N − n is the number of singular values equal to zero, and U, V are (N × N) unitary matrices, i.e. U T = U −1 and V T = V −1 . By substituting (14) for E in (1) and pre-multiplying both sides of the equation by U T , it follows that where z(t ) ∈ n is a candidate state vector and w(t ) ∈ p is a vector of additional variables. Since V T = V −1 , (16) can alternatively be written as From (15)-(17), it follows that By definingÃ = U T AV ,B = U T B and choosing appropriate partitions for these matrices, one may write If the (p × p) matrixÃ 22 is invertible, an algebraic equation for w(t) can be obtained from (21) as After substituting (22) for w(t) in (20) and premultiplying both sides by −1 1 , one arrives at a state equation of the form (12) with Moreover, from (2), (17) and (22), with an appropriate partition for matrices C and V, it follows that which is an output equation of the form (13) with Finally, in view of (17) and (22), the descriptor vector x(t) can be calculated from z(t) and u(t) as Remark 2.1: The use of SVD to convert the original descriptor equation (1) into (20), (21) is a wellknown form of restricted equivalent transformation (see Appendix). Indeed, decomposition (14) is a particular case of the transformation (A3) in Appendix, with Q −1 = U and P −1 = V T . The transformed model obtained through this SVD procedure is said to be in the dynamics decomposition form (Duan, 2010). Remark 2.2: A problem of practical relevance is to determine a state transformation T such that where = [ 1 2 ] is some matrix of interest and adequate dimension. Substituting (28) for x(t) in (29) results in Thus, if = [ 1 2 ]satisfies the condition it follows that T will be given by which can be used as transformation matrix if it is nonsingular.

Illustrative examples
This section presents two examples involving simple electrical circuits. The first example illustrates a case in which the descriptor model can be converted into state-space form by using either the shuffle algorithm or the SVDbased method described in Section 2. In the second example, this will not be possible. Such a problem will serve as a motivation for the proposed method, which will be introduced in Section 4. In what follows, resistances (R), inductances (L) and capacitances (C) will be written with upright letters, in order to avoid confusion with the matrices defined in Section 2. Moreover, the notation d/dt for time derivative will be employed to facilitate reading. All the numerical calculations were carried out by using the MATLAB R software. Figure 1 present the RC circuit considered in this first example. The input signal u(t) corresponds to the source voltage and the output signal y(t) will be taken as the current i S (t). By applying Kirchoff 's current and voltage laws and noting that the current through the capacitors is given by C 1 dv C1 (t )/dt and C 2 dv C2 (t )/dt, it follows that

First example: RC circuit
These equations can be cast into the descriptor form (1) with and Since y(t) = i S (t), the output equation is of the form (2) with In this case, the descriptor equation is already in the form (3) assumed in the shuffle algorithm, with As can be seen, [Ê T 1Â T 2 ] T is non-singular for any positive values of R S , C 1 , C 2 . Therefore, from (7) and (40), it follows that the state-space model can be written with a single state variable In what follows, the numerical procedures for conversion to state-space form will be illustrated with R S = C 1 = C 2 = 1 (normalised units). With these component values, the expression for the state variable z(t) becomes The Q, R matrices defined in (10) are calculated as and the resulting model matrices are given by A similar result can be obtained by using the SVD procedure. For this purpose, matrix E is decomposed as where 1 = √ 2. Since n = 1 (only one singular value different from zero), the state z(t) is a single variable. From (16) and (45), it can be seen that , which corresponds to (43) up to a scale factor.
As can be seen in (46),Ã 22 is a non-singular (3 × 3) matrix. By using (23), (24), (26), (27) and (46), the statespace model matrices are calculated as F = −1/2, G = √ 2/2, H = − √ 2/2, L = 1. These values correspond to those obtained above, up to scale factors in the input gain G and output gain H, which result from the definition of the state variable z(t) in this case. In fact, the two statespace models are different realisations of the same transfer function (s), which is given by where U(s) and Y(s) denote the Laplace transforms of the input signal u(t) and the output signal y(t), respectively. Since u(t) and y(t) are the source voltage and current in Figure 1, the transfer function (47) corresponds to the admittance of the circuit (including the source resistance R S ) with the normalised values for R S , C 1 , C 2 .
Remark 3.1: It is clear that a first-order state-space model and the associated transfer function can also be obtained by using an equivalent capacitance (C 1 + C 2 ) in Figure 1. However, simplifications of this form are not always straightforward in the development of models for networks with a large number of components. In such a case, the systematic methods presented in Section 2 may be of value to obtain a state-space model from the descriptor equation that arises from the use of Kirchoff 's laws together with an incidence matrix description of the network topology. This approach was employed, for instance, to model large 3-D RC networks as described in Galvao, Hadjiloucas et al. (2013).

Second example: RL circuit
This second example is concerned with the RL circuit presented in Figure 2. As in the previous case, the input signal u(t) corresponds to the source voltage and the output signal is taken as y(t) = i S (t). By applying Kirchoff 's current and voltage laws and noting that the voltage across the inductors is given by These equations can be cast into the descriptor form (1), (2) with and Again, the descriptor equation is already in the form assumed in Equation (3), witĥ However, in this case [Ê T 1Â T 2 ] T is singular for any values of R S , L 1 , L 2 , as can be seen from the final column of zeros in bothÊ 1 andÂ 2 . Therefore, the standard shuffle algorithm described in Section 2 cannot be employed to obtain a model in state-space form. A similar problem arises if one attempts to use the SVD-based approach. Indeed, let R S = L 1 = L 2 = 1 (normalised units) for illustration. The SVD operations then lead to As can be seenÃ 22 is not invertible, which prevents the use of this approach for converting the model to statespace form.
However, even though the methods described in Section 2 cannot be applied in this case, (48)-(52) can actually be used to obtain a state-space model. In fact, from (48)-(51 ), it follows that At this point, the problem lies in the absence of an explicit equation relating di L1 (t )/dt to di L2 (t )/dt. However, this problem can be circumvented by noting that (52) holds not only for a specific time t, but for every time. Therefore, it can be concluded that From (61) and (62), one arrives at which is a state equation of the form (12) with z(t) = i L1 (t). In view of (51), the output equation becomes which is in the form (13). Finally, by using (48)-(52) together with (63), the descriptor vector x(t) can be recovered from z(t) and u(t) as The method proposed in the next section generalises the reasoning employed in this example, thus providing an alternative to obtain a state-space model in cases where the procedures described in Section 2 cannot be applied.

Proposed method
Let be an ((N + m) × N) matrix defined as with E, B from the descriptor equation (1). Assume that is rank-deficient (i.e. the N columns of are linearly dependent) and let be an (N × q) matrix with columns forming a basis for the null space of , where q = N − rank( ). The matrix is such that E T = 0 N×q and B T = 0 m×q , which is equivalent to writing Pre-multiplying both sides of (1) by T leads to From (67) and (68), it follows that T Eẋ(t ) = 0 q×1 and T Bu(t ) = 0 q×1 , ∀t. Therefore, (69) can be rewritten as Now, the key point lies in the fact that (70) is valid not only for a specific time t, but for every t. Therefore, d/dt[ T Ax(t )] = 0 q×1 holds as well. Since , A are constant matrices, it follows that Equations (1) and (71) can then be merged into an extended descriptor equation of the form The use of (71) to obtain the new descriptor representation (72) resembles the procedure employed in the standard shuffle algorithm to arrive at (5). However, by using the matrix with the properties (67) and (68), one avoids the introduction of derivative inputu(t ) terms in the model. In what follows, the SVD procedure described in Section 2 will be applied to the new representation (72), with suitable modifiations to account for the fact thatĒ andĀ are no longer square matrices.
By using an SVD, the ((N + q) × N) matrixĒ can be expressed asĒ wheren is the number of singular values different from zero,p = N −n,¯ 1 is an (n ×n) diagonal matrix andŪ , V are unitary matrices with dimensions ((N + q) × (N + q)) and (N × N), respectively. Now, let or, alternatively, wherez(t ) ∈ n is a candidate state vector andw(t ) ∈ p is a vector of additional variables. After substituting (74) forĒ and (76) for x(t) in (72) and choosing appropriate partitions for the matrices˜Ā =Ū TĀV ,˜B =Ū TB , one may write Equations (78) and (79) are similar to (20) and (21) Moreover, from (2), (76) and (80), with an appropriate partition for matrices C andV , the output equation can be rewritten as which is of the form y(t ) = Hz(t ) + Lu(t ) with Finally, in view of (76) and (80), the descriptor vector x(t) can be calculated fromz(t ) and u(t) as Remark 4.1: Unlike the SVD procedure presented in Section 2, the proposed method to convert the original descriptor Equation (1) into the form (78), (79) is not a restricted equivalence transformation, because it involves an augmentation of the E, A, B matrices as indicated in (73).

The RL example revisited
The proposed method will now be illustrated by using the RL example presented in Section 3.2, with normalised component values R S = L 1 = L 2 = 1. In view of the E, B matrices presented in (54), (55), the matrix defined in (66) is given by In this case, N = 5 (number of descriptor variables) and rank( ) = 3. Therefore, the dimension of the null space of is q = 2. An (N × q) matrix that satisfies (67) and (68) TheĒ,Ā,B matrices defined in (73) are then given by The SVD procedure leads tō As can be seen,˜Ā 22 is a (4 × 2) matrix with linearly independent columns. By proceeding with the conversion to state-space form, the resulting model matrices are These matrices can be used to obtain the following transfer function (s): which corresponds to the admittance of the RL network in Figure 2 (including the source resistance R S ) with the normalised values for R S , L 1 , L 2 . It is worth noting that the cancellation of the s 2 factors in the numerator and denominator of the transfer function indicates the presence of two non-controllable and/or non-observable modes (Chen, 1984). By using a Kalman decomposition (Chen, 1984), these modes can be removed in order to reduce the order of the state-space model from three to one. As a result, it can be shown that the simplified model becomesż ( (which corresponds to the current through the inductors, up to a scale factor).

Remark 4.2:
It is clear that a first-order state-space model and the associated transfer function can also be obtained by using an equivalent inductance (L 1 + L 2 ) in Figure 2. However, such simplifications may not be straightforward in networks with a large number of components, as in the 3-D RLC case study which will now be presented.

Case study: descriptor modelling of three-dimensional RLC networks
The RLC network topology under consideration comprises N Z layers, as depicted in Figure 3, each one consisting of a matrix of nodes with N X rows and N Y columns. A previous work involving this 3-D RLC structure was reported in a conference (Paiva, Duarte, Galvão, & Hadjiloucas, 2013) within an investigation concerning the detection of changes in the composition of amorphous dielectric materials. However, that paper was not concerned with the conversion from descriptor to state-space form. Relevant applications associated with networks of this type may include the modelling of performance characteristics in oxide-based, metallic or Febased superconductors and high-temperature superconductors in wire, tape, ribbon, thin-film or thick-film form and micro-strip resonator designs. For example, one can see applications within the context of superconducting surface plasmon interfaces (Ghamsari & Majedi, 2011), and electromagnetic modelling of meta-materials, with application to the design of THz filters or studies of effective resistance in mesh structures (Zhongyang Li, & Ding, 2013) or other multilayer structures (Gutin, Ytterdal, Kachorovskii, Muraviev, & Shur, 2013). Such network models may also be of value to study the broadband response of photonic band-gap structures and subwavelength gratings (Halir et al., 2014). As in the illustrative examples presented in Section 3, the input signal u(t) in Figure 3 corresponds to the source voltage and the output signal y(t) is the current i S (t) entering the network. In what follows, the notation (t) for the time dependence will be omitted for brevity.
The 3-D RLC network topology can be expressed in terms of an incidence matrix W N N ×N B , as described in Galvao, Hadjiloucas et al. (2013) for the case of an RC network. A +1 entry in matrix W indicates that the branch at the corresponding column leaves the node at the corresponding row, a −1 entry indicates that the branch arrives at the node, and a 0 entry indicates that the branch is not connected to the node. The number of nodes N N includes the input electrode, but not the ground electrode, since the latter is used as the electric potential reference in Kirchoff 's laws. The number of branches N B includes the source branch plus the branches in all layers and the branches connecting consecutive layers. The number of resistors (not including the source resistor), inductors and capacitors is denoted by N R , N L and N C , respectively, so that N B = N R + N L + N C + 1. The fractions of resistors, inductors and capacitors in the network are defined as Kirchoff 's current and voltage laws (KCL and KVL) can be written as (Chua, Desoer, & Kuh, 1987) where i N B ×1 , v N B ×1 and e N N ×1 are the vectors of branch currents, branch voltages and node potentials, respectively. The i vector can be written as where i S is a scalar value corresponding to the source current and i R , i L and i C are the vectors of currents at the branches corresponding to the N R resistors, N L inductors and N C capacitors, respectively. P is a (N B × N B ) permutation matrix (Galvao, Hadjiloucas et al., 2013) that can be used to randomise the position of the R, L, C components in the network. It is worth noting that P is unitary, i.e. P −1 = P T . In a similar manner, the v vector can be written as where v S , v R , v L and v C denote the branch voltages at the source, resistors, inductors and capacitors, respectively. The relations between the branch currents and voltages can be expressed as Equations (96)- (103) can be combined into a descriptor representation of the network dynamics. In fact, from (96), (98) and (103), it follows that whereas from (97) and (99)- (102), it follows that which can be rewritten as or where M i S , M i R , M i L , M dv C /dt , M e correspond to the column blocks of the M (N N +N B )×(N N +N B ) matrix in (106) associated to i S , i R , i L , dv C /dt, e, respectively. By moving the i S , i R , i L , e terms to the right-hand side, and the di L /dt term to the left-hand side, (107) can be rewritten as T . An output equation of the form (2) for y = i S can be written with C = [1 0 1×N R 0 1×N L 0 1×N C 0 1×N N ] and D = 0. After conversion to the state-space form (12), (13), the transfer function (s) = Y (s)/U (s) = H(sI − F ) −1 G + L can be obtained and the admittance at a given frequency ω (rad/s) can thus be calculated as (jω). It is worth noting that (jω) includes the contribution of the source resistance R S . However, simple calculations for the series association of R S and the 3-D network can be used to obtain the admittance of the latter.

Results and discussion
The MATLAB R software was employed to implement the RLC network equations, obtain the corresponding statespace model, calculate (jω) and separate the contribution of the source resistance from the 3-D network admittance. The network dimensions were set to N X = N Y = 5, Section 2 were both singular, thus preventing the use of the standard shuffle algorithm and SVD-based approach for conversion to state-space form. In contrast, by using the method proposed in Section 4, the˜Ā T 22Ã 22 matrix was invertible in all the 100 cases, thus allowing the corresponding state-space models to be obtained. The computational effort involved in the conversion process was not substantial, requiring an average of 0.8 second for each network in a computer with i7-3687U processor (2.10GHz) and 8GB RAM.
For illustration, Figure 4(a) presents the admittance responses (amplitude and phase) obtained for two RLC networks, after removing the contribution of the source resistance R S . As can be seen, the admittance plots are slightly different for the two networks (solid and dotted lines), as a result of the random allocation of the R, L, C elements. However, overall both behave as band-reject filters with similar stop bands. It is worth noting that the admittance grows large at low and high frequencies due to the increase in the inductive and capacitive susceptances, respectively. As shown in the inset, the change from inductive to capacitive behaviour is reflected in a phase change from −90°to +90°as the frequency increases.
The contribution of the source resistance R S = 0.001 to the overall admittance can be visualised by comparing Figure 4(a) with the corresponding plot in Figure 4(b). In Figure 4(b), the admittance becomes upper-bounded by 1/R S = 10 +3 (60 dB), which establishes pass-bands with no ripple at low and high frequencies. The levelling of the amplitude plot within the pass bands is reflected in a change of the phase values towards zero, as can be seen in the inset.
It is interesting to note that, apart from the ripple in the stop band, the behaviour of these large 3-D networks resembles that of a simple network with single R, L, C components in a parallel configuration. Such an analogy can be further illustrated by investigating the effect of changes in the component values. In the case of a simple parallel RLC network, the admittance is given by which has a second-order numerator with natural frequency ω n = 1/ √ LC and damping factor ξ = (2R) −1 √ L/C. Therefore, an increase in L reduces the central frequency of the stop band and leads to larger damping. For illustration, Figure 5(a-c) present the admittance response of a 3-D RLC network with the same configuration used in Figure 4(a) and three different values of L, namely 0.002, 0.02 and 0.2. As can be seen, an increase in L reduces the central frequency of the stop band, as expected. In addition, the ripple within the stop band is also reduced. Although such an effect does not have a direct counterpart in the simple parallel RLC network, it may be argued that a reduction in the ripple can be associated to a larger overall damping of the network response. A similar analysis can be carried out with respect to changes in the capacitance C, as shown in Figure 5(d-f). An increase in C shifts the stop band to lower frequencies and increases the ripple, which is consistent with a reduction in both the natural frequency and overall damping of the network response. Finally, Figure 5(g-i) illustrates the result of changing the resistance R. The central frequency of the stop band remains unaltered, but the ripple increases with growing R. These effects are consistent with the expected results of keeping the natural frequency constant and reducing damping, which is again in line with the expressions for the simple parallel RLC network.

Conclusion
This paper has shown that standard methods for conversion of descriptor representations to state-space form are not always applicable to electrical circuit models derived in a systematic manner from Kirchoff 's laws and constitutive element equations. In particular, two well-known methods (shuffle algorithm and SVD) failed to provide a state-space model from descriptor representations of a simple RL circuit and complex 3-D RLC networks. A novel procedure to address this problem was proposed as an extension of the classic shuffle algorithm combined with the SVD approach. In a case study involving 100 different 3-D RLC networks with random allocation of components, the proposed method always succeeded in converting the descriptor representation to state-space form, which illustrates its advantage over the classic approaches.
Although the examples presented herein were concerned with electrical circuits, the relevance of the proposed method for applications involving different types of systems can be inferred from the existing analogies between circuit models and other physical domains (Hinaje, Raël, Noiying, Nguyen, & Davat, 2012;Schwarz, 2000;Smith, 2002). The case study investigated in the present paper also bridges a gap with respect to previous work concerning the modelling of 3-D RC networks (Galvao, Hadjiloucas et al., 2013;Galvão et al., 2013;Jacyntho et al., 2016) and opens up the possibility of investigating emergent properties of 3-D RLC networks in greater detail. Such investigations may be of value for several applications, including the study of percolation processes (Stauffer & Aharony, 1992), neuroscience (Hill, Wang, Riachi, Schürmann, & Markram, 2012) and critical phenomena (Hatef, Sadeghi, Fortin-Deschênes, Boulais, & Meunier, 2013). Simulations using 3-D RLC networks may also be used to study a range of complex or emergent behaviours encountered across all parts of the electromagnetic spectrum and are of much relevance to the optoelectronics, terahertz, photonics and plasmonics communities where complex dynamics across a wide range of frequencies are observed (Avouris & Freitag, 2014). Such simulations can be particularly valuable for the interpretation of results at grain boundaries, e.g. charge transport in graphene-based electronics (Tsen et al., 2012) and interconnects (Nishad & Sharma, 2014) in experiments in compounds displaying complex (e.g. stretched multi-exponential mixed with resonant) responses (Scheller, 2011).
Moreover, studies concerning the use of fractional order models to approximate the dynamic behaviour of 3-D RLC networks may be of particular interest to complement the studies involving the 3-D RC case (Galvao, Hadjiloucas et al., 2013;Jacyntho et al., 2016). This is a topic of significant interest within the systems and control literature (Elwakil, 2010;Mesbahi, Haeri, Nazari, & Butcher, 2015;Radwan, 2013), which we shall pursue in future work.

Disclosure statement
No conflict of interest is reported by the authors.