Next Article in Journal
A Prediction Method for the Average Winding Temperature of a Transformer Based on the Fully Connected Neural Network
Previous Article in Journal
Integration of Relation Filtering and Multi-Task Learning in GlobalPointer for Entity and Relation Extraction
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

INVSID 1.0: An Inverse System Identification Toolbox for MATLAB

by
Runzhe Han
1,*,
Christian Bohn
2 and
Georg Bauer
2
1
College of Optical, Mechanical and Electrical Engineering, Zhejiang A&F University, Hangzhou 311300, China
2
Institut für Elektrische Informationstechnik, Technische Universität Clausthal, 38678 Clausthal-Zellerfeld, Germany
*
Author to whom correspondence should be addressed.
Appl. Sci. 2024, 14(15), 6838; https://doi.org/10.3390/app14156838
Submission received: 21 May 2024 / Revised: 10 July 2024 / Accepted: 14 July 2024 / Published: 5 August 2024
(This article belongs to the Section Electrical, Electronics and Communications Engineering)

Abstract

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

1. 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
G ( z ) = C ( z I n A ) 1 B + D = A B C D ,
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 ) = A B D C B D D C D .
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-preview-based approaches [4,6]:
(a)
Preview-based inversion approaches can be further categorized into infinite preview-based 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 time-domain 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.

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

2.1. 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:
(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.
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 ( A d , B d , C d , D d ) represents a minimal realization of G d ( z ) , the corresponding state-space model of the system G d can be represented as
{ x d ( k + 1 ) = A d x d ( k ) + B d u ( k ) , y ( k ) = C d x d ( k ) + D d u ( k ) ,
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:
S = { u m ( k ) , m = 1 , 2 , , N } ,
where
u m ( k ) = sin ( 2 π f m k T s ) ,
where f m denotes the frequency in Hz, and the sequence ( f m ) m = 1 N is an arithmetic sequence, each sine wave in the set S can be represented as the output of a state-space model G m , i.e.,
{ x m ( k + 1 ) = A m x m ( k ) , u m ( k ) = C m x m ( k ) ,
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
A m = cos ( 2 π f m T s ) sin ( 2 π f m T s ) sin ( 2 π f m T s ) cos ( 2 π f m T s )
and
C m = 1 0 ,
respectively.
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:
f m = f b + ( m 1 ) d ,
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
{ x d ( k + 1 ) = A d x d ( k ) + B d u m ( k ) , y m ( k ) = C d x d ( k ) + D d u m ( k ) ,
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
{ x a , m ( k + 1 ) = A a , m x a , m ( k ) , y m ( k ) = C a , m x a , m ( k ) ,
where the state vector x a , m ( k ) can be denoted as
x a , m ( k ) = x d ( k ) x m ( k ) ,
and the matrices A a , m and C a , m can be denoted as
A a , m = A d B d C m 0 A m
and
C a , m = C d D d C m ,
respectively.
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
x ^ a , m ( k + 1 ) = ( A a , m L a , m C a , m ) x ^ a , m ( k ) + L a , m y m ( k ) ,
where L a , m denotes the observer gain, and x ^ a , m ( k ) denotes the reconstructed value of x a , m ( k ) using the observer.
With the reconstructed value x ^ a , m ( k ) , we can obtain the reconstructed value of u m ( k ) , which can be calculated by the following equation:
u ^ m ( k ) = C r , m x ^ a , m ( k ) ,
where
C r , m = 0 × C d C m .
By combining (15) with (16), we can obtain a state-space model
{ x ^ a , m ( k + 1 ) = ( A a , m L a , m C a , m ) x ^ a , m ( k ) + L a , m y m ( k ) , u ^ m ( k ) = C r , m x ^ a , m ( k ) ,
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 G inv , m ( e j Ω m T s ) where Ω m = 2 π f m .
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 G ^ 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 G ^ inv in series, and the resulting model can be represented as
G s = G ^ inv G d
of which the frequency response function is written as
G s ( e j Ω T s ) = | G s ( e j Ω T s ) | e j G s ( e j Ω T s ) ,
where Ω = 2 π f with f the ordinary frequency in Hz, so if the identified inverse model G ^ inv is perfect, the frequency response of the model G s should satisfy:
(a)
| G s ( e j Ω T s ) | = 1 .
(b)
G s ( e j Ω T s ) { θ θ = 2 k π with k = 0 , 1 , 2 , } .
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 G ^ inv can be validated.

2.2. 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.
>> systf=tf(numerator,denominator,Ts);
>> sysss=ss(systf);
>> figure;
>> opts=bodeoptions;
>> opts.FreqUnits=‘Hz’;
>> 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.,
F N = f 1 f 2 f N .
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.,
A N = A 1 A 2 A N ,
and
C N = C 1 C 2 C N .
We can obtain the following MATLAB codes for the above.
>> dim=N;
>> FN=zeros(dim,1);
>> for m=1:dim
    FN(m,:)=fb+(m-1)*d;
   end
>> AN=zeros(2*dim,2);
>> for m=1:dim
    AN(((2*m-1):2*m),:)=[cos(2*pi*FN(m,:)*Ts) sin(2*pi*FN(m,:)*Ts)
                         -sin(2*pi*FN(m,:)*Ts) cos(2*pi*FN(m,:)*Ts)];
   end
>> CN=zeros(dim,2);
>> for m=1:dim
    CN(m,:)=[1 0];
   end
					   
Step 3:
With the created matrices A m and C m for m = 1 , 2 , , N , create the matrices A a , m , C a , m , and C r , m for m = 1 , 2 , , N by using Equation (13), Equation (14), and Equation (17), respectively. Then stack the matrices A a , m and C a , m , m = 1 , 2 , , N , into the matrices A N a , C N a , and C N r , respectively, i.e.,
A N a = A a , 1 A a , 2 A a , N ,
C N a = C a , 1 C a , 2 C a , N ,
and
C N r = C r , 1 C r , 2 C r , N .
Based on the expression of the model (18), calculate the transfer function of the reconstructors G inv , m by using
G inv , m ( z ) = C r , m [ z I ( A a , m L a , m C a , m ) ] 1 L a , m ,
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]:
L a , m = A a , m P m C a , m T ( C a , m P m C a , m T + R ) 1 ,
where the value of P m can be derived as the unique solution of the following algebraic Riccati equation:
P m = A a , m [ P m P m C a , m T ( C a , m P m C a , m T + R ) 1 C a , m P m ] A a , m T + Q ,
under the following conditions:
(a)
( A a , m , C a , m ) is detectable. (A system is detectable if all the unobservable states are stable)
(b)
( A a , m , Q ) is controllable.
Then, stack all the calculated observer gains L a , m for m = 1 , 2 , , N into the matrix L N , i.e.,
L N = L a , 1 L a , 2 L a , N .
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
E N = G inv , 1 ( z ) G inv , 2 ( z ) G inv , N ( z ) .
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.,
G N = G inv , 1 ( e j Ω 1 T s ) G inv , 2 ( e j Ω 2 T s ) G inv , N ( e j Ω N T s ) .
The above process can be realized by the following MATLAB codes.
>> r1=size(sysss.a,1)+2;
>> AaN=zeros(r1*dim,r1);
>> for m=1:dim
    r2=1+(m-1)*r1;
    AaN(r2:r1*m,:)=[sysss.a sysss.b*CN(m,:)
                    zeros(2,size(sysss.a,1)) AN(((2*m-1):2*m),:)];
   end
>> CaN=zeros(dim,r1);
>> for m=1:dim
    CaN(m,:)=[sysss.c sysss.d*CN(m,:)];
   end
>> CrN=zeros(dim,r1);
>> for m=1:dim
    CrN(m,:)=[zeros(1,size(sysss.a,1)) CN(m,:)];
   end
>> LN=zeros(r1,dim);
>> for m=1:dim
    Q=pc*eye(r1);
    R=mc;
    r2=1+(m-1)*r1;
    LN(:,m)=dlqr(AaN(r2:r1*m,:)’,CaN(m,:)’,Q,R)’;
   end
>> g=cell(dim,1);
>> GN=zeros(dim,1);
>> for m=1:dim
    r2=1+(m-1)*r1;
    sysr=ss(AaN(r2:r1*m,:)-LN(:,m)*CaN(m,:),LN(:,m),CrN(m,:),0,Ts);
    if isstable(sysr)==1
     g{m,:}=sysr;
     GN(m,:)=frd(g{m,:},FN(m,:),‘Hz’).ResponseData;
    else
     break
   end
Remark 3.
The values of Q and R can be tuned by changing the values of pc and mc.
Step 4:
With the calculated 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 MATLAB function n4sid, then the Bode plot of the identified inverse model can be obtained.
>> fdata=idfrd(GN,2*pi*FN,Ts);
>> opt=n4sidOptions("EnforceStability",1);
>> Ginv=n4sid(fdata,nx,’Ts’,Ts,opt);
>> figure;
>> opts=bodeoptions;
>> opts.FreqUnits=‘Hz’;
>> bode(Ginv,opts);
Remark 4.
The values of nx denote the inverse model order which can be specified.
Step 5:
The following MATLAB command can be used for the series connection of the models G d and G ^ inv .
>> Gs=series(sysss,Ginv);
Then, the Bode plot of the combined model can be displayed using the following codes.
>> figure;
>> opts=bodeoptions;
>> opts.FreqUnits=‘Hz’;
>> bode(Gs,opts);

2.3. Inverse System Identification Toolbox Creation

Building the inverse system identification toolbox consists of two parts.
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.
Listing 1. MATLAB function of inverse system identification.
1
function Ginv = INVSIDToolbox(numerator,denominator,Ts,fb,
    d,N,pc,mc,nx)
2
                          % numerator and denominator: The numerator and denominator
     coefficients of the transfer function of the nomial
    model G_d.
                        
3
                          % Ts: The sampling period of the nominal model G_d.
                        
4
                          % fb: The smallest frequency among the frequency
    components for inverse system identification.
                        
5
                          % d: The common difference.
                        
6
                          % N: The number of the frequency components for inverse
    system identification.
                        
7
                          % pc: The covariance of the process noise.
                        
8
                          % mc: The covariance of the measurement noise.
                        
9
                          % nx: Vector of model orders to scan.
                        
10
                          % Ginv: The identified inverse model.
                        
11
12
                          %% Step I
                        
13
systf=tf(numerator,denominator,Ts);
14
sysss=ss(systf);
15
figure;
16
opts=bodeoptions;
17
opts.FreqUnits=‘Hz’;
18
bode(sysss,opts);
19
20
                          %% Step II
                        
21
dim=N;
22
FN=zeros(dim,1);
23
for m=1:dim
24
FN(m,:)=fb+(m-1)*d;
25
                          end
                        
26
AN=zeros(2*dim,2);
27
for m=1:dim
28
AN(((2*m-1):2*m),:)=[cos(2*pi*FN(m,:)*Ts) sin(2*pi*FN(m,:)*Ts)
29
-sin(2*pi*FN(m,:)*Ts) cos(2*pi*FN(m,:)*Ts)];
30
                          end
                        
31
CN=zeros(dim,2);
32
for m=1:dim
33
CN(m,:)=[1 0];
34
                          end
                        
35
36
                          %% Step III
                        
37
r1=size(sysss.a,1)+2;
38
AaN=zeros(r1*dim,r1);
39
for m=1:dim
40
r2=1+(m-1)*r1;
41
AaN(r2:r1*m,:)=[sysss.a sysss.b*CN(m,:)
42
zeros(2,size(sysss.a,1)) AN(((2*m-1):2*m),:)];
43
                          end
                        
44
CaN=zeros(dim,r1);
45
for m=1:dim
46
CaN(m,:)=[sysss.c sysss.d*CN(m,:)];
47
                          end
                        
48
CrN=zeros(dim,r1);
49
for m=1:dim
50
CrN(m,:)=[zeros(1,size(sysss.a,1)) CN(m,:)];
51
                          end
                        
52
LN=zeros(r1,dim);
53
for m=1:dim
54
Q=pc*eye(r1);
55
R=mc;
56
r2=1+(m-1)*r1;
57
LN(:,m)=dlqr(AaN(r2:r1*m,:)’,CaN(m,:)’,Q,R)’;
58
                          end
                        
59
g=cell(dim,1);
60
GN=zeros(dim,1);
61
for m=1:dim
62
r2=1+(m-1)*r1;
63
sysr=ss(AaN(r2:r1*m,:)-LN(:,m)*CaN(m,:),LN(:,m),CrN(m,:),0,Ts);
64
if isstable(sysr)==1
65
g{m,:}=sysr;
66
GN(m,:)=frd(g{m,:},FN(m,:),‘Hz’).ResponseData;
67
                          else
                        
68
                          break
                        
69
                          end
                        
70
                          end
                        
71
72
                          %% Step IV
                        
73
fdata=idfrd(GN,2*pi*FN,Ts);
74
opt=ssestOptions("EnforceStability",1);
75
Ginv=ssest(fdata,nx,‘Ts’,Ts,opt);
76
figure
77
opts=bodeoptions;
78
opts.FreqUnits=‘Hz’;
79
bode(Ginv,opts);
80
81
                          %% Step V
                        
82
Gs=series(sysss,Ginv);
83
figure;
84
opts=bodeoptions;
85
opts.FreqUnits=‘Hz’;
86
bode(Gs,opts);
87
                          end
                        
Part 2:
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/matlabprog/create-and-share-custom-matlab-toolboxes.html).

3. Numerical Studies

3.1. 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
G d ( z ) = z z 2 1.5 z + 0.7
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:
(a)
G d is stable.
(b)
G d is proper.
(c)
G d is minimal-realized.
(d)
G d has one nonminimum-phase zero.
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 G ^ inv is recommended to be 4. The frequency response properties of the model G ^ inv with fourth order are demonstrated in Figure 4.
In addition, the inversion G ^ inv is identified using the MATLAB function n4sid with stability enforcement, so the identified model G ^ inv is stable and causal.
By connecting the model G d and the model G ^ 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 G ^ inv in the frequency range of interest. The frequency response of the obtained model G ^ 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
G d * ( z ) = z z 2 5 z + 6
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:
(a)
G d * is unstable.
(b)
G d * is proper.
(c)
G d * is minimal-realized.
(d)
G 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 inverse model G ^ inv * is recommended to be 4. The frequency response properties of the model G ^ inv * with fourth order are demonstrated in Figure 7.
In addition, the inversion G ^ inv * is identified using the MATLAB function n4sid with stability enforcement, so the identified model G ^ inv * is stable and causal.
By connecting the model G d * and the model G ^ inv * in series using Equation (19), the model G s * can be obtained. The frequency response of the obtained model G ^ 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.

3.2. 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
Y ( z ) R ( z ) = F ( z ) G ( z ) = G ^ inv ( z ) G ( z ) = G s ( z ) ,
where F ( z ) = G ^ 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.978354 z + 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:
(a)
G is stable.
(b)
G is proper.
(c)
G is minimal-realized.
(d)
G has one nonminimum-phase zero at around 2.5 .
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 G ^ inv with tenth order are demonstrated in Figure 11.
The inversion G ^ inv is identified using the MATLAB function n4sid with stability enforcement, so the identified model G ^ inv is stable and causal.
As shown in Figure 9, by connecting the model G and the model G ^ inv in series using Equation (19), the model G s can be obtained. The frequency response of the obtained model G ^ 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.
Furthermore, the well-known zero phase error tracking algorithm [11] is also implemented into the feedforward controller design, and according to the comparison between INVSID 1.0 and zero phase error tracking algorithm. For this example, it can be concluded that within the specified frequency range, INVSID 1.0 performs better than the zero phase error tracking algorithm in terms of phase, and in addition, INVSID 1.0 can be implemented without any preview.

4. Conclusions and Perspectives

In this paper, an alternative system inversion approach is proposed, based on which the toolbox named INVSID 1.0 is developed. The advantages of the toolbox INVSID 1.0 can be concluded as follows:
(a)
The proposed inverse identification toolbox can be used for stable or unstable systems, minimum-phase or nonminimum-phase systems.
(b)
Preview is not needed, i.e., the causality of the identified inverse model can be guaranteed.
(c)
Stability of the identified inverse model can be guaranteed.
(d)
The frequency range of interest can be specified.
(e)
Subspace identification is used such that there is no non-convex problem.
Furthermore, according to the theoretical derivation of the proposed system inversion approach, it can be indicated that the proposed approach can be used for systems with noise, because an observer is involved in the approach.
Currently, the inverse identification toolbox INVSID 1.0 is used for single-input single-output systems, while in the future, the proposed inverse system identification approach will be extended to identify the inverse models of general multiple-input multiple-output systems such that more advanced versions of the INVSID toolbox can be created, and squaring down approaches [28] have a potential to solve the extension problem.

Author Contributions

Conceptualization, R.H.; methodology, R.H.; software, R.H.; validation, R.H.; formal analysis, R.H.; investigation, R.H.; resources, R.H.; data curation, R.H.; writing—original draft preparation, R.H.; writing—review and editing, R.H., C.B. and G.B.; visualization, R.H.; supervision, R.H.; project administration, R.H.; funding acquisition, R.H. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by the Natural Science Foundation of Zhejiang province, P.R. China (Grant No. LQ23F030005) and the University Research Development Foundation in Zhejiang A&F University, P.R. China (Grant No. 203402000601).

Data Availability Statement

The original contributions presented in the study are included in the article, further inquiries can be directed to the corresponding author.

Conflicts of Interest

The authors declare no conflicts of interest.

References

  1. Goodwin, G.C.; Graebe, S.F.; Salgado, M.E. Control System Design; Prentice Hall: Upper Saddle River, NJ, USA, 2001. [Google Scholar]
  2. Golub, G.H.; Van Loan, C.F. Matrix Computations; John Hopkins University Press: Baltimore, MD, USA, 2012. [Google Scholar]
  3. Jung, Y. Inverse System Identification with Applications in Predistortion; Linköping University: Linköping, Sweden, 2018. [Google Scholar]
  4. van Zundert, J.; Oomen, T. On inversion-based approaches for feedforward and ILC. Mechatronics 2018, 50, 282–291. [Google Scholar] [CrossRef]
  5. Zhou, K.; Doyle, J.C.; Glover, K. Robust and Optimal Control; Prentice Hall: Englewood Cliffs, NJ, USA, 1996. [Google Scholar]
  6. Rigney, B.P.; Pao, L.Y.; Lawrence, D.A. Nonminimum phase dynamic inversion for settle time applications. IEEE Trans. Control Syst. Technol. 2009, 17, 989–1005. [Google Scholar] [CrossRef]
  7. Devasia, S.; Chen, D.; Paden, B. Nonlinear inversion-based output tracking. IEEE Trans. Autom. Control 1996, 41, 930–942. [Google Scholar] [CrossRef]
  8. Hunt, L.; Meyer, G.; Su, R. Noncausal inverses for linear systems. IEEE Trans. Autom. Control 1996, 41, 608–611. [Google Scholar] [CrossRef]
  9. Widrow, B.; Walach, E. Adaptive Inverse Control: A Signal Processing Approach; John Wiley & Sons: Hoboken, NJ, USA, 2008. [Google Scholar]
  10. Sogo, T. On the equivalence between stable inversion for nonminimum phase systems and reciprocal transfer functions defined by the two-sided Laplace transform. Automatica 2010, 46, 122–126. [Google Scholar] [CrossRef]
  11. Tomizuka, M. Zero phase error tracking algorithm for digital control. J. Dyn. Syst. Meas. Control 1987, 109, 65–68. [Google Scholar] [CrossRef]
  12. Gross, E.; Tomizuka, M.; Messner, W. Cancellation of discrete time unstable zeros by feedforward control. J. Dyn. Syst. Meas. Control 1994, 116, 33–38. [Google Scholar] [CrossRef]
  13. Roover, D.D.; Bosgra, O.H. Synthesis of robust multivariable iterative learning controllers with application to a wafer state motion system. Int. J. Control 2000, 73, 968–979. [Google Scholar] [CrossRef]
  14. Mirkin, L. On the H fixed-lag smoothing: How to exploit the information preview. Automatica 2003, 39, 1495–1504. [Google Scholar] [CrossRef]
  15. Hazell, A.; Limebeer, D.J.N. An efficient algorithm for discrete-time H preview control. Automatica 2008, 44, 2441–2448. [Google Scholar] [CrossRef]
  16. Zhang, Y.; Liu, S. Pre-actuation and optimal state to state transition based precise tracking for maximum phase system. Asian J. Control 2016, 18, 1728–1738. [Google Scholar] [CrossRef]
  17. Zhu, Q.; Zhang, Y.; Xiong, R. Stable inversion based precise tracking for periodic systems. Asian J. Control 2020, 22, 217–227. [Google Scholar] [CrossRef]
  18. Zou, Q. Optimal preview-based stable-inversion for output tracking of nonminimum-phase linear systems. Automatica 2020, 45, 230–237. [Google Scholar] [CrossRef]
  19. Jetto, L.; Orsini, V.; Romagnoli, R. Accurate output tracking for nonminimum phase nonhyperbolic and near nonhyperbolic systems. Eur. J. Control 2014, 20, 292–300. [Google Scholar] [CrossRef]
  20. Romagnoli, R.; Garone, E. A general framework for approximated model stable inversion. Automatica 2019, 101, 182–189. [Google Scholar] [CrossRef]
  21. Bohn, C.; Cortabarria, A.; Härtel, V.; Kowalczyk, K. Active control of engine-induced vibrations in automotive vehicles using disturbance observer gain scheduling. Control Eng. Pract. 2004, 12, 1029–1039. [Google Scholar] [CrossRef]
  22. Han, R.; Bohn, C.; Bauer, G. Comparison between calculated speed-based and sensor-fused speed-based engine in-cylinder pressure estimation method. IFAC-PapersOnLine 2022, 55, 223–228. [Google Scholar] [CrossRef]
  23. McKelvey, T.; Akçay, H.; Ljung, L. Subspace-based multivariable system identification from frequency response data. IEEE Trans. Autom. Control 1996, 41, 960–979. [Google Scholar] [CrossRef]
  24. O’Reilly, J. Observers for Linear Systems; Academic Press: London, UK, 1983. [Google Scholar]
  25. Lewis, F.L.; Xie, L.; Popa, D. Optimal and Robust Estimation: With an Introduction to Stochastic Control Theory; CRC Press: Boca Raton, FL, USA, 2008. [Google Scholar]
  26. Chen, X.; Tomizuka, M. A minimum parameter adaptive approach for rejecting multiple narrow-band disturbances with application to hard disk drives. IEEE Trans. Control Syst. Technol. 2011, 20, 408–415. [Google Scholar] [CrossRef]
  27. Wang, D.; Chen, X. H-based selective inversion of nonminimum-phase systems for feedback controls. IEEE/CAA J. Autom. Sin. 2020, 7, 702–710. [Google Scholar] [CrossRef]
  28. Saberi, A.; Sannuti, P. Squaring down by static and dynamic compensators. IEEE Trans. Autom. Control 1988, 33, 358–365. [Google Scholar] [CrossRef]
Figure 1. Inverse model-based feedforward–feedback control (r: reference signal; F: feedforward controller; C: feedback controller; G: plant model; G c : closed-loop system model; e: error; f: feedforward controller output; y: control system output).
Figure 1. Inverse model-based feedforward–feedback control (r: reference signal; F: feedforward controller; C: feedback controller; G: plant model; G c : closed-loop system model; e: error; f: feedforward controller output; y: control system output).
Applsci 14 06838 g001
Figure 2. Flowchart of the proposed inverse system identification approach. (The symbol “ ^ ” denotes the reconstructed or the estimated value).
Figure 2. Flowchart of the proposed inverse system identification approach. (The symbol “ ^ ” denotes the reconstructed or the estimated value).
Applsci 14 06838 g002
Figure 3. Bode plot of G d .
Figure 3. Bode plot of G d .
Applsci 14 06838 g003
Figure 4. Bode plot of G ^ inv .
Figure 4. Bode plot of G ^ inv .
Applsci 14 06838 g004
Figure 5. Bode plot of G s .
Figure 5. Bode plot of G s .
Applsci 14 06838 g005
Figure 6. Bode plot of G d * .
Figure 6. Bode plot of G d * .
Applsci 14 06838 g006
Figure 7. Bode plot of G ^ inv * .
Figure 7. Bode plot of G ^ inv * .
Applsci 14 06838 g007
Figure 8. Bode plot of G s * .
Figure 8. Bode plot of G s * .
Applsci 14 06838 g008
Figure 9. Block diagram to illustrate the goal of the system inversion design. Note that F can be implemented as a feedforward controller.
Figure 9. Block diagram to illustrate the goal of the system inversion design. Note that F can be implemented as a feedforward controller.
Applsci 14 06838 g009
Figure 10. Bode plot of G.
Figure 10. Bode plot of G.
Applsci 14 06838 g010
Figure 11. Bode plot of G ^ inv .
Figure 11. Bode plot of G ^ inv .
Applsci 14 06838 g011
Figure 12. Bode plot of G s . ZPETC: zero phase error tracking algorithm.
Figure 12. Bode plot of G s . ZPETC: zero phase error tracking algorithm.
Applsci 14 06838 g012
Table 1. Parameters for inverse identification.
Table 1. Parameters for inverse identification.
ParameterValue in MATLAB
numerator [ 0 , 1 , 0 ]
denominator [ 1 , 1.5 , 0.7 ]
T s 1 × 10 5
f b 10
d10
N50
pc 1 × 10 3
mc 1 × 10 3
nx2:10
Table 2. Parameters for inverse identification.
Table 2. Parameters for inverse identification.
ParameterValue in MATLAB
numerator [ 0 , 1 , 0 ]
denominator [ 1 , 5 , 6 ]
T s 1 × 10 5
f b 10
d10
N50
pc 1 × 10 3
mc 1 × 10 3
nx2:10
Table 3. Parameters for inverse identification.
Table 3. Parameters for inverse identification.
ParameterValue in MATLAB
numerator [ 0 , 0 , 0 , 1.4477 , 3.6845 , 0.1836 ]
denominator [ 1 , 1.9784 , 0.9788 , 0 , 0 , 0 ]
T s 1 / 26 , 400
f b 0.1
d 0.1
N7000
pc 1 × 10 3
mc 1 × 10 3
nx2:10
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Han, R.; Bohn, C.; Bauer, G. INVSID 1.0: An Inverse System Identification Toolbox for MATLAB. Appl. Sci. 2024, 14, 6838. https://doi.org/10.3390/app14156838

AMA Style

Han R, Bohn C, Bauer G. INVSID 1.0: An Inverse System Identification Toolbox for MATLAB. Applied Sciences. 2024; 14(15):6838. https://doi.org/10.3390/app14156838

Chicago/Turabian Style

Han, Runzhe, Christian Bohn, and Georg Bauer. 2024. "INVSID 1.0: An Inverse System Identification Toolbox for MATLAB" Applied Sciences 14, no. 15: 6838. https://doi.org/10.3390/app14156838

APA Style

Han, R., Bohn, C., & Bauer, G. (2024). INVSID 1.0: An Inverse System Identification Toolbox for MATLAB. Applied Sciences, 14(15), 6838. https://doi.org/10.3390/app14156838

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop