INVSID 1.0: An Inverse System Identification Toolbox for MATLAB

: The inverse system identification toolbox named INVSID 1.0 for MATLAB, which is used to identify the inversion of single-input single-output systems, is developed. The complete process from theoretical derivation to toolbox creation of developing the toolbox is demonstrated. Afterwards, numerical examples are illustrated to describe how the toolbox can be used to solve inverse identification problems. Simulation results demonstrate the effectiveness of the toolbox.


Introduction
Nowadays, in many motion control systems, there are requirements of high performance, such as short motion times and small settling times.To fulfill these demands, combining feedback with feedforward control is normally implemented [1]. Figure 1 displays how a feedforward controller can be involved in feedback control systems.In total, there are two kinds of modes.The feedback controller guarantees stability and improves disturbance rejection, while the feedforward controller enhances tracking performance such that the feedforward controller should be designed as F = G † (in the first mode) and F = G † c (in the second mode), where the symbol " † " denotes the Moore-Penrose pseudo-inverse [2].So system inversion is the key to the problem of feedforward control.Actually, in addition to applications in control systems, system inversion is frequently used in the areas of sensor calibration, loudspeaker linearization, digital predistortion for radio frequency communications, and so on [3].So system inversion plays an important role in various research areas.System inversion can be conducted by direct inversion and indirect inversion [4].There exist several kinds of intrinsic limitations of direct inversion approaches; here, an example is used to illustrate this, denoting the transfer function of a finite-dimensional, discrete-time, single-input single-output, linear, constant dynamical system as where A ∈ R n×n , B ∈ R n , C ∈ R 1×n , D ∈ R, I n denotes an n-dimensional identity matrix, and (A, B, C, D) is a state-space realization of G(z).Based on (1), the inversion of G(z) can be represented as [5] G † (z) There are at least two challenges in the direct inversion (2).(i) When D † does not exist, the direct inversion cannot be conducted; (ii) when there exist nonminimum-phase zeros (For discrete-time systems, nonminimum-phase zeros are zeros that lie outside the unit disk) in G(z), the inversion G † (z) will be unstable.
Due to the limitations of direct inversion approaches, various indirect inversion approaches have thus been proposed.A possible classification of existing indirect system inversion approaches consists in distinguishing between preview-based and non-previewbased approaches [4,6]: (a) Preview-based inversion approaches can be further categorized into infinite previewbased approaches and finite preview-based approaches.Infinite preview-based approaches admit an exact stable inversion solution; however, such a solution may require an infinite pre-actuation [7][8][9][10].Because the length of the pre-actuation is proportional to the length of the desired output preview, infinite preview-based approaches are not applicable from a practical point of view.To handle the problem of applicability, finite preview-based approaches have been proposed [11][12][13][14][15][16][17][18].(b) Non-preview-based inversion approaches are preferred in practice.A family of approaches called pseudo-inversion, which can be conducted without preview, has been proposed [19,20].However, such approaches will encounter other problems such as the difficulty of choosing a suitable basis function.Direct system identification-based inversion approaches use the input-output data from the system to be inverted to identify the inverse system directly; however, the system identification cannot be conducted when the system to be inverted is not stable [3].For signal modeling-based inversion approaches, the input signal, which is to be reconstructed, must be a periodic signal under stationary operating conditions [21,22].
It should be noted that guaranteeing the stability of obtained solutions is a priority of both direct and indirect inversion approaches.So for some indirect inversion approaches, stability is ensured, but infinite or finite pre-actuation is needed, which cannot be applied well in practice because sometimes the desired output in unknown.
In this paper, an entirely different system inversion approach by combining timedomain observer design and frequency-domain subspace identification is proposed.The presented approach can guarantee the stability of obtained system inversion, and simultaneously the proposed approach does not need any pre-actuation.Moreover, the approach can be applied to stable or unstable systems/minimum-phase or nonminimum-phase systems to be inverted, and there is also no requirement for the type of input and output signals.Furthermore, it does not suffer from non-convex or input noise problems.
Even though the proposed approach has the above advantages by comparing with other approaches, the current version of the proposed approach can just used for solving the system inversion problem of systems of which the inputs do not have overlapping frequency ranges.In this paper, only single-input single-output systems are considered.
To facilitate the use of the proposed system inversion approach by third parties, a MATLAB toolbox implementing the approach is created in this paper after theoretical derivation of the approach.The full name of the MATLAB toolbox is INVerse System IDentification, with the first version abbreviated as INVSID 1.0.
The remainder of the paper is organized as follows.In Section 2, the inversion approach is first proposed, and then the corresponding MATLAB codes are generated, based on which the toolbox NIVSID 1.0 is created, followed by Section 3, in which the usage of the toolbox NIVSID 1.0 is validated by both simulation and practical examples.Finally, conclusions and perspectives are given in Section 4.

Creation of INVSID 1.0
This section is illustrated by using three subsections with a progressive relationship.The proposed inverse identification approach is first presented, then on the basis of which the corresponding MATLAB codes and the toolbox INVSID 1.0 are finally created.It should be noted that the toolbox INVSID 1.0 is only used for the inverse identification of single-input single-out systems.

Inverse Identification Approach
A finite-dimensional, discrete-time, single-input single-output, linear, constant dynamical system G d which is minimal-realized (The system is minimal-realized if and only if it is both controllable and observable) and proper (The system is proper when the degree of the numerator does not exceed the degree of the denominator of its transfer function; otherwise, the system is improper) is given, and the given system G d can be either stable or unstable, and the sampling period of system G d is T s in seconds.
Figure 2 is used to demonstrate the basic idea behind the proposed inverse system identification approach, based on which the inverse model of the nominal model G d can be derived.As can be seen in Figure 2, the proposed approach mainly consists of five steps: The above steps are discussed in detail as follows.
Step 1: If the transfer function of the system G d is denoted as G d (z), and represents a minimal realization of G d (z), the corresponding state-space model of the system G d can be represented as where x d (k) ∈ R n d is the state vector of the model (3), and u(k) ∈ R and y(k) ∈ R are the input and output of the system G d , respectively.Step 2: A set S containing N discrete-time sine signals is given as: where where f m denotes the frequency in Hz, and the sequence ( f m ) N m=1 is an arithmetic sequence, each sine wave in the set S can be represented as the output of a statespace model G m , i.e., where x m (k) ∈ R 2 denotes the state vector of the model ( 6), and the matrices A m and C m can be denoted as and Remark 1.The frequencies of the signals u m (k) for m = 1, 2, . . ., N in the set S can be specified by the following rule: where f b is a non-negative value, and d is a positive value.The rule in Equation ( 9) is not the only way to specify the frequencies.The values of f m for m = 1, 2, . . ., N belong to the range (0, f s 2 ) with f s = 1 T s .
Step 3: If the signal u m (k) is used as the input signal of the model (3), we can obtain where y m (k) denotes the output of the model (3) when the input signal is u m (k).
By augmenting the model (10) with the state vector of the model ( 6), an augmented model G a,m can be obtained, and it can be represented as where the state vector x a,m (k) can be denoted as and the matrices A a,m and C a,m can be denoted as and Based on the state-space model representation (11) for the augmented models G a,m , m = 1, 2, . . ., N, we can totally build N full-order observers corresponding to the augmented models G a,m , m = 1, 2, . . ., N, and we denote the m-th observer for the m-th augmented model G a,m as E a,m , which can be described by where L a,m denotes the observer gain, and xa,m (k) denotes the reconstructed value of x a,m (k) using the observer.
With the reconstructed value xa,m (k), we can obtain the reconstructed value of u m (k), which can be calculated by the following equation: where By combining (15) with ( 16), we can obtain a state-space model which can be regarded as a reconstructor of the input u m (k) of the model (10).
Let the model (18) be denoted as G inv,m , and the value of the frequency response function of the model G inv,m at the frequency f m can then be represented as Step 4: Let the inverse model of the model G d be denoted as G inv , and with the frequency response function values G inv,m (e −jΩ m T s ), m = 1, 2, . . ., N, the inverse model G inv in state-space representation can be identified by using the subspace-based system identification method in the frequency domain [23].After system identification, the identified inverse model is denoted as Ĝinv .
Remark 2. The effective frequency range of the inverse model G inv can be specified by selecting the range of the frequencies of the signals u m (k), m = 1, 2, . . ., N, in the set S.
Step 5: Connect the models G d and Ĝinv in series, and the resulting model can be represented as of which the frequency response function is written as where Ω = 2π f with f the ordinary frequency in Hz, so if the identified inverse model Ĝinv is perfect, the frequency response of the model G s should satisfy: By observing whether the frequency response function of the model G s within the frequency range specified by inverse system identification satisfies the above conditions (a) and (b) and to what extent it satisfies them, the goodness of the identified inverse model Ĝinv can be validated.

MATLAB Codes
The corresponding MATLAB (R2020b) commands to realize the above mentioned inverse identification approach are illustrated step by step.
Step 1: Firstly the transfer function of a nominal model G d is given, then a state-space realization of the transfer function is made, afterwards the Bode plot is plotted.The process can be realized by following MATLAB codes.» bode(sysss,opts); Step 2: According to the bandwidth of the nominal model G d , specify the frequencies f m , m = 1, 2, . . ., N, by using the rule stated in Remarks 1 and 2, and then stack all the frequencies into the vector F N , i.e., Based on Equations ( 7) and ( 8), create the matrices A m and C m for m = 1, 2, . . ., N in the models of the N discrete-time sine signals in the set S, and then stack them into the matrices A N and C N , respectively, i.e., and We can obtain the following MATLAB codes for the above. and Based on the expression of the model (18), calculate the transfer function of the reconstructors G inv,m by using where the observer gain L a,m can be chosen in a linear least-squares sense for stochastic systems, i.e., the observer gain L a,m can be obtained as the gain of the steady-state Kalman filter for the model (11) with process noise and measurement noise or can be obtained in a minimum mean-integral squared error sense [24].Denote the process noise and measurement noise as w(k) ∈ R n a with n a = n d + 2 and v(k) ∈ R, respectively, and assume that both {w(k), k = 1, 2, . ..} and {v(k), k = 1, 2, . ..} are white Gaussian sequences, w(k) ∼ N(0, Q) with Q > 0, v(k) ∼ N(0, R) with R > 0, and assume that the distribution of x a,m (0) is Gaussian, and assume that {w(k), k = 1, 2, . ..} and {v(k), k = 1, 2, . ..} are uncorrelated with x a,m (0) and with each other.Then, derive the gain of the steady-state Kalman filter for the model (11) by using the following equation [25]: where the value of P m can be derived as the unique solution of the following algebraic Riccati equation: under the following conditions: Then, stack all the calculated observer gains L a,m for m = 1, 2, . . ., N into the matrix L N , i.e., After calculating the reconstructor transfer functions G inv,m (z) for m = 1, 2, . . ., N based on (27), stack all the obtained transfer functions into the vector function E N , which can be denoted as Obtain the frequency response function values G inv,m (e −jΩ m T s ) for m = 1, 2, . . ., N by replacing z in Equation ( 27) with e −jΩ m T s for m = 1, 2, . . ., N, then stack all the frequency response function values into the vector G N , i.e., The above process can be realized by the following MATLAB codes.Step 5: The following MATLAB command can be used for the series connection of the models G d and Ĝinv .
» Gs=series(sysss,Ginv); Then, the Bode plot of the combined model can be displayed using the following codes.

Inverse System Identification Toolbox Creation
Building the inverse system identification toolbox consists of two parts.With the MATLAB m-file created in the first part, the inverse system identification toolbox INVSID 1.0 can be installed.The complete installation procedure contains five steps, from the first step about the selection of the item Package Toolbox from the Add-Ons menu to the end step about saving the created toolbox (Installation procedure of MATLAB toolboxes can be referred to the following link: https://www.mathworks.com/help/matlab/matlab_prog/create-and-share-custommatlab-toolboxes.html).

Pure Simulation Examples
In this section, two simulation examples are used to validate the effectiveness of the toolbox INVSID 1.0, i.e., check the effectiveness of the proposed inverse system identification approach.
Firstly, consider a discrete-time, single-input single-output, linear, constant dynamical system G d , of which the transfer function is described by with the sampling period T s = 1 × 10 −5 seconds.The Bode plot of the system G d is displayed in Figure 3.According to the transfer function (33), the following observations can be made: According to the above observations, it can be known that the direct inversion of the system G is a challenging problem.So we now turn to using the proposed toolbox INVSID 1.0 to identify the inverse model of the system G.The proposed inverse identification approach can specify the frequency range of interest, i.e., by selecting the values of f b , m, and d in Equation ( 9), the frequency range to be identified can be determined.The parameters shown in Table 1 are used as the inputs of the inverse identification toolbox.With the above inputs, the final output of the inverse identification toolbox is the identified inverse model which is the best model corresponding to the recommended singular value.The model order of the identified inverse model Ĝinv is recommended to be 4.The frequency response properties of the model Ĝinv with fourth order are demonstrated in Figure 4.
In addition, the inversion Ĝinv is identified using the MATLAB function n4sid with stability enforcement, so the identified model Ĝinv is stable and causal.By connecting the model G d and the model Ĝinv in series using Equation (19), the model G s can be obtained.The obtained model G s can then be used for validating the goodness of the identified inverse model Ĝinv in the frequency range of interest.The frequency response of the obtained model Ĝinv is shown in Figure 5; in the specified frequency range from 10 Hz to 500 Hz, the magnitude is nearly a constant near 0 dB, and the phase is nearly a constant around 0 degrees.The values of magnitude and phase can indicate the effectiveness of the proposed inverse identification toolbox for stable systems to be inverted.In practice, unstable systems are also frequently encountered.So the second numerical example is about using the toolbox INVSID 1.0 to solve the system inversion problem of an unstable system.Consider a discrete-time, single-input single-output, linear, constant dynamical system G * d , of which the transfer function is described by with the sampling period T s = 1 × 10 −5 s.The Bode plot of the system G * d is displayed in Figure 6.According to the transfer function (34), the following observations can be made: d has one nonminimum-phase zero.The parameters shown in Table 2 are used as the inputs of the inverse identification toolbox.With the above inputs, the final output of the inverse identification toolbox is the identified inverse model which is the best model corresponding to the recommended singular value.The model order of the identified model Ĝ * inv is recommended to be 4.The frequency response properties of the model Ĝ * inv with fourth order are demonstrated in Figure 7.In addition, the inversion Ĝ * inv is identified using the MATLAB function with stability enforcement, so the identified model Ĝ * inv is stable and causal.By connecting the model G * d and the model Ĝ * inv in series using Equation (19), the model G * s can be obtained.The frequency response of the obtained model Ĝ * inv is shown in Figure 8; in the specified frequency range from 10 Hz to 500 Hz, the magnitude is nearly a constant near 0 dB, and the phase is nearly a constant around 0 degrees.The values of magnitude and phase can indicate the effectiveness of the proposed inverse identification toolbox for unstable systems to be inverted.

Practical Hard Disk Drive System
Figure 9 shows a block diagram to illustrate the goal of the system inversion design, where r(k) ∈ R, u(k) ∈ R, and y(k) ∈ R represent the discrete-time reference, the input, and the output signals, respectively.Note that F can be implemented as a block in the feedforward controller design.In Figure 9, the overall transfer function from the reference signal r(k) to the output signal y(k) is where F(z) = Ĝinv denotes the inversion of the system G(z).Equation ( 35) can reflect the accuracy of the inversion F(z).We take the hard disk drive system as an illustrative example [26,27].The transfer function of the system with a sampling frequency of 26.4 kHz is G(z) = z −3 1.447663(z + 0.050852)(z + 2.494311) z 2 − 1.978354z + 0.978808 .
The Bode plot of the system G is shown in Figure 10.Based on the transfer function (36), the following observations can be made for the system G: The parameters shown in Table 3 are used as the inputs of the inverse identification toolbox.With the above inputs, the final output of the inverse identification toolbox is the identified inverse model, of which the model order is chosen as 10.The frequency response properties of the model Ĝinv with tenth order are demonstrated in Figure 11.The inversion Ĝinv is identified using the MATLAB function n4sid with stability enforcement, so the identified model Ĝinv is stable and causal.
As shown in Figure 9, by connecting the model G and the model Ĝinv in series using Equation (19), the model G s can be obtained.The frequency response of the obtained model Ĝ * inv is shown in Figure 12; in the specified frequency range from 0.1 Hz to 700 Hz, the magnitude is nearly a constant near 0 dB, and the phase is nearly a constant around 0 degrees.The above results can indicate that the proposed system inversion approach is effective in the practical application in the feedforward controller design of the hard disk drive system.
(a) Obtain a state-space representation of the nominal system G d .(b) Generate a number of sine signals u m (k), m = 1, 2, . . ., N, corresponding to a number of N specified frequencies, model them in state-space representation such that N signal models G m , m = 1, 2, . . ., N, can be obtained.(c) By combining the signal models G m for m = 1, 2, . . ., N with the model G d , respectively, the augmented models G a,m , m = 1, 2, . . ., N, can be obtained.Then, based on using the observers for the augmented models G a,m , m = 1, 2, . . ., N, inverse models G inv,m with m = 1, 2, . . ., N, which can be used for reconstructing the input signals u m (k), m = 1, 2, . . ., N, can be obtained.(d) Use frequency-domain system identification approaches to identify the inverse model of G d based on the frequency response function values of the models G inv,m with m = 1, 2, . . ., N at specified frequencies.(e) Choose the best inverse model by validating the identified inverse models.

Figure 2 .
Figure 2. Flowchart of the proposed inverse system identification approach.(The symbol " ˆ" denotes the reconstructed or the estimated value).

Figure 3 .
Figure 3. Bode plot of G d .

Figure 5 .
Figure 5. Bode plot of G s .

Figure 9 .
Figure 9. Block diagram to illustrate the goal of the system inversion design.Note that F can be implemented as a feedforward controller.
(a) G is stable.(b) G is proper.(c) G is minimal-realized.(d)G has one nonminimum-phase zero at around −2.5.

Part 1 :
Based on the MATLAB codes of inverse identification obtained in Section 2.2, a MATLAB function file, which is an m-file, can be created.The specific content of the m-file is given in Listing 1.

Table 1 .
Parameters for inverse identification.

Table 2 .
Parameters for inverse identification.

Table 3 .
Parameters for inverse identification.