INVSID 1.0: An Inverse System Identification Toolbox for MATLAB

2.1. Inverse Identification Approach

2.2. MATLAB Codes

Step 1

Given the transfer function of a nominal model $G_{\mathrm{d}}$, make a state-space realization of the transfer function:

>> systf=tf(numerator,denominator,Ts);

>> sysss=ss(systf);

>> Figure (1);

>> opts=bodeoptions;

>> opts.FreqUnits=’Hz’;

>> bode(sysss,opts);

>> set(title(’Bode plot of $G_{\mathrm{d}}$’),’interpreter’,’latex’);

>> grid;

Step 2

According to the bandwidth of the nominal model $G_{\mathrm{d}}$, specify the frequencies $f_m$, $m=1,2,\dots ,N$, using the rule stated in Remark 2 and Remark 3, and stack all the frequencies into the vector ${\mathit{F}}_N$, i.e.,

$${\mathit{F}}_N=\left(\begin{array}{c}{\mathit{f}}_1\\ {\mathit{f}}_2\\ ⋮\\ {\mathit{f}}_N\end{array}\right).$$

(21)

Based on Equation (7) and Equation (8), create the matrices ${\mathit{A}}_m$ and ${\mathit{C}}_m$ for $m=1,2,\dots ,N$ in the models of the N discrete-time sine signals in the set S, and stack them into the matrices ${\mathit{A}}_N$ and ${\mathit{C}}_N$, respectively, i.e.,

$${\mathit{A}}_N=\left(\begin{array}{c}{\mathit{A}}_1\\ {\mathit{A}}_2\\ ⋮\\ {\mathit{A}}_N\end{array}\right),$$

(22)

and

$${\mathit{C}}_N=\left(\begin{array}{c}{\mathit{C}}_1\\ {\mathit{C}}_2\\ ⋮\\ {\mathit{C}}_N\end{array}\right),$$

(23)

so we can get the following MATLAB codes correspondingly:

>> 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)

     AN(((2*m-1):2*m),:)=[-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 ${\mathit{A}}_m$ and ${\mathit{C}}_m$ for $m=1,2,\dots ,N$, create the matrices ${\mathit{A}}_{\mathrm{a},m}$, ${\mathit{C}}_{\mathrm{a},m}$, and ${\mathit{C}}_{\mathrm{r},m}$ for $m=1,2,\dots ,N$ by using Equation (13), Equation (14), and Equation (17), respectively. Then stack the matrices ${\mathit{A}}_{\mathrm{a},m}$ and ${\mathit{C}}_{\mathrm{a},m}$, $m=1,2,\dots ,N$, into the matrices ${\mathit{A}}_N^{\mathrm{a}}$, ${\mathit{C}}_N^{\mathrm{a}}$, and ${\mathit{C}}_N^{\mathrm{r}}$, respectively, i.e.,

$${\mathit{A}}_N^{\mathrm{a}}=\left(\begin{array}{c}{\mathit{A}}_{\mathrm{a},1}\\ {\mathit{A}}_{\mathrm{a},2}\\ ⋮\\ {\mathit{A}}_{\mathrm{a},N}\end{array}\right),$$

(24)

$${\mathit{C}}_N^{\mathrm{a}}=\left(\begin{array}{c}{\mathit{C}}_{\mathrm{a},1}\\ {\mathit{C}}_{\mathrm{a},2}\\ ⋮\\ {\mathit{C}}_{\mathrm{a},N}\end{array}\right),$$

(25)

and

$${\mathit{C}}_N^{\mathrm{r}}=\left(\begin{array}{c}{\mathit{C}}_{\mathrm{r},1}\\ {\mathit{C}}_{\mathrm{r},2}\\ ⋮\\ {\mathit{C}}_{\mathrm{r},N}\end{array}\right).$$

(26)

Based on the expression of the model (18), calculate the transfer function of the reconstructors $G_{\mathrm{inv},m}$ by using

$$G_{\mathrm{inv},m}\left(z\right)={\mathit{C}}_{\mathrm{r},m}{[\mathit{z}\mathit{I}-({\mathit{A}}_{\mathrm{a},m}-{\mathit{L}}_{\mathrm{a},m}{\mathit{C}}_{\mathrm{a},m})]}^{-1}{\mathit{L}}_{\mathrm{a},m},$$

(27)

where the observer gain ${\mathit{L}}_{\mathrm{a},m}$ can be chosen in a linear least-squares sense for stochastic systems, i.e., the observer gain ${\mathit{L}}_{\mathrm{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 $\mathit{w}\left(k\right)\in {\mathbb{R}}^{n_{\mathrm{a}}}$ with $n_{\mathrm{a}}=n_{\mathrm{d}}+2$ and $v\left(k\right)\in \mathbb{R}$, respectively, and assume that both $\left\{\mathit{w}\right(k),k=1,2,\dots \}$ and $\left\{v\right(k),k=1,2,\dots \}$ are white Gaussian sequences, $\mathit{w}\left(k\right)\sim N(0,\mathit{Q})$ with Q > 0, $v\left(k\right)\sim N(0,R)$ with $R>0$, and assume that the distribution of ${\mathit{x}}_{\mathrm{a},m}\left(0\right)$ is Gaussian, and assume that $\left\{\mathit{w}\right(k),k=1,2,\dots \}$ and $\left\{v\right(k),k=1,2,\dots \}$ are uncorrelated with ${\mathit{x}}_{\mathrm{a},m}\left(0\right)$ and with each other. Then derive the gain of the steady-state Kalman filter for the model (11) by using the following equation [25]

$${\mathit{L}}_{\mathrm{a},m}={\mathit{A}}_{\mathrm{a},m}{\mathit{P}}_m{\mathit{C}}_{\mathrm{a},m}^{\mathrm{T}}{({\mathit{C}}_{\mathrm{a},m}{\mathit{P}}_m{\mathit{C}}_{\mathrm{a},m}^{\mathrm{T}}+\mathit{R})}^{-1},$$

(28)

where the value of ${\mathit{P}}_m$ can be derived as the unique solution of the following algebraic Riccati equation

$${\mathit{P}}_m={\mathit{A}}_{\mathrm{a},m}[{\mathit{P}}_m-{\mathit{P}}_m{\mathit{C}}_{\mathrm{a},m}^{\mathrm{T}}{({\mathit{C}}_{\mathrm{a},m}{\mathit{P}}_m{\mathit{C}}_{\mathrm{a},m}^{\mathrm{T}}+\mathit{R})}^{-1}{\mathit{C}}_{\mathrm{a},m}{\mathit{P}}_m]{\mathit{A}}_{\mathrm{a},m}^{\mathrm{T}}+\mathit{Q},$$

(29)

under the following conditions:

(a)

$({\mathit{A}}_{\mathrm{a},m},{\mathit{C}}_{\mathrm{a},m})$ is detectable4.

(b)

$({\mathit{A}}_{\mathrm{a},m},\mathit{Q})$ is controllable.

Then stack all the calculated observer gains ${\mathit{L}}_{\mathrm{a},m}$ for $m=1,2,\dots ,N$ into the matrix ${\mathit{L}}_N$, i.e.,

$${\mathit{L}}_N=\left(\begin{array}{cccc}{\mathit{L}}_{\mathrm{a},1}& {\mathit{L}}_{\mathrm{a},2}& \dots & {\mathit{L}}_{\mathrm{a},N}\end{array}\right).$$

(30)

After calculating the reconstructor transfer functions $G_{\mathrm{inv},m}\left(z\right)$ for $m=1,2,\dots ,N$ based on (27), stack all the obtained transfer functions into the vector function ${\mathit{E}}_N$ which can be denoted as

$${\mathit{E}}_N=\left(\begin{array}{c}{G}_{\mathrm{inv},1}\left(z\right)\\ G_{\mathrm{inv},2}\left(z\right)\\ ⋮\\ G_{\mathrm{inv},N}\left(z\right)\end{array}\right).$$

(31)

Obtain the frequency response function values $G_{\mathrm{inv},m}\left(e^{-j{\Omega }_m{T}_{\mathrm{s}}}\right)$ for $m=1,2,\dots ,N$ by replacing z in Equation (27) with $e^{-j{\Omega }_m{T}_{\mathrm{s}}}$ for $m=1,2,\dots ,N$, then stack all the frequency response function values into the vector ${\mathit{G}}_N$ , i.e.,

$${\mathit{G}}_N=\left(\begin{array}{c}{G}_{\mathrm{inv},1}\left(e^{-j{\Omega }_1{T}_{\mathrm{s}}}\right)\\ G_{\mathrm{inv},2}\left(e^{-j{\Omega }_2{T}_{\mathrm{s}}}\right)\\ ⋮\\ G_{\mathrm{inv},N}\left(e^{-j{\Omega }_N{T}_{\mathrm{s}}}\right)\end{array}\right).$$

(32)

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,:)

     AaN(r2:r1*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=qc*eye(r1);

     R=mc;

     sysa=ss(AaN(r2:r1*m,:),zeros(r1,r1),CaN(m,:),zeros(1,r1),Ts);

     [∼,LN(:,m),∼]=kalman(sysa,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 4. 

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_{\mathrm{inv},m}\left(e^{-j{\Omega }_m{T}_{\mathrm{s}}}\right)$, $m=1,2,\dots ,N$, the inverse model $G_{\mathrm{inv}}$ in state-space representation can be identified by using the MATLAB function n4sid:

>> 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);

>> set(title(’Bode plot of ${\widehat{G}}_{\mathrm{inv}}$’),’interpreter’,’latex’);

>> grid;

Remark 5. 

The values of nx denotes 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_{\mathrm{d}}$ and ${\widehat{G}}_{\mathrm{inv}}$.

>> Gs=series(sysss,Ginv);

Then the Bode plot of the combined model can be displayed using:

>> figure;

>> opts=bodeoptions;

>> opts.FreqUnits=’Hz’;

>> bode(Gs,opts);

>> set(title(’Bode plot of $G_{\mathrm{s}}$’),’interpreter’,’latex’);

>> grid;

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 as follows.

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