Offline Least-Square System Identification

In certain applications where tight specifications are needed, simple model-free control such as PID may not be adequate. More sophisticated schemes can be customized to the system under control. Linear controllers such as H_2/H_{\infty} are classified as model-based because the design process requires a mathematical model of the plant. Though achieving a perfect model of a real system is not feasible, to be a good enough representation, the model should capture all dominant dynamics, especially the “troublemakers” that might be present.

While physics laws can be used to form a theoretical model of a plant, they often disregard a few effects such as delay time, component aging, etc. Some complicated process such as reactive ion etching has interaction very difficult, if not impossible, to be modeled from physics.

With nowadays technology, several methods to practically obtain a good math model, commonly termed system identification, are developed and applied effectively in the control engineering field. One may suit a particular application better than the others, quite often depending on the input choice. For example, trying to model a loudspeaker from its step response is a bad idea. A strong impulse from a hammer works well for large structure such as a bridge, but not for a delicate hard drive.

In this article, we focus on an offline* system identification scheme based on least-square estimation, shortly addressed as LSID. This method, along with others, was discussed in our previous article System Modeling 101, with MATLAB examples. Here we just show how to compute and simulate LSID using Scilab and Xcos. While we create a different example, the theory, excitation input creation, and LSID algorithm structure, remain the same.

* By offline we mean the input and output data are collected first, then the LISD algorithm is executed to achieve a discrete-time transfer function of the plant.

PRBS Input Generation

Most system identification methods require a rich input to excite plant dynamics as much as it could. A convenient input signal for such purpose is Pseudo-Random Binary Sequence (PRBS). It is easy to generate in both simulation and embedded application software. With its magnitude varied between two constant levels, PRBS is particularly useful for motor control applications when we can select a proper magnitude to drive a particular motor (in contrast to a random magnitude signal such as white noise). It is the pattern of PRBS that does not repeat itself within the length chosen by number of delay bits; i.e., the maximum sequence length is 2^N -1 . Figure 1 shows how to generate a PRBS sequence using unit delay and XOR blocks in Xcos. The resulting sequence is shown in Figure 2.

Figure 1 PRBS generator using XCOS delay and Xor blocks

Figure 2 PRBS sequence generated from Xcos model in Figure 1

To elaborate the detail in Figure 1, the PRBS generator consists of 13 unit delay blocks connected in series. The output of block 1,3,4, and 13 are Xored and fed back to the input of block 1. It is essential to initialize the delays by setting initial conditions of some blocks to 1; otherwise, the sequence is not generated. Note that the output of unit delay 13 switches between 0 and 1. So it needs to be amplified and biased to achieve the signal in Figure 2.

The diagram in Figure 1 has an output port. So it is ready to be contained in an Xcos superblock later on.

LSID Algorithm

The offline least-square algorithm can be classified as an parametric identification problem. Consider a discrete-time transfer function of order n

(1)   \begin{equation*}  P(z) = \frac{b_1z^{n-1}+b_2z^{n-2}+\ldots + b_{n-1}z+b_n}{z^n+a_1z^{n-1}+\ldots+a_{n-1}z+a_n} \end{equation*}

with coefficient vector

(2)   \begin{equation*}  \theta = \left[ a_1, \ldots a_n, b_1, \ldots , b_n \right] ^T \end{equation*}

The goal is to find \theta_{LS}, the best estimate of \theta in the least-square sense. Let u[k] and y[k] be the input and output data sequence measured from a real plant, respectively. Suppose N_p samples are gathered. Form a matrix

(3)   \begin{eqnarray*} \mathbf{X} = \left( \begin{array}{cccccc} -y[n-1] & \ldots & -y[0] & u[n-1] & \ldots & u[0] \\ \nonumber -y[n] & \ldots & -y[1] & u[n] & \ldots & u[1] \\ \nonumber  : & : & : & : & : & : \\ \nonumber -y[N_p] & \ldots & -y[N_p - n +1] & u[N_p] & \ldots & u[N_p - n + 1]   \end{array} \right) \end{eqnarray*}

Using this matrix together with column vector Y=[y[n], \ldots , y[N_p]]^T containing N_p-n+1 samples of output data, we can compute \theta_{LS} from

(4)   \begin{equation*}  \theta_{LS} = (\mathbf{X}^T\mathbf{X})^{-1}\mathbf{X}^TY \end{equation*}

Obviously, \theta_{LS} exists only when \mathbf{X}^T\mathbf{X} is non-singular. The condition is true when the input signal is persistently exciting, which is the case for PRBS input

Example : We want to apply the LSID algorithm to a DC motor equipped with harmonic drive, whose linear model developed in our article HIL Implementation of Harmonic Drive Motor (Part I). Assume that the real plant is represented by

(5)   \begin{equation*}  P(s) = \frac{100000}{0.3s^5+3.008s^4+3400.08s^3+4056s^2+1000060s} \end{equation*}

This transfer function is unknown to us. Construct an Xcos setup hdm_sysid.zcos as shown in Figure 3.

Figure 3 hdm_sysid.zcos Xcos model for LSID

This is the way we would gather input and output data in practice. The real plant transfer function (5) is created by Scilab script hdminit.sce . The poles are at

 ans  = 
  - 4.5241679 + 104.88968i  
  - 4.5241679 - 104.88968i  
  - 0.4891654 + 17.383778i  
  - 0.4891654 - 17.383778i  

The PRBS generator in Figure 1, encapsulated in a superblock, is used as excitation input to the plant. Input and output data is then gathered with sampling period selected as T_s = 0.01 sec. Also notice the addition of random signal to the output labeled as measurement noise. In practice we cannot measure the output with 100% accuracy, so this makes the problem setup more realistic.

Figure 4 shows the PRBS input and plant output captured to workspace. The input magnitude is set to +/- 500. The output swings within 25 degree range. In a real motor setup, you should adjust the input magnitude not too large to avoid saturation, while not too small to combat with friction and motor deadband.

Figure 4 PRBS input and plant response from simulation

To apply LSID algorithm, first we need to decide on the order of estimated transfer function. The order must be high enough to capture all the essential dynamics of the real plant. The drawback of too-high the order is the unnecessary complexity of the model. Model reduction scheme may be applied afterwards to get rid of the insignificant dynamics.

Recall that we do not have information on the order of real plant. So in this experiment we try a model of order 7; i.e., the estimated discrete-time model is in the form

(6)   \begin{equation*}  P(z) = \frac{b_1z^6+b_2z^5+b_3z^4+b_4z^3+b_5z^2+b_6z+b_7}{z^7+a_1z^6+a_2z^5+a_3z^4+a_4z^3+a_5z^2+a_6z+a_7} \end{equation*}

with coefficient vector

(7)   \begin{equation*}  \theta = \left[ a_1, a_2, a_3, a_4, a_5, a_6, a_7, b_1, b_2, b_3, b_4, b_5, b_6, b_7 \right] ^T \end{equation*}

Form matrix \mathbf{X} in (3) and vector Y with Scilab commands as follows

uv = ut.values;  // PRBS input vector
qv = qt.values;  // plant output vector
Y = qv(8:490,1);
X = [];
for j=1:483,
    X = [
   -qv(j+6), -qv(j+5), -qv(j+4), -qv(j+3), -qv(j+2), -qv(j+1), -qv(j), uv(j+6), uv(j+5), uv(j+4), uv(j+3), uv(j+2), uv(j+1), uv(j)];

resulting in \mathbf{X} and Y of size 483×14 and 483×1, respectively. Then use (4) to compute the coefficient vector (7)

theta_ls = inv(X'*X)*X'*Y;

Now both the numerator and denominator coefficients of the estimated model (6) are contained in theta_ls. For clarity, we separate them to 2 vectors

bn = [theta_ls(8) theta_ls(9) theta_ls(10) theta_ls(11) theta_ls(12) theta_ls(13) theta_ls(14) ];
an = [theta_ls(1) theta_ls(2) theta_ls(3) theta_ls(4) theta_ls(5) theta_ls(6) theta_ls(7) ];

The numerator and denominator coefficients have values

 bn  =
    0.0000016    0.0000271    0.0000338    0.0000542    0.0000786      
             0.0000143  - 0.0000072  
 an  =
  - 1.4119431  - 0.2107028    0.3533239    0.2627608    0.1407983    
             0.1995993  - 0.3335632

Then use syslin command to construct the estimated plant model from these coefficients

z=poly(0,'z'); Ts = 0.01;
numz_id = bn(1)*z^6 + bn(2)*z^5 + bn(3)*z^4 + bn(4)*z^3 + bn(5)*z^2 + bn(6)*z + bn(7);
denz_id = z^7 + an(1)*z^6 + an(2)*z^5 + an(3)*z^4 + an(4)*z^3 + an(5)*z^2 + an(6)*z + an(7);
Pzsys_id = syslin(Ts,numz_id,denz_id);

To verify this model, we want to somehow compare it to the real plant (5). So our next step is to convert the discrete-time LSID model to continuous-time. Unfortunately, Scilab does not have a specific command to do so. We have written an article Convert Transfer Function from Discrete to Continuous-time on how to write a custom function for the conversion process using bilinear transform method.

For the problem at hand, we can perform the conversion using the following Scilab code

fs = -(s*Ts+2);
gs = (s*Ts-2);
nums_id = 0;
dens_id = fs^nj;
for j=1:nj,
   nums_id = nums_id + bn(j)*(fs^(nj-j))*(gs^j);
   dens_id = dens_id + an(j)*(fs^(nj-j))*(gs^j);

Now that we have the numerator and denominator polynomials, the continuous-time LSID model can be constructed

P_id = syslin('c',nums_id,dens_id)
 P_id  =
                                       2            3            4            
  - 0.0259087 + 0.0001201s + 0.0000005s - 8.126D-09s + 4.770D-11s - 
          5            6            7  
6.095D-14s - 1.988D-16s + 1.117D-19s   
                                       2            3            4            
  - 0.0349645 - 0.2005056s - 0.0037309s - 0.0006416s - 0.0000025s - 
          5            6            7  
1.697D-08s - 1.980D-11s - 2.503D-14s

In reality it is not possible to compare numerical values of the transfer functions since the real plant is not known, so the verification process is often done qualitatively, such as feeding some rich input (different from the sequence used in LSID) to both models and compare the responses. Another effective and convenient method is to compare the frequency response of the LSID model to the real plant. We do so in this example.

bode([P; P_id]);

The result is shown in Figure 5. The frequency response of LSID model resembles the real plant in the 0.1 – 10 Hz range. The integrator in real plant is estimated by a low frequency pole, while the responses diverge in the high frequency region. This result is expected in a real setup, since most sensors lose their accuracy at higher frequency.

Figure 5 Bode plot comparison between the real plant and LSID model

For your convenience, the whole LSID process discussed above is contained in a script file hdm_sysid.sce .


In this article we address a specific system identification based on least-square estimation, and show how to create a setup and simulate using Scilab/Xcos. Once the process is understood it can be applied to a real plant with confidence. The LSID takes less steps than some others such as swept-sin or DFT-based identification methods, and the result is in a transfer function form ready for analysis and control design.

The PRBS input is easy to generate by writing a C function, and is quite suitable to drive a motor since its magnitude switches between two selected values instead of varying randomly. We have implemented a PRBS generator in C to excite a DC motor, captured input and output data to Scilab, run the LSID algorithm, and used the estimated model to synthesize H-infinity controllers for position control. The whole procedure works smoothly and yields satisfactory results.

Supplement : Scilab script and Xcos model used in the article.



Comments are closed.