VLabs logo

Virtual Labs
IIT Kharagpur

IIT KGP Logo

Simulation of Control Systems


Theory

In this experiment, the behaviour of some control systems will be analysed and simulation of their responses will be observed with the help of four problem statements.

Problem 1
Observe the response of a unity feedback system (Fig. 1) and its rise time, settling time and percentage maximum over shoot for different values of ζ (damping ratio).

Fig. 1. A simple second order system

Unity feedback system
If the output or some part of the output is returned to the input side and utilized as part of the system input, then it is known as feedback. Feedback plays an important role in order to improve the performance of the control systems. Negative feedback reduces the error between the reference input R(s) and system output. The above figure (Fig. 1) shows the block diagram of the negative feedback control system. Since H(s) = 1, hence it is a unity feedback system.
Transfer function of negative feedback control system is, $$T = \frac{G}{1 + GH}$$ where,
T is the transfer function or overall gain of negative feedback control system.
G is the open loop gain, which is function of frequency.
H is the gain of feedback path, which is function of frequency.

Second order system
The open loop transfer function of the second order unity feedback system is given by, $$G(s) = \frac{\omega_{n}^2}{s ( s + 2 \zeta \omega_n )}$$ The closed loop transfer function of the second order unity feedback system is given by, $$\frac{Y(s)}{R(s)} = \frac{G(s)}{1 + G(s)} = \frac{\omega_{n}^2}{s^2 + 2 \zeta \omega_n s + \omega_{n}^2}$$ where,
R(s) = Laplace transform of the input signal r(t),
Y(s) = Laplace transform of the output signal y(t),
ζ = Damping ratio,
ωn = Natrural frequency of oscillations.

The characteristic equation of the second order system is given by equating the denominator of the closed loop transfer function to zero. $$( s^2 + 2 \zeta \omega_n s + \omega_{n}^2 ) = 0$$ The expression for the response of the second order system can be written as, $$Y(s) = \frac{ \omega_{n}^2 }{ s^2 + 2 \zeta \omega_n s + \omega_{n}^2 } R(s) \tag{1}$$ When ζ = 0, the system is undamped.
When ζ = 1, the system is critically damped.
When 0 < ζ < 1, the system is under damped.

Unit step response of second order system
Apply unit step signal at the input of the second order system,
$$r(t) = u(t)$$ Taking laplace transform on the both sides
$$R(s) = \frac{1}{s}$$ Case 1 : when ζ = 0 i.e. system is undamped
The equation 1 becomes, $$Y(s) = \frac{\omega_{n}^2}{s^2 + \omega_{n}^2} R(s)$$ Substituting $$R(s) = \frac{1}{s}$$ $$Y(s) = \frac{\omega_{n}^2}{s ( s^2 + \omega_{n}^2 )} $$ Taking the inverse laplace transform on both the sides, we have, $$Y(t) = ( 1 - cos ( \omega_n t )) \ u(t) \tag{2}$$ Equation (2), shows that the unit step response of an undamped system is a continuous time signal of constant amplitude and frequency.

Case 2 : when ζ = 1 i.e. system is critically damped
The equation 1 becomes, $$Y(s) = \frac{\omega_{n}^2}{s^2 + 2 \omega_n s + \omega_{n}^2} R(s) = \frac{\omega_{n}^2}{( s + \omega_{n} )^2} R(s)$$ Substituting $$R(s) = \frac{1}{s}$$ $$Y(s) = \frac{\omega_{n}^2}{ s (s + \omega_{n} )^2 } $$ After the partial fraction, taking the inverse laplace transform on both the sides, we have, $$Y(t) = ( 1 - e^{- \omega_n t} - \omega_n t e^{- \omega_n t} ) \ u(t) \tag{3}$$ When the system is critically damped then, the equation (3) shows, that the unit step response of the second order system would try to reach the steady state step input.

Case 3 : when 0 < ζ < 1 i.e. system is under damped
The equation 1 can be written as, $$Y(s) = \frac{\omega_{n}^2}{(( s + \zeta \omega_{n} )^2 + \omega_{n}^2 ( 1 - \zeta^2 ))} R(s)$$ Substituting $$R(s) = \frac{1}{s}$$ $$Y(s) = \frac{\omega_{n}^2}{s (( s + \zeta \omega_{n} )^2 + \omega_{n}^2 ( 1 - \zeta^2 ))} $$ After the partial fraction, taking the inverse laplace transform on both the sides, we have, $$Y(t) = ( 1 - \frac{e^{- \zeta \omega_n t}}{\sqrt{(1 - \zeta^2)}} sin ( \omega_n \sqrt {( 1 - \zeta^2 )} \ t + cos^{-1}\zeta)) \ u(t) \tag{4}$$ Equation (4) shows, when the system is under damped then, the unit step response of the system is having damped oscillations i.e. response of decreasing amplitude.

Unit impulse response of second order system
Apply unit impulse signal at the input of the second order system,
$$r(t) = \delta(t)$$ Taking laplace transform on the both sides
$$R(s) = 1$$ The expression for the response of the second order system for unit impulse input can be written as, $$Y(s) = \frac{ \omega_{n}^2 }{ s^2 + 2 \zeta \omega_n s + \omega_{n}^2 } \tag{5}$$ After the necessary calculations, taking inverse laplace on both the sides, we get,
Case 1 : when ζ = 0 i.e. system is undamped
$$Y(t) = \omega_n sin ( \omega_n t) \ for \ t \ \geq \ 0 \tag{6}$$ Case 2 : when ζ = 1 i.e. system is critically damped
$$Y(t) = \omega_{n}^2 t e^{- \omega_n t} \ for \ t \ \geq \ 0 \tag{7}$$ Case 3 : when 0 < ζ < 1 i.e. system is under damped
$$Y(t) = \frac{ \omega_n e^{- \zeta \omega_n t}}{\sqrt{(1 - \zeta^2)}} sin ( \omega_n \sqrt {( 1 - \zeta^2 )} \ t ) \ for \ t \ \geq \ 0 \tag{8}$$

Problem 2
Plot the root loci for the unity feedback system given in Fig. 1. where, $$G( s ) = K G_1( s ) = \frac{K}{s ( s + 1 )( s + 2 )}$$ Assume that the amplifier gain (K) is varied from 0 to 50. Indicate the value of the gain K for which the root locus crosses the imaginary axis. Plot the output response for K = 0.4, 2, 6 and 12.
The characteristic equation of the system (considered in problem-2) is $$1 + K \ G_1(s) H(s) = 0 \tag{9}$$ The root locus is the path of the roots of characteristic equation traced out in s-plane as gain K is changed from 0 to ∞ and it is symmetrical about the real axis.
Construction of Root Locus
Branch
The root locus branches start at the open loop poles and end at open loop zeros. So, the number of root locus branches N is equal to the number of finite open loop poles P or the number of finite open loop zeros Z, whichever is greater. Mathematically, the number of root locus branches N can be written as $$N = P \tag{10}$$ when $$P\geq Z$$ $$N = Z \tag{11}$$ when $$Z\gt P$$ If odd number of open loop poles and zeros exist to the left side of a point on the real axis, then that point is on the root locus branch.

Centroid (α)
$$\alpha = \frac{\sum{Real \ part \ of \ finite \ open \ loop \ poles} - \sum{Real \ part \ of \ finite \ open \ loop \ zeros}}{P - Z} \tag{12}$$ Angle of asymptotes (θ) $$\theta = \frac{(2q + 1)180^\circ}{P - Z} \tag{13}$$ where q = 0, 1, 2,...(P-Z)-1

Break-away and Break-in points
If there exists a real axis root locus branch between two open loop poles, then there will be a break-away point in between these two open loop poles. If there exists a real axis root locus branch between two open loop zeros, then there will be a break-in point in between these two open loop zeros.
Note : Break-away and break-in points exist only on the real axis root locus branches.

Steps to find break-away and break-in points
1) Write K in terms of s from the characteristic equation (equation 9).
2) Differentiate K with respect to s and make it equal to zero. Substitute these values of s in the above equation (found from step 1).
3) The values of s for which the K value is positive are the break points.

Problem 3
Check the stability of the unity feedback system given in Fig.1. where, $$G ( s ) = \frac{2}{s ( s + 1 )( s + 2 )}, \ H ( s ) = 1$$ by drawing the Bode and Nyquist diagrams and indicate gain margin, phase margin, gain crossover frequency, phase crossover frequency.

Nyquist plots are used to draw the complete frequency response of the open loop transfer function.
If P = Number of open loop poles in the the right half of the s-plane and Z = Number of closed loop poles in the the right half of the s-plane then
the number of encirclements N can be written as, $$N = P - Z \tag{14}$$ The open loop control system is stable if there is no open loop pole in the the right half of the s-plane. $$P = 0 ; N = -Z$$ The closed loop control system is stable if there is no closed loop pole in the right half of the s-plane. $$Z = 0 ; N = P$$ General steps for drawing Nyquist plots
i) Locate the poles and zeros of open loop transfer function G(s)H(s) in s plane.
ii) Draw the polar plot by varying ω (angular frequency) from zero to infinity. If pole or zero present at s = 0, then varying ω from 0+ to infinity for drawing polar plot.
iii) Draw the mirror image of above polar plot for values of ω ranging from −∞ to zero (0− if any pole or zero present at s = 0).
iv) The infinite radius half circle will start at the point where the mirror image of the polar plot ends and this infinite radius half circle will end at the point where the polar plot starts.
v) After drawing the Nyquist plot, we can find the stability of the closed loop control system using the Nyquist stability criterion. If the critical point (-1+j0) lies outside the encirclement, then the closed loop control system is absolutely stable.

Stability Analysis using Nyquist Plots
From the Nyquist plots, we can identify whether the control system is stable, marginally stable or unstable based on the values of the following parameters.
Gain cross over frequency (ωgc) and phase cross over frequency (ωpc).
Gain margin (GM) and phase margin (PM)

Phase cross over frequency (ωpc)
The frequency at which the Nyquist plot intersects the negative real axis (phase angle is 180°) is known as the phase cross over frequency. It is denoted by ωpc.

Gain cross over frequency (ωgc)
The frequency at which the Nyquist plot is having the magnitude of one is known as the gain cross over frequency. It is denoted by ωgc.

The stability of the control system based on the relation between phase cross over frequency and gain cross over frequency is discussed below.

i) If the phase cross over frequency (ωpc) is greater than the gain cross over frequency (ωgc), then the system is stable.
ii) If the phase cross over frequency (ωpc) is equal to the gain cross over frequency (ωgc), then the system is marginally stable.
iii) If phase cross over frequency (ωpc) is less than gain cross over frequency (ωgc), then the system is unstable.

Gain Margin
The gain margin (GM) is equal to the reciprocal of the magnitude of the Nyquist plot at the phase cross over frequency. $$GM = \frac{1}{M_{pc}} \tag{15}$$ where, Mpc is the magnitude in normal scale at the phase cross over frequency.

Phase Margin
The phase margin (PM) is equal to the sum of 180° and the phase angle at the gain cross over frequency. $$PM = 180^\circ + \phi_{gc} \tag{16}$$ where, φgc is the phase angle at the gain cross over frequency.

The stability of the control system based on the relation between the gain margin and the phase margin is discussed below.
i) If the gain margin (GM) is greater than one and the phase margin (PM) is positive, then the control system is stable.
ii) If the gain margin (GM) is equal to one and the phase margin (PM) is zero degrees, then the control system is marginally stable.
iii) If the gain margin (GM) is less than one and / or the phase margin (PM) is negative, then the control system is unstable.

Bode Plot
The Bode plot or the Bode diagram consists of two plots :-
1. Magnitude plot
2. Phase plot
In both the plots, x-axis represents angular frequency (logarithmic scale). whereas, yaxis represents the magnitude (linear scale) of open loop transfer function in the magnitude plot and the phase angle (linear scale) of the open loop transfer function in the phase plot. The magnitude of the open loop transfer function in dB is - $$M = 20 \ log \ |G ( j \omega ) H ( j \omega )| \tag{17}$$ The phase angle of the open loop transfer function in degrees is - $$\phi = \angle G ( j \omega ) H ( j \omega ) \tag{18}$$ Note − The base of logarithm is 10.

Problem 4
Obtain the system response of a permanent magnet dc motor (given in Fig. 2)

Fig. 2. A shunt dc motor diagram

where, ω : Speed (rad/sec) of the motor.
Ra : Armature resistance (ohms)
La : Armature inductance (henry)
T : Load torque (newton-m)
V : Supply voltage (volts)
DC Motor model
Let us assume a general schematic diagram of a dc motor, shown in Fig. 3. Assume the following notations are used.

Fig. 3. Schematic diagram of dc motor

ea : Armature voltage (volts)
ia : Armature current (amp.)
Ra : Armature resistance (ohms)
La : Armature inductance (henry)
eb : Back emf (volts)
if : Field current (amp.)
TM : Motor torque (newton-m)
TL : Load torque (newton-m)
ω : Angular velocity (rad/sec)
J : Moment of inertia of the rotor (newton-m/rad/sec2)
B : Viscous friction coefficient (newton-m/rad/sec)
KT : Torque constant
Kb : Back emf constant

Upper case notations Ea, Ia, Eb, If are used for steady state values of the respective variables ea, ia, eb and if
In the present setup a permanent magnet dc motor is used, the field winding is thus absent and the air gap flux is constant. The input drive may therefore be applied to the armature only, i.e. only armature controlled operation is possible. The mathematical equations in this operating mode are, $$T_{M} = K_{T} I_a \tag{19}$$ $$e_{b} = K_{b} \omega \tag{20}$$ Armature circuit model $$L_a\frac{di_a}{dt} + R_a i_a + e_b = e_a \tag{21}$$ Mechanical model $$J\frac{d\omega}{dt} + B\omega + T_L = T_M \tag{22}$$ Taking Laplace Transform of (21) and (22), $$\frac{\omega(s)}{E_a(s)} = \frac{K_T}{(sL_a + R_a)(sJ + B) + K_T K_b} \tag{23}$$ Assuming the inductance of the armature circuit to be very small, the motor transfer function may be written as, $$G_M(s) = \frac{\omega(s)}{E_a(s)} = \frac{K_T/R_a}{Js + B + \frac{K_T K_b}{R_a}} = \frac{K_M}{s\tau_m + 1} \tag{24}$$ Motor gain constant (KM) $$K_M = \frac{K_T}{R_a B + K_T K_b}$$ Motor time constant (τm) $$\tau_m = \frac{R_a J}{R_a B + K_T K_b}$$ The armature controlled motor therefore has a first order type-0 transfer function and the two constant KM and τm depend upon motor parameters.