This article is contained in *Scilab Control Engineering Basics* study module, which is used as course material for International Undergraduate Program in Electrical-Mechanical Manufacturing Engineering, Department of Mechanical Engineering, Kasetsart University.

### Module Key Study Points

- How to create a transfer function in Scilab
- Xcos CTR block usage
- Frequency response concept
- Scilab commands for plotting frequency responses

In general, the first step for control system analysis and design is to acquire a model that represents the actual plant to be controlled. Then a feedback diagram is constructed with this plant model and a controller described as transfer functions, either in continuous or discrete time domain. For analysis and design in frequency domain such as the so-called classical method, loopshaping, or Quantitative Feedback Theory (QFT), some form frequency response data is needed. Hence, in this module we show how to formulate a transfer function in Scilab and plot its frequency response.

### Transfer Function Formulation

To be concrete, we consider in Figure 1 a simple diagram of robot joint driven by DC motor through a gear transmission with ratio r:1 [1].

Figure 1 robot joint connected to DC motor via a gear transmission

Let J_{m} = J_{a} + J_{g} be the sum of motor and gear inertia. By simple calculation, it is easy to show that the rotational motion in terms of θ_{m} is described by

where k_{t} represents torque constant. We want to describe a model in transfer function form so that a block diagram can be drawn. To develop the electrical side of DC motor, consider the model shown in Figure 2.

Figure 2 a model of permanent magnet DC motor

By Kirschoff’s voltage law, we have

where k_{e} is back emf constant. From now on we omit the a subscript in the armature inductance and resistance. It is left to the reader to verify that, in Laplace domain, the joint dynamics in Figure 1 can be described by

This can be drawn as a block diagram in Figure 3.

Figure 3 block diagram of the robot joint dynamics in Figure 1

The transfer function from V(s) to θ_{m} can be derived by setting τ_{l} = 0 , which gives

Similarly, the transfer function from τ_{l} to θ_{m} is found by setting V=0.

To simplify the equation further, we can assume that the electrical constant L/R is much smaller than the mechanical constant J_{m}/B_{m}. So the transfer functions in (5) and (6) reduce to

respectively. These two equations correspond to second order differential equation in time domain

By omitting parameter subscripts, (9) can be rewritten as

with B = B_{m} + k_{e}k_{t}/R represents effective damping, u = (K_t/R)V control input, and d = τ_{l}(t)/r disturbance input. The reduced block diagram of (10) can be drawn as in Figure 4.

Figure 4 reduced block diagram of robot joint dynamics

So, the transfer function for a robot joint driven by DC motor we will be using in our study modules is in the form

Let’s put some values to the parameters, say, J = 10, B = 0.1. Hence the resulting transfer function becomes

Now we demonstrate how to construct a transfer function such as (12) in Scilab. One method begins by creating the Laplace variable s

-->s=poly(0,'s') s = s

and then use it to form P as described by (12)

-->P = 1/(10*s^2+0.1*s) P = 1 --------- 2 0.1s + 10s

So far so good. It already looks like (12). However, this data format still lacks some inside information necessary for further processing such as frequency response plot. So one more step is needed to convert it to a continuous-time linear transfer function, by using the syslin() command

-->P = syslin('c',P) P = 1 --------- 2 0.1s + 10s

The interactive response shown in Scilab console does not look any different than before. But the first argument ‘c’ passed to syslin tells Scilab that this is a continuous-time system. Now P is ready to be used by other commands such as bode, freqresp etc.

Another method that eventually yields the same P is shown below.

-->num = 1; -->den = poly([0 0.1 10],'s','c'); -->P = syslin('c',num,den) P = 1 --------- 2 0.1s + 10s

In this method, the numerator and denominator of P are created first, as a number 1 for the former and a polynomial of ‘s’ for the latter. Note that the row vector passed to poly(), representing the polynomial coefficients, starts with the lowest order of s.

### Transfer Function Block in an Xcos Model

Xcos is a simulation engine that comes with Scilab. Here we assume that the reader has some familiarity with basic functionality of Xcos to the level that she could create some simple diagram using blocks from standard palettes.

The standard block for transfer function can be located in Continuous time systems palette, under the name CLR. When a user clicks on this block, a parameter dialog window emerges. For this block, numerator and denominator in terms of polynomial of s can be put into the input field directly. For our transfer function (12), the data is typed in as shown in Figure 5.

Figure 5 A dialog window for setting transfer function parameters

In the next section, we will connect some input to this plant and measure its output by an oscilloscope.

### Frequency Response

Figure 6 depicts frequency response concept in a nutshell. In words, for a Linear Time-Invariant (LTI) system driven by a sinusoid input, the output is a sinusoid with same frequency, only its magnitude A and phase φ might change. When the input frequency varies, this results in new values for A and φ. This pair of data through out a range of frequency, actually a vector of complex numbers, constitutes a frequency response for an LTI system.

Figure 6 An LTI system driven by sinusoid input

To experiment with this, build a model shown in Figure 7, or download sinio_compare.zcos, and run the simulation. Note that the transfer function used is the DC motor model from the voltage input to * shaft velocity * output; i.e., the integrator 1/s is removed.

Figure 7 sinusoid input and output comparison

See the input and output comparison from the scope. Change the input frequency by clicking on and put in new value. Do you see the change in amplitude and phase of the output?

### Relationship between transfer function and frequency response

You may remember from linear systems course that, for a continuous-time transfer function described in terms of Laplace variable s, frequency response can be achieved by letting s = jω.

By this relationship, a frequency response of a transfer function can be plotted “the hard way.” Using (12) as an example, we solve for P(jω) manually by substituting s = jω, which gives

Now we can plot the magnitude and phase versus frequency by using the following set of commands

f=linspace(0.01,1,1000); // frequency vector (Hz) w=2*%pi*f; // convert to rad/s Pmag=zeros(1,1000);Pph=zeros(1,1000); for k=1:1000, // compute P(jw) at each frequency point P=1/(-10*w(k)^2+0.1*%i*w(k)); [Pmag(k),Pph(k)]=polar(P); end Pph = (180/%pi)*Pph; // convert to degree figure(1) subplot(211),plot2d(f,Pmag); // magnitude plot xlabel("Frequency (Hz)"); ylabel("Magnitude"); subplot(212),plot2d(f,Pph); // phase plot xlabel("Frequency (Hz)"); ylabel("Phase (degree)")

that yields Figure 8. First a linear frequency vector is created in Hz and converted to radian/sec. Vectors for keeping magnitude and phase information are initialized

Figure 8 a simple frequency response plot

Then, the frequency response of P is computed at each frequency point , converted to polar form, and put into the magnitude and phase vectors. The phase vector is converted from radian to degree. Last, the magnitude and phase are plotted in the same figure using subplot() command.

Some drawback of this plot is obvious. The frequency range is quite limited and the plots are concentrated only on narrow low frequency region. This is the reason that a frequency response is plotted on a semi-logarithm scale and the magnitude is expressed in decibel, where A_{dB} = 20log_{10}(A). Slightly adjust the above Scilab code to yield a much nicer plot range as in Figure 9.

Figure 9 a frequency response plot that covers broader frequency range

f=logspace(-5,2,1000); w=2*%pi*f; // convert to rad/s Pmag=zeros(1,1000);Pph=zeros(1,1000); for k=1:1000, // compute P(jw) at each frequency point P=1/(-10*w(k)^2+0.1*%i*w(k)); [Pmag(k),Pph(k)]=polar(P); end Pmag = 20*log10(Pmag); // convert Pmag to dB Pph = (180/%pi)*Pph; // convert Pph to degree figure(2) subplot(211),plot2d("ln",f,Pmag); // magnitude plot xlabel("Frequency (Hz)"); ylabel("Magnitude"); subplot(212),plot2d("ln",f,Pph); // phase plot xlabel("Frequency (Hz)"); ylabel("Phase (degree)");

This type of frequency response is known as Bode plot. An easier way to plot from a transfer function created by syslin() is by the command bode().

-->s=poly(0,'s'); -->P = 1/(10*s^2+0.1*s); -->P = syslin('c',P); -->bode(P)

This yields the plot in Figure 10. Type help bode to see options such as how to adjust frequency range.

Figure 10 frequency response using Scilab bode command

*Keyword: *the frequency range that a plant is responsive is called bandwidth. Strictly speaking, bandwidth covers frequency region such that the gain is above 0.707 or -3 dB. When trying to identify bandwidth from a Bode plot, we can roughly indicate the frequency point where the magnitude curve touches 0 dB line. For example, from Figure 10, the bandwidth is about 0.05 Hz.

Another type of frequency response useful for control design is called a Nyquist plot. Actually, it is the same data expressed in different format. Recall that each point of frequency response is just a complex number. When it is described in polar form, we get a Bode plot. A Nyquist plot, on the other hand, is the frequency response described in rectangular form x + jy. So it constitutes a graph in complex plane. A Nyquist plot is easily produced with command nyquist, say,

-->nyquist(P)

yields the plot in Figure 11. We will discuss this type of data more in later modules.

Figure 11 frequency response from Scilab nyquist command

For the last part of this study module, we discuss an application of using frequency response that may be useful beyond control engineering field. A filter is an electronic circuit, either analog or digital, that can be used to alter the frequency response of a system to suit some particular needs. The most common one is a low pass filter (LPF) used to attenuate high frequency noise.

Let’s make our plant more realistic. Suppose in our DC motor robot joint with transfer function described by (12), the joint angle is measured by a potentiometer and the resulting voltage is fed to a 12-bit A/D input of a microcontroller. The interface circuit are properly designed to use the full range of A/D; i.e., joint angle at 0 degree corresponds to A/D value 0, and at 360 degree, it reads as 4095. Adding this conversion factor in series with output of (12) in effect raises the gain by 4095/360 or 11.375. The new plant is

The bandwidth of (15) is about 0.2 Hz, which you can verify by yourself from a Bode plot, or search it more accurately from data as follows

-->s=poly(0, 's'); -->P=syslin('c',11.375/(10*s^2+0.1*s)); -->f=linspace(0.01,1,1000); -->fres=repfreq(P,f); // compute frequency respons -->ibw=find(abs(fres)fbw=f(ibw(1)) // locate BW frequency with index fbw = 0.2022523

In this chunk of code, we use another Scilab frequency response function repfreq() , which returns a vector of frequency response without plotting anything. Then use function find() to locate the vector index with gain less than 0.707, and use that index to find the bandwidth frequency.

Now, suppose that during operation, the potentiometer circuit that measures plant output is somehow contaminated by some noise with higher frequency spectrum than the plant bandwidth. To make it simple, suppose the noise is a sin wave of magnitude 0.5, and its frequency is 5 Hz. We want to get rid of this noise using a simple passive LPF circuit as shown in Figure 12.

Figure 12 a passive LPF circuit

consisting of only a resistor and a capacitor. We want to design the LPF with cutoff frequency at the plant bandwidth 0.2 Hz. Consult some analog filter design cookbook to get that the cutoff frequency of this filter can be selected from

Simple calculation shows that practical values of R = 80KΩ and C = 10 μF yields the cutoff frequency 1.25 rad/s or approximately 0.2 Hz.

A transfer function for this LPF can be derived easily. It is left to the reader to verify that it equals

Substituting the chosen component values yields

Plot the frequency response to verify the cutoff frequency. Since we are now interested in only the magnitude plot, it is convenient to use Scilab function gainplot(), which does as its name says

-->H=syslin('c',1/(1+0.8*s)); -->gainplot(H)

The plot result in Figure 13 shows that our LPF has cutoff frequency about 0.2 Hz as desired.

Figure 13 gain plot of LPF (18)

Exercise: Construct an Xcos diagram and simulate to verify that this LPF filter works. Use (15) as plant transfer function and connect some command input to it, say, a narrow pulse. Then add a sinusoid noise signal with magnitude 0.5 and frequency 5 Hz to the plant output. Use CMSCOPE block that has two channels. On the first channel, measure the output contaminated by the noise, and on the second channel, measure the plant output after passing through the LPF. Do you see the noise attenuated? Try varying the noise magnitude and frequency. What happens when the noise frequency lies within the plant (and LPF) bandwidth?

### Summary

In this study module, we discuss transfer function and frequency response basics, concentrating on Scilab commands and Xcos block usages. This serves as a good foundation for later development, such as feedback control analysis and design that relies on block diagram manipulation and frequency domain. Bode plot is a powerful tool used in classical control design since WWII. In the past, plotting manually on a sheet of paper could be tedious. Nowadays, software such as MATLAB or Scilab can generate magnitude and phase plots easily and accurately. Some basic knowledge taught in a undergrad control course is still necessary for sanity check.

### Reference

- M.W.Spong, S. Hutchinson and M. Vidyasagar, Robot Modeling and Control. John Wiley & Sons. 2006.

### Supplement

Exercise solution: LPF_exercise.zcos