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 GenerationMost 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 . 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.
LSID AlgorithmThe offline least-square algorithm can be classified as an parametric identification problem. Consider a discrete-time transfer function of order
(1)with coefficient vector
(2)The goal is to find , the best estimate of in the least-square sense. Let and be the input and output data sequence measured from a real plant, respectively. Suppose samples are gathered. Form a matrix
(3)Using this matrix together with column vector containing samples of output data, we can compute from
(4)Obviously, exists only when 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)This transfer function is unknown to us. Construct an Xcos setup hdm_sysid.zcos as shown in Figure 3.
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 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.
-->roots(P.den) ans = - 4.5241679 + 104.88968i - 4.5241679 - 104.88968i - 0.4891654 + 17.383778i - 0.4891654 - 17.383778i 0
(6)with coefficient vector
(7)Form matrix in (3) and vector with Scilab commands as follows
resulting in and of size 483×14 and 483×1, respectively. Then use (4) to compute the coefficient vector (7)
uv = ut.values; // PRBS input vector qv = qt.values; // plant output vector Y = qv(8:490,1); X = ; for j=1:483, X = [ 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)]; end
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
theta_ls = inv(X'*X)*X'*Y;
The numerator and denominator coefficients have values
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) ];
Then use syslin command to construct the estimated plant model from these coefficients
-->bn bn = 0.0000016 0.0000271 0.0000338 0.0000542 0.0000786 0.0000143 - 0.0000072 -->an an = - 1.4119431 - 0.2107028 0.3533239 0.2627608 0.1407983 0.1995993 - 0.3335632
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
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);
Now that we have the numerator and denominator polynomials, the continuous-time LSID model can be constructed
fs = -(s*Ts+2); gs = (s*Ts-2); nj=7; 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); end
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.
P_id = syslin('c',nums_id,dens_id) -->P_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
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.