INVSID 1.0: An Inverse System Identification Toolbox for MATLAB
Abstract
:1. Introduction
- (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].
2. Creation of INVSID 1.0
2.1. Inverse Identification Approach
- (a)
- Obtain a state-space representation of the nominal system .
- (b)
- Generate a number of sine signals , , corresponding to a number of N specified frequencies, model them in state-space representation such that N signal models , , can be obtained.
- (c)
- By combining the signal models for with the model , respectively, the augmented models , , can be obtained. Then, based on using the observers for the augmented models , , inverse models with , which can be used for reconstructing the input signals , , can be obtained.
- (d)
- Use frequency-domain system identification approaches to identify the inverse model of based on the frequency response function values of the models with at specified frequencies.
- (e)
- Choose the best inverse model by validating the identified inverse models.
- Step 1:
- If the transfer function of the system is denoted as , and represents a minimal realization of , the corresponding state-space model of the system can be represented as
- Step 2:
- A set S containing N discrete-time sine signals is given as:
- Step 3:
- If the signal is used as the input signal of the model (3), we can obtainBy augmenting the model (10) with the state vector of the model (6), an augmented model can be obtained, and it can be represented asBased on the state-space model representation (11) for the augmented models , , we can totally build N full-order observers corresponding to the augmented models , , and we denote the m-th observer for the m-th augmented model as , which can be described byWith the reconstructed value , we can obtain the reconstructed value of , which can be calculated by the following equation:By combining (15) with (16), we can obtain a state-space modelLet the model (18) be denoted as , and the value of the frequency response function of the model at the frequency can then be represented as where .
- Step 4:
- Let the inverse model of the model be denoted as , and with the frequency response function values , , the inverse model 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 .
- Step 5:
- Connect the models and in series, and the resulting model can be represented as
- (a)
- .
- (b)
- .
By observing whether the frequency response function of the model 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 can be validated.
2.2. MATLAB Codes
- Step 1:
- Firstly the transfer function of a nominal model 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 , specify the frequencies , , by using the rule stated in Remarks 1 and 2, and then stack all the frequencies into the vector , i.e.,Based on Equations (7) and (8), create the matrices and for in the models of the N discrete-time sine signals in the set S, and then stack them into the matrices and , respectively, i.e.,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 and for , create the matrices , , and for by using Equation (13), Equation (14), and Equation (17), respectively. Then stack the matrices and , , into the matrices , , and , respectively, i.e.,Based on the expression of the model (18), calculate the transfer function of the reconstructors by using
- (a)
- is detectable. (A system is detectable if all the unobservable states are stable)
- (b)
- is controllable.
Then, stack all the calculated observer gains for into the matrix , i.e.,After calculating the reconstructor transfer functions for based on (27), stack all the obtained transfer functions into the vector function , which can be denoted asObtain the frequency response function values for by replacing z in Equation (27) with for , then stack all the frequency response function values into the vector , i.e.,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
- Step 4:
- With the calculated frequency response function values , , the inverse model 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);
- Step 5:
- The following MATLAB command can be used for the series connection of the models and .>> 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
- 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
- (a)
- is stable.
- (b)
- is proper.
- (c)
- is minimal-realized.
- (d)
- has one nonminimum-phase zero.
- (a)
- is unstable.
- (b)
- is proper.
- (c)
- is minimal-realized.
- (d)
- has one nonminimum-phase zero.
3.2. Practical Hard Disk Drive System
- (a)
- G is stable.
- (b)
- G is proper.
- (c)
- G is minimal-realized.
- (d)
- G has one nonminimum-phase zero at around .
4. Conclusions and Perspectives
- (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.
Author Contributions
Funding
Data Availability Statement
Conflicts of Interest
References
- Goodwin, G.C.; Graebe, S.F.; Salgado, M.E. Control System Design; Prentice Hall: Upper Saddle River, NJ, USA, 2001. [Google Scholar]
- Golub, G.H.; Van Loan, C.F. Matrix Computations; John Hopkins University Press: Baltimore, MD, USA, 2012. [Google Scholar]
- Jung, Y. Inverse System Identification with Applications in Predistortion; Linköping University: Linköping, Sweden, 2018. [Google Scholar]
- van Zundert, J.; Oomen, T. On inversion-based approaches for feedforward and ILC. Mechatronics 2018, 50, 282–291. [Google Scholar] [CrossRef]
- Zhou, K.; Doyle, J.C.; Glover, K. Robust and Optimal Control; Prentice Hall: Englewood Cliffs, NJ, USA, 1996. [Google Scholar]
- 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]
- Devasia, S.; Chen, D.; Paden, B. Nonlinear inversion-based output tracking. IEEE Trans. Autom. Control 1996, 41, 930–942. [Google Scholar] [CrossRef]
- Hunt, L.; Meyer, G.; Su, R. Noncausal inverses for linear systems. IEEE Trans. Autom. Control 1996, 41, 608–611. [Google Scholar] [CrossRef]
- Widrow, B.; Walach, E. Adaptive Inverse Control: A Signal Processing Approach; John Wiley & Sons: Hoboken, NJ, USA, 2008. [Google Scholar]
- 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]
- Tomizuka, M. Zero phase error tracking algorithm for digital control. J. Dyn. Syst. Meas. Control 1987, 109, 65–68. [Google Scholar] [CrossRef]
- 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]
- 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]
- Mirkin, L. On the H∞ fixed-lag smoothing: How to exploit the information preview. Automatica 2003, 39, 1495–1504. [Google Scholar] [CrossRef]
- Hazell, A.; Limebeer, D.J.N. An efficient algorithm for discrete-time H∞ preview control. Automatica 2008, 44, 2441–2448. [Google Scholar] [CrossRef]
- 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]
- Zhu, Q.; Zhang, Y.; Xiong, R. Stable inversion based precise tracking for periodic systems. Asian J. Control 2020, 22, 217–227. [Google Scholar] [CrossRef]
- Zou, Q. Optimal preview-based stable-inversion for output tracking of nonminimum-phase linear systems. Automatica 2020, 45, 230–237. [Google Scholar] [CrossRef]
- 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]
- Romagnoli, R.; Garone, E. A general framework for approximated model stable inversion. Automatica 2019, 101, 182–189. [Google Scholar] [CrossRef]
- 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]
- 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]
- 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]
- O’Reilly, J. Observers for Linear Systems; Academic Press: London, UK, 1983. [Google Scholar]
- 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]
- 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]
- 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]
- Saberi, A.; Sannuti, P. Squaring down by static and dynamic compensators. IEEE Trans. Autom. Control 1988, 33, 358–365. [Google Scholar] [CrossRef]
Parameter | Value in MATLAB |
---|---|
numerator | |
denominator | |
10 | |
d | 10 |
N | 50 |
pc | |
mc | |
nx | 2:10 |
Parameter | Value in MATLAB |
---|---|
numerator | |
denominator | |
10 | |
d | 10 |
N | 50 |
pc | |
mc | |
nx | 2:10 |
Parameter | Value in MATLAB |
---|---|
numerator | |
denominator | |
d | |
N | 7000 |
pc | |
mc | |
nx | 2: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. |
© 2024 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).
Share and Cite
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
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 StyleHan, 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 StyleHan, 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