APPROXIMATE INVERSION OF MIMO SYSTEMS FOR PRECISION TRACKING

The problem considered is the design of a digital control system for precision tracking control of a fully-coupled non-minimum phase MIMO plant. The first step is to design a 2-degree-of-freedom feedback tracking system using standard pole-placement or linear quadratic regulator techniques. The result is a stable closed-loop system having zero steady-state error to step inputs. In order to obtain precision tracking for other types of inputs, some kind of feedforward control is needed. Two different tracking architectures are considered in this thesis, both feature a feedforward inverse filter to enable precision tracking. The filters are the inverse of a closed-loop system. The derivations of the filters are given; for non-minimum phase systems the result is approximately a decoupled system of delays over a certain bandwidth. Two different possibilities to design the inverse filters are considered: the first one is based on a novel frequency approximation, while the other one relies on the addition of feedback. Several options to design the feedback are discussed. The tracking performances of the resulting precision tracking architectures are demonstrated and evaluated for several non-minimum phase example systems via simulations.


Introduction
Tracking is often used to achieve zero steady-state error for control systems.
This thesis, however, focuses on the transient response of control systems, and does not only care about the steady-state error. For this reason, the approach is called precision tracking because it is designed to precisely follow the reference trajectory during all time. In order to reach this goal, some kind of system inversion is required.
If the system under consideration is non-minimum phase, the realization of an exact inversion of the system is not used since this scenario leads to unstable eigenvalues in the inverted system. Therefore, approximations are needed. In this thesis, two architectures for the approximate inversion of discrete-time linear systems will be considered. Even if the continuous-time system is minimum phase, the corresponding discrete-time system can be non-minimum phase, due to the presence of sampling zeros (for example there will always be zeros outside the unit circle if the relative degree of the continuous-time system is greater than or equal to three, and the sampling period is sufficiently short [2]).
Classic feedback control systems with integral control may face the problem of having only a small bandwidth for precision tracking, since they only guarantee zero steady-state error. This work will address this problem by using and designing a feedforward filter to improve the precision tracking bandwidth of the control system.
Precision tracking is important in many applications such as atomic force microscopes (AFM) [3,4,5], scanning tunneling microscopes (STM) [6], piezoelectricstack actuated nanopositioning platforms [7,8] and hypersonic vehicles [9]. Two of the applications will be briefly explained below to demonstrate the importance of precision tracking.
The AFM is a versatile instrument which is able to image nanoscale structures with an increasing importance in molecular biology [10]. It is particularly interesting for control engineers since the imaging depends utterly on the feedback control loop. Advantages of the AFM over optical methods comprise a higher resolution and the fact that the AFM produces a 3D surface map. The applications include material sciences (photovoltaic cells, crystallization, semiconductor properties) and the scan of biologically relevant materials (membrane stability, cell motility) [10].
The AFM can also be used to manipulate material, making it a useful actuator for nanotechnology. The quality and speed of AFM images is dependent on the overall dynamics of the AFM system [4], hence a good tracking of the reference trajectory is desired.
The STM is an important tool in nanofabrication, but its problem is that it cannot compete with more established techniques due to its limited operating speed (throughput) [6]. The throughput, however, is limited because of positioning errors in the STM system, and therefore a high-speed precision positioning system is required.

Overview of Thesis
The thesis is structured as follows: Chapter 2 gives an overview of the existing literature.
Next, in Chapter 3, the considered approaches to invert a given system are explained: two different ways to derive a stable feedforward filter which approximately inverts the system are presented. Starting point is the exact, unstable feedforward filter. The first design possibility relies on considering a certain number of advances for the system outputs to come up with a stable filter, while the other technique focuses on designing feedback to stabilize the filter. In the context of these inversion techniques, two different precision tracking architectures will be presented.
In Chapter 4, several alternatives to calculate the mentioned feedback controller are presented. Among these alternatives, especially an idea to formulate an optimization problem to obtain an ideal feedback gain matrix which optimizes the system's frequency response excels.
Chapter 5 introduces various example systems which are to be inverted for precision tracking and provides simulation results. Finally, Chapter 6 gives the conclusion.

Review of Literature
Inversion of non-minimum phase systems has been studied previously by various authors. First, some basic thoughts on the inversion of non-minimum phase systems are presented. In the following sections, a few approaches and ideas shall be discussed, together with a (personal) evaluation of the advantages and disadvantages between these approaches and the method derived in this thesis.

Different kinds of inversion approaches and architectures can be distinguished
and exist in the literature. In general, inversion approaches can be split up into closed-loop inversion feedforward (CLIF) and plant-inversion feedforward (PIF) architecures, as shown in Fig 2.1 [11,12]. The inverse filter is denoted by F, G c is a feedback controller and G p a plant model. The goal of the architectures is to track the desired output y d with the actual plant output y.
Moreover, different approaches to find the inverse filter can be distinguished.
On the one hand, various stable approximate model-inversion techniques exist [11,13] (e.g. Tomizuka's zero phase error tracking control [14]), especially for discrete-time, single-input, single-output (SISO) plants. These approaches try to exactly invert the system model, but replace the unstable part of the zero dynamics with a stable approximation [13]. On the other hand, techniques have been published which use the exact, unstable inversion and ensure stability of the system by using noncausal plant inputs [13] (e.g. in [15,6,16]). Some of the noncausal approaches use some kind of preview time to make the inversion more applicable to applications where the desired output trajectory is not completely pre-specified [15].
An aspect that many inversion approaches (both continuous-and discrete-  time) have in common is that they specifically exclude systems with a transmission zero on the imaginary axis or on the unit circle, respectively (e.g. [15,6,17,16,18]).
The inverse filter that will be derived in Chapter 3 is an approximate inverse for a discret-time non-minimum phase system which uses an advanced version of the reference input (i.e. advanced by a finite number of samples) to drive the inverse filter.

Continuous-time approaches
One continuous-time approach can be found in [15]. It calculates the feedforward input that leads to precision output tracking. The method will be discussed in more detail here.
The starting point is the square state-space systeṁ y(t) = Cx(t) with x(t) ∈ R n and u(t), y(t) ∈ R p . First, the relation where ξ d (t) is the desired ξ(t), which is known when the desired output and its time derivatives are specified.
Finally, the internal dynamics η(t) have to be calculated. If the system is hyperbolic (i.e. none of the zeros of the system is located on the imaginary axis), the internal dynamics can be decoupled with a transformation matrix U into a stable (σ s (t)) and an unstable (σ u (t)) subsysteṁ . In this hyperbolic case, the solutions for the internal dynamics are given by (2.9) i.e. the internal dynamics η(t) can be calculated (if the system is non-hyperbolic, the author developed an idea to deal with this case in [19]). It is important to notice that the desired output must be completely specified (including future information) in order to compute the solution for the unstable part σ u (t). For online computation, it is assumed that y d is known for a preview time of T p seconds, i.e. Y d (τ ) is known for all t ≤ τ ≤ t + T p , and is the approximate solution for the internal dynamics. It can be shown that it is possible to make σ u (t) −σ u (t) arbitrarily small by having a large enough preview time T p . The proposed method was tested on an experimental flexible structure consisting of two discs, which are connected by a thin freely rotating shaft. The input of the system is the voltage applied to a DC motor, and the output is the angular rotation of the second disc. The method was able to track the desired output signal, and increasing the preview time improved the tracking performance.
It was also possible to specify the trajectory online.
A similar approach by the same authors is presented in [6]. The starting point is the same as above, i.e. a square continuous-time state-space system, with the Laplace-domain representation given by y(s) = C (sI − A) −1 Bu(s) = G(s)u(s). (2.11) The optimal inversion problem is formulated as an optimization problem. The goal is to minimize the cost function with u opt (jω) = G opt (jω)y d (jω). (2.14) The problem with this filter is that it tends to be unstable if the system under consideration is non-minimum phase. Due to the non-causality of the filter, a preview-based implementation approach for the optimal inverse filter is developed.
First, the filter has to be rewritten u opt (s) = G opt (s)y d (s) =Ĝ opt (s)ŷ d (s) (2.15) such thatĜ opt (s) is proper andŷ d (s) is the Laplace transformation of a linear combination of the desired output and its time derivatives. Next, the filter is decoupled into a stable (G s opt (s)) and an unstable (G u opt (s)) part G opt (s) = G s opt (s) + G u opt (s) (2.16) by partial fraction expansion. Let the state-space representations be given bẏ x s (t) = A s x s (t) + B sŷ d (t) u s opt (t) = C s x s (t) + D sŷ d (t) (2.17) andẋ u u opt (t) = C u x u (t) + D uŷ d (t). (2.18) If the desired output and its time derivativesŷ d (t) are bounded in time, the bounded solution to the optimal-inverse input u opt (t) is found as u opt (t) = u s opt (t) + u u opt (t).
(2. 19) The computation of the optimal-inverse input, at any time t, requires the knowledge of all future values of the desired output trajectory y d (t). The input can be approximated by using a finite preview time T p , so that the values are known within the time interval [t, t + T p ]. The approximation for the unstable part becomes is the result for the finite preview based optimal-inverse input. It is shown that the tracking error can be made arbitrarily small by choosing a sufficiently large preview time. A rule of thumb is introduced as well, the preview time T p should be "greater than four times the time constant of of the dominant unstable pole" of the optimal inversion filter G opt (s). This method was applied to the STM system, and compared to a "dc-gain approach", where the desired trajectory is scaled by the dc gain of the STM model. The optimal-inverse can greatly improve the tracking performance of the STM system in comparison to the dc-gain approach, when using a large enough preview time. An insufficient preview time, however, can lead to substantial tracking errors. An advantage of the method considered in this thesis in comparison to this continuous-time approach may be that neither in [15] nor in [6] the author says anything about the digital implementation of the approach, e.g. the sampling rate necessary to perform the integral calculations with sufficient accuracy.
Another continuous-time approach can be found in [9]. It deals with the control of an unstable, non-minimum phase hypersonic vehicle model. The basic idea of the approach is to add a stabilizing term to a standard dynamic inversion method to move the unstable zero into the left-half plane. The proposed method is valid for multiple-input, multiple-output (MIMO) systems. One disadvantage of the approach is that it involves a transformation to Jordan form, which can be difficult in practice, especially when the matrix which is to be converted possesses multiple eigenvalues.
In [20], an approximate-inverse method for SISO systems (that can also include a term G 0 (s) = e −sT 0 ) is presented, which uses no preview time and a causal inverse control law. For this method, the desired output trajectory must be known at any time instant t ≥ 0. The perfect inversion control law for a plant G(s) with m transmission zeros is found by partitioning the plant: The structure can be seen in Fig. 2.2, and the ideal feedforward inversion control law follows as where the coefficients c k , k = 1, . . . , n ν are determined by the coefficients in G 10 (s) (n ν depends on the system order and the order of the Taylor approximation of G 0 (s) = e −sT 0 ) [20]. The problem is that (2.23) cannot be used, since signals y (k) 1d are not available for a designer. The paper focuses on finding good approximations for this control law. An advantage of this approach is that it specifically includes , which is not the case for the method in Chapter 3, but the author does not generalize the approach to MIMO systems.

Discrete-time approaches
The theoretical problem of inverting non-minimum phase discrete-time MIMO systems is solved in [16]. It is shown, however, that the perfect inversion requires an infinite number of "preaction" samples. Almost perfect tracking can be achieved if a preview of the reference trajectory is available that is "significantly greater than the maximum time constant associated to the inverses of the controlled system invariant zeros".
One of the standard approaches for discrete-time non-minimum phase SISO systems, the zero phase error tracking control (ZPETC), was presented in [14].
The method designs a feedforward controller for a closed loop system, where a combination of pole/zero cancellation and phase cancellation (applied to uncancellable zeros) is used in order to ensure that the desired and the actual output are in phase for all frequencies. In [1], a filter is proposed to ensure unity gain, instead of zero phase. The ZPETC could be successfully applied to motion control of a robot arm, excelling both in terms of tracking error and smoothness of velocity.
Unfortunately, the original ZPETC, as proposed in [14], is sensitive to modeling errors and plant uncertainties, whereas the method from Chapter 3 explicitly deals with stability robustness. Therefore, an adaptive ZPETC was presented in [21].
In [17], a trajectory tracking control for non-minimum phase SISO systems is introduced by factorizing the system into a minimum phase and a zero phase system. By inverting both systems individually, a feedforward controller can be constructed. For the zero phase system, a discrete zero-phase FIR filter was de- signed. The method was tested successfully on the mock model of a XY gantry stage.
A different approach for perfect tracking for discrete-time systems was proposed in [22]. The authors eliminate the unstable zero problem in the design of the inverse discrete-time system by using a multirate feedforward controller. A highly robust performance is achieved by the controller. Advantages of the method in comparison to ZPETC could be shown through simulation and experiments, in the context of a position control using a dc servomotor.
An approach for a filter design using optimization can be found in [1]. It is called model matching and the basic idea is to find a filter F which minimizes the worst case frequency domain gain deviation, i.e.
where G is the non-minimum phase plant and W a weight chosen based on the frequencies contained in the desired output. A choice of W = 1 for all frequencies would cause difficulties, since the minimization will try to approach F = G −1 .
Instead, W should be chosen close to unity for the frequency range of interest. An idea from the same authors to improve robustness can be found in [1,23]. The first step is to design an inverse dynamic filter F for a feedback control system (the papers propose several approaches for that filter, e.g. ZPETC or model matching). The authors' concern, however, derives from an unsatisfactory filter performance due to a mismatch between the model and the physical system. Therefore, they propose an iterative refinement approach to update the input of the control system based on the output tracking error. An initial input sequence is given, which is called the nominal command input. It derives from the filter F .
A gradient descent approach was used to develop an algorithm for the update process. The algorithm takes the input of the proceeding loop and updates it in the following way: Assume a SISO system with its transfer function G (from input u to output y), where the feedback-loop is already included. Given is a desired output sequence The starting point of the algorithm is to set u = u * , and it works as follows: 1. Apply u to physical system and obtain output sequence y: 3. Iterate until y − y * or ∆u becomes sufficiently small The architecture can be seen in Fig. 2.3. For the update process, not the ideal desired output trajectory y d is used, since the use of it turned out to be too aggressive and leading to saturation. Instead, filtering is applied to y d (e.g y * = Gu * could be used).
For this algorithm, a complete run of the control system is needed in order to make one update. For a real-time implementation, it is proposed by the authors to use the results from the iterative refinement to train a FIR filter mapping y d to the corrective input ∆u: with n = n 1 + n 2 as the order of the filter and n 1 as the look-ahead horizon. The filter coefficients can be obtained through a least square fit to the data obtained from the iterative refinement. The method was tested on a high speed positioning system. The system model was obtained by identification. The fixed-coefficient FIR filter was able to improve the tracking performance. The tracking error was significantly and uniformly reduced in comparison to the case where just the inverse dynamic filter was used. An advantage of the proposed inversion method in Chapter 3 could be that, as mentioned earlier, the stability robustness is taken into account while designing the entire inversion-based tracking system, without running an algorithm to update the ideal inversion input.

CHAPTER 3 Stable MIMO System Inversion
In the following chapter, the tracking architectures that were examined and derived in the course of this thesis will be presented. They were first introduced in [24,25]. The design of digital control systems for precision tracking will be discussed. The first step is to design a feedback control system using classic poleplacement or linear-quadratic regulator techniques, together with integral control, so a stable closed-loop system having zero steady-state error to step inputs is guaranteed. In order to have precision tracking for other types of inputs, different kinds of feedforward filters for the closed-loop system will be derived, based on a novel approximation. If the closed-loop system is non-minimum phase, the filter and the closed-loop system are approximately a decoupled system of delays over a certain bandwidth. Moreover, an alternative feedback design technique to stabilize the inverse filter is considered and presented.

Feedback Tracking System
The first step is to design a 2-degree-of-freedom feedback control system in order to achieve a desired bandwidth with good stability robustness.
Consider an nth-order linear, m-input, p-output plant with state-space model (3.1) For digital control with sampling time T seconds the zero-order hold (ZOH) equivalent plant model is In order to compute the feedback gain K 1 and the integrator gain matrix K 2 in Fig. 3.1, a design model is needed. The equations (3.6) are used to form the system (3.7) With the definitions  [27]), which finds the feedback gain matrix that maximizes a combination of the input-multiplicative (δ 1 ) and the input-feedback (δ 2 ) stability robustness bounds.
The algorithm described in [27] returns the optimal δ 1 , but it was slightly modified in this thesis so that it optimizes a combination of both bounds (a weighted measure for δ 2 was added to the cost function).
The two robustness bounds will be shortly introduced. For δ 1 , an inputmultiplicative plant pertubation model, shown in Fig 3.2a, is considered [28]. Let the system from w[k] to v[k], when all external inputs are set to zero, be denoted by H 1 (z). The result is shown is Fig. 3.2b. Note that H 1 (z) is stable, since it has the same poles as the nominal closed-loop control system that is stable by design [28]. If the plant is a p-input system, then ∆ 1 is a p-input, p-output system which is assumed to be stable. It represents a perturbation to the nominal plant model.
The small gain theorem is used to examine the stability of the system shown in Fig. 3.2. It gives a sufficient condition for stability for the feedback interconnection of two stable systems, as shown in Fig. 3.2b. It states states that the closed-loop system is guaranteed to be stable if holds, where · ∞ denotes the system infinity norm [28]. This inequality may be rewritten to obtain a bound on the size of the unknown pertubation system: The right-hand size of this inequality is defined as the input-multiplicative stability robustness bound δ 1 [28]. For the feedback tracking system in Fig. 3.1, the system H 1 (z) is given as the following state-space system (note that w[k] refers to the notation in Fig. 3.2 here, and not to Fig. 3.1): (3.11) With this system description, δ 1 can be computed as the reciprocal of the system infinity norm.
The input-feedback stability robustness bound (δ 2 ) is defined in a similar way.
The difference to δ 1 is that a input-feedback plant perturbation model, shown in Fig. 3.3, is considered [28]. The procedure to determine δ 2 is the same, i.e. find Pertubed Plant Model (a) Input-muliplicative plant pertubarion model (b) Closed-loop system for input-muliplicative plant pertubarion model when all external inputs are set to zero Pertubed Plant Model (3.12) The system H 2 (z) for the feedback tracking system in Fig. 3.1 is given by: (3.13) The larger the bound δ 1 is, the more tolerant the control system is to errors in the plant model, the same goes for δ 2 . As a rule of thumb, a control system with δ 1 > 0.5 is desired [28].
Recommendations for the pole locations of feedback tracking systems can also be found in [27]. Another possibility for the feedback design is to use the discretetime linear quadratic regulator (dlqr ) formulas.
In the next sections, the derivations of filters is given which, by cascading them with the closed-loop system in Fig. 3.1, increase the precision-tracking bandwidth.

Algorithm for Stable System Inversion
In this section, the general procedure to invert a stable digital system, with input w[k] and output y[k], is presented [24,25]. Let the n c th order, p input and p output, system be given by its state-space representation The goal of the following filter design is to invert the model of the system (3.14).
Consider r advances of the plant output, y[k + r]. Let r be the smallest integer for which C c Φ r−1 c Γ c is a nonzero matrix. It follows that holds. This result can be easily shown: . . .
with m < r. Define and rearrange (3.15) to obtain the ideal inversion control law w[k]. Substituting this equation into (3.14a) yields Equations (3.19) and (3.18) show that a system where D f is given by (3.17) and where holds may be used as a feedforward filter to invert the system (3.14), since it produces the signal w[k], which inverts the system model (3.14) with a delay of r samples, from the advanced plant output y[k + r]. That is, the cascade of the filter (3.20) and the system (3.14) is a pure delay of r samples on each output signal.
Thus, the filter in (3.20) is theoretically able to achieve perfect tracking with a delay of r samples.
However, the eigenvalues of the filter will include the zeros of the system under consideration (see Section 3.2.1). Therefore, if the system (3.14) is non-minimum phase (either because the continuous-time system is already non-minimum phase or due to the presence of unstable "sampling zeros"), the filter in (3.20) cannot be used as a feedforward filter, since it is an unstable system. Additional calculations are needed to obtain a stable approximate inverse of the system.

Eigenvalues of the Inverse Filter
This section shows why the system (3.20), which inverts the system (3.14) with a delay of r samples, cannot be used as a feedforward filter for non-minimum phase systems (see also [9]). As mentioned before, this is due to the fact that it contains the zeros of system (3.14) as eigenvalues.
First, we will make an assumptions about system (3.14): we assume that it is minimal, because then the transmission zeros of the system coincide with its invariant zeros [29].
The transmission zeros of a system with transmission matrix G(z) are defined as the values η i for which holds [30]. For systems where the number of inputs is equal to the number of outputs, this condition reduces to It is noted that other definitions for the transmission zeros exist in the literature, they can also be defined with the help of the Smith-McMillan form [31,32]. Particularly in [32] it is mentioned that (3.23) "can not in general be used to find or define the zeros and poles of a square matrix G(s)". This is due to the fact that in det (G (s)) = α z(s) p(s) (3.24) the polynomials z(s) and p(s) "are not necessarily relatively prime" [32]. In this thesis, however, it is assumed that the definition for the transmission zeros in (3.23) and the definition with the Smith-McMillan form coincide because we consider only minimal state-space models.
The invariant zeros are defined with the help of the Rosenbrock matrix P(z) [30,31,32]. The Rosenbrock matrix for a system (A, B, C, D) is defined as The invariant zeros are the values η i for which the Rosenbrock matrix is rankdeficient: rank (P (η i )) < max z rank (P (z)) . Again, for a system with the same number of inputs and outputs, this reduces to det (P (η i )) = 0. (3.27) Let the transmission matrix of system (3.14) be denoted by G(z): Therefore, its transmission zeros can be found according to (3.23). We now consider the same system, but we assume that the output is y[k + r], and not y[k]. The transfer function of this system is given by where Z{·} denotes the z-transform and the shifting property of the z-transform was used (see for example [26,33]). Thus, the transmission zeros of system (3.14), where the output is advanced by r samples, can be found by det G(z) = det (z r G(z)) = z p·r det (G(z)) = 0. (3.30) For this computation, the fact that for a square n × n matrix A det (λA) = λ n det (A) (3.31) holds, was used [26]. It follows from (3.30) that the system with the advanced output has the same transmission zeros as the original system, but additionally it has p · r transmission zeros in 0.
The state-space representation of the system, when the output is advanced by r samples, is given by (3.32) Due to our assumption, the transmission zeros of the system coincide with its invariant zeros, so the Rosenbrock matrix can be used as well to find the transmission zeros. The Rosenbrock matrix for the system (3.32) is defined as and for every invariant zero η i of system (3.32), this matrix is rank-deficient. Each invariant zero η i is associated with an invariant-zero direction z 0 = z x0 z w0 which lies in the kernel or null space of P (η i ) [32]: From the definition of the zero direction, the two equations follow. Equation (3.35b) can be rearranged where the definitions (3.18) and (3.21a) have been used. This result is then plugged into (3.35a): so it becomes clear that the invariant zeros η i of system (3.32) are the eigenvalues of , but this is also the state-transition matrix of the inverse filter. Together with the previous result from this section, it can be concluded that the eigenvalues of the inverse filter in (3.20) consist of the invariant zeros of system (3.14) (this includes possible non-minimum phase zeros) and p · r eigenvalues at 0. Consequently, more steps have to be taken in order to stabilize the inverse filter (3.20), which will be discussed in the following section.

Adding Additional Advances
Now, the stabilization of the filter (3.20) by adding more advances is shown.
Consider adding s more advances to the r advances in (3.15). Let d = r + s. Then it follows: With the approximation which will hold true for low frequency signals, equation (3.38) becomes Now the same procedure as before is used, i.e. solving holds. The advantage of adding more delays is that the new filter (3.42) will be stable for a sufficiently large d. The state-transition matrix is Φ f = Φ c + Γ c C f and the entries of C f = −D f C c Φ d c can be made arbitrarily small by increasing d, so the eigenvalues of Φ f will move towards the eigenvalues of Φ c , from which we assume that it is stable.
This result can be verified using the Bauer/Fike Theorem [34]. Let A ∈ C n×n be a diagonalizable matrix with eigenvalues λ i , i = 1, . . . , n, so follows. Since any matrix norm induced by a vector norm is submultiplicative [34,35], we can further write, together with (3.42a): where K 2 has been defined as K 2 = Γ c 2 C c 2 . The 2-norm of a matrix can also be expressed as its maximum singular value [34], and it can be shown that holds [36]. Here, σ min (·) denotes the minimum singular value of a matrix. In the light of this result, it becomes clear that has to be shown, where the definition in (3.41) was used. Since we assume that Φ c is stable, we have |λ i | < 1, i = 1, . . . , n c for every eigenvalue of Φ c , thus Φ d c converges to the zero matrix as d tends towards infinity [26]. This means that the numerator of (3.48) tends towards zero. It is further well known that σ min > 0 holds for a full-rank matrix [26]. Because D f is defined as the inverse of the full rank (because otherwise the entire method is not applicable), and therefore the denominator in (3.48) is greater than zero, so Φ d c 2 D f 2 does converge to zero as d tends towards infinity.
From (3.46), this implies that Using (3.42c), this results leads to the recognition that, as d approaches infinity, the eigenvalues of Φ f and Φ c coincide. Due to the assumption that Φ c is stable, it follows that the inverse filter will eventually be stabilized when more advances d are added.
Consequently, the value of d is chosen to be smallest integer for which the state-transition matrix (3.42c) of the inverse filter is stabilized. That is where R is a user-defined maximum pole radius. If the assumption that Φ c is stable is fair, the choice of R can be based on the eigenvalues of Φ c . It is recommended to chose R to be exactly in between the largest absolute value of the eigenvalues of Φ c and 1. That is, if λ max denotes the eigenvalue with the largest absolute value, R is calculated according to

Extension to Feedback Approach
In this section, an alternative approach to stabilize the inverse filter (3.20), which inverts the system (3.14) with a delay of r samples, is presented. In contrast to the algorithm shown in Section 3.2.2, this new approach does not necessarily use more advances, but uses feedback to stabilize the filter. The idea for this approach was taken from [19].
The goal is still to find an inverse filter for the system (3.14), and the first steps are exactly the same as presented in Section 3.2, i.e. consider r advances of the output and rearrange (3.15) to obtain the ideal inversion control law (3.18).
Next, a possibility to stabilize (3.20) is discussed. Similar to [19], an additive term v[k] is included in the optimal inversion control law (3.18) is obtained as the "inverse" system. This system, however, is not the exact inverse anymore, since the ideal inversion control law (3.18)  may be introduced to move the eigenvalues of the inverse system into the unit circle. Due to the feedback (3.54), the approximate inverse filter becomes From (3.56d) it becomes evident that F must be chosen such that the matrix Φ f +Γ c F is stable. This can be done by any state-space controller design technique (e.g. any pole placement algorithm).
The filter (3.55) may be used as a feedforward filter to invert system (3.14).
Due to the addition of v[k] to the ideal inversion control law in (3.52), the cascade of the filter and the system is not a pure delay of r samples. The precision tracking bandwidth has to be evaluated with the help of the Bode plot of the precision tracking system.
shall be shortly demonstrated. Consider the state-transition matrix of the inverse filter in (3.56d). With the definiton in (3.21c), it can be rewritten as The eigenvalues of the inverse filter are therefore determined by the eigenvalues of the matrix It becomes clear now that the feedback gain matrixF can place the eigenvalues arbitrarily if the system (Φ c , Γ c ) is controllable. Even if not every eigenvalue of (Φ c , Γ c ) is controllable, the feedbackF can stabilize the inverse filter if the noncontrollable eigenvalues of Φ c are stable, which we know is true since we assume Once the matrixF has been found, it is easy to obtain the feedback gain matrix F: The matrix C f , which is needed for this computation, only depends on the matrices of system (3.14), i.e. it is known at any point. The rewriting of the inverse filter state-transition matrix Φ f was only carried out here to demonstrate that F is able to stabilize Φ f over Γ c ; the design of F in this thesis will always be based on the The idea of adding delays from Section 3.2.2 and the feedback approach presented in this section may also be combined. Assume that s max denotes the number of advances needed in the algorithm described in Section 3.2.2 to stabilize the inverse filter. Now consider that s ≤ s max advances are added to the r advances shown in (3.15). Again, let d = r + s, so (3.38) is obtained and the approxima- is added to the resulting control law, so follows. The resulting inverse filter is still given by (3.55) and (3.56), only the definitions of C f and D f in (3.21a) and (3.17) change to (3.62) The remaining definitions Φ f = Φ c + Γ c C f and Γ f = Γ c D f stay the same, but the new definitions from (3.62) have to be plugged in. Note that the design of F is the eigenvalues of Φ f will be different. The resulting inverse filter approximately inverts system (3.14) with a delay of d = r + s ≤ d max samples.
The advantage of combining the two methods is that fewer advances than in Section 3.2.2 may be needed, and that the additive term v[k] may be kept "smaller" (and thus does not "disturb" the ideal inversion control law (3.18) as much) if the feedback has to stabilize a filter that inverts the system with d = r + s delays, instead of r delays.

Command Shaping Filter
First, a filter that will be called the command shaping filter (CSF) is presented.
It follows the inversion approach shown in Section 3.2. The system (3.14), which will be considered for the inversion, is the feedback tracking system in Fig. 3.1.
Let the closed-loop system from w[k] to y[k] be denoted by so for the CSF, the system considered in Section 3.2 is given by (3.63), and the notation in (3.63) coincides with the notation in (3.14).
Next, a short argument shall be given that the system to be inverted contains the zeros of the discrete-time plant model (3.2). In order to do so, it is assumed that the order of the additional dynamics is equal to the number of plant outputs, n a = p, and that both Γ c and K 2 are not rank deficient, i.e. they are p×p matrices with a nonzero determinant. The invariant zeros of system (3.63) are found via its Rosenbrock matrix P c (z) (see also Section 3.2.1): It is known that so that we can as well find the values η i for which P c (η i ) is rank deficient in order to find the invariant zeros of (3.63). Next, the rule is used [36], and since B = 0 in our case and it is assumed that Γ a is a square matrix with full rank follows. This equation can be further rewritten and due to our assumption that K 2 and Γ a have full rank, we finally obtain with det (Γ a ) = 0 and det (K 2 ) = 0. The matrix P(z), however, is the Rosenbrock matrix of the discrete-time plant (3.2). Thus the invariant zeros of (3.63) and (3.2) coincide if the order of the additional dynamics match the number of plant outputs and both Γ a and K 2 have full rank (Γ a has full rank if integral additional dynamics are used). If these conditions are not met, the argumentation has to be modified.
Since the CSF is the inverse of a closed-loop system, it belongs to the class of closed-loop-inversion feedforward (CLIF) architectures [11,12], and we know derived by using either of the algorithms shown in Section 3.2. If the feedback approach from Section 3.2.3 is applied, it should be checked whether the system (3.63) is controllable, but it is assumed that this should be the case if no pole/zero cancellations occur during the design process of the feedback tracking system. If they do occur and the system is not controllable, a model reduction and inverting the reduced system may solve the problem. It is noted that, according to Section 3.2, the feedback is always able to stabilize the inverse filter; however, it may not be possible to place its eigenvalues arbitrarily.
An advantage of the CSF is that the feedback design in Section 3.1 determines the eigenvalues of the state-transition matrix Φ c , which is inverted in Section 3.2.
Thus, it is ensured that Φ c is stable. Additionally, if all eigenvalues of Φ c are distinct (which can be achieved by design), the matrix T which diagonalizes Φ c is known to have full rank [37], so its condition number in (3.45) will be finite.
A precision tracking system is obtained by inserting a feedforward CSF into

Inverse Modified Plant
The second inverse filter that is considered in this thesis is the so-called inverse modified plant (IMP). The difference between the IMP and the CSF presented in Section 3.3 is that not the model of the entire feedback tracking system in Fig. 3.1 is inverted, but the model for the system from v[k] to y[k], which will be called the modified plant. It is given as follows: (3.73) Thus, for the IMP, the system (3.14), which is inverted in Section 3.2, corresponds to (3.73). The system (3.73) has the same invariant zeros as the discrete-time plant (3.2). This can be shown with its Rosenbrock matrix: where P(z) denotes the Rosenbrock matrix of the plant (3.2).
possesses the desired pole locations, so the feedback design is not based on the system (Φ, Γ)). For that reason, the IMP is classified as PIF. However, since the algorithm in Section 3.2 assumes that the system to be inverted is stable, it should be checked whether the eigenvalues of Φ c = Φ−ΓK 1 are all located inside the unit circle or not. If this is not the case, the feedback matrix K d may be recalculated (e.g. by selecting new pole locations or using a different design technique), or the usage of the CSF from Section 3.3 may be considered.

Inclusion of a State Observer
The preceding precision tracking architectures assumed that all the state variables of the plant are available in the design of the feedback tracking system. In practice, however, this will not always be the case. Therefore, this section focuses on the question whether and how the inclusion of a state observer will affect the CSF and IMP precision tracking systems in Sections 3.3 and 3.4, respectively.
Consider system (3.2). A state observer provides an estimatex[k] for the actual state vector x[k] [26]: (3.76) follows. A feedback tracking system which includes a plant, additional dynamics and state observer, can therefore be written as  (3.77) Since the state-transition matrix of the feedback tracking system in (3.77) is a block diagonal matrix, the eigenvalues of the tracking system are determined by the eigenvalues of the matrices of the observer dynamics can be found in [27].
When an observer is used in the feedback tracking system in Fig. 3.1, the system (Φ c , Γ c , C c ) that is inverted in Section 3.2 has to be adjusted for both CSF and IMP. For the CSF, the system from w[k] to y[k] is inverted, which is given by (3.78) The IMP, in contrast, inverts the system from v[k] to y[k]: (3.79) Furthermore, the systems (3.11) and (3.13), which are used to calculate the robustness bounds δ 1 and δ 2 in Section 3.1, have to be adjusted if an observer is added to the feedback tracking system. System H 1 (z) is now given by 

CHAPTER 4 Feedback Inversion Design Techniques
This chapter deals with possibilities to design the feedback matrix F, which was introduced in Section 3.2.3 to stabilize the inverse filter. Clearly, F has to move the eigenvalues of Φ f = Φ f +Γ c F inside the unit circle, but there are multiple ways how this can be done. Since there are no strict rules for the exact pole locations of the inverse filter, this degree of freedom can be used to try out different design approaches.
A first idea is to use "classical" feedback design techniques, such as poleplacement or dlqr. For pole-placement, the inverse filter's pole location could be based on the influence of pole and zero locations of a dynamic system on its frequency response, and the weights for a dlqr design on the demand that the control effort deriving from the addition of the feedback is low. Thoughts on these design techniques are presented in Section 4.1.
As another step, a design was developed in this thesis which focuses on finding a matrix F which stabilizes the system while having a minimal norm. It is assumed that the influence of the feedback on the ideal inversion control law is reduced if the norm of F is minimal. This design is explicated in Section 4.2.
Finally, inspired by the optimization in Section 4.2, another optimization design is introduced in Section 4.3. The main idea is to design F so that it "compensates" the deviation of the tracking system's frequency response from the ideal frequency response (due to the approximation in (3.39)), while a stable inverse filter is guaranteed. From all design techniques that were considered in order to calculate F, this proved to be the most successful, which is why it will be considered as the main design method for the inversion via feedback in the discussion in Chapter 5.

Classical Feedback Design
The easiest way to design the feedback matrix F is to choose the desired locations for the filter with the feedback included. For the discussion of the "closedloop" poles, we assume that the state-transition matrix Φ f is given as where Φ f is defined according to (3.21c), i.e. Φ f is the state-transition matrix of the (unstable) filter which inverts system (3.14) with a pure delay of r samples (see Section 3.2.3 for details). Consequently, the eigenvalues of Φ f are moved with the help of F to obtain the eigenvalues of Φ f . Assume that system (3.14) has n z zeros.
Then, n z eigenvalues of Φ f coincide with the zeros of (3.14). The remaining n c −n z eigenvalues are located in zero. Since the idea of the feedback approach is to simply stabilize the filter, without interfering too much with the ideal inversion control law (3.18), it is proposed here to keep all the stable eigenvalues in their respective spot and just move the unstable eigenvalues to obtain a stable filter. Sometimes, problems might arise if the attempt is made to place multiple eigenvalues in the same spot (e.g. the n c − n z eigenvalues in zero). In this thesis, the eigenvalues were not placed exactly in zero, but n c − n z equally spaced eigenvalues at radius Another possibility is to reflect the unstable eigenvalues into the unit circle, so that the magnitude is maintained in the Bode plot. Assume that Φ f has n z unstable eigenvalues, that are given as the roots of the polynomial Then, the roots of the polynomial should be chosen as the eigenvalues of Φ f [13].
If the system under consideration is a MIMO system, it is proposed to select the eigenvalues according to this section, and to simultaneously try to minimize the ) in the inversion control law (3.52). A reasonable approach is to minimize a norm F , which was achieved in this thesis by adapating the algorithm described in [27].
Another possibility to calculate F is via a dlqr design. For a state-space system with n variables and p inputs, a dlqr design minimizes the following quadratic function of the states and inputs: where Q ∈ R n×n and S ∈ R p×p are symmetric, positive-definite weighting matrices.
If S is large with respect to Q, the resulting regulator will stabilize the plant without using much control effort [26]. In the light of this discussion, a dlqr design for F may be used with S being much larger than Q. As a reminder: the term was added to the inversion control law (3.18) to stabilize the inverse filter, so it can be argued that the tracking performance will be best if the control effort v[k] is as small as possible, while it is still large enough to stabilize the filter.

Norm Minimization Approach
One of the feedback design approaches that is considered and developed in this thesis is to stabilize the inverse filter and simultaneously attempt to minimize the Frobenius norm of the resulting gain matrix F. The idea behind this is that if the norm of F is minimal, the influence of the additive term v[k] = Fx c [k] (see Section 3.2.3) will probably also be limited.
In order to design the matrix F, an idea from [9] is used. The author designs feedback to stabilize the inverse of a continuous-time hypersonic vehicle plant, which has one non-minimum phase zero. The basic idea is generalized here to the general case of multiple (and possibly complex) non-minimum phase zeros for a digital system.
The result of the design process is a set of linear equality constraints and nonlinear inequality constraints, which could be solved "by hand" and, by themselves, do not ensure a minimum norm. An advantage of this design is that an analytical relation between the eigenvalues of the stabilized inverse filter and the entries of F is obtained. In this thesis, they were included in an optimization problem which minimizes the norm of the feedback gain matrix. One of the (theoretic) advantages of minimizing the norm is the Lagrange dualism could be applied to the problem, and the Lagrange dual problem is known to be convex [38].
This section is organized as follows: first, the constraints for the entries of the gain matrix F are derived. After this, they are included in an optimization problem, which is solved using the Lagrange dualism. A discussion of the approach is given in the last subsection.

Constraints for the Feedback Gain Matrix
Following the idea in [9], it is desired that the matrix F moves the unstable eigenvalues of the inverse filter (resulting from the plant's non-minimum phase zeros) into the stability region, while it keeps the remaining (stable) eigenvalues in their respective spot. In order to do so, system (3.53) is transformed to canonical modal form. A state transformation is the result. The transformation matrix V is the right eigenvector matrix of is replace by the real part of the jth eigenvector, and the (j + 1)th column by the imaginary part of the jth eigenvalue. One realization aspect that was not discussed in [9], however, is that Φ f possesses a multiple eigenvalue in 0 (with multiplicity p · r, see Section 3.2.1), and usually there are less than p · r linearly independent eigenvectors associated with the eigenvalue in 0 (i.e. V would be rank-deficient). For that reason, generalzed eigenvectors have to be used [37]. Let can be calculated according to [37]: The usage of the generalized eigenvectors results in Jordan blocks for the multiple eigenvalue in 0, i.e. in ones above the main diagonal.
With the definitons where Λ is in canonical modal form. Next, the partitioningŝ are introduced. Furthermore, it is assumed (in this section) that the first n z eigenvalues of Λ are located outside the unit circle (without loss of generality).
These are the eigenvalues that correspond to the non-minimum phase zeros of (3.14). It is also assumed that the first 2n i eigenvalues within the first n z unstable eigenvalues are complex (i.e. there are n i pairs of of conjugate-complex unstable eigenvalues). The conjugate complex eigenvalues of Λ are denoted by δ j ± ω j and the real eigenvalues by λ j . With these definitons and assumptions, we have for the first 2n i modal states (which are affected by the complex unstable eigenvalues of Φ f ) and for the next n z − 2n i modal states (which are affected by the real-valued unstable eigenvalues of Φ f ). The remaining n c − n z modal states are only affected by the stable eigenvalues of Φ f , which we do not wish to move.
Following the idea in [9], we want the modal states associated with the unstable eigenvalues (i.e. the first n z modal states) to be only excited by them- are made, so that the last n c − n z modal states do not affect the first n z modal states. Next, it is desired that holds, so that the modal states associated with the complex eigenvalues are only excited by themselves. Similarly, we want γ T ψ,jf l = 0 j = 2n i + 1, 2n i + 2, . . . , n z l = 1, 2, . . . , n z , l = j (4. 15) for the modal states associated with the real unstable eigenvalues. If we achieve (4.14) and (4.15) and use (4.13), we have for j = 1, 3, . . . , 2n i − 1, and for the first n z states. Note that the first n z states of the closed-loop system Λ − Γ ψF are in canonical modal form as well, so that the first n z columns ofF allow us to place the corresponding eigenvalues. The complex eigenvalues of the closed-loop system have to be conjugate-complex again, so we want (from (4.16)) for j = 1, 3, . . . , 2n i − 1. This means that for the ith conjugate-complex eigenvalue pair (i = 1, 2, . . . , n i ) in modal form we need a block Moreover, we want the closed-loop system to be stable, so we demand and Finally, it has to be examined how the addition ofF has influenced the remaining n c − n z eigenvalues (with (4.13) we only ensured that the last n c − n z modal states do not affect the first n z modal states; the first n z modal states, however, do have an influence on the last n c − n z modal states). In order to do so, we first state that we can partition where Λ 1 and Λ 2 are matrices in canonical modal form containing the unstable and stable eigenvalues of Φ f , respectively (if generalized eigenvectors have been used, Λ 2 additionaly has ones above the main diagonal). Together with the other results from this section, we can write for the transformed system witĥ This summarizes the whole approach. BecauseΛ 2 = 0 in (4.24e), the eigenvalues of the closed-loop system are the eigenvalues ofΛ 1 andΛ 4 (Λ is a block tridiagonal matrix). We can place the eigenvalues inΛ 1 by the choice of the first n z columns of F, so we can make sure that we move the unstable eigenvalues (which derive from the non-minimum phase zeros) into the unit circle. Moreover, because of (4.24g), the remaining n c − n z eigenvalues are the stable eigenvalues of Φ f (hence we have not changed them by introducing v[k]).
In conclusion, by introducing feedback of the form v[k] = Fx c [k] to the ideal control law (3.18) and by designing the feedback gain matrix according to this section, we are able to move all unstable eigenvalues of the filter into the unit circle and keep the already stable eigenvalues in their respective spot.

Design by Optimization
The idea of this section is to include the equality demands (4.14), (4.15), (4.18) and the inequality demands (4.21), (4.20) (which ensure that the the inverse filter is stabilized) for the design of the feedback gain matrix F into an optimization problem so that the norm of the resulting F is minimized to limit its influence on the ideal inversion control law. Therefore, an optimization problem was formulated to design F and which includes the stability demands (4.21), (4.20) as inequality constraints. Further, the demands (4.14), (4.15), (4.18) are included as equality constraints.
The gain matrix F is obtained after designingF (according to (4.8e)): As design parameters, we have the p·n z entries of the first n z columns ofF (because we set the last n c − n z columns ofF to zero in (4.13)), so we define as the decision variable (please note that x in this case does not refer to any state variable). Thus we define as the optimization problem, together with the inequality constraints (4.21) and (4.20) and the equality constraints (4.14), (4.15) and (4.18). The squared norm of F is minimized because it simplifies the calculations later on.
This optimization problem with nonlinear constraints can be either solved with Matlab's fmincon function, or some further analysis can be done to solve it with the Lagrange dualism, which is what will be done here.
The Lagrange dualism is very briefly introduced here, see for example [38] for details. Consider an optimization problem in standard form with decision variable x ∈ R n , and it is assumed that its domain D is nonempty (please note that the notation here is independent of the usual notation in this thesis: n, m and p in (4.28) do not refer to the plant order or to the number of plant inputs or outputs). By defining the Lagrangian L : R n × R m × R p → R the constraints are taken into consideration: (4.30) We refer to µ i as the Lagrange multiplier associated with the i-th inequality constraint f i (x) ≤ 0; similarly we refer to ν i as the Lagrange multiplier associated with the i-th equality constraint h i (x) = 0. Furthermore, we define the Lagrange dual function g : R m × R p → R as the minimum value of the Lagrangian over x: The dual function yields lower bounds on the optimal value p * of the primal problem (4.28): so the question is what the best lower bound is that can be obtained by the parameters µ, ν. This leads to the optimization problem maximize g(µ, nu) subjet to µ ≥ 0 where µ ≥ 0 means that all elements of the vector µ should be ≥ 0. This problem is called the Lagrange dual problem associated with the primal problem (4.28). It is a convex optimization problem because the objective to be maximized over is concave and the constraint is convex. This is the case whether the primal problem (4.28) is convex or not [38].
Next, we want to find an expression for the Lagrange dual function for our problem at hand. The details of the calculations and the definitions of the following matrices are described in Appendix A. It is first noted that the primal problem (4.27) can be expressed as: where the details can be found in Appendix A. After some calculations (see Appendix A for a detailed discussion) the Lagrangian can be written as The notation Q (µ) denotes that the matrix Q is not constant, but depends on the Lagrangian multiplier µ. It can be shown that Q is positive definite. Thus, the dual function becomes and a minimum is obtained when (4.36) is differentiated and the derivative is set to zero, since Q is positive definite. It follows: where Q = Q + Q T was defined. Plugging this result back into (4.36) yields for the dual function, where the two facts were used [34]. Thus we have as the (convex) dual optimization problem. Since we can express any maximization problem as a minimization problem so we finally obtain as the dual minimization problem that we want to solve using Matlab's fmincon function. We have to use the fmincon function as we still have to consider a constrained optimization problem.
Once the optimization problem (4.43) is solved, and we have obtained µ * and ν * as its solutions, we can conclude by using (4.37), hence the entries ofF (and therefore F as well) are known.  vanishes and this approach yields better results than before, but they are still not a real improvement to the classic design of adding advances until the inverse filter is stabilized, and the derivation and implementation are much more complicated.

Discussion of the Approach
It is concluded that the success of this approach does not only depend on the non-minimum phase zeros of the plant, but also on how "well" the inverse filter can be transformed to canonical modal form. Furthermore, the chosen objective function (i.e. the Frobenius norm of the resulting controller F) may not be ideal to obtain a "good" precision tracking system. Despite this disappointing conclusion, this approach played an important role in this thesis, since the work on this initial idea to involve optimization in the design process of F led to more thoughts on how optimization can improve tracking and how a corresponding optimization problem has to be formulated. This resulted in the approach presented in Section 4.3.

Frequency Optimization Approach
This section focuses on the idea of optimizing the frequency response of the cascade of the inverse filter and the system whose output is desired to be tracked.
In theory, the filter-based tracking system should be a pure delay of d samples (d ≥ r), i.e. the ideal frequency response of the cascade is given by Let H c e jωT denote the frequency response of the system, and H f e jωT the frequency response of the inverse filter with feedback matrix, F (defined in (3.56)).
Assume that a grid of frequency points ω k , k = 1, . . . , N is given, based on the frequencies included in the system. In the simulation results presented in Chapter 5, ω N was chosen to be the maximum frequency for which none of the magnitudes of the main diagonal elements of the feedback tracking system shown in Fig. 3.1 has dropped by 3 dB, and ω 1 = 0 rad s (dc gain).
According to (4.45), an initial idea for the optimization could be where · F denotes the Frobenius norm. As a quick reminder, Φ f denotes the unstable state-transition matrix of the inverse filter that approximately inverts the system with a delay of d samples, and to which the feedback gain matrix F has not yet been added (inversion is exact for d = r). This initial optimization problem, even though it is a formulation of what we want to achieve, can prove to be quite wasteful in practice. The reason for that is that the optimization, for a MIMO system, attempts to optimize the off-diagonal elements of H c e jωT H f e jωT as well (i.e. zero magnitude and phase over all frequencies), but this is not necessarily required for a good inversion-based tracking performance. The phase of the offdiagonal elements is not important for the tracking performance if the attenuation is small enough, so the demand for the off-diagonal elements is to "simply" provide a sufficient attenuation over the frequency range of interest.
The following sections deal with different formulation for the objective function of the optimization problem that were considered in this thesis. For better readability, a few definitions are introduced first. Define z k = e jω k T . Introduce and let the elements of G (z k ) be denoted by g ij (z k ), i, j = 1, . . . , p: (4.48)

Main Diagonal Elements
The main diagonal elements are the main concern of the optimization. For a good tracking performance, the frequency responses g ii (z k ), i = 1, . . . , p have to be "close" to the ideal frequency response g ideal (z k ) = e −jωdT . Therefore, the following cost function is introduced for the main diagonal elements where g * ii (z k ) denotes the conjugate-complex of g ii (z k ). It optimizes the distance between the actual frequency response g ii (z k ) and the desired value in the complex plant over the frequencies of interest.

Off-Diagonal Elements
Two approaches were considered for the off-diagonal elements. Since we only care about the attenuation, the first cost function tries to minimize the squared magnitudes of the off-diagonal elements: The factor α can be used to weight the cost function for the off-diagonal elements in comparison to the cost function for the main diagonal elements. Naturally, the weighting factor could also be different for every k, i and j, so a factor α kij would be possible, if weighting for different frequencies and elements is desired.
Another possibility is to define a desired attenuation for the off-diagonal elements. First, let the element g ij (z k ) be given by For i = j, an attenuation of is demanded. It can bee seen that g ij (z k ) has to lie within a circle with radius R a in the complex plane, centered in the origin. An idea for a penalty function (see [39]) is given by The argument of the exponential function was chosen such that it is negative if g ij (z k ) lies within the desired circle, and positive if it is not. Together with an adequate choice of the weighting factor α, g ij (z k ) has a very high contribution to the penalty function if it is outside the desired circle (i.e. if the attenuation at this frequency is higher than desired), and it has basically no contribution if the attenuation demand is satisfied. It goes without saying that α and R a could be chosen differently for every k, i and j, if so desired.

Stability Constraint
As indicated in (4.46), the optimization has to be constrained. Namely, the feedback gain matrix F has to stabilize Φ f , i.e. the eigenvalues of the matrix Φ f + Γ c F must be all located inside the unit circle. Usually, constraint functions f (x) must be defined such that they return a negative value if the constraint is satisfied and a positive value otherwise. In our case, one possibility for the constraint function is to consider the spectral radius ρ of the matrix It is defined as the largest absolute value of the eigenvalues of the matrix [34].
Hence, the constraint function could be formulated as Alternatively, a constraint function for every eigenvalue of the filter could be defined: where λ i denotes the ith eigenvalue of Φ f + Γ c F and n c is the order of the filter.
Another possibility is to not explicitly constrain the optimization, but to include the stability constraint in the objective function. In order to do so, a penalty function is designed, similar to Section 4.3.2. The desired area in the complex plane, in which the eigenvalues must be kept in, is the unit circle. If is a reasonable penalty function for the stability constraint. As before, the eigenvalue λ i has a very high contribution if it is outside the unit circle, and hardly any contribution otherwise. To ensure a successful design, the parameter β should be assigned a pretty high value (e.g. β = 10 4 ).

Discussion of the Approach
In theory, this frequency optimization approach can be used in an attempt to get a sufficient tracking performance within the frequency range of interest, while only using the minimum number of r advances to invert the system, i.e. the cascade of the filter and the closed-loop system is an (approximate) delay of r samples. For most cases, however, the optimization will not be able to find a F which stabilizes the filter while also maintaining a respectable tracking performance over the desired frequency range, if only r advances are used. That being the case, it is recommended to use a "hybrid approach" between adding delays and designing the feedback matrix F (via optimization) to stabilize the filter, discussed at the end of Section 3.2.3. Accordingly, the system matrices for the inverse filters are given in (3.56), with C f and D f being defined in (3.62). A reasonable procedure to find the optimal inverse filter design is to let s run from 0 to s max and use the number of advances s that produced the smallest value of the objective function.
As an initial value F 0 , a result from the classical feedback design techniques, presented in Section 4.1, can be used. In particular, it proved to be useful to determine F 0 via a dlqr design.

CHAPTER 5 Example Systems and Simulation Results
In this chapter, the example systems that were examined in the course of this thesis will be presented, together with the tracking results obtained from simulations in Matlab. The desired reference trajectory will be denoted by r[k], the actual output by y[k].
The feedback matrix K d = K 1 −K 2 for the tracking system in Fig. 3.1 was designed according to the guidelines in [27]. These rules use normalized Bessel poles, which can be found in Table 5.1 [25,27].  architectures are able to perfectly track this reference. Let A ij (jω) = |g ij e jωT | be the magnitude and ϕ ij (jω) = ∠g ij e jωT be the phase of the frequency response of the transfer function from the jth input to the ith output of the entire filterbased tracking system (shown in Fig. 3.5 and Fig. 3.6, respectively). Assume that we choose a pure sinusoid with frequency ω as the reference signal for every output, i = 1, . . . , p. The difference between the ideal output and the actual output will be used to formulate a condition for the precision tracking bandwidth.
As a first step, we derive the output of the precision tracking system for the chosen reference input. With the presented notation, we obtain as the steady-state response of the ith ouput [26]. According to Sections 3. A ij (jω) sin (ωkT + ϕ ij (jω)) . (5.7) follows. Obviously, δ[k] is a (discrete) function of time, t = kT , but we want to derive a measure for the tracking error at a specific frequency ω, independent of the time t, which we will call (jω). For this reason, we define as the performance measure for every output. Next, Euler's formula e jx = cos(x) + j sin(x) (5.9) will be used to find an expression for (jω) [40] 1 . This is done by rewriting (5.7) (for ease of notation and better readability kT is replaced by t and the ω-dependency of the magnitudes and the phase shifts is omitted): Im{A ij e j(ωt+ϕ ij ) } = Im e j(ωt+ϕ 0 ) − p j=1 A ij e j(ωt+ϕ ij ) = Im e jωt e jϕ 0 − p j=1 A ij e jϕ ij :=Â i e jφ i = Im{Â i e jϕ i e jωt } = Im{Â i e j(ωt+ϕ i ) } =Â i sin(ωt +φ i ). (5.10) In the light of this observation, it can be concluded that δ i [k] is a sinusoid with amplitudeÂ i and phaseφ i , so i (jω) is defined as the amplitude of this sinusoid.
Finally, we define the precision tracking bandwidth ω b to be the largest value ω for which i (jω) < K b holds for every output, i = 1, . . . , p. For the following simulation results we chose K b = 10 −2 , i.e.
The precision tracking bandwidth will for example be used to evaluate how much of the reference trajectory's energy is included in this bandwidth. In order to do so, we introduce (with Ω k = ω k T ) . Therefore, the amount of the signal's energy that is contained up to the frequency ω k ≤ π T is calculated as in the following sections.
Moreover, two other performance measures were established in the thesis and will be taken into account for the example systems in the following sections. These performance measures were taken from [4], where the tracking of continuous-time SISO systems is discussed. One performance measure (J e ) is introduced to quantify the energy in the tracking error, another one (J m ) to quantify the peak deviation from the desired trajectory. Assume that y i [k] denotes the ith output and r i [k] the corresponding desired trajectory. Then, we can define and as the respective performance measures for one output, where w i [k] is a weighting function defined by and r i,t > 0 denotes a threshold value for the ith output. In the following, it is set to i.e. 10% of the maximum absolute value of the reference trajectory for the ith output. If the reference for the ith output is r i [k] = 0 for all k, it is set to r i,t = 1.
The weighting function was introduced so that the tracking error y i [k] − r i [k] is normalized and becomes comparable for different systems. However, if the reference input is zero or close to zero, the performance measures would become very large if they were divided by r i [k], and therefore the threshold r i,t was introduced.
Matlab's trapz function was used to carry out the numerical integration in (5.17).
For MIMO systems, the performance measures for the single outputs will be combined, so that and will be used as the performance measures in the following sections.
Another measure that will be discussed is the (absolute) tracking error . This is not to be confused with the stability robustness bounds δ 1 and δ 2 introduced in Section 3.1, as they will also be mentioned in the discussions and evaluations of the example systems.
The "standard" approach to calculate the CSF and IMP in the following sections will be the design presented in Section 3.2.2, i.e. advances will be added until all the eigenvalues of the inverse filter are contained in the chosen pole radius R. This result will then be contrasted and compared with a design that involves a feedback gain matrix F in order to stabilize the inverse filter. Most of the time a hybrid approach between adding advances and designing a feedback matrix will yield the best result.
Please note that all the Bode plots shown in the following sections are for the closed-loop system shown in Fig. 3.5 and Fig. 3.6, respectively. Thus, they are not indicative of gain or phase margin of the control system.

H-Frame System
In this section, an H-Frame XY positioning system, which consists of two stationary motors, eight pulleys and a single drive belt, is considered. A detailed description of the system modeling and an 8th-order state-space model can be found in [42]. (see (5.1)). According to the rules in [27], damping was added to the complex plant poles s 1 to s 4 , and the remaining six poles were chosen to be s 6 T S (see Table  5.1). These poles were mapped using the ZOH pole-mapping formula, λ i = e s i T , and the algorithm from [27] was used to calculate the (discrete-time) feedback gain matrix K d . The attained input-multiplicative and input-feedback stability robustness bounds were δ 1 = 0.6490 and δ 2 = 0.7395, respectively. After K d is calculated, the (digital) feedback tracking system (Fig. 3.1) can be formed. The

Command Shaping Filter
The CSF can invert the closed-loop system in          is available for controlling the system (see Section 3.5 for details). The desired observer pole locations were calculated via a dlqr design on the system Φ T , C T , with weighting matrices Q = 100 · I n and S = I p (since the continuous-time plant does not have any zeros, the guidelines in [27] were not applied). The locations returned by this design were then used in the algorithm from [27], in order to find As a next step, to validate the robustness of the tracking architecture, it is assumed that the system description is subject to uncertainties. The inverse filter design was based on the nominal plant model, but when the simulation was carried out, every parameter of the H-Frame system (see [42] for a description of the parameters) was changed by 10%. The same observer was used as in the previous result. As to be expected, adding uncertainties substantially reduced the precision tracking bandwidth, we now have ω b = 3.1145 rad s . It has to be evaluated via a simulation if this bandwidth is sufficient to track the desired reference without visible tracking errors. The plot of the tracking performance and the corresponding tracking errors is shown in Fig. 5.6. Especially for the second output, the tracking error gets worse, as deviations between reference input r 2 [k] and plant output  The results for an uncertain plant in combination with an observer are comparable to the CSF results.
All the results can be seen in Table 5.2.

Atomic Force Microscope
The model for an Atomic Force Microscope (AFM) was taken from [4], where also details on the AFM can be found. As a versatile instrument, the AFM is able to image nanoscale structures, and it is of particular interest for control engineers since the imaging depends utterly on the feedback control loop. In this section, a model for motion in the X direction will be considered [4]. It is given as the discrete-time transfer function model which uses a sampling rate of f = 20.833 kHz, so the sampling time is T = 1 f = 0.048 ms.
In order to use the rules given in [27] for the design of the feedback track-    The reference trajectory is a triangular wave, whose fundamental frequency is f 0 = 100 Hz and ranges from −9 µm to 9 µm (because the microscope shall move in a "back-and-forth motion" [4]). For the simulation, the Fourier series of this signal was formed and the first N harmonics were included into the reference trajectory.

IMP FTS
The result for the Fourier series is b k = 3.6 × 10 −5 k 2 π 2 sin k π 2 − sin k 3π 2 For the simulations in the following sections, the maximum frequency included in the reference trajectory r(t) was f max = 1900 Hz.

Command Shaping Filter
The CSF can invert the AFM model exactly with a delay of r = 2 samples.
Since     Ultimately, an uncertain plant is considered, and it is assumed that the statevariables of the plant are not measurable, so that an observer has to be designed. If  Table 5.1) and the observer gain is calculated with place on the system Φ T , C T . The resulting stability robustness bounds are δ 1 = 1 and δ 2 = 0.5372.
After that, the inverse filter is formed with the nominal plant model (

Inverse Modified Plant
In the case of the AFM, it was examined whether the IMP delivers more satisfying results if not only integral additional dynamics are used, but the dynamics of the reference trajectory are included as additional dynamics. It can be seen in (5.30) that the reference trajectory consists of a sum of sine waves with frequency ω k . For the following simulations, the first harmonic ω 1 of the reference was included in the additional dynamics, i.e. the matrix Φ f is ought to have the two (discrete-time) eigenvalues λ 1 = e jω 1 T and λ 2 = e −jω 1 T . In order to guarantee that the design model (3.8) is controllable, the additional dynamics are given in controllable canonical form:

(5.31)
This design did improve the overall performance of the IMP precision tracking system, especially the precision tracking bandwidth.
Therefore, the design model is a 9th-order system, which is designed in the same way as before (i.e. the real parts of the eigenvalues are moved to −2000).
After the feedback tracking system design, the modified plant has the following (digital) poles   Fig. 5.14. The precision tracking bandwidth of ω b = 909.6986 rad s is larger than for the first presented CSF result, but smaller than for the hybrid CSF approach.
In Fig. 5.15 the tracking performance (Fig. 5.15a) and error (Fig. 5.15b) can be seen. Similar to the tracking performance of the CSF, the tracking error has a peak during the first rising edge of the reference trajectory and at the first maximum. The maximum tracking error is |δ max | = 1.3 µm for the IMP, this is an improvement of 28% in comparison to the CSF, and also a slight improvement in comparison to the hybrid CSF approach. After this initial peak, the tracking error stays in between ±3.5 µm, just as for the CSF. While J e has about the same value The best frequency optimization result was obtained when s = 8 additional advances were added, and an unconstrained optimization problem was solved with (4.56) as penalty function for the stability constraint (β = 10 4 , ω N = 1.6897 × 10 3 rad s ) The tracking results for the AFM can be found in Table 5.3.

Comparison to Source
In this section, the tracking performance achieved with the CSF and IMP will be compared to the tracking results in [4]. approaches were considered (e.g H ∞ or 1 control) to design a feedback controller C and a feedforward controller F . Due to the fact that several approaches were considered, multiple tracking performances were presented in [4]. In general, the results presented in [4] have a larger maximal tracking error |δ max |. This large tracking error is not only observed at the first turnaround point, but at every minimum and maximum of the reference trajectory.
The first approach presented is a H ∞ design. During the rising/falling edges of the reference trajectory, a good tracking performance is achieved. At the turnaround points, however, the tracking error is at almost |δ max | = 2 µm. In contrast to this result, the best CSF performance improves the maximum tracking error by 30%; the IMP by 35%. Furthermore, the maximum tracking error for the CSF and the IMP occurred during the first rising edge, the tracking error at the later turnaround points is much smaller for both CSF and IMP: the maximum tracking error for both architectures at the later turnaround points is 0.35 µm, which improves the error by 82.5% in contrast to the H ∞ design.
If C and F are designed via an approximate model-inversion presented in [1,23], the tracking error at the turnaround points could be improved to be at around |δmax| = 1 µm [4], but it achieves a worse tracking performance away from the turnaround points. Nevertheless, CSF and IMP improve the tracking error by 65% at the turnaround points.

Scanning Tunneling Microscope
Models for the x-and y-Dynamics for a Scanning Tunneling Microscope (STM) were taken from [6]. In this application, a piezo scanner moves the STM probe across a sample surface. During this movement, the distance between the  holds for the y-Dynamics (from the input voltage u y in V to the piezo-position p y inÅ) [6]. It is important to mention that the Laplace variable s is in rad ms . In this thesis, the two models for x-and y-Dynamics were combined to one single state-space system (i.e. the overall systems consists of two decoupled subsytems) with input u = u x u y T and output y = p x p y T , so that the algorithms from Section 3.2 can be applied. Naturally, the state-space models contains Together with the additional dynamics, the design model is a 14th-order system. For the design of the feedback tracking system, the desired settling time is set to T S = 0.5 ms. According to the rules in [27], damping was added to the complex poles, and the remaining six real-valued poles (s 5 , s 6 , s 11 , s 12 and the two integrator poles from the additional dynamics) were moved to the 6th order Bessel poles (see Table 5.1). Based on the desired settling time, the sampling time was chosen to be T = 2 µs. The discrete-time modified plant has four non-minimum phase zeros, as well as the feedback tracking system. The tracking system (without inverse filter) has a precision tracking bandwidth of ω b = 0.1 rad ms . The reference trajectory causes the STM probe to move in a raster pattern [6].
First, the probe is moved from the center to the top-left point of the image area.
During the forward (left-right) scan, the y-position is fixed, and it is incremented while the x-position is returned back to the left. This procedure is repeated until the entire desired are is scanned. In this case, the scan rate is s r = 1 Tr = 250 Hz, i.e. T r is the time to complete one back-and-forth motion. This is illustrated in  The tracking error is always smaller than |δ max | = 0.26Å, while the maximum tracking error in the simulation results presented in [6] is at approximately |δ max | = 0.6Å, i.e. the tracking error shrank by 57% in comparison to the result in [6]. It can be stated that the tracking error for the first output, δ 1 , oscillates with a relatively high frequency around 0Å. This also translates to the input voltage u x , so it would be necessary to evaluate on the physical system if this control input was realizable.
The performance measures obtained by the CSF are J e = 0.0162 and J m = 0.0357.

Inverse Modified Plant
The modified plant needs r = 1 delay to be perfectly inverted.  tracking bandwidth is ω b = 2.6906 rad ms , thus the bandwidth was increased by 40% in comparison to the CSF result. However, the maximum absolute tracking error is a bit higher than for CSF tracking, namely |δ max | = 0.34Å. The performance measures for the IMP are J e = 0.0039 and J m = 0.0366. Especially J e is much smaller than for the CSF case, we have J e,CSF J e,IMP ≈ 4. A summary of the obtained results is given in Table 5.4.

Bell 205 Helicopter
In [18,19], the linearized model of a Bell 205 helicopter is considered. The model is a non-minimum-phase near non-hyperbolic system (i.e. the non-minimum phase zeros are close to the imaginary axis), and the model represents the helicopter at a nominal 5 • pitch attitude, with mid-range weight, a mid-position center of gravity, and operating in-ground effect at near sea level [18].
The (A, B, C) continous-time model is given as [19]:   The input vector is u = δ C δ B δ A δ P T where δ C is collective, δ B longitudinal, δ A lateral cyclic and δ P tail rotor collective, the output vector is A hybrid approach between frequency optimization and adding advances has been carried out for the helicopter as well. It proved to be most successful to solve an unconstrained minimization problem with (4.56) as penalty function for the stability constraint and (4.50) as cost function for the off-diagonal elements, where α = 10 and β = 10 4 were used as weights. The best result was obtained when s = 4 additional advances were added to the inversion control law , and then F was designed to minimize the frequency error between the tracking system and the ideal tracking behaviour, so that the cascade of inverse filter and feedback tracking system is an approximate delay of d = 6 samples. After the optimization (with ω N = 17.829 rad s ), the CSF has a spectral radius of Φ f = 0.9998. With this approach, the precision tracking bandwidth is ω b = 0.6889 rad s , which is more than three times higher than the prior precision tracking bandwidth. Now, the closest frequency to ω b is ω = 1.2556 rad s , and 97.86% of the reference trajectory's energy is included within this frequency range.
The absolute tracking error of the second output could be decreased significantly by this approach, however, the third output now shows the worst All the results can be found in Table 5.5.

Comparison to Source
In [19], an adaption for the method from [15] is presented (see Section 2.1 for details). This adaption is introduced in order to be able to deal with nonhyperbolic systems. In the framework presented in [19], non-minimum phase zeros close to the imaginary axis result in a very large preactuation time. Therefore, the major idea presented in [19] is to first introduce feedback to move the nonminimum phase zeros further away from the imaginary axis, and then to apply the stable inversion for non-minimum phase system from [15]. It is a trade-off between precision tracking and required preactuation time.
The trade-off approach in [19] especially suffers from a (comparatively) poor tracking performance for the yaw rate R (output y 4 ), for it has a maximum value of |δ 4,max | ≈ 0.015 rad s . The maximum tracking error for the forward velocity (y 1 ) is |δ 1,max | ≈ 2 × 10 −3 m s . It has to be noticed that the tracking errors are so small that they cannot be distinguished from zero if the approach from [15] is used without moving the near non-hyperbolic zeros first, so this result cannot be compared with the CSF and IMP performance. The apprehension that this results in a large preactuation time, however, proved to be justified, i.e. the practicality of this result is limited.
On the other hand, the tracking errors of the vertical and lateral velocities are not shown in detail in [19], so that a discussion and comparison of the results is not possible.

Multilink Flexible Manipulator
The model for a multilink flexible manipulator was taken from [43]. They are used in fields like assembling of electronic hardware, space exploration or precision welding, but they suffer from vibrating links at high operation speeds, which delays the precise positioning of the end effector [43]. Since conditioning of the For the design of the feedback tracking system, the desired settling time is set to T S = 2.5 s, which leads to a sampling time of T = 6.25 ms. According to the rules in [27], ten plant eigenvalues were kept at their respective spot, damping was added to four of them and the remaining six were chosen as the 6th-order Bessel poles. The feedback gain was calculated with the algorithm in [27], δ 1 = 0.4713 and δ 2 = 0.6574 were attained as the stability robustness bounds. The resulting feedback tracking system has four non-minimum phase zeros on the unit circle.
A step with an amplitude of 20 • lasting for 10 s after which the manipulator goes back to its vertical position for ten more seconds is the reference trajectory.
Since a step would cause infinite joint velocities during the rising and falling edges, leading to a potential mechanical breakdown of the manipulator, the reference is filtered with with the following lowpass-filter [43]: The order of the filter is set to n = 2, and with the adjustable paramter κ the filter roll-off can be changed, which determines the speed of the response [43]. In this case, it was chosen to be κ = 0.2.

Command Shaping Filter
The spectral radius for the feedback tracking system is ρ (Φ c ) = 0.9895 so that the pole radius becomes R = 0.9948. The exact inversion is possible with r = 2 delays, the CSF, however, needs s = 17 additional delays to attain the pole radius R, which makes the filter-based tracking system an approximate delay of d = 19 samples.
With this pole radius, the CSF achieves a precision tracking bandwidth of ω b = 3.2099 rad s , and the tracking error is never larger than |δ max | = 0.49 • . Furthermore, the performances measures J e = 0.001 and J m = 0.2120 were obtained.
The CSF needs a relatively large number of delays, given the fact that the non-minimum phase zeros are located right on the stability border. Therefore, the pole radius was adjusted and set to R = 0.9999, i.e. the CSF is supposed to just For this system, stabilization of the inverse filter with the help of feedback via a dlqr design proved to be quite successful (see Section 4.1). As weighting matrices, Q = 10 −10 · I nc and S = 100 · I p were chosen. Since r = 2 delays are still needed to invert the system, the cascade of CSF and feedback tracking system is an approximate delay of r = 2 samples. With the dlqr design, perfect tracking should be achieved over the whole considered frequency grid, ω b ≥ 500 rad s . The performance measures are much smaller as well, they become J e = 1.4294 × 10 −7 and J m = 2.2870 × 10 −5 .
Clearly, this result is due to the ideal circumstances of the simulation. Therefore, an observer was added to the system and uncertainties were simulated. A dlqr design on the system Φ T , C T was used to calculate the desired observer pole locations, with Q = 10 3 · I n and S = I p as weighting matrices. Then, the algorithm from [27] was applied to obtain the observer gain matrix L with optimized stability robustness. The stability robustness bounds reduced to δ 1 = 0.3981 and δ 2 = 0.5011 in comparison to the non observer-based tracking system. As for the uncertainty, the inverse filter was based on the nominal plant model, but for the simulation, every entry of the plant's state-transition matrix was changed by 0.5%, i.e. Φ = (1 + ∆) Φ and ∆ = 0.005.
As to be expected, adding an observer and uncertainties aggravated the tracking performance. The best result was obtained with a frequency optimization approach, where s = 7 additional delays were added before the design of the feedback gain matrix F, i.e. the tracking system is an approximate delay of d = 9 samples. The cost function (4.50) was used for the off-diagonal elements and (4.56) as penalty function for the stability constraint, with α = 10 and β = 10 4 . The precision tracking bandwidth dropped down to ω b = 0.0439 rad s , the performance measures became J e = 1.5939 and J m = 3.8149. Even though the precision tracking bandwidth reduced significantly, a sufficient tracking performance may still have been achieved for the desired reference trajectory. The tracking errors at the end of the rising and falling edges are way larger than before, |δ max | ≈ 3.5 • , and it has to be evaluated on the physical system if that error is acceptable or not.

Inverse Modified Plant
The spectral radius of the modified plant is ρ (Φ c ) = 0.9888, the corresponding pole radius is R = 0.9944, and the IMP needs s = 14 additional advances. The result achieved with this pole radius is slightly worse than the initial CSF result.
The precision tracking bandwidth is ω b = 0.8475 rad s , while J e = 0.0446 and J m = 0.2820 holds.
Again, changing the pole radius to R = 0.9999 yields better results (only s = 1 additional advance is required), but the best result is obtained with a feedback approach via a dlqr design. As in the CSF case, the entire considered frequency grid belongs to the precision tracking bandwidth (ω b ≥ 500 rad s ), and the performance measures almost completely disappear: J e = 3.9575 × 10 −21 and J m = 9.3296×10 −21 . Then, the same changes were made to the simulation settings as in the previous section, i.e. an observer and plant uncertainties were added. In the IMP case, the best result is achieved when s = 3 additional advances are added before the frequency optimization. The precison tracking bandwidth and the performance measures are comparable to the CSF case.
The results are summarized in Table 5.6.

Comparison to Source
In [43], a continuous-time feedforward filter is proposed, in a closed-loop inversion feedforward architecture (see Fig. 2.1). Similar to the feedback approach in Section 3.2.3, an additive term is added to the inverse control law to move the eigenvalues of the inverse filter into the stability region. Several filter parameters κ were considered (see (5.41)). For high values of κ (e.g. κ = 0.7), the authors detected that the inverse filter could not improve the tracking in comparison to non-filter based tracking. For small values (e.g. κ = 0.2), the output trajectories settled faster to the desired angles (first a step to 20 • , and then return back to 0 • ) and are also closer to the desired trajectory. For these filter parameters, however, the tracking system produces overshoot. For κ = 0.2, the overshoot is around 35% on the worst channel (27 • instead of the desired 20 • ), and when the manipulator is supposed to return to the initial 0 • angles, the worst output "shoots over" up to almost −10 • . In this regard, a clear enhancement can be seen in the simulation results obtained by the CSF and IMP, as no overshooting is noticeable when the nominal plant is used in the simulations. When uncertainties were added to the plant model, the absolute tracking error was still smaller than 10 • . Furthermore, the control effort is comparable for κ = 0.2 for the inverse filter proposed in [43] and the CSF and IMP filters discussed in the previous subsection. Clearly, it has to be mentioned that in [43], actual experimental results are considered while in this thesis, only simulations were carried out.

Overhead Crane
In [44], a linearized model for an overhead crane is presented. It consists of a cart and a load, which is connected to the cart via a rope.
where the input u is the force applied to the cart and the output y is the horizontal position of the suspended load. The parameters of the system are given in Table   5.7, see [44] for more details. For the design of the feedback tracking system, a settling time of T S = 0.5 s is desired, which leads to T = 5 ms as the choice for the sampling time. Taking the guidelines in [27] into account, damping is added to the complex pole pair s 1,2 , and the third order Bessel poles (see Table 5 The reference input should cause the center of gravity of the suspended crane load to move by 300 mm, in a time interval τ = 2 s. In [44], a τ -parameterized transition polynomial is proposed as the desired output function. The normalized reference input is given as and it allows "an arbitrarily smooth transition between 0 and 1", and it can be shown thatr(t; τ ) ∈ C (n) [44]. In this thesis, n was set to n = 2, the resulting trajectory is shown in Fig. 5.19.

Command Shaping Filter
The The best optimization result is attained when (4.56) is used as penalty function for the stability constraint and s = s max = 5 additional advances are considered, with a weight of β = 10 4 and ω N = 5.7309 rad s . This approach, however, does not deliver significantly better results.

Inverse Modified Plant
The modified plant has a spectral radius of ρ (Φ c ) = 0.9717, resulting in  samples, while ω b = 57.7365 rad s holds for the actual precision tracking bandwidth, which is somewhat worse than the CSF result, but this bandwidth still contains the same energy as the CSF bandwidth.
Just as the CSF result, the IMP tracking error has a peak at the beginning of the reference trajectory, with |δ max | = 1.2 × 10 −6 m, and it has another (smaller) peak at the end of the reference with |δ| = 1.2×10 −7 m. The performance measures become J e = 2.8911 × 10 −11 and J m = 1.6488 × 10 −9 .
Again, the best optimization result can be observed if the maximum number of additional advances is considered, s = s max = 4, but no noticeable improvement is obtained.
Finally, all the tracking results are shown in Table 5.8.

Comparison to Source
In [44], a dynamic inversion technique is presented by the authors to design suitable position (and velocity) set-point feedforward signals. In their approach, they apply a continuous-time input-output inversion and are able to get rid of the postaction (postaction means that the invertng signal reaches its steady-state after the transition time), which usually occurs as a problem in motion control inversion problems and interferes with the practicability of these approaches [44].
Besides the desired output trajectory, [44] specifies limitations for the control input u. The cart is actuated with a brushless servomotor through a pulley and toothed belt system, and the relationship between the motor torque T and the force applied to the cart u (the input signal of the plant model) is given by where i is the reduction ratio of an epicycloidal speed reducer, η its mechanical efficiency and r the pulley radius [44]. The maximum continuous motor torque is T max = 3 Nm, and the peak is at T peak = 6 Nm. It is noted, however, that there has to be a typing error in (5.47) (for example because the units are not consistent, u is in Newtons, r in meters and T in Newton-meters, so [rT ] = Nm 2 holds for the unit of the product of torque and radius). It is assumed that at least the reciprocal of r has to be used in (5.47) (the units are consistent then), but it is not known whether the reciprocals of η and i have to be taken as well. Anyway, the motor torque following the CSF and IMP design is larger than T peak , but it cannot be evaluated if this would be the case for an experiment on the real physical system as well, since it is unknown how (5.47) has to be modified to adequately describe the system's dynamics.
A similar behaviour as in the CSF and IMP cases for the tracking error can be seen in the results presented in [44] for the position control, i.e. it has a relatively sharp peak at the beginning of the reference trajectory and converges towards zero in the course of the tracking process. Its maximum value is at around |δ max | = 0.8 mm, which is worse than the CSF results of |δ max | = 4.8 µm. Naturally, the tracking results are hard to compare since they were obtained from an experimental setup in [44], while simulations were used in this thesis.

Two Discs
In [15], a flexible structure consisting of two discs which are connected by a thin freely rotating shaft is considered for tracking. The input is the voltage U (t) applied to a DC motor, while the output is the angular rotation (in degrees) of the second disc, θ 2 (the disc further away from the motor). The plant is given as the This values was chosen because a settling time of T S = 30 s is desired for the feedback tracking system (Fig. 3.1). The feedback tracking system was designed as follows (according to the guidelines in [27]): the second eigenvalue pair s 3,4 is kept, and the remaining poles in the continuous-time domain are chosen as the third order Bessel poles (see Table 5.1). After these pole locations were mapped into the discrete-time domain using the ZOH pole-mapping formula, and place was used to calculate the feedback gain K d , the (digital) eigenvalues of the feedback tracking system are The CSF reaches a precision tracking bandwidth of ω b = 0.6842 rad s . After simulations have been carried out, it can be concluded that this bandwidth is sufficient: the tracking error is never higher than |δ max | = 0.01 • for the desired reference. The performance measures are J e = 1.0266 × 10 −5 and J m = 2.0780 × 10 −6 .
The optimization approach to design a stabilizing feedback matrix F can be used to enhance the tracking performance. When an unconstrained optimization with (4.56) as penalty function (β = 10 4 ) and ω N = 0.143 rad s is used, where F 0 is calculated with a dlqr design (Q = 10 −10 · I nc , S = 100 · I p ), the precision tracking bandwidth becomes ω b = 1.5239 rad s which more than doubles the previous bandwidth. The maximum tracking error shrinks down to |δ max | = 0.002 • , and the performance measures become J e = 1.5010 × 10 −6 and J m = 1.2895 × 10 −7 .

Inverse Modified Plant
The spectral radius of the modifed plant is ρ (Φ c ) = 0.9480, correspondingly the pole radius becomes R = 0.9740. While the modified plant can be inverted exactly with a delay of r = 1 samples, it takes the algorithm s = 2 more advances to move the filter eigenvalues inside the desired pole radius, which is why the filter-based tracking system is an approximate delay of d = 3 samples.
The achieved precision tracking bandwidth is the same as for the initial CSF design, ω b = 0.6842 rad s , though the other performances measures differ from the CSF result. The maximum tracking error is slightly worse (|δ max | = 0.016 • ), but since the maximum value of the reference input is r max = 60 • , this deviation is negligible. The performance measures J e and J m stay roughly the same.
As in the CSF case, slightly better results can be obtained if the frequency optimization approach is used to design a feedback matrix F to stabilize this filter.
For the IMP, however, it proved to be the best choice to add s = 1 additional delay and then design the feedback F, i.e. the hybrid approach between adding delays and optimizing F is used (Section 3.2.3). The initial value is the same as in the CSF case (but naturally, the considered system is different). The most significant refinement in comparison to the standard IMP design is the precision tracking bandwidth, it could be expanded to ω b = 1.0943 rad s (this is faintly worse than the CSF optimization result though).
The results are summarized in Table 5

Comparison to the Source
As mentioned earlier, a real comparison is difficult since the reference trajectory was created on-line in [15], so it is not exactly the same as the one shown in Fig. 5.20. Further, the results presented in [15] are experimental results, and not simulations.
In [15], the internal (or zero) dynamics are decoupled into a stable and unstable part, and bounded solutions are found for both parts . For the unstable dynamics, however, an infinite preview time is required (i.e. the desired output must be completely specified, which the author does not want to assume). Therefore, the solution for the unstable dynamics are approximated using a finite preview time T p (details can also be found in Section 2). Two experiments were carried out and presented, the first one uses a preview time of T p = 20 s, the second one uses T p = 50 s.
The tracking performance for the first preview time is rather bad, as large tracking errors can be observed during the rising edges of the reference trajectory.
For the first rising edge, the output first goes down to θ 2 ≈ −20 • , before it approaches the reference again. This behaviour vanishes for T p = 50 s, but still small tracking erros are present, this time during the falling edges (albeit they are neglectably small in comparison to the other preview time). Even though the preview time in [15] and the advances needed by the inversion algorithm in this thesis are hardly comparable, it can be noted that both the best CSF and MP results only need d = 2 delays (i.e. a "preview" of T p = 0.6 s).

CHAPTER 6 Conclusion
Methods to approximately invert linear non-minimum phase MIMO systems were presented. The success of the inversion techniques was confirmed with simulation results.
Two different tracking architectures were considered. The first one is the command shaping filter (CSF), shown in Fig. 3.5, which inverts a feedback-based tracking system to extend the precision tracking bandwidth. The second one is the inverse modified plant (IMP, see Fig. 3.6), which inverts a modified plant that is part of a feedback tracking system.
First, an approach to add advances in order to design a stable inverse was presented. This makes the cascade of the inverse filter and the system to be inverted an approximate decoupled system of pure delays. The inverse filter design was extended to a feedback approach to stabilize the filter, which led in general to a fewer number of required delays. Several possibilities to calculate the corresponding feedback gain matrix were discussed, among which especially a frequency optimization approach excelled.
An advantage of the IMP is that the additional dynamics of the feedback tracking system can be based on the reference trajectory's dynamics, which showed to be an improvement for some systems. However, this could also be interpreted as a disadvantage, since the CSF always got its best tracking results with "standard" integrator additional dynamics, thus the CSF performance apparently depends less on the chosen reference input and the architecture does not have to be changed if a different signal is desired to be tracked. It can be noted, though, that the IMP tends to achieve the larger precision tracking bandwidth. The CSF, however, has some practical advantages, e.g. the modified plant (which is inverted by the IMP algorithm) is not guaranteed to be stable.
An advantage of the "standard" way to compute the CSF and IMP (i.e. adding advances until the filter is stabilized) is a very easy implementation, and the derivation of this filter is based on a reasonable approximation. A disadvantage may be that the designer has a rather limited influence on the resulting tracking performance, as the only parameter he can change is the pole radius R. The design of the feedback tracking system has an influence on the tracking performance as well, but it is not clear how changes in (Φ a , Γ a ) or K d will translate to the tracking performance. On the other hand, if the frequency optimization approach is used, the designer has some more influence on the inverse filter, e.g. through the choice of the considered frequency grid, which allows him to a certain degree to affect the precision tracking bandwidth. Moreover, it is possible to test several values for the optimization parameters and different formulations for the objective function, as presented in Section 4.3. Clearly, the design of the inverse filter with this approach becomes an iterative task, which is more complicated than the standard design.
As a conclusion, it can be stated that the choice of the tracking architecture and the filter design technique depends on the tracking problem and the system under consideration; in the ideal case both architectures are tested with different design techniques in simulations and the best performance is selected.

Areas of Future Work
A problem that has to be further examined is the proof that the eigenvalues of the inverse filter (3.21c) contain the zeros of the discrete-time plant (3.2). In Section 3.2.1 it was assumed that (3.23) holds for the transmission zeros of a system, but as mentioned in Section 3.2.1, other definitions exist in the literature.
It has to be evaluated when these two different definitions coincide, since it is thought that they do for the class of systems that were considered during this thesis. The results from Section 3.2.1 were confirmed by Matlab computations and simulations, i.e. the inverse filter (3.21c) always contained the zeros of the plant (3.2). Furthermore, an argument has to be given that the invariant zeros Further, the inversion methods may be tested on more example systems and experiments should be carried out to confirm the achieved success in the simulations.

Derivation of the Dual Problem
In this appendix, detailed calculations for the derivation of the Lagrange dual function are presented (see Section 4.2.2).
First, the inequality (4.21), (4.20) and equality (4.14), (4.15), (4.18) constraints that were included in the optimization problem shall be expressed in terms of the decision variable x. As a reminder, x is defined as i.e. the first n z columns ofF. For better readability, the bigger part of these calculations is omitted (because they just involve a lot of indexing) and only the results are presented. Note that γ T ψ,j denotes the jth row of matrix Γ ψ defined in (4.8c), n i the number of complex non-minimum-phase zeros and n z the number of all non-minimum phase zeros (see Section 4.2.1).
First, the constraint (4.14) is considered. It can be expressed as    with (blkdiag refers to the Matlab function of the same name):  As a next step, we consider the inequality constraints. The calculations are presented in more detail here, since the absolute value is a nonlinear function and is therefore harder to find an expression in terms of x. First, we take a look at (4.21). This is the constraint that the real-valued unstable eigenvalues have to be moved into the unit circle. Since (4.21) is the absolute value of a real number, we can rewrite it as − 1 < λ j − γ ψ,jf j < 1 j = 2n i + 1, 2n i + 2, . . . , n z (A.8) which gives us the two demands λ j − γ T ψ,jf j − 1 < 0 −λ j + γ T ψ,jf j − 1 < 0 (A.9) for each j which can be expressed as . . .
Finally, the constraints for the complex eigenvalues (4.20) have to be considered.
Next, we want to introduce the Lagrangian with parameters µ ∈ R 2nz−3n i and ν ∈ R nz(nz−1) (we have 2n z − 3n i inequality constraints and n z (n z − 1) equality constraints in (A.28)). Let µ be partitioned as follows for the Lagrangian. We take a closer look at the sum: (A.32) With the definitions we can simplify the sum to (A. 34) Please note that H depends on the first n i entries of the Lagrangian multiplier µ.
Additionally, we define holds, so S j is an upper triangular matrix with elements ≥ 0 on the main diagonal.
As a consequence, H is an upper triangular matrix with elements ≥ 0 on the main diagonal as well. This leads to the result that Q = Y + H is an upper triangular matrix with elements > 0 on the main diagonal, so it is positive definite.