A general framework for approximated model stable inversion

Automatica 101 (2019) 182–189

Contents lists available at ScienceDirect

Automatica journal homepage: www.elsevier.com/locate/automatica

Brief paper

A general framework for approximated model stable inversion✩ Raffaele Romagnoli *, Emanuele Garone Université Libre de Bruxelles, B-1050 Bruxelles, Belgium

article

info

Article history: Received 31 July 2017 Received in revised form 28 September 2018 Accepted 13 November 2018

Keywords: Feedforward control Stable inversion Least squares Basis functions

a b s t r a c t In this paper we propose an approximated method for the design of feedforward actions based on model stable inversion. In the proposed framework the solution of the model stable inversion is approximated in the least squares sense by a suitable combination of any arbitrary set of basis functions. The paper makes use of operator theory to derive a closed-form formula to compute such approximation and to provide boundedness properties for the computed solution. Simulations are provided to show the effectiveness of the proposed method with different families of basis functions. © 2018 Elsevier Ltd. All rights reserved.

1. Introduction Model inversion (Horowitz, 2013) represents a classic approach for the design of feedforward inputs capable of making the output follow a pre-defined trajectory. As well known for nonminimum phase (NMP) systems, inversion generates a diverging feedforward input (Moylan, 1977). For this kind of systems, inversion is addressed by the so-called model stable inversion methods. The aim of stable inversion methods is to find a stable solution for the model inversion problem for NMP systems. A large amount of methods have been proposed in the literature, many of which arose to solve specific practical problems. A possible classification of the existing methods consists in distinguishing between exact and approximated stable inversion methods (Rigney, Pao, & Lawrence, 2009). In line of principle, using pre-actuation, the model stable inversion problem admits an exact solution. However such a solution may require an infinite pre-actuation (Devasia, Chen, & Paden, 1996; Hunt, Meyer, & Su, 1996), which makes it not applicable from the practical viewpoint. ✩ This work is performed in the framework of the BATWAL project financed by the Walloon region (Belgium). This research has been funded by the Mandats d’Impulsion Scientific ‘‘Optimization-free Control of Nonlinear Systems subject to Constraints’’ of the Fonds de la Recherche Scientifique (FNRS), Ref. F452617F. The material in this paper was not presented at any conference. This paper was recommended for publication in revised form by Associate Editor Denis Arzelier under the direction of Editor Richard Middleton. Corresponding author. E-mail addresses: [email protected] (R. Romagnoli), [email protected] (E. Garone).

*

https://doi.org/10.1016/j.automatica.2018.11.044 0005-1098/© 2018 Elsevier Ltd. All rights reserved.

To face this problem, methods based on a finite preview of the desired output that generate a finite pre-actuation have been proposed in Zou and Devasia (1999, 2007). Note that, the length of the desired output preview is proportional to the length of the pre-actuation. In many applications a long preview may not be acceptable. For this reason the model inversion problem has been addressed proposing approximated methods that modify the system to be inverted or the desired output in order to reduce or avoid the preview (Benosman & Le Vey, 2003; Devasia, 1997; Graichen, Hagenmeyer, & Zeitz, 2005). Redesigning or imposing a particular desired output is also important to provide optimal stable inputs in presence of constraints or uncertainties (Clayton, Tien, Leang, Zou, & Devasia, 2009; Piazzi & Visioli, 2001). In the last few years a new family of approximated approaches called pseudo-inversion has been proposed which addresses the problem from a different point of view. In the first paper on the subject (Jetto, Orsini, & Romagnoli, 2013) the idea was to assign a piece-wise polynomial function to the unknown input, and estimate the polynomial’s coefficients using a least square algorithm. Interestingly enough, pseudo-inversion does not need pre-actuation, and can be applied also in the cases where methods using pre-actuation do not work, such as nonhyperbolic systems (Jetto, Orsini, & Romagnoli, 2014). Moreover, it does not require to modify the desired output or the system. Recently, a generalization of the pseudo-inversion has been introduced in Ramani, Duan, Okwudire, and Ulsoy (2017) for SISO discrete time linear systems. This method, called Filtered Basis Function (FBF) method solves the inversion problem over a finite time interval using a general independent basis functions.

R. Romagnoli and E. Garone / Automatica 101 (2019) 182–189

Other approaches based on the least square pseudo-inversion and B-splines basis functions are presented in Duan, Yoon, and Okwudire (2017) and Jetto, Orsini, and Romagnoli (2017). It is also interesting to remark that the use of basis function for control is common in the context of iterative learning control (ILC) (Bolder & Oomen, 2015) and that inversion represents a basis for many ILC algorithms (van Zundert & Oomen, 2017; Wang, Kim, & Zou, 2013). An approach that merges iterative control, classical inversion approach, least squares and Fourier basis functions has been proposed in Zhang and Liu (2014). In this paper we propose a novel and general framework able to solve the pseudo-inversion problem for MIMO continuous-time system using any arbitrary set of bounded L2 basis functions. This framework includes as particular cases all existing methods based on basis functions and least squares, such as pseudo-inversion (Jetto, Orsini, & Romagnoli, 2015a, c; Jetto et al., 2017) and FBF (Duan et al., 2017; Ramani et al., 2017). Moreover the presented method does not need any pre-actuation, can be applied to systems with zeros on the imaginary axis, and does not need any output or system redefinition. Furthermore, the proposed method allows also to define the conditions for obtaining a solution with a prescribed degree of smoothness. From a practical viewpoint this last condition is extremely important in order to provide smooth control signals to the actuators. The effectiveness of the method is demonstrated in simulation. To facilitate its use by third parties, a MATLAB toolbox implementing the proposed method is made available in ToolboxMatlab (2018). The paper is organized as follows: after the problem statement in Section 1, the proposed solution is given in Section 2 where the main theorem is stated. In Section 3 the introduction of smoothness constraints is discussed. The effectiveness of the proposed method is demonstrated by means of numerical examples in Section 4. In Appendix the proof of the main theorem stated in Section 2 is reported. 2. Problem statement Consider the continuous-time linear time invariant (LTI) system x˙ (t) = Ax(t) + Bu(t) y(t) = Cx(t) + Du(t) q

(1) p

n

with y(·) ∈ R , u(·) ∈ R and x(·) ∈ R denoting the output, the input, and the state, respectively. Let yd (·) be a desired output function that we would like to track. We can define the tracking error as e(·) ≜ yd (·) − y(·).

(2)

183

Problem 1 admits a solution only under very precise conditions. The first condition is that yd (t) must be a signal that (at least for some suitable initial condition) can be followed perfectly without violating (3). This condition can be overcome by assuming that yd (t) is sufficiently smooth. The second (and more relevant) condition is that, for a given yd (t), requirements (3) and (4) can be jointly satisfied only for some suitable initial states x(0). To overcome this problem, the classic approach used in the literature (Seifried, 2014) is to relax the problem and to assume that it is possible to use a pre-actuation, i.e. that it is possible to control the system from time t = −τ so as to ensure that x(0) is so that Problem 1 admits a solution. However the theory (Hirschorn, 1979; Moylan, 1977; Silverman, 1969) ensures that this can be done using a finite τ and a bounded input u(t) ∈ [−τ , ∞) only for minimum phase systems. For nonminimum phase systems, the previously mentioned inversion algorithms produce an unbounded input u(t). As shown in Devasia et al. (1996) and Hunt et al. (1996) in the general case a bounded u(t) can be obtained using infinite horizon pre-actuation (i.e. τ → ∞) and only in the case the system has no zeros on the imaginary axis. The above discussion shows that, despite the theoretical importance of Problem 1, from a practical viewpoint (and especially in the case of nonminimum phase systems) only approximated solutions can be used. In this paper we will consider one of the possible approximated formulations of Problem 1 based on the following standing assumptions A1 the matrix A of (1) is Hurwitz; A2 yd (t) is such that there exists a known input signal us (t) satisfying condition (3) such that

∫

+∞

eTs (τ )es (τ )dτ < +∞.

(5)

0

where es is defined as es (·) ≜ yd (·) − ys (·) and ys (·) is the response of the system for the known input us (·). Remark 2. Assumptions A1 is a standard assumption for stable inversion problems. Assumption A2 is reasonable for a wide range of problems. For instance, for all the cases where the asymptotic behavior of the desired output yd (t) can be represented as a rational transfer function with a finite number of poles/zeros in the nonnegative half plane and the system is right-invertible, a possible us (t) can be a computed in closed form. A complete report on how to compute us (t) in these cases, can be found in Romagnoli and Garone (2018). In this paper the starting idea is to rewrite the input u(t) as the sum of two contributions

In this paper we focus on the model stable inversion problem, which can be defined as follows.

u(·) = us (·) + ut (·).

Problem 1 (Model Stable Inversion Problem). Let an initial state x(0) and a desired output signal yd (t), t ∈ [0, +∞), bounded for any bounded t . Compute an input signal u(t), t ∈ [0, +∞) that

where us (·) is as per Assumption 2, and ut (·) is a signal to be computed. By the superposition property of (1), the system output can be written as y(·) = ys (·) + yt (·), where yt can be seen as the transient response given by ut (·). Hence, we can define the transient output tracking error

• it is bounded with respect to yd (t) ∃ k∗ , M ∗ : ∥u(t)∥ ≤ M ∗ ∥yd (t)∥ + k∗ , ∀t ≥ 0;

(3)

(7)

which allows us to define an Approximated Model Stable Inversion Problem as follows

• makes the error (2) equal to zero. e(t) = 0.

et (·) ≜ es (·) − yt (·).

(6)

(4)

Note that, this definition also includes the case that yd (t) does not admit a limit for t → ∞ (e.g. if it tends to infinity), as far as u(t) does not ‘‘grow faster’’ than yd (t).

Problem 3 (Approximated Model Stable Inversion Problem). Let an initial state x(0) and a desired output signal yd (t), and let us (t), t ∈ [0, +∞) as per Assumption 2, compute a finite-energy transient input signal uˆ t (t), t ∈ [0, +∞) such that

184

R. Romagnoli and E. Garone / Automatica 101 (2019) 182–189

functions ϕj (t) ∈ L2 ([0, +∞))p , j = 1, . . . , nb , i.e.

• there exists ¯ : ∥ˆut (t)∥ ≤ M ¯ , ∀t ≥ 0; M

(8)

• the L2 norm of the transient tracking error is minimized, i.e. ∫ ∞ (9) min eTt (τ )et (τ )dτ . uˆ t (t)

0

The solution of Problem 3 allows to build an approximated solution of Problem 1 as uˆ (·) = us (·) + uˆ t (·).

(10)

Remark 4. Problem 3, is one of the many possible approximated reformulations of Problem 1. As will be clear in the next sections, the main advantage of this reformulations is that it allows us to use L2 arguments and operators theory, and eventually leads to a simple closed-form approximated solution. It is also worth to note that this reformulation is the one (implicitly) used for all the socalled pseudo-inversion techniques proposed in the literature (see e.g. Duan et al., 2017; Jetto et al., 2014, 2015c, 2017; Ramani et al., 2017).

nb ∑

uˆ t (t) =

λj ϕj (t) = ϕ (t)λ,

where the coefficients λj are the unknown parameters to be computed, λ = [λ1 , . . . , λnb ], and ϕ (t) = [ϕ1 (t), . . . , ϕnb ]. The aim of this paper is to provide a closed form solution to compute the coefficients (15) which give the best approximation of Problem 3 for the selected set of basis functions. To this end, let us first define the kernel of (14) as K (t , τ ) ≜ CeA(t −τ ) B. The following theorem summarizes the main result of this paper Theorem 5. Let uˆ t (t) as in (15) and define

Φ≜

+∞

∫

τ

(∫

K (τ , θ )ϕ (θ )dθ

0

0

The approximated model stable inversion problem is based on the transient output tracking error defined in (7). By assumption A2, us exists and is known, hence es can be computed as follows es (t) = yd (t) − C At x0 +

t

∫

CeA(t −τ ) Bus (τ )dτ .

(11)

0







ys (t)

The transient response is t

∫

CeA(t −τ ) Buˆ t (τ )dτ .

yt (t) =

(12)

0

At this point, the minimization problem (9) can be rewritten as

∫

∞

(

min uˆ t (·)

e s (τ ) −

∫

τ

CeA(τ −θ ) uˆ t (θ )dθ

0

0

(

es (τ ) −

τ

∫

·

CeA(τ −θ ) uˆ t (θ )dθ

)

dτ .

(13)

It is interesting to note that problem (13) is equivalent (Beutler, 1965) to seeking for the ’’best approximated solution’’ (in the operator theory sense) of the linear integral equation

∫

t

CeA(t −τ ) Buˆ t (τ )dτ

· K (τ , θ )ϕ (θ )dθ

)

dτ ,

(16)

0

Ψ (es (·)) ≜

τ

∫ +∞ (∫

K (τ , θ )ϕ (θ )dθ

)T

es (τ )dτ ,

(17)

0

where Φ is a matrix in Rnb ×nb and Ψ is a vector in Rnb . The least square solution of the linear integral equation (14) is the vector λ ∈ Rnb which can be computed in closed-form as

λ = Φ † Ψ (es )

(18)

†

where denotes the usual Moore–Penrose pseudo-inverse. Under assumptions A1 and A2, the above equation always admits a finite solution. Furthermore, the obtained uˆ t (t) is bounded on the time p interval [0, +∞), is of class L2 , and the overall control input uˆ (t) = us (t) + uˆ t (t) is an approximated solution for Problem 1. Proof. The main idea of the proof is to resort to the operator theory (Kammerer & Nashed, 1972) which not only allows to compute the above closed-form, but also gives guarantee of solvability and boundedness. For details, please refer to Appendix.

)T

0

es (t) =

)T

τ

(∫

0

3. Proposed solution

(15)

j=1

t ∈ [0, +∞].

(14)

0

Since both es and uˆ t have finite energy (and so they belong to the Hilbert spaces L2 ([0, +∞))q and L2 ([0, +∞))p , respectively) Eq. (14) is a well-known linear integral equation called Fredholm integral equation of the first kind (Kress, 2014). It is known that, in the general case, the exact solution of this linear integral equation may not exist (Kress, 2014). However, even if (14) does not admit a solution, it is always possible to find an approximation (Kammerer & Nashed, 1972) which is the ‘‘best approximation’’ in the operator theory sense. Note that in general this solution belongs to the span of the infinite dimensional space L2 ([0, +∞))p and so its exact computation is prohibitive from the computational viewpoint. In this paper we propose a finite dimensional approximation of Problem 3 obtained by representing uˆ t (t) as the linear combination of a finite number of arbitrary linearly independent bounded basis

Remark 6. Note that the result of Theorem 5 is very general as it works for whatever finite set of basis functions φj , j = 1, . . . , nb such that φj is bounded and is of class L2 . In this view, many results presented in the literature (e.g. pseudo-inversion approaches (Jetto et al., 2014, 2015a; Jetto, Orsini, & Romagnoli, 2015b) and prefiltered basis function (FBF) methods (Ramani et al., 2017)) can be seen as a special case of the method proposed in this paper using different basis functions (B-splines, radial basis functions, piecewise polynomial functions, etc.). It is also worth noticing that if instead of working on an infinite horizon we consider the case of finite time intervals, assumptions A1 and A2 can be dropped and we can always choose us = 0. Remark 7. Note that, since the matrix A is Hurwitz, the integrals (16) and (17) can be numerically computed with arbitrary precision. Moreover, both integrals depend on the computation of the following convolution integral τ

∫

K (τ , θ )ϕ (θ )dθ

0

which is independent of the particular desired output and thus for any given set of basis functions can be computed once and for all. Remark 8. Note that, clearly, the choice of the ‘‘starting’’ steadystate solution us (t) may affect the overall performances of the proposed inversion method. As a rule of thumb, the more one

R. Romagnoli and E. Garone / Automatica 101 (2019) 182–189

uses a priori knowledge to select a us (t) which makes es (t) small, the better the performance. The other aspects that influence the performance are the choice of the form and the number of basis functions. Clearly the more the basis function has ‘‘a similar nature’’ to the signal one wants to follow and the larger the number of basis functions chosen, the better. It is also interesting to remark that in line of principle an iterative approach could be used: one may start with a ‘‘basic’’ us (t) and a certain set of basis functions, get a solution uˆ r , impose as a new us (t) = us (t) + uˆ r and use a different set of basis function with . Note also that when us (t) does not admit a steady-state does not represent a problem for the proposed method as far as an us (t) ensuring assumption A2 can be found.

4. Continuity and smoothness of the solution The solution uˆ (t) obtained through the proposed method does not guarantee any property of continuity and smoothness. This can be a problem as, in several real applications, possible discontinuities of uˆ (t) can generate undesired behavior. More in general, a certain degree of smoothness of the signal uˆ (t) could be required. If us (t) is of class C k , the following Lemma provides necessary and sufficient conditions to ensure that uˆ t (t) (and in turn uˆ (t)) is of class C k . Lemma 9. Consider a function uˆ t (·) ∈ L2 ([0, +∞))p defined as in (15). Define the set D of time instants where at least one of the basis functions or at least one of their first k derivatives is discontinuous, i.e. ti : ∃m ∈ {0, . . . , k}, j ∈ {1, . . . , nb }

i

The function uˆ t (·) is continuous of class C if and only if

λ ∈ ker(Ck ),

(19)

+

⎤

Ck (t1 ) − Ck (t1 )

⎥ .. ⎦, . − + Ck (tℓ ) − Ck (tℓ )

(20)

⎢

dϕ (t)T

ϕ (t)T

t →ti−

dt dϕ (t)T

[

+

Ck (ti ) ≜ lim

t →ti+

ϕ (t)

T

dt

... ...

dk−1 ϕ (t)T

]T

dt k dk−1 ϕ (t)T

]T

dt k

,

(21)

,

(22)

⏐ ⏐ ⏐

t =ti−

=

dm (uˆ t (t)) ⏐ dt

, m = 0, . . . , k.

(24)

dm (ϕ (t)) ⏐ dt

⏐ ⏐ ⏐

−

]

dm (ϕ (t)) ⏐ dt

t =ti−

⏐ ⏐ ⏐

λ = 0, m = 0, . . . , k

t =ti+

(25)

and finally as Ck (ti− ) − Ck (ti+ ) λ = 0.

[

]

(26)

By imposing the above constraints for all ti ∈ D, it follows that uˆ t (·) is of class C k if and only if λ is such that Ck λ = 0, i.e. λ ∈ Ker{Ck }. Condition (19) can be exploited to reduce the searching space of the stable inversion problem to all and only the parameters that ensure that the resulting uˆ (·) is of class C k . The main idea to do so is to use a coordinate transformation of the vector λ into the null space of Ck . More precisely, if the rank of Ck is nc and nc < nb , it is possible to find nb − nc linearly independent vectors bi ∈ Rnb such that ker(Ck ) = span b1 , . . . , bnb −nc .

{

}

(27)

Defining the mapping matrix M = [b1 , . . . , bnb −nc ] from ker(Ck ) to the space of λ ∈ Rnb , all and only the λ ∈ ker{Ck } can be generated as

¯ λ = M λ,

(28)

(29)

¯ † Ψ¯ (es ), λ¯ = Φ

(30)

+∞

∫

τ

(∫

0

K (τ , θ )ϕ¯ (θ )dθ

0

τ

(∫

)T ·

K (τ , θ )ϕ¯ (θ )dθ

)

⏐ ⏐ ⏐

t =ti+

, m = 0, . . . , k.

τ

+∞ (∫

∫ 0

K (τ , θ )ϕ¯ (θ )dθ

)T

dτ ,

(31)

es (τ )dτ .

(32)

0

5. Numerical example

Proof. For a generic point t ̸ ∈ D, all the functions ϕj (·) are of class C k in a neighborhood of t. Since uˆ t (·) is a linear combination of the basis functions ϕj (·), it follows that also uˆ t (·) is of class C k in a neighborhood of t for whatever choice of the vector λ. As a consequence, the continuity of the function and of the first k derivatives around a point ti ∈ D, is ensured if and only if the following continuity constraints hold true dt

[

Ψ¯ (es (·)) ≜

and ℓ = card{D}.

dm (uˆ t (t)) ⏐

t =ti+

0

Ck (ti ) and Ck (ti ) are Ck (ti− ) ≜ lim

⏐ ⏐ ⏐

dt

t =ti−

+

[

dm (ϕ (t)λ) ⏐

which can be rewritten as

¯ ≜ Φ

Ck ≜ ⎣ −

dt

=

where

where Ck is the matrix of the k-continuity constraints defined as −

⏐ ⏐ ⏐

At this point, define as a new set of basis functions ϕ¯ (t) ≜ ϕ (t)M. The inversion problem with continuity/smoothness constraints is solved applying Theorem 5 and finding the best approximation to our problem in terms of the new variables λ¯ as follows

k

⎡

dm (ϕ (t)λ) ⏐

¯ uˆ t (·) = ϕ (t)M λ,

⏐ ⏐ } dm ϕj (t) ⏐ dm ϕj (t) ⏐ ⏐ ̸= ⏐ . such that dt ⏐t − dt ⏐t + i

Using the definition of uˆ t in (15), Eq. (23) becomes

for any λ¯ ∈ Rnb −nc . This allows to rewrite uˆ t in terms of λ¯ as

{ D=

185

(23)

A flexible robot arm represents a classical example of nonminimum phase system (De Luca, Lanari, Lucibello, Panzieri, & Ulivi, 1990). Several methods of model stable inversion and ILC have been proposed to solve the output tracking problem of this kind of systems, e.g. in Benosman and Le Vey (2003), Jetto et al. (2014), Vakil, Fotouhi, and Nikiforuk (2009), and Zhang and Liu (2014). Consider a two-link robot with a flexible forearm as depicted in Fig. 1. The inputs of the system are the torque τ1 at the joint j1 and the torque τ2 at the second joint j2 . The outputs are the position of the joint j2 that we denote as yj2 , and the terminal of the flexible forearm ytip . Controlling the ytip , the considered system is non-minimum phase. The standard nonlinear model of the system is described in De Luca et al. (1990). For the purpose of this paper

186

R. Romagnoli and E. Garone / Automatica 101 (2019) 182–189

Fig. 1. Two-link robot with flexible forearm. Fig. 2. Desired trajectories (ydj , ydtip ) and steady state response (ysj , ystip ). 2

2

we consider the linearized model proposed in Hussein and Nemah (2015) which is described by x˙¯ (t) = Ax¯ (t) + Bu

(33)

y(t) = C x¯ (t) where A =

(34)

[

0

I

−M −1 × K

−M −1 × H

]

,B =

[

0

]

M −1 × G

[

, C = C1

]

0 .

Please refer to Hussein and Nemah (2015) for further details about the above matrices and the state space vector x¯ (t). The considered system has a nonminimum phase zero at 7.41. The applied control law is the same proposed in De Luca et al. (1990) and Hussein and Nemah (2015) uc (t) = KP θ (t) + KD θ˙ (t) where KP = diag(11.5, 6), KD = diag[2, 0.8], and θ (t) = [θ1 (t), θ2 (t)]T , θ˙ (t) = [θ˙1 (t), θ˙2 (t)]T . The objective is to track the desired trajectories depicted in Fig. 2. The blue dashed line represents the desired trajectory of the position of yj2 , and the red dashed line is the terminal of flexible link ytip . Since the system is asymptotically stable, it is possible to compute the static gain G, which is an invertible matrix, then the tracking of the steady state is achieved using us (t) = G−1 [25◦ , 60◦ ]T . The system response for us (t) is also reported in Fig. 2, where the blue solid and the red solid lines refer to yj2 and ytip respectively. The theory developed in this paper will be applied using three different basis functions:

• B-splines (Boor, 1978). We assume that the input signal uˆ t (t) is represented as a B-spline, i.e. nb

uˆ t (t) =

∑

λi Bi,d (t) t ∈ [tˆ1 , tˆnb +d+1 ]

(35)

i=1 n

n +d+1

where, (λi )i=b 1 is a set of nb control points and (tˆ )i=b 1 is a set of non decreasing knot points. The basis functions Bi,d (t) associate to a B-spline can be obtained using the Cox–de Boor recursion formula and are Bi,d (t) =

t − tˆi tˆi+d − tˆi

Bi,d−1 (t) +

tˆi+1+d − t tˆi+1+d − tˆi+1

Bi+1,d−1 (t)

(36)

and

{ Bi,0 (t) =

1, 0,

tˆi ≤ t ≤ tˆi+1 otherwise

(37)

In this example we consider each knot point fixed at a particular tj . The λi are scalar values, and the choice of the values tj is made partitioning the time axis in a finite number of subintervals. For more details see Boor (1978) and Jetto et al. (2017).

Fig. 3. Position of joint 2 yj2 and position of terminal flexible link ytip . The blue solid line represents the output obtained by B-splines, the red dashed line sinc functions, the green dotted dashed line piecewise polynomial functions, and the black dotted line is the desired output.

• Sinc (Gearhart & Schultz, 1990). For the sinc function, we consider uˆ t (t) as uˆ t (t) =

nb ∑

λi

sin(π (t − ti ))

π (t − ti )

i=1

(38)

where ti are assigned values that divide the time axis in a finite number of subintervals. • Polynomials. For piece-wise polynomial functions, we consider polynomials of order 3. The j−th basis function has the following form

ϕj (t) =

⎧ 3 ⎪ ⎨∑

λ j,k t k ,

t ∈ [tj , tj+1 )

⎪ ⎩ k=0

0

elsewhere,

where, accordingly with (15), λj ≜ [λj,0 , . . . , λj,3 ]T . To compare the three different basis functions, we used in total nb = 320 basis functions, 160 for each input signal, for all the three cases. Fig. 3 shows the system responses for the three different cases compared with the desired outputs. In Table 1 the tracking errors are reported for each output and each used basis function. In this particular application, the use of B-splines shows the best results. The computed inputs are reported in Figs. 4 and 5. To carry out the computation, a small MATLAB toolbox has been realized and is

R. Romagnoli and E. Garone / Automatica 101 (2019) 182–189

187

Table 1 Tracking errors for the different used basis functions.

||e||2

B-spline

Sinc

Polynomial

yj2 ytip

0.047893 0.083262

0.0926 0.109268

0.07797 0.127574

Fig. 4. The solid line represents the input obtained by B-splines, the dashed line the sinc functions. Blue and red lines correspond to uˆ 1 and uˆ 2 respectively. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

now available at ToolboxMatlab (2018) together with the code of this example. Remark 10. The best basis function cannot be determined a priori. In fact one can obtain different results depending on the signal to be tracked. Remark 11. The presence of unstable low frequency zeros represents one of the main limitation of classical inversion methods as, in the case of zeros on the imaginary axis, these methods cannot be applied (Devasia et al., 1996; Hunt et al., 1996; Zou & Devasia, 1999). To overcome this limitation, several solutions have been proposed in the literature which require to modify the system dynamics (Devasia, 1997) and/or the desired output to obtain acceptable results (Benosman & Le Vey, 2003; Graichen et al., 2005). In our case, if assumption A2 is satisfied, the method can be applied independently to the position of the zeros. Clearly the presence of low frequency zeros may slow down the response of the system and so to follow sufficiently well fast signals one has to ensure that the selected family of the basis function is sufficiently rich to provide the requested dynamics, e.g. by selecting a sufficiently high number of basis functions. 6. Conclusions In this paper we have addressed the design of a feedforward controller for the output tracking problem. In order to generalize existing methods based on least squares and basis functions, a general pseudo-inversion method has been proposed. The proposed method always provides a solution to the stable inversion problem, even for nonminimum phase systems. The few assumptions considered, that can be even dropped in the case of finite horizons, make the proposed stable inversion a method that generalizes all the methods presented in literature based on basis functions and least squares. One additional feature of the proposed method is the possibility to choose the degree of smoothness of the solution in case of discontinuous or non smooth basis functions.

Fig. 5. The inputs obtained from piecewise polynomial functions, the blue dashed line is obtained without considering continuity constraints, and the red solid line is the results of the application of Lemma 9.

Appendix For the reader’s convenience some basic notion of operator theory and preliminary results needed for the proof of Theorem 5 are reported. We focus on the linear integral equation described in (14), which can rewritten as t

∫

K (t , τ )uˆ t (τ )dτ = es (t)

t ∈ I ≜ [0, +∞)

0

where K (t , τ ) ≜ CeA(t −τ ) B is defined as the kernel of the linear integral equation. The above linear integral equation can be rewritten in the operator notation as T uˆ t = es

(39)

where T is a linear operator from a Hilbert space Uˆ = L2 ([0, +∞))p to the Hilbert space Y = L2 ([0, +∞))q defined as

∫ (T uˆ t )(t) ≜

t

K (t , τ )uˆ t (τ )dτ

t ∈ I ≜ [0, +∞).

(40)

0

Thanks to assumption A1 (i.e. A is Hurwitz), the linear operator (40) is bounded in the L2 norm, i.e. ∃M ≥ 0 : ∥T ∥2 ≤ M. The objective is to approximate (39) in the least square sense. A best approximation can be always provided, even when (39) is not solvable (Kammerer & Nashed, 1972). A function uˆ t is the least square solution of (39) if it satisfies the normal equation T ∗ T u¯ = T ∗ es

(41)

188

R. Romagnoli and E. Garone / Automatica 101 (2019) 182–189

where T ∗ : Y → Uˆ represents the adjoint operator (Kammerer & Nashed, 1972), i.e. an operator such that

⟨Th, g ⟩ = ⟨h, T ∗ g ⟩

(42)

for all h ∈ Uˆ and g ∈ Y . The reader can refer to Balakrishnan (1976) for the definition of the inner product used in (42). The key idea of this paper is to represent uˆ t as a linear combination of basis functions (15). Substituting (15) into (40) we define a new operator Tϕ as follows (Tϕ λ)(t) =

t

∫

K (t , τ )ϕ (τ )dτ · λ

t ∈ I.

(43)

0

The new domain of the linear operator Tϕ is Rnb since the basis functions ϕ (t) are known. Hence, (Tϕ ) : Rnb → L2 ([0, +∞))q .

(44)

The following proposition shows how to compute the adjoint matrix of Tϕ . Proposition 12. The linear bounded operator Tϕ : Rnb → L2 ([0, +∞))q defined in (43), admits as adjoint operator Tϕ∗ : L2 ([0, +∞))q → Rnb the following operator (Tϕ∗ es )(t) =

+∞

∫

τ

∫

ϕ T (θ )K T (τ , θ )dθ es (τ )dτ

(45)

0

0

The operator Tϕ∗ is a bounded operator over the interval I = [0, +∞). Proof. Using the property (42) of the adjoint operator

⟨Tϕ λ, es ⟩ = ⟨λ, Tϕ∗ es ⟩

(46)

we compute the first term of the above equation

⟨Tϕ λ, es ⟩ =

+∞

∫

eTs (τ )

τ

∫

0

K (τ , θ )ϕ (θ )λdθ dτ

0

using Fubini’s theorem (Beard, 2002; Kress, 2014; Priestley, 1997) we get +∞

∫

τ) ∫ =

eTs ( 0

τ

∫

K (τ , θ )ϕ (θ )λdθ dτ

0

+∞

+∞

∫

eTs (τ )K (τ , θ )ϕ (θ )dτ dθλ

θ

0

+∞

(∫

+∞

∫

ϕ T (θ )K T (τ , θ )es (τ )dτ dθ

= θ

0

)T

= ⟨λ, T ∗ es ⟩,

λ (47)

nb

∗

where ⟨λ, T y⟩ is the scalar product in R . Hence, the adjoint operator is (Tϕ∗ es )(t) =

∫

+∞

θ

0

∫

+∞

∫

+∞

∫

= 0

ϕ T (θ )K T (τ , θ )es (τ )dτ dθ τ

ϕ T (θ )K T (τ , θ )dθ es (τ )dτ .

(48)

0

By recalling the definition of K (τ , θ ) ≜ CeA(τ −θ ) B and the fact that A is Hurwitz, and by remembering that the basis functions ϕ (θ ) are bounded and belong to L2 ([0, +∞]), we can conclude that Tϕ∗ is a bounded operator. Proof of Theorem 5. As mentioned, the least square solution of the integral equation (39) with uˆ t represented by the linear combination of the basis function ϕ (t) (15) must satisfy the following normal equation (Tϕ∗ Tϕ )λ = Tϕ∗ es . For any es , the above equation always admits at least one solution (Kammerer & Nashed, 1972). Writing the above normal equation using the matrix Φ and

the vector Ψ , λ is the solution of the linear system Φ λ = Ψ which is always solvable and the solution satisfies the normal equation Φ ∗ Φ λ = Φ ∗ Ψ . Among the possible solutions of this normal equation, the minimal norm solution can be obtained as λ = Φ † Ψ where Φ † is the usual Moore–Penrose pseudo-inverse of Φ . From Proposition 12, Tϕ∗ is bounded then ∥Ψ ∥2 ≤ L. The linear operator Φ is bounded, then it has a bounded generalized inverse (Kammerer & Nashed, 1972), ∥λ∥2 = ∥Φ † Ψ ∥2 ≤ ∥Φ † ∥2 ∥Ψ ∥2 ≤ N. Since ϕ (t) consists in a finite number of bounded and L2 basis functions, and since λ is bounded, it results that also uˆ t (t) = ϕ (t)λ is bounded and belongs to L2 . Since the input u is uˆ (t) = us (t) + uˆ t (t) and us satisfies the assumption A2, the resulting uˆ is an approximated solution for Problem 1. References Balakrishnan, A. (1976). Applications of mathematics: Applied functional analysis. Springer-Verlag. Beard, R. W. (2002). Linear operator equations with applications in control and signal processing. IEEE Control Systems, 22(2), 69–79. Benosman, M., & Le Vey, G. (2003). Stable inversion of SISO nonminimum phase linear systems through output planning: an experimental application to the one-link flexible manipulator. IEEE Transactions on Control Systems Technology, 11(4), 588–597. Beutler, F. J. (1965). The operator theory of the pseudo-inverse I. Bounded operators. Journal of Mathematical Analysis and Applications, 10(3), 451–470. Bolder, J., & Oomen, T. (2015). Rational basis functions in iterative learning control — With experimental verification on a motion system. IEEE Transactions on Control Systems Technology, 23(2), 722–729. Boor, D. (1978). A practical guide to splines, vol. 27. Springer-Verlag New York. Clayton, G. M., Tien, S., Leang, K. K., Zou, Q., & Devasia, S. (2009). A review of feedforward control approaches in nanopositioning for high-speed SPM. Journal of Dynamic Systems, Measurement, and Control, 131(6), 061101. De Luca, A., Lanari, L., Lucibello, P., Panzieri, S., & Ulivi, G. (1990). Control experiments on a two-link robot with a flexible forearm. In Proceedings of the 29th IEEE conference on decision and control, 1990 (pp. 520–527). IEEE. Devasia, S. (1997). Output tracking with nonhyperbolic and near nonhyperbolic internal dynamics. Journal of Guidance, Control and Dynamics, 20(3), 573–580. Devasia, S., Chen, D., & Paden, B. (1996). Nonlinear inversion-based output tracking. IEEE Transactions on Automatic Control, 41(7), 930–942. Duan, M., Yoon, D., & Okwudire, C. E. (2017). A limited-preview filtered B-spline approach to tracking control — With application to vibration-induced error compensation of a 3D printer. Mechatronics. Gearhart, W. B., & Schultz, H. (1990). The function sin(x)x. The College Mathematics Journal, 2(2), 90–99. Graichen, K., Hagenmeyer, V., & Zeitz, M. (2005). A new approach to inversionbased feedforward control design for nonlinear systems. Automatica, 41(12), 2033–2041. Hirschorn, R. M. (1979). Invertibility of nonlinear control systems. SIAM Journal on Control and Optimization, 17(2), 289–297. Horowitz, I. M. (2013). Synthesis of feedback systems. Elsevier. Hunt, L., Meyer, G., & Su, R. (1996). Noncausal inverses for linear systems. IEEE Transactions on Automatic Control, 41(4), 608–611. Hussein, M. T., & Nemah, M. N. (2015). Control of a two-link (rigid-flexible) manipulator. In 2015 3rd RSI international conference on robotics and mechatronics, ICROM (pp. 720–724). IEEE. Jetto, L., Orsini, V., & Romagnoli, R. (2013). Almost perfect tracking through mixed numerical-analytical stable pseudo-inversion of non minimum phase plants. In 52nd IEEE conference on decision and control (pp. 1453–1460). IEEE. Jetto, L., Orsini, V., & Romagnoli, R. (2014). Accurate output tracking for nonminimum phase nonhyperbolic and near nonhyperbolic systems. European Journal of Control, 20(6), 292–300. Jetto, L., Orsini, V., & Romagnoli, R. (2015a). A mixed numerical — analytical stable pseudo-inversion method aimed at attaining an almost exact tracking. International Journal of Robust and Nonlinear Control, 25(6), 809–823. Jetto, L., Orsini, V., & Romagnoli, R. (2015b). Optimal transient performance under output set-point reset. International Journal of Robust and Nonlinear Control. Jetto, L., Orsini, V., & Romagnoli, R. (2015c). Spline based pseudo-inversion of sampled data non-minimum phase systems for an almost exact output tracking. Asian Journal of Control, 17(5), 1866–1879. Jetto, L., Orsini, V., & Romagnoli, R. (2017). B-splines and pseudo-inversion as tools for handling saturation constraints in the optimal set-point regulation. In American control conference (ACC), 2017 (pp. 1041–1048). IEEE. Kammerer, W., & Nashed, M. (1972). Iterative methods for best approximate solutions of linear integral equations of the first and second kinds. Journal of Mathematical Analysis and Applications, 40(3), 547–573.

R. Romagnoli and E. Garone / Automatica 101 (2019) 182–189 Kress, R. (2014). Linear integral equations, vol. 82. Springer-Verlag. Moylan, P. (1977). Stable inversion of linear systems. IEEE Transactions on Automatic Control, 22(1), 74–78. Piazzi, A., & Visioli, A. (2001). Optimal inversion-based control for the set-point regulation of nonminimum-phase uncertain scalar systems. IEEE Transactions on Automatic Control, 46(10), 1654–1659. Priestley, H. A. (1997). Introduction to integration. Oxford University Press. Ramani, K. S., Duan, M., Okwudire, C. E., & Ulsoy, A. G. (2017). Tracking control of linear time-invariant nonminimum phase systems using filtered basis functions. Journal of Dynamic Systems, Measurement, and Control, 139(1), 011001. Rigney, B. P., Pao, L. Y., & Lawrence, D. A. (2009). Nonminimum phase dynamic inversion for settle time applications. IEEE Transactions on Control Systems Technology, 17(5), 989–1005. Romagnoli, R., & Garone, E. (2018). Steady-state input calculation for achieving a desired steady-state output of a linear systems. arXiv preprint arXiv:1804. 11185. Seifried, R. (2014). Dynamics of underactuated multibody systems. Springer. Silverman, L. (1969). Inversion of multivariable linear systems. IEEE Transactions on Automatic Control, 14(3), 270–276. ToolboxMatlab, (2018). General least square inversion toolbox. https://approxmod elinversion.wordpress.com/ . Vakil, M., Fotouhi, R., & Nikiforuk, P. (2009). Piece-wise causal inversion by output redefinition for a flexible link manipulator. Transactions of the Canadian Society for Mechanical Engineering, 33(2), 57–78. van Zundert, J., & Oomen, T. (2017). On inversion-based approaches for feedforward and ILC. Mechatronics. Wang, H., Kim, K., & Zou, Q. (2013). B-spline-decomposition-based output tracking with preview for nonminimum-phase linear systems. Automatica, 49(5), 1295–1303. Zhang, L., & Liu, S. (2014). Basis function based adaptive iterative learning control for flexible manipulator. In The 11th world congress on intelligent control and automation (pp. 828–833). IEEE. Zou, Q., & Devasia, S. (1999). Preview-based stable-inversion for output tracking of linear systems. Journal of Dynamic Systems, Measurement, and Control, 121(4), 625–630.

189

Zou, Q., & Devasia, S. (2007). Precision preview-based stable-inversion for nonlinear nonminimum-phase systems: The VTOL example. Automatica, 43(1), 117–127.

Romagnoli Raffaele completed the Master Engineering degree in Industrial Automation Engineering from the Universitá Politecnica delle Marche (UNIVPM) — Italy in 2011. In 2015 he received the Ph.D. in Control System and Automation specialized in Optimal and Robust Control System Theory at the same institute. In the same year he is a postdoctoral researcher at the Department of Control Engineering and System Analysis (SAAS) at the Universitè libre de Bruxelles (ULB)-Belgium. Currently he is a postdoctoral researcher at the Electrical and Computer Engineering of the Carnegie Mellon University, PA, USA. He is currently working on secure control of cyberphysical systems and the control of systems subject to software rejuvenation. Other research interests are: model stable inversion, output tracking problems, reference governor, set-theoretic control, optimal and robust control, control of Li-ion battery cell, space applications in micro-gravity conditions.

Emanuele Garone received the master’s (Laurea) and Ph.D. degrees in Systems Engineering from the University of Calabria, Rende, Italy, in 2005 and 2009, respectively. Since 2010, he has been an Associate Professor in the Control Design and System Analysis Department, Universitè Libre de Bruxelles, Brussels, Belgium. He has authored more than 100 papers in peer-reviewed conferences and journals. His interests focus on constrained control, reference governors, networked control systems, and aerial robotics.