Abstract
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 Taylor models 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 continuous-time linear equalizer.
1. 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–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 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].
2. Problem Statement
Given a set of independent input variables , with associated uncertainty bounds , 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 is the center of the domain , is a th-order multivariate polynomial, and is an interval remainder such thatThe 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 , 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.
3. Basic Scalar Operations
With the previously introduced definitions, basic scalar operations like addition and subtraction between TMs are readily computed aswhere the operation over the remainders is again performed according to the IA, with the subtraction being defined as
The multiplication between two TMs yieldswhere is the part of the product 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 aswhere 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].
4. 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 . Each variable is represented as a structure with the field triplet :(i) is an array that collects the coefficients of the multivariate polynomial;(ii) is an auxiliary two-element vector with the lower and upper bound of the polynomial;(iii) is a two-element vector containing the lower and upper bound of the interval remainder.To ease the calculation of the polynomial bounds (see Section 5), the coefficients of the polynomial are given in terms of rescaled variables , , instead of the original variables .
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 and and calculates , , and . >> d = 2; % number of variables d_i=1,2 >> d_1 = 1; % index of variable x1 >> d_2 = 2; % index of variable x2 >> n = 3; % expansion order >> >> % Taylor model representation >> initialize_TMsim(n,d); >> x1=TaylorModel(2,3,d_1,’bounds’); >> x2=TaylorModel(-1,1,d_2,’bounds’); >> y1=x12+x2; >> y2=x1*x22; >> z=x1*y2+x2*y1;The results areIt is worth noting that the third-order representation is sufficient to represent and without error, but not whose remainder is non-null (this would have been the case by setting ). The estimated bound on is , which is a conservative estimation of the true bound .
5. Polynomial Bounds
The exact calculation of the bounds of a multivariate polynomials can be performed exactly only for order . 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–16].
By means of a change of basis, a standard univariate polynomial in the variable is converted into the Bernstein formwhere the Bernstein basis polynomials are defined asThe corresponding Bernstein coefficients arefrom whichThe 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 formwhere and is a multi-index defining the exponents of each monomial, with and . The vector of multivariate Bernstein coefficients is calculated aswhere is a vector with entries
In order to conveniently calculate the Bernstein coefficients through (16), the independent variables are rescaled into the domain . It is important to remark that, for a given order , the vectors up to for the conversion in (16) are predetermined and stored into a look-up table.
6. Nonlinear Functions
The TM resulting from the application of a nonlinear function is calculated by considering the functional Taylor expansion:centered at . A shifted version of the original TM has been introduced. Moreover, the function derivatives are defined asand 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 of a function on is defined as [17–19]where , with and .
An alternative formulation of the remainder is the Cauchy form, which reads [19]
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–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 . According to [7], the resulting TM is computed aswhere 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 ifHowever, 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 equalitywhere , the Cauchy remainder is expressed asSince and under the (quite practical) assumption (i.e., the variation of a variable is smaller than its central value) [19],which leads toFrom the above expression,which implies that the convergence region is as opposed to (27).
As in the case of the Lagrange remainder, the bound of the Cauchy remainder is calculated using IA as
6.2. Square Root
Consider now the square root function . The resulting TM is [7]where the remainder in Lagrange form is
The bound of the above remainder is computed by mean of IA as
However, also in this case
The use of the Cauchy remainder, which reads [19]whereallows again mitigating the convergence issue. Under the assumptions that lead to (30), the corresponding bound is obtained as
6.3. Multiplicative Inverse
Consider the function . The resulting TM is [7]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 up to the order is considered
The exact remainder for is obtained via the identity
The above equation is recast to obtain the remainder of aswhere in this case . The bound of the above remainder is computed asfor .
6.4. Sine and Cosine
For the case of sine and cosine functions, the coefficients in (18) are given by [7]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 or 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 . 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]with .
6.5. Illustrative Example
As an illustrative example, the function with is considered. The corresponding TM is computed by means of the MATLAB code >> initialize_TMsim(9,1); >> x=TaylorModel(3,9,1,’bounds’); >> 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 ninth-order 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 and , respectively. It is appreciated that increasing the order yields tighter bounds. However, both bounds are conservative over the entire domain of the input variable .
7. 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 . 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.
7.1. Matrix Inversion
The SM formula provides an exact solution for the computation of the inversion when has unitary rank. Under this assumption, the resulting inverse matrix iswhere and only the inversion of matrix is required.
For its application to a matrix TM , this is first split into a constant matrix and a shifted TM, as was done in Section 6 for the scalar case:where is a standard matrix and . If the nonconstant part has unitary rank, then (52) is applied directly and the inverse TM iswith . 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 , 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 , which is carried out as described in Section 6.3.It is worth noting that in practice the matrix 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 is first split into the sum of rank-one matrices:where is defined as a matrix TM with null entries except for the th column, which coincides with the th column of the original TM , thus being rank-1 by construction. The following iterative procedure is then put forward:(1)the inversion is computed via (50) by considering and ;(2)the inversion is calculated by letting and . The inverse of , 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.
7.2. 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 .
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–23]. It is important to remark that .
The following MATLAB script calculates the inverse of by means of TMsim:
>> initialize_TMsim(3,1);
>> x_TM = TaylorModel(-1,1,1,’bounds’)
>> a_TM = exp(1+0.3*x_TM);
>> b_TM = log(3+x_TM);
>> c_TM = sin(pi+pi/4*x_TM);
>> M_TM = [a_TM,b_TM,c_TM;
b_TM,a_TM,b_TM;
c_TM,b_TM,a_TM];
>> invM_TM = inv(M_TM);
The resulting matrix TM provides a parametric representation of the nonlinear matrix over the entire domain and a conservative estimation of its upper and lower bounds. The minimum and maximum values of the entries of , estimated based on 10’000 MC samples, arerespectively.
These bounds are also estimated from the TM by considering that . This, for an expansion order , leads to
Furthermore, Figure 2 compares the behavior of (black line) with the polynomial model provided by a TM of order (dashed red line). The bounds for expansion orders of and 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 and its inverse. In classical algebra, this would result in the identity matrix. With the proposed TM framework and an expansion order , the result isIt 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 yields
8. Application Example
In this section, the accuracy and the strength of TMsim are further investigated by considering two real-life application examples, namely, a quarter-car model and a continuous-time linear equalizer (CTLE). The code for these and other application examples is available online [12].
8.1. Quarter-Car Model
The first example considers the quarter-car model depicted in 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 readswhere the column vector collects the state variables of the system and is a vector defining the sources of the system. The matrix and the column vector are defined aswhere and denote the identity matrix and the null matrix of dimensions , respectively. The remaining terms are defined as follows:
The frequency-domain response of the system, denoted as , is obtained via the Fourier transform of (3) and by considering as excitation vector :where .
For the simulation, it is assumed kg, kg, and Ns/m. In addition, three independent uncertainty variables are considered; that is, Ns/m, N/m, and 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.
(a)
(b)
8.2. 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), , ;(b), , ;(c), .
(a)
(b)
(c)
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 “”, 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 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.
(a)
(b)
(c)
9. 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.
Acknowledgments
The authors would like to thank Luca Castellazzi, Politecnico di Torino, for providing the mechanical example on the quarter-car model. This work was supported in part by the Research Foundation Flanders (FWO), of which Dr. Manfredi is a Postdoctoral Research Fellow.