Mathematical Problems in Engineering

Volume 2017, Article ID 6739857, 12 pages

https://doi.org/10.1155/2017/6739857

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

^{1}Department of Electronics and Telecommunications, Politecnico di Torino, Corso Duca degli Abruzzi 24, 10129 Torino, Italy^{2}Department of Information Technology, IDLab, Ghent University-IMEC, iGent Tower, Technologiepark-Zwijnaarde 15, 9052 Ghent, Belgium

Correspondence should be addressed to Riccardo Trinchero; ti.otilop@orehcnirt.odraccir

Received 3 November 2016; Accepted 8 March 2017; Published 30 March 2017

Academic Editor: J.-C. Cortés

Copyright © 2017 Riccardo Trinchero et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

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