TMsim : An Algorithmic Tool for the Parametric and Worst-Case Simulation of Systems with Uncertainties

This paper presents a general purpose, algebraic tool—named TMsim—for the combined parametric and worst-case analysis of systems with bounded uncertain parameters.The tool is based on the theory of Taylormodels and represents uncertain variables on a bounded domain in terms of a Taylor polynomial plus an interval remainder accounting for truncation and round-off errors.This representation is propagated from inputs to outputs by means of a suitable redefinition of the involved calculations, in both scalar and matrix form. The polynomial provides a parametric approximation of the variable, while the remainder gives a conservative bound of the associated error. The combination between the bound of the polynomial and the interval remainder provides an estimation of the overall (worst-case) bound of the variable. After a preliminary theoretical background, the tool (freely available online) is introduced step by step along with the necessary theoretical notions. As a validation, it is applied to illustrative examples as well as to real-life problems of relevance in electrical engineering applications, specifically a quarter-car model and a continuoustime linear equalizer.


Introduction
With the ever-growing complexity of modern design scenarios in every domain of engineering, today's applications are constantly facing the need for techniques that are capable of taking parameter uncertainty into account. Researchers in the field of uncertainty quantification have striven to find alternatives to state-of-the-art approaches, which often turn out to be computationally intensive or inaccurate. In this framework, the strategies can be classified into two groups: probabilistic approaches take the actual distribution of uncertain parameters into account and seek for a statistical characterization of the system variables of interest. A classical example that has been utilized in virtually every engineering problem is the Monte Carlo (MC) method [1]. An alternative and more efficient probabilistic approach that was proposed in relatively recent times is polynomial chaos [2,3].
On the contrary, interval or worst-case (WC) approaches assume bounded uncertain input parameters and provide a conservative estimation of the lower and upper bounds of the system variables without explicit information on the actual distribution. The calculations involved in the analysis are suitably modified to provide a conservative estimation of the bound of the result. These methods include interval arithmetic or interval analysis (IA) [4], which is however rather basic and often results in large overestimation since it is unable to track variable dependency [5]. Affine arithmetic (AA) is an improved approach that takes variable dependence into account by means of a linear parametrization [6]. Finally, the Taylor model (TM) utilizes a higher-order polynomial representation in conjunction with a residual interval remainder that takes truncation and round-off errors into account, yielding a well-defined and robust theoretical framework [7][8][9]. TM calculations are carried out by using standard operations over polynomials or, whenever a nonlinear operator occurs, by the rules of Taylor expansion. The interval remainder is propagated and updated by means of Taylor remainder theory and IA. The TM algebra combines 2 Mathematical Problems in Engineering the strength of a high-order parametric representation and of IA, while limiting the overestimation issue to the interval remainder. By adding the interval remainder to the bound of the polynomial part, a rigorous and conservative estimation of the overall (WC) bounds of the variable is obtained.
This paper presents an open source and general purpose algebraic tool, named TMsim, that implements TM algebra and proposes specific improvements in the calculation of remainders and an extension to some matrix operations. Specifically, the remainder of nonlinear operations, which is usually computed in Lagrange form [7], is now calculated also according to the Cauchy formulation. The tightness and convergence differ between the Lagrange and Cauchy forms; thus the best of the two is chosen. Moreover, for the special case of the multiplicative inverse, an alternative and exact formulation of the remainder is implemented. The tool also handles basic matrix operations and implements a recently proposed algorithm to implement in TM algebra the matrix inversion [10], which is a critical operation in the solution of many engineering problems.
A preliminary framework was introduced with specific application to the frequency-domain analysis of circuits [10] and transmission lines [11]. In the present paper, the tool is further extended and presented in a more systematic and comprehensive way. It is implemented in MATLAB and made available online [12].

Problem Statement
Given a set of independent input variables = [ 1 , . . ., ] ∈ , with associated uncertainty bounds = [ 1 , 1 ] × ⋅ ⋅ ⋅ × [ , ] ⊂ R , the purpose of the proposed TMsim tool is to algebrically and parametrically calculate a set of output variables and their corresponding uncertainty bounds. Each variable that depends on the inputs is understood to be represented according to the TM formulation [7] where 0 is the center of the domain , is a th-order multivariate polynomial, and = [ , ] is an interval remainder such that The TM of ( ) is therefore guaranteed to lie between the upper and lower bound defined by the value of the polynomial and the corresponding remainder . Ideally, if = [0, 0], the polynomial provides an exact parametric representation.
Moreover, the information on the overall (WC) bounds of ( ) is readily obtained as ( ) + , where (⋅) denotes the bound operator; that is, and the addition between intervals is intended in the IA-sense [4]; that is, The determination of the bound of a multivariate polynomial is a nontrivial task that is handled as detailed in Section 5. In practice, a common maximum order is assumed for each TM variable. Hence, for the ease of notation, the superscript is dropped from now on when denoting the polynomial part of a TM.

Basic Scalar Operations
With the previously introduced definitions, basic scalar operations like addition and subtraction between TMs are readily computed as where the operation over the remainders is again performed according to the IA, with the subtraction being defined as The multiplication between two TMs yields where ⋅ is the part of the product ( − 0 ) ⋅ ( − 0 ) up to order , whereas is the remaining higher-order contribution. The higher-order part and the remaining terms are encompassed into an interval remainder calculated as where the product of two intervals is defined in IA as Complex calculations are handled by suitably operating on the real and imaginary part with real-valued computations [10].

TMsim Overview
TMsim is a MATLAB tool that implements TM algebra and defines a specific class of variables for which the standard basic operations (see Section 3) as well as other available nonlinear functions or matrix operations (see Sections 6 and 7, respectively) are replaced by the corresponding TM calculations.
The tool requires an initialization that defines the order of the TM representation and the number of independent variables and creates some global auxiliary variables. An independent input variable is defined by its upper and lower bounds and the index ∈ {1, . . . , }. Each variable is represented as a structure with the field triplet ( , , ): With the above definitions, the basic scalar operations are readily implemented as described in Section 3. These also include the raise to an integer power, which is implemented by means of iterative multiplications.
As a brief example, the code below defines two variables It is worth noting that the third-order representation is sufficient to represent 1 and 2 without error, but not whose remainder is non-null (

Polynomial Bounds
The exact calculation of the bounds of a multivariate polynomials can be performed exactly only for order ≤ 2. No rigorous approach exists for higher-order polynomials. TMsim implements a strategy, based on Bernstein polynomials, that allows for a conservative estimation of the bounds [7,10]. Bernstein polynomials have the notable property of being bounded by their largest and smallest coefficient [13][14][15][16]. By means of a change of basis, a standard univariate where the Bernstein basis polynomials are defined as The corresponding Bernstein coefficients are from which The above expression emphasizes that the obtained bounds are not necessarily exact, but rather a conservative estimation of the actual bounds. For the multivariate case, consider a polynomial in the form wherẽ∈ [−1, 1] ∀ and k = [ 1 , . . . , ] is a multi-index defining the exponents of each monomial, with 0 ≤ ≤ and ∑ =1 ≤ . The vector of multivariate Bernstein coefficients is calculated as where is a vector with entries In order to conveniently calculate the Bernstein coefficients through (16), the independent variables are rescaled into the domain = [−1, 1] . It is important to remark that, for a given order , the vectors 0 up to for the conversion in (16) are predetermined and stored into a look-up table.

Nonlinear Functions
The TM resulting from the application of a nonlinear function (⋅) is calculated by considering the functional Taylor expansion: centered at = (0). A shifted version = − of the original TM has been introduced. Moreover, the function derivatives are defined as and ( ( )) is the pertinent remainder in Lagrange or Cauchy formulation. It should be noted that the first term ( ) + in the r.h.s. of (18) merely involves sums and products of TMs, which are computed with the rules outlined in Section 3. The remainder , resulting from these calculations, is added to the bound of the Lagrange/Cauchy remainder to obtain the overall remainder .
The Lagrange remainder of the th-order Taylor expansion centered in 0 of a function ( ) ∈ +1 on ∈ is defined as [17][18][19] , ( ( )) fl where * ∈ [ 0 , ] = 0 + , with = − 0 and ∈ [0, 1]. An alternative formulation of the remainder is the Cauchy form, which reads [19] , ( ( )) fl It is important to remark that both the Lagrange and the Cauchy remainder provide a conservative estimation of the truncation error. However, as described in the following, they have in general different bounds, and one of the two may even diverge as the order is increased. Therefore, a single definition of the remainder does not guarantee a tight estimation of the error. For this reason, the TMsim tool calculates the remainder in both Lagrange and Cauchy form and selects the one with the smaller bound.
The remaining part of this section details the implementation of a representative set of general purpose nonlinear functions, including the logarithm, the square root, the multiplicative inverse, and the sin and cosine functions. The implementation of the exponential function is well documented in the literature [7][8][9]. All the aforementioned functions are available in the current version of the tool. Other nonlinear functions can be implemented by following a similar procedure.
6.1. Logarithm. Consider the logarithmic function ( ( )) = log ( ( )). According to [7], the resulting TM ( ) is computed as where the remainder is given in the Lagrange formulation (20) as The bound of the remainder is obtained by means of the IA as follows: which is equivalent to Unfortunately, the Lagrange remainder does not always converge to zero when the expansion order is increased [19]. In particular, from (24), which implies that the Lagrange remainder converges only if However, since ( ) represents the variation of a parameter around its means value, its codomain can include negative values. The Cauchy remainder (21) allows extending the convergence region, thus alleviating the limitations of the Lagrange remainder. Starting with the simple equality log ( ( )) = log ( ) + log (1 +̃( )) , Mathematical Problems in Engineering 5 wherẽ( ) = ( )/ , the Cauchy remainder is expressed as Since 0 < < 1 and under the (quite practical) assumption |̃| < 1 (i.e., the variation of a variable is smaller than its central value) [19], which leads to From the above expression, which implies that the convergence region is −1 ≤ ( )/ ≤ 1 as opposed to (27). As in the case of the Lagrange remainder, the bound of the Cauchy remainder is calculated using IA as where the remainder in Lagrange form is The bound of the above remainder is computed by mean of IA as ( , ( ( ))) = ( ( )) +1 / +1 (1 + ( ( )) / ) +1/2 .
However, also in this case The use of the Cauchy remainder, which reads [19] , ( ( )) =̃+ where in principle the remainder can be computed via either the Lagrange or Cauchy formulation. However, in this special case, an exact remainder is calculated based on the properties of geometric series [17].
To begin with, the Taylor expansion of 1/(1 − ) up to the order is considered The exact remainder for | | < 1 is obtained via the identity The above equation is recast to obtain the remainder of 1/ ( ) as where in this casẽ= − / . The bound of the above remainder is computed as for |̃| < 1.
respectively. The bounds of the Lagrange remainder are calculated as By substituting (46) or (47) in (48), the determination of the bounds requires the calculation of the bounds (sin( )) or (cos( )) in the pertinent domain of ( ). The estimation of these bounds is complicated by the periodic behavior of the sine and cosine functions. Since the sine and cosine are bounded functions, a naive solution would be to simply assume (sin( )) = (cos( )) = [−1, 1]. However, tighter bounds can be obtained by considering the function as locally monotonic and by including the effect of a possible slope change only when a stationary point falls into the domain of ∈ ( ( )). This leads to [11] min with ∈ Z.
6.5. Illustrative Example. As an illustrative example, the function ( ) = sin(1/ log(√ )) with ∈ [3,9] is considered. The corresponding TM is computed by means of the MATLAB code >> g=sin(inv(log(sqrt(x)))); The result is shown in Figure 1. The black curve is the actual value of the function ( ). The dashed red curve is the approximation provided by the polynomial part of a ninthorder TM. Finally, the green and magenta lines are the upper and lower bounds obtained with the inclusion of the interval remainder for an expansion order of = 6 and = 9, respectively. It is appreciated that increasing the order yields tighter bounds. However, both bounds are conservative over the entire domain of the input variable .

Matrix Operations
A matrix TM is a matrix where all the entries are defined via the standard TM representation (1) and is here denoted by means of bold quantities as T ( ) = P ( ) + I . The basic scalar calculations introduced in Section 3 are readily generalized to the corresponding matrix computations by operating element-wise. For instance, the product of two matrix TMs merely involves additions and multiplications of scalar TMs.
On the other hand, despite a generalization of the Taylor expansion to matrices being defined, the implementation of nonlinear matrix operations in the TM algebra is hindered by the fact that no corresponding definition of the remainder is available. A rigorous solution to handle the matrix inversion was proposed in [10] based on the Sherman-Morrison (SM) formula [20]. The implementation of other nonlinear matrix operations, such as the exponentiation, is still under investigation. An envisaged approach is to seek for an eigenvalue decomposition in TM algebra in order to reduce the problem to a scalar nonlinear operation over the eigenvalues.

Matrix Inversion.
The SM formula provides an exact solution for the computation of the inversion (A + B) −1 when B has unitary rank. Under this assumption, the resulting inverse matrix is where = trace (A −1 B) and only the inversion of matrix A is required.
For its application to a matrix TM T , this is first split into a constant matrix and a shifted TM, as was done in Section 6 for the scalar case: where c = P (0) is a standard matrix and P ( − 0 ) = P ( − 0 ) − c . If the nonconstant part T has unitary rank, then (52) is applied directly and the inverse TM is Hence, the SM formula reduces the inversion of a matrix TM to standard and available calculations, namely, (1) sums and multiplications of matrix TMs; (2) the inversion of the constant matrix c , which is carried out according to classical algebra; (3) the calculation of the trace operator, which merely involves the scalar sum of the TMs on the matrix diagonal; (4) the calculation of the multiplicative inverse 1/(1 + ), which is carried out as described in Section 6.3.
It is worth noting that in practice the matrix c is always invertible as it corresponds to the matrix for the central value of the variable . The generalization of (52) to the case of full-rank matrices requires an iterative application of the algorithm [10]. The matrix T is first split into the sum of rank-one matrices: where T , is defined as a matrix TM with null entries except for the th column, which coincides with the th column of the original TM T , thus being rank-1 by construction. The following iterative procedure is then put forward: (1) the inversion (c + T ,1 ) −1 is computed via (50) by considering A = c and B = T ,1 ; (2) the inversion (c +T ,1 +T ,2 ) −1 is calculated by letting A = c +T ,1 and B = T ,2 . The inverse of A, as needed in (50), is available from the previous step.
The above procedure is iterated until all terms in (53) are accounted for, thus yielding the end result.

Illustrative Example.
The outlined procedure is general and allows inverting an arbitrary matrix whose elements are nonlinear functions. As an example, consider the following matrix: where ∈ = [−1, 1].
The above matrix belongs to the class of Toeplitz matrices, which are rather common in many fields of engineering as they result from Fourier expansions [21][22][23]. It is important to remark that det ] .
(56) Furthermore, Figure 2 compares the behavior of [M −1 ( )] 11 (black line) with the polynomial model provided by a TM of order = 6 (dashed red line). The bounds for expansion orders of = 3 and = 6 are indicated by the green and magenta curves, respectively. This example highlights the capability of TMsim to provide a parametric approximation of the entries of the inverse of a matrix with nonlinear elements as well as a conservative estimation of their upper and lower bounds.
It is interesting to calculate the product between T and its inverse. In classical algebra, this would result in the identity matrix. With the proposed TM framework and an expansion order = 3, the result is It is remarkable that the polynomial part still defines an identity matrix. This property stems from the exactness of the SM formula. All the approximations arising from the TM operations are included in the interval remainder. By increasing the expansion order, the interval remainder is reduced. For instance, using = 6 yields  Figure 3 which is widely used for the analysis of the dynamical behavior of vehicle suspensions under the assumption that the front and the rear part of the vehicle are uncoupled. The model has two degrees of freedom, specifically, the vertical displacements of the chassis and of the tire from the reference frame. The road profile ℎ is considered as the excitation of the system. The proposed quarter-car model is defined by the following parameters, the masses of the chassis and of the tire , the spring coefficients of the suspension and of the tire , and the damping coefficients of the suspension and of the tire . According to [24], the quarter-car model is described by the following matrix equation:  where and are the displacements from the static equilibrium position with respect to an inertial frame.
After some straightforward manipulation, the above equation is rewritten in its state-space controllable form, which readsẋ where the column vector x( ) = [V ( ), V ( ), ( ), ( )] ∈ R 4 collects the state variables of the system and ( ) = [ℎ( ),ḣ ( )] ∈ R 2 is a vector defining the sources of the system. The matrix A ∈ R 4×4 and the column vector B ∈ R 4×2 are defined as where I and 0 denote the identity matrix and the null matrix of dimensions 2 × 2, respectively. The remaining terms are defined as follows: The frequency-domain response of the system, denoted as X( ), is obtained via the Fourier transform of (3) and by considering as excitation vector ( ) = 1/(2 ) exp( )[1, ] : where U( ) = [1, ] . For the simulation, it is assumed = 332.5 kg, = 38 kg, and = 0 Ns/m. In addition, three independent uncertainty variables are considered; that is, = (6000 ± 20%) Ns/m, = (22000 ± 20%) N/m, and = (20000 ± 20%) N/m. Figure 4 shows the variation (gray area) of the real (top panel) and imaginary (bottom panel) part of the transfer function ( ) = ( )/ ( ), related to the state variable ( ), computed with 10'000 MC samples over 100 frequency points. The WC bounds obtained with a fifth-order analysis in TMsim (magenta lines) compare very well with the spread given by the MC samples. Figure 5 provides an analogous comparison for the magnitude of the transfer function. The excellent agreement validates the capability of TMsim to accurately deal with large parameter variations in solving a matrix system of equations.

Continuous-Time Linear
Equalizer. The second proposed application example deals with the analysis of three CTLEs whose electrical circuits are depicted in Figure 6. CTLEs are commonly used in telecommunications to compensate for high-frequency channel losses with a "high-pass" type of response. The impact of parameter uncertainty in these CTLE configurations was analyzed in [25]. The circuit components are as follows: (a) = 35 Ω ± 20%, = 1 pF ± 20%, = 3 nH ± 20%; (b) = 100 Ω ± 20%, = 1 pF ± 20%, = 3 nH ± 20%; (c) = 100 Ω ± 20%, = 1 pF ± 20%.  The impact of the uncertainty in the circuit components on the transmission from the input port (left side) to the output port (right side), commonly denoted in electrical engineering as " 21 ", is investigated by analyzing the three circuits of Figure 6 by means of the modified nodal analysis (MNA) [10,26]. The results for the real part, imaginary part, and magnitude of 21 for the three CTLE configurations are collected in the first, second, and third row of Figure 7, respectively. The gray areas indicate the spread of the response resulting from a MC analysis with 10'000 samples, whereas the magenta lines are the WC bounds computed from a fourth-order TM. This second example confirms the excellent accuracy provided by the TM analysis despite the large variability of the circuit components.

Conclusions
This paper presents an open source and freely available algebraic tool, named TMsim, that allows for the parametric and WC analysis of systems with bounded uncertain parameters. TMsim implements TM algebra, which represents uncertain variables as Taylor expansions complemented by an interval remainder accounting for approximation errors. While the Taylor expansion provides a parametric model of the variable, the information on the remainder allows obtaining a conservative estimation of the WC bounds. The TM representation is propagated from input to output uncertainties via a suitable redefinition of the linear and nonlinear operations involved in the system analysis. TMsim implements special improvements in the estimation of remainders and extensions to matrix operations. Specifically, a more accurate calculation of the remainders is obtained by choosing the tightest option between the Lagrange and Cauchy form. Moreover, an exact determination of the remainder is accomplished for the relevant case of the multiplicative inverse function. The tool also implements an effective algorithm for the inversion of a matrix TM.
The tool is validated by means of ad hoc illustrative examples as well as real-life application examples, that is, a quarter-car model and three CLTE electrical circuits. Excellent accuracy in the conservative estimation of the WC bounds of the output variables is achieved.

Conflicts of Interest
The authors declare that they have no conflicts of interest.