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

Authors: Riccardo Trinchero; Paolo Manfredi; Igor S. Stievano
Journal: Mathematical Problems in Engineering
 (2017)

---

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


---

## Body

## 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 andexact 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 ofd independent input variables x=[x1,…,xd]∈D, with associated uncertainty bounds D=[a1,b1]×⋯×[ad,bd]⊂Rd, 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 f that depends on the inputs is understood to be represented according to the TM formulation [7](1)fx∈Tf≔Pfnx-x0+If,where x0 is the center of the domain D, Pfn is a nth-order multivariate polynomial, and If=[a,b] is an interval remainder such that(2)Pfnx-x0+a≤fx≤Pfnx-x0+b∀x∈D.The TM of f(x) is therefore guaranteed to lie between the upper and lower bound defined by the value of the polynomial Pfn and the corresponding remainder If. Ideally, if If=[0,0], the polynomial Pfn provides an exact parametric representation.Moreover, the information on the overall (WC) bounds off(x) is readily obtained as B(Pf)+If, where B(·) denotes the bound operator; that is,(3)Bfx=minx⁡fx,maxx⁡fxand the addition between intervals is intended in the IA-sense [4]; that is,(4)a,b+c,d=a+c,b+d.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 ordern is assumed for each TM variable. Hence, for the ease of notation, the superscript n 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 as(5)Tf±Tg≔Pfx-x0±Pgx-x0+If±Ig,where the operation over the remainders is again performed according to the IA, with the subtraction being defined as(6)a,b-c,d≔a-d,b-c.The multiplication between two TMs yields(7)Tf·Tg≔Pf·gx-x0+Pex-x0+Pfx-x0Ig+Pgx-x0If+IfIg,where Pf·g is the part of the product Pf(x-x0)·Pg(x-x0) up to order n, whereas Pe is the remaining higher-order contribution. The higher-order part and the remaining terms are encompassed into an interval remainder calculated as(8)If·g=BPe+BPfIg+BPgIf+IfIg,where the product of two intervals is defined in IA as(9)a,b·c,d≔min⁡ac,ad,bc,bd,max⁡ac,ad,bc,bd.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 Section3) 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 ordern of the TM representation and the number d of independent variables and creates some global auxiliary variables. An independent input variable is defined by its upper and lower bounds and the index i∈{1,…,d}. Each variable is represented as a structure with the field triplet (P,B,I):(i)
P is an array that collects the coefficients of the multivariate polynomial;(ii)
B is an auxiliary two-element vector with the lower and upper bound of the polynomial;(iii)
I 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 x~i=2xi-ai-bi/(bi-ai)∈[-1,1], ∀i=1,…,d, instead of the original variables xi∈[ai,bi].With the above definitions, the basic scalar operations are readily implemented as described in Section3. 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 variablesx1=[2,3] and x2=[-1,1] and calculates y1=x12+x2, y2=x1x22, and z=x2·y1+x1·y2.>> 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=x1∧2+x2;>> y2=x1*x2∧2;>> z=x1*y2+x2*y1;The results are(10)y1=P=6.25+2.5x~1+x~2+0.25x~12,B=3,10,I=0,0y2=P=2.5x~22+0.5x~1x~22,B=-1,3,I=0,0z=P=6.25x~2+2.5x~1x~2+7.25x~22+0.25x~12x~2+2.5x~1x~22,B=-6.25,18.75,I=-0.05,0.25.It is worth noting that the third-order representation is sufficient to represent y1 and y2 without error, but not z whose remainder is non-null (this would have been the case by setting n=4). The estimated bound on z is [-6.25,18.75]+[-0.05,0.25]=[-6.3,19], which is a conservative estimation of the true bound [-2.025,19].


## 5. Polynomial Bounds

The exact calculation of the bounds of a multivariate polynomials can be performed exactly only for ordern≤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–16].By means of a change of basis, a standard univariate polynomialP(x)=∑i=0npixi in the variable x∈[a,b] is converted into the Bernstein form(11)PBx=∑i=0nbiBix,where the Bernstein basis polynomials are defined as(12)Bix=nib-anx-aib-xn-i,0≤i≤n.The corresponding Bernstein coefficients are(13)bi=∑j=0iijnjb-aj∑k=jnkjak-jpk,from which(14)mini⁡bi≤minx⁡Px≤Px≤maxx⁡Px≤maxi⁡bi.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(15)Px=∑kpk∏i=1dx~iki,where x~i∈[-1,1]  ∀i and k=[k1,…,kd] is a multi-index defining the exponents of each monomial, with 0≤ki≤n and ∑i=1d‍ki≤n. The vector of multivariate Bernstein coefficients is calculated as(16)b=∑kpk·βk1⊗βk2⊗⋯⊗βkd,where βk is a vector with entries(17)βi,k=∑j=0min⁡i,kijnj2jkj-1k-j.In order to conveniently calculate the Bernstein coefficients through (16), the independent variables are rescaled into the domain D=-1,1d. It is important to remark that, for a given order n, the vectors β0 up to βn 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 functiong(·) is calculated by considering thefunctional Taylor expansion:(18)gTfx=Pgx+Ig=∑k=0ngkcfk!T¯fxk︸Pgx+Ip+RnTfxcentered at cf=Pf(0). A shifted version T¯f=Tf-cf of the original TM has been introduced. Moreover, the function derivatives are defined as(19)gkcf≔dgTfxdTfxx=cfand Rn(Tf(x)) is the pertinent remainder in Lagrange or Cauchy formulation. It should be noted that the first term Pg(x)+Ip 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 Ip, resulting from these calculations, is added to the bound of the Lagrange/Cauchy remainder Rn to obtain the overall remainder Ig.The Lagrange remainder of thenth-order Taylor expansion centered in x0 of a function g(x)∈Cn+1 on x∈D is defined as [17–19](20)RL,ngx≔gn+1x∗n+1!x-x0n+1≡gn+1x0+θx¯n+1!x¯n+1,where x∗∈[x0,x]=x0+θx¯, with x¯=x-x0 and θ∈[0,1].An alternative formulation of the remainder is the Cauchy form, which reads [19](21)RC,ngx≔x-x∗nn!x-x0gn+1x∗≡1-θnn!x¯n+1gn+1x0+θx¯.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 ordern 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 functiong(Tf(x))=log⁡(Tf(x)). According to [7], the resulting TM  Tg(x) is computed as(22)log⁡Tfx=log⁡cf+∑k=1n-1k+11kT¯fkxcfk+RnTfx,where the remainder is given in the Lagrange formulation (20) as(23)RL,nTfx=-1n+2n+1T¯fn+1x/cfn+11+θT¯fx/cfn+1.The bound of the remainder is obtained by means of the IA as follows:(24)BRL,nTfx=-1n+2n+1BT¯fxn+1/cfn+1B1+T¯fx/cfn+1,which is equivalent to(25)BRL,nTfx=-1n+2n+1BT¯fn+1BTfn+1.Unfortunately, the Lagrange remainder does not always converge to zero when the expansion ordern is increased [19]. In particular, from (24),(26)RL,nTfx⟶n→∞0⟺BT¯f/cf1+T¯fx/cf<1which implies that the Lagrange remainder converges only if(27)0≤T¯fxcf≤1.However, since T¯f(x) 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(28)log⁡Tfx=log⁡cf+log⁡1+T~fx,where T~f(x)=T¯f(x)/cf, the Cauchy remainder is expressed as(29)RC,nTfx=-1n+11-θ1+θT~fxnT~fn+1x1+θT~fx.Since 0<θ<1 and under the (quite practical) assumption T~f<1 (i.e., the variation of a variable is smaller than its central value) [19],(30)1-θ1+θT~fx<1which leads to(31)RC,nTfx≤-1n+1T~fn+1x1+θT~fx.From the above expression,(32)RC,nTfx⟶n→∞0⟺T~fx≤1which implies that the convergence region is  -1≤T¯f(x)/cf≤1 as opposed to (27).As in the case of the Lagrange remainder, the bound of the Cauchy remainder is calculated using IA as(33)BRC,n=-1n+1BT~fxn+1B1+T~fx.


### 6.2. Square Root

Consider now the square root functiong(Tf(x))=Tf(x). The resulting TM Tg(x) is [7](34)Tfx=cf+∑k=1n-1k-1cf2k-3‼k!2kT¯fkxcfk+RnTfx,where the remainder in Lagrange form is(35)RL,nTfx=CLT¯fn+1x/cfn+11+θT¯fx/cfn+1/2,CL=-1ncf2k-1‼n+1!2n+1.The bound of the above remainder is computed by mean of IA as(36)BRL,nTfx=CLBT¯fxn+1/cfn+11+BT¯fx/cfn+1/2.However, also in this case(37)RL,nTfx⟶n→∞0⟺0≤T¯fxcf≤1.The use of the Cauchy remainder, which reads [19](38)RC,nTfx=CCT~fn+1x1+θT~fx1-θ1+θT~fxn,where(39)CC=-1ncf2n-1‼n+1n+1!2n+1,allows again mitigating the convergence issue. Under the assumptions that lead to (30), the corresponding bound is obtained as(40)RC,nTfx=CCBT~fxn+11+BT~fx.


### 6.3. Multiplicative Inverse

Consider the functiong(Tf(x))=1/Tf(x). The resulting TM Tg(x) is [7](41)1Tfx=∑k=0n-1kT¯fkcfk+1+RnTfx,where in principle the remainder can be computed via either the Lagrange or Cauchy formulation. However, in this special case, anexact remainder is calculated based on the properties of geometric series [17].To begin with, the Taylor expansion of1/(1-x) up to the order n is considered(42)11-x=∑k=0nxk+Rnx.The exact remainder forx<1 is obtained via the identity(43)Rnx=∑k=n+1∞xk=limN→∞⁡xn+1-xN+11-x=xn+11-x.The above equation is recast to obtain the remainder of1/Tf(x) as(44)RnTfx=1cfT~fn+1x1-T~fx,where in this case T~f=-T¯f/cf. The bound of the above remainder is computed as(45)BRnTfx=1cfBT~fxn+1B1-T~fxfor |T~f|<1.


### 6.4. Sine and Cosine

For the case of sine and cosine functions, the coefficientsg(k)(cf) in (18) are given by [7](46)gkcf=-1k/2sin⁡cfkeven-1k-1/2cos⁡cfotherwise,(47)gkcf=-1k/2cos⁡cfkeven-1k-1/2+1sin⁡cfotherwise,respectively.The bounds of the Lagrange remainder are calculated as(48)BRnTfx=Bgn+1Tfxn+1!BT¯fxn+1.By substituting (46) or (47) in (48), the determination of the bounds requires the calculation of the bounds B(sin⁡(Tf)) or B(cos⁡(Tf)) in the pertinent domain of Tf(x). 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 B(sin⁡Tf)=B(cos⁡(Tf))=[-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 ζ∈B(Tf(x)). This leads to [11](49a)min⁡sin⁡ζ=-13π2+2πk∈BTfmin⁡sin⁡BTfxotherwise(49b)max⁡sin⁡ζ=+1π2+2πk∈BTfmax⁡sin⁡BTfxotherwise(49c)min⁡cos⁡ζ=-1π+2πk∈BTfmin⁡cos⁡BTfxotherwise(49d)max⁡cos⁡ζ=+12πk∈BTfmax⁡cos⁡BTfxotherwise,with k∈Z.


### 6.5. Illustrative Example

As an illustrative example, the functiong(x)=sin⁡1/log⁡x with x∈[3,9] 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 g(x). 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 n=6 and n=9, respectively. It is appreciated that increasing the order yields tighter bounds. However, both bounds are conservative over the entire domain D of the input variable x.Figure 1
TM approximation of the functiong(x). Black line: actual value; dashed red line: polynomial approximation provided by a ninth-order TM; green and magenta lines: bounds obtained with orders 6 and 9, respectively.


## 6.1. Logarithm

Consider the logarithmic functiong(Tf(x))=log⁡(Tf(x)). According to [7], the resulting TM  Tg(x) is computed as(22)log⁡Tfx=log⁡cf+∑k=1n-1k+11kT¯fkxcfk+RnTfx,where the remainder is given in the Lagrange formulation (20) as(23)RL,nTfx=-1n+2n+1T¯fn+1x/cfn+11+θT¯fx/cfn+1.The bound of the remainder is obtained by means of the IA as follows:(24)BRL,nTfx=-1n+2n+1BT¯fxn+1/cfn+1B1+T¯fx/cfn+1,which is equivalent to(25)BRL,nTfx=-1n+2n+1BT¯fn+1BTfn+1.Unfortunately, the Lagrange remainder does not always converge to zero when the expansion ordern is increased [19]. In particular, from (24),(26)RL,nTfx⟶n→∞0⟺BT¯f/cf1+T¯fx/cf<1which implies that the Lagrange remainder converges only if(27)0≤T¯fxcf≤1.However, since T¯f(x) 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(28)log⁡Tfx=log⁡cf+log⁡1+T~fx,where T~f(x)=T¯f(x)/cf, the Cauchy remainder is expressed as(29)RC,nTfx=-1n+11-θ1+θT~fxnT~fn+1x1+θT~fx.Since 0<θ<1 and under the (quite practical) assumption T~f<1 (i.e., the variation of a variable is smaller than its central value) [19],(30)1-θ1+θT~fx<1which leads to(31)RC,nTfx≤-1n+1T~fn+1x1+θT~fx.From the above expression,(32)RC,nTfx⟶n→∞0⟺T~fx≤1which implies that the convergence region is  -1≤T¯f(x)/cf≤1 as opposed to (27).As in the case of the Lagrange remainder, the bound of the Cauchy remainder is calculated using IA as(33)BRC,n=-1n+1BT~fxn+1B1+T~fx.


## 6.2. Square Root

Consider now the square root functiong(Tf(x))=Tf(x). The resulting TM Tg(x) is [7](34)Tfx=cf+∑k=1n-1k-1cf2k-3‼k!2kT¯fkxcfk+RnTfx,where the remainder in Lagrange form is(35)RL,nTfx=CLT¯fn+1x/cfn+11+θT¯fx/cfn+1/2,CL=-1ncf2k-1‼n+1!2n+1.The bound of the above remainder is computed by mean of IA as(36)BRL,nTfx=CLBT¯fxn+1/cfn+11+BT¯fx/cfn+1/2.However, also in this case(37)RL,nTfx⟶n→∞0⟺0≤T¯fxcf≤1.The use of the Cauchy remainder, which reads [19](38)RC,nTfx=CCT~fn+1x1+θT~fx1-θ1+θT~fxn,where(39)CC=-1ncf2n-1‼n+1n+1!2n+1,allows again mitigating the convergence issue. Under the assumptions that lead to (30), the corresponding bound is obtained as(40)RC,nTfx=CCBT~fxn+11+BT~fx.


## 6.3. Multiplicative Inverse

Consider the functiong(Tf(x))=1/Tf(x). The resulting TM Tg(x) is [7](41)1Tfx=∑k=0n-1kT¯fkcfk+1+RnTfx,where in principle the remainder can be computed via either the Lagrange or Cauchy formulation. However, in this special case, anexact remainder is calculated based on the properties of geometric series [17].To begin with, the Taylor expansion of1/(1-x) up to the order n is considered(42)11-x=∑k=0nxk+Rnx.The exact remainder forx<1 is obtained via the identity(43)Rnx=∑k=n+1∞xk=limN→∞⁡xn+1-xN+11-x=xn+11-x.The above equation is recast to obtain the remainder of1/Tf(x) as(44)RnTfx=1cfT~fn+1x1-T~fx,where in this case T~f=-T¯f/cf. The bound of the above remainder is computed as(45)BRnTfx=1cfBT~fxn+1B1-T~fxfor |T~f|<1.


## 6.4. Sine and Cosine

For the case of sine and cosine functions, the coefficientsg(k)(cf) in (18) are given by [7](46)gkcf=-1k/2sin⁡cfkeven-1k-1/2cos⁡cfotherwise,(47)gkcf=-1k/2cos⁡cfkeven-1k-1/2+1sin⁡cfotherwise,respectively.The bounds of the Lagrange remainder are calculated as(48)BRnTfx=Bgn+1Tfxn+1!BT¯fxn+1.By substituting (46) or (47) in (48), the determination of the bounds requires the calculation of the bounds B(sin⁡(Tf)) or B(cos⁡(Tf)) in the pertinent domain of Tf(x). 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 B(sin⁡Tf)=B(cos⁡(Tf))=[-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 ζ∈B(Tf(x)). This leads to [11](49a)min⁡sin⁡ζ=-13π2+2πk∈BTfmin⁡sin⁡BTfxotherwise(49b)max⁡sin⁡ζ=+1π2+2πk∈BTfmax⁡sin⁡BTfxotherwise(49c)min⁡cos⁡ζ=-1π+2πk∈BTfmin⁡cos⁡BTfxotherwise(49d)max⁡cos⁡ζ=+12πk∈BTfmax⁡cos⁡BTfxotherwise,with k∈Z.


## 6.5. Illustrative Example

As an illustrative example, the functiong(x)=sin⁡1/log⁡x with x∈[3,9] 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 g(x). 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 n=6 and n=9, respectively. It is appreciated that increasing the order yields tighter bounds. However, both bounds are conservative over the entire domain D of the input variable x.Figure 1
TM approximation of the functiong(x). Black line: actual value; dashed red line: polynomial approximation provided by a ninth-order TM; green and magenta lines: bounds obtained with orders 6 and 9, respectively.


## 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 Tf(x)=Pf(x)+If. 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 anexact solution for the computation of the inversion A+B-1 when B has unitary rank. Under this assumption, the resulting inverse matrix is(50)A+B-1=A-1-11+αA-1BA-1,where α=trace(A-1B) and only the inversion of matrix A is required.For its application to a matrix TMTf, this is first split into a constant matrix and a shifted TM, as was done in Section 6 for the scalar case:(51)Tf=cf+P¯fx-x0+If=cf+T¯f,where cf=Pf(0) is a standard matrix and P¯f(x-x0)=Pf(x-x0)-cf. If the nonconstant part T¯f has unitary rank, then (52) is applied directly and the inverse TM is(52)Tg=Tf-1=cf+T¯f-1=cf-1-11+αcf-1T¯fcf-1,with α=trace(cf-1T¯f). 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 matrixcf, 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 inverse1/(1+α), which is carried out as described in Section 6.3.It is worth noting that in practice the matrix cf is always invertible as it corresponds to the matrix for the central value of the variable x.The generalization of (52) to the case of full-rank matrices requires an iterative application of the algorithm [10]. The matrix T¯f is first split into the sum of rank-one matrices:(53)T¯f=T¯f,1+T¯f,2+T¯f,3+⋯,where T¯f,i is defined as a matrix TM with null entries except for the ith column, which coincides with the ith column of the original TM T¯f, thus being rank-1 by construction. The following iterative procedure is then put forward:(1)
the inversioncf+T¯f,1-1 is computed via (50) by considering A=cf and B=T¯f,1;(2)
the inversioncf+T¯f,1+T¯f,2-1 is calculated by letting A=cf+T¯f,1 and B=T¯f,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.


### 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:(54)Mx=exp⁡1+0.3xlog⁡3+xsin⁡π+π4xlog⁡3+xexp⁡1+0.3xlog⁡π+π4xsin⁡π+π4xlog⁡π+π4xexp⁡1+0.3x,where x∈D=[-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–23]. It is important to remark that detMx≠0  ∀x∈[-1,1].The following MATLAB script calculates the inverse ofM(x) 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 TMTM-1(x) provides a parametric representation of the nonlinear matrix M-1(x) over the entire domain D and a conservative estimation of its upper and lower bounds. The minimum and maximum values of the entries of M-1(x), estimated based on 10’000 MC samples, are(55)min⁡M-1x=0.3754-0.2252-0.1598-0.22520.4216-0.2252-0.1598-0.22520.3754,max⁡M-1x=0.6055-0.15340.1487-0.15340.6022-0.1534-0.1487-0.15340.6055,respectively.These bounds are also estimated from the TM by considering that[Mx)-1ij∈BTM-1xij. This, for an expansion order n=6, leads to(56)min⁡BTM-1x=0.3229-0.2831-0.1720-0.28550.3967-0.2820-0.1795-0.28630.3220,max⁡BTM-1x=0.6230-0.11270.2542-0.11030.6614-0.11370.2617-0.10950.6238.Furthermore, Figure2 compares the behavior of M-1x11 (black line) with the polynomial model provided by a TM of order n=6 (dashed red line). The bounds for expansion orders of n=3 and n=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.Figure 2
TM approximation ofM-1x11. Black line: actual value; dashed red line: polynomial approximation provided by a sixth-order TM; green and magenta lines: bounds obtained with orders 3 and 6, respectively.It is interesting to calculate the product betweenTM and its inverse. In classical algebra, this would result in the identity matrix. With the proposed TM framework and an expansion order n=3, the result is(57)TMx·TM-1x=100010001︸Px+±1.262±1.226±1.190±1.306±1.305±1.274±1.289±1.268±1.252︸I.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 n=6 yields(58)I=±0.095±0.086±0.075±0.108±0.096±0.086±0.106±0.095±0.087.


## 7.1. Matrix Inversion

The SM formula provides anexact solution for the computation of the inversion A+B-1 when B has unitary rank. Under this assumption, the resulting inverse matrix is(50)A+B-1=A-1-11+αA-1BA-1,where α=trace(A-1B) and only the inversion of matrix A is required.For its application to a matrix TMTf, this is first split into a constant matrix and a shifted TM, as was done in Section 6 for the scalar case:(51)Tf=cf+P¯fx-x0+If=cf+T¯f,where cf=Pf(0) is a standard matrix and P¯f(x-x0)=Pf(x-x0)-cf. If the nonconstant part T¯f has unitary rank, then (52) is applied directly and the inverse TM is(52)Tg=Tf-1=cf+T¯f-1=cf-1-11+αcf-1T¯fcf-1,with α=trace(cf-1T¯f). 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 matrixcf, 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 inverse1/(1+α), which is carried out as described in Section 6.3.It is worth noting that in practice the matrix cf is always invertible as it corresponds to the matrix for the central value of the variable x.The generalization of (52) to the case of full-rank matrices requires an iterative application of the algorithm [10]. The matrix T¯f is first split into the sum of rank-one matrices:(53)T¯f=T¯f,1+T¯f,2+T¯f,3+⋯,where T¯f,i is defined as a matrix TM with null entries except for the ith column, which coincides with the ith column of the original TM T¯f, thus being rank-1 by construction. The following iterative procedure is then put forward:(1)
the inversioncf+T¯f,1-1 is computed via (50) by considering A=cf and B=T¯f,1;(2)
the inversioncf+T¯f,1+T¯f,2-1 is calculated by letting A=cf+T¯f,1 and B=T¯f,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.


## 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:(54)Mx=exp⁡1+0.3xlog⁡3+xsin⁡π+π4xlog⁡3+xexp⁡1+0.3xlog⁡π+π4xsin⁡π+π4xlog⁡π+π4xexp⁡1+0.3x,where x∈D=[-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–23]. It is important to remark that detMx≠0  ∀x∈[-1,1].The following MATLAB script calculates the inverse ofM(x) 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 TMTM-1(x) provides a parametric representation of the nonlinear matrix M-1(x) over the entire domain D and a conservative estimation of its upper and lower bounds. The minimum and maximum values of the entries of M-1(x), estimated based on 10’000 MC samples, are(55)min⁡M-1x=0.3754-0.2252-0.1598-0.22520.4216-0.2252-0.1598-0.22520.3754,max⁡M-1x=0.6055-0.15340.1487-0.15340.6022-0.1534-0.1487-0.15340.6055,respectively.These bounds are also estimated from the TM by considering that[Mx)-1ij∈BTM-1xij. This, for an expansion order n=6, leads to(56)min⁡BTM-1x=0.3229-0.2831-0.1720-0.28550.3967-0.2820-0.1795-0.28630.3220,max⁡BTM-1x=0.6230-0.11270.2542-0.11030.6614-0.11370.2617-0.10950.6238.Furthermore, Figure2 compares the behavior of M-1x11 (black line) with the polynomial model provided by a TM of order n=6 (dashed red line). The bounds for expansion orders of n=3 and n=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.Figure 2
TM approximation ofM-1x11. Black line: actual value; dashed red line: polynomial approximation provided by a sixth-order TM; green and magenta lines: bounds obtained with orders 3 and 6, respectively.It is interesting to calculate the product betweenTM and its inverse. In classical algebra, this would result in the identity matrix. With the proposed TM framework and an expansion order n=3, the result is(57)TMx·TM-1x=100010001︸Px+±1.262±1.226±1.190±1.306±1.305±1.274±1.289±1.268±1.252︸I.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 n=6 yields(58)I=±0.095±0.086±0.075±0.108±0.096±0.086±0.106±0.095±0.087.


## 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 Figure3 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 zs and of the tire zu from the reference frame. The road profile h is considered as the excitation of the system.Figure 3
Schematic of the quarter-car model considered in the first application, with definition of parameters and state variables.The proposed quarter-car model is defined by the following parameters, the masses of the chassisms and of the tire mu, the spring coefficients of the suspension K and of the tire P, and the damping coefficients of the suspension C and of the tire Ct.According to [24], the quarter-car model is described by the following matrix equation:(59)ms00muz¨stz¨ut+C-C-CCz˙stz˙ut+K-K-KK+Pzstzut=0Cth˙t+Pht,where zs and zu 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(60)x˙=Ax+Bu,where the column vector x(t)=vst,vut,zst,zutT∈R4 collects the state variables of the system and u(t)=ht,h˙tT∈R2 is a vector defining the sources of the system. The matrix A∈R4×4 and the column vector B∈R4×2 are defined as(61)A=-M-1C-M-1KI0,B=M-1F0,where I and 0 denote the identity matrix and the null matrix of dimensions 2×2, respectively. The remaining terms are defined as follows:(62a)M=ms00mu(62b)C=C-C-CC+Ct(62c)K=K-K-KK(62d)F=00PCt.The frequency-domain response of the system, denoted asX(ω), is obtained via the Fourier transform of (3) and by considering as excitation vector u(t)=1/2πexp⁡(jωt)1,jωT:(63)Xω=jωI-A-1BUω,where U(ω)=[1,jω]T.For the simulation, it is assumedms=332.5 kg, mu=38 kg, and Ct=0 Ns/m. In addition, three independent uncertainty variables are considered; that is, C=(6000±20%) Ns/m, K=(22000±20%) N/m, and P=(20000±20%) N/m.Figure4 shows the variation (gray area) of the real (top panel) and imaginary (bottom panel) part of the transfer function T(ω)=Zs(ω)/H(ω), related to the state variable Zs(ω), 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.Figure 4
Real and imaginary parts of the transfer functionT(ω)=Zs(ω)/H(ω). Gray area: spread of the response resulting from parameter uncertainty and obtained with MC; magenta lines: WC bounds estimated from a fifth-order TM.
(a)
(b)Figure 5
Magnitude of the transfer functionT(ω)=Zs(ω)/H(ω). Line identification is as in Figure 4.


### 8.2. Continuous-Time Linear Equalizer

The second proposed application example deals with the analysis of three CTLEs whose electrical circuits are depicted in Figure6. 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)
R=35Ω±20%, C=1pF±20%, L=3nH±20%;(b)
R=100Ω±20%, C=1pF±20%, L=3nH±20%;(c)
R=100Ω±20%, C=1pF±20%.Figure 6
Schematics of the three different CTLE configurations considered.
(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 “S21”, 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 S21 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.Figure 7
Real part, imaginary part, and magnitude of the transmission coefficientS21 in the three CTLE configurations of Figure 6. Gray area: spread of the response from a MC analysis. Magenta lines: WC bounds obtained with a fourth-order TM.
(a)
(b)
(c)


## 8.1. Quarter-Car Model

The first example considers the quarter-car model depicted in Figure3 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 zs and of the tire zu from the reference frame. The road profile h is considered as the excitation of the system.Figure 3
Schematic of the quarter-car model considered in the first application, with definition of parameters and state variables.The proposed quarter-car model is defined by the following parameters, the masses of the chassisms and of the tire mu, the spring coefficients of the suspension K and of the tire P, and the damping coefficients of the suspension C and of the tire Ct.According to [24], the quarter-car model is described by the following matrix equation:(59)ms00muz¨stz¨ut+C-C-CCz˙stz˙ut+K-K-KK+Pzstzut=0Cth˙t+Pht,where zs and zu 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(60)x˙=Ax+Bu,where the column vector x(t)=vst,vut,zst,zutT∈R4 collects the state variables of the system and u(t)=ht,h˙tT∈R2 is a vector defining the sources of the system. The matrix A∈R4×4 and the column vector B∈R4×2 are defined as(61)A=-M-1C-M-1KI0,B=M-1F0,where I and 0 denote the identity matrix and the null matrix of dimensions 2×2, respectively. The remaining terms are defined as follows:(62a)M=ms00mu(62b)C=C-C-CC+Ct(62c)K=K-K-KK(62d)F=00PCt.The frequency-domain response of the system, denoted asX(ω), is obtained via the Fourier transform of (3) and by considering as excitation vector u(t)=1/2πexp⁡(jωt)1,jωT:(63)Xω=jωI-A-1BUω,where U(ω)=[1,jω]T.For the simulation, it is assumedms=332.5 kg, mu=38 kg, and Ct=0 Ns/m. In addition, three independent uncertainty variables are considered; that is, C=(6000±20%) Ns/m, K=(22000±20%) N/m, and P=(20000±20%) N/m.Figure4 shows the variation (gray area) of the real (top panel) and imaginary (bottom panel) part of the transfer function T(ω)=Zs(ω)/H(ω), related to the state variable Zs(ω), 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.Figure 4
Real and imaginary parts of the transfer functionT(ω)=Zs(ω)/H(ω). Gray area: spread of the response resulting from parameter uncertainty and obtained with MC; magenta lines: WC bounds estimated from a fifth-order TM.
(a)
(b)Figure 5
Magnitude of the transfer functionT(ω)=Zs(ω)/H(ω). Line identification is as in Figure 4.


## 8.2. Continuous-Time Linear Equalizer

The second proposed application example deals with the analysis of three CTLEs whose electrical circuits are depicted in Figure6. 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)
R=35Ω±20%, C=1pF±20%, L=3nH±20%;(b)
R=100Ω±20%, C=1pF±20%, L=3nH±20%;(c)
R=100Ω±20%, C=1pF±20%.Figure 6
Schematics of the three different CTLE configurations considered.
(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 “S21”, 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 S21 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.Figure 7
Real part, imaginary part, and magnitude of the transmission coefficientS21 in the three CTLE configurations of Figure 6. Gray area: spread of the response from a MC analysis. Magenta lines: WC bounds obtained with a fourth-order TM.
(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.



---

*Source: 6739857-2017-03-30.xml*