Experiment No: 01
Experiment Name: Introduction of MATLAB software and determining the roots of a polynomial
equation.
Required Software: MATLAB 2014a or Later
MATLAB:
Polynomial Equation:
i. Express the following equation into the MATLAB software and hence find out the roots of
this equation
P ( s )=s 4 +6 s 3 +5 s 2+ 8 s +20
Result:
clc -5.1906 + 0.0000i
clear -1.6881 + 0.0000i
P=[1 6 5 8 20]; 0.4394 + 1.4455i
rootsP=roots(P) 0.4394 - 1.4455i
pause
ii. Express the following equation into the MATLAB software and hence find out the roots of
this equation
P ( s )=(s ¿ ¿ 5+6 s +5 s +8 s +20)( s ¿ ¿ 4+26 s +7 s +8 s+20)¿ ¿
3 2 3 2
Result:
clc
clear -25.7389 + 0.0000i -0.9031 + 0.0000i
P=conv([1 0 6 5 8 20],[1 26 7 8 20]);
rootsP=roots(P) -0.0000 + 2.0000i 0.3210 + 0.8703i
pause
-0.0000 - 2.0000i 0.3210 - 0.8703i
0.6641 + 1.8230i
0.6641 - 1.8230i
-1.3283 + 0.0000i
Experiment No: 02
Experiment Name: Consider the following transfer function, determine partial
fraction expansion using MATLAB software.
C( s) 2 s3 +5 s2 +3 s+6
=
R( s) s 3 +6 s 2 +11s +6
Result:
clc
clear r=
num=[2 5 3 5];
den=[1 6 11 6];
[r,p,k]=residue(num,den)
-6.5000
-3.0000
2.5000
C( s) 2 s3 +5 s2 +3 s+6
=
R( s) s 3 +6 s 2 +11s +6
−6 −3 2.5
¿ + + +2
s +3 s+2 s+ 1 p=
-3.0000
-2.0000
-1.0000
k=
2
Experiment No: 03 a
Experiment Name: Demonstrate the ration expression in frequency domain using
MATLAB software where given
Numerator = [1 2 7]
Result:
Denominator = [1 5 4 0]
T=
Program:
clc
s^2 + 2 s + 7
clear
-----------------
num=[1 2 7];
s^3 + 5 s^2 + 4 s
den=[1 5 4 0];
T=tf(num,den)
Experiment No: 03 b
Experiment Name: Demonstrate the state space expression using MATLAB software
where given
[ ][]
0 1 0 0 Result:
ẋ= 0 0 1 x+ 0 r
−3 −2 −5 10 StaeSpace =
a= b=
y= [ 1 0 0 ] x + 0.r
x1 x2 x3 u1
x1 0 1 0 x1 0
Program:
clc
x2 0 0 1 x2 0
clear
x3 -3 -2 -5 x3 10
A=[0 1 0;0 0 1;-3 -2 -5];
B=[0;0;10];
C=[1 0 0];
D=[];
StaeSpace=ss(A,B,C,D)
c=
d=
x1 x2 x3
u1
y1 1 0 0
y1 0
Experiment No: 04 a
Experiment Name: Represent the following transfer function in state space
Result:
C( s) 25
=
R( s) s 3 +13 s 2+36 s+ 25 A= B= C=
Program: -13 -36 -25 1 0 0 25
clc 1 0 0 0
clear
num=[25]; 0 1 0 0
den=[1 13 36 25];
[A,B,C,D]=tf2ss(num,den) D=
0
% StateSpace=ss(A,B,C,D)
Experiment No: 04 b
Experiment Name: Convert the following state space equation into transfer function,
G(s)=Y(s)/R(s)
Result:
[ ][]
0 1 0 0 num = T=
ẋ= 0 0 1 x+ 0 r
−3 −2 −5 10
0 0 0 10 10
---------------------
y= [ 1 0 0 ] x + 0.r
s^3 + 5 s^2 + 2 s + 3
Program:
den =
clc
clear
A=[0 1 0;0 0 1;-3 -2 -5];
B=[0;0;10]; 1.0000 5.0000 2.0000 3.0000
C=[1 0 0];
D=[];
[num,den]=ss2tf(A,B,C,D)
T=tf(num,den)
Experiment No: 05
Experiment Name: Simulink model of DC motor transfer function
Software: MATLAB
Theory: A common actuator in control systems is the DC motor. It directly provides rotary
motion and, coupled with wheels or drums and cables, can provide translational motion. The
electric circuit of the armature and the free-body diagram of the rotor are shown in the
following figure
we will assume that the input of the system is the voltage source ( ) applied to the motor's
armature, while the output is the rotational speed of the shaft . The rotor and shaft are
assumed to be rigid. We further assume a viscous friction model, that is, the friction torque is
proportional to shaft angular velocity.
The physical parameters for our example are:
(J) moment of inertia of the rotor 0.01 kg.m^2
(b) motor viscous friction constant 0.1 N.m.s
(Ke) electromotive force constant 0.01 V/rad/sec
(Kt) motor torque constant 0.01 N.m/Amp
(R) electric resistance 1 Ohm
(L) electric inductance 0.5 H
In general, the torque generated by a DC motor is proportional to the armature current and the
strength of the magnetic field. In this example we will assume that the magnetic field is
constant and, therefore, that the motor torque is proportional to only the armature current i by
a constant factor Kt as shown in the equation below. This is referred to as an armature
controlled motor.
T= Kt i
The back emf e, is proportional to the angular velocity of the shaft by a constant factor K e.
e= Keϴ’
In SI units, the motor torque and back emf constants are equal, that is K t= Ke; therefore, we
will use to represent both the motor torque constant and the back emf constant.
Applying Newton's law and Kirchoff's law to the motor system to generate the following
equations:
Simulink Model:
Experiment No: 06
Experiment Name: Observation of unit step response of overdamped, undamped,
underdamped and critically damped of second order systems.
Software: MATLAB
Theory:
Overdamped responses:
Poles: Two real at σ1; σ2
Natural response: Two exponentials with time constants equal to the reciprocal of the
pole locations, or
c t K1eσ1t K2eσ2t
Underdamped responses:
Poles: Two complex at σd jωd
Natural response: Damped sinusoid with an exponential envelope whose time
constant is equal to the reciprocal of the pole’s real part. The radian frequency of
the sinusoid, the damped frequency of oscillation, is equal to the imaginary part of the
poles, or
c t Aeσdt cos ωdt ϕ
Undamped responses
Poles: Two imaginary at jω1
Natural response: Undamped sinusoid with radian frequency equal to the imaginary part
of the poles, or
c t A cos ω1t ϕ
Critically damped responses
Poles: Two real at σ1
Natural response: One term is an exponential whose time constant is equal to the
reciprocal of the pole location. Another term is the product of time, t, and an exponential
with time constant equal to the reciprocal of the pole location, or
c t K1eσ1t K2teσ1t
The step responses for the four cases of damping has been determined by MATLAB
simulation software as follows
b
The general second order equation can be written as 2
s + as+ b
b
2
1 s + a s +b
R ( s )= C (s)
s
MATLAB Simulation:
1. Overdamped:
2. Undamped:
3. Underdamped:
4. Critically Damped
Four of this conditions are in one diagram
Discussion:
Experiment No: 07
Experiment Name: Calculate characteristics of a second-order system, such as damping
ratio, ζ; natural frequency, ωn; percent overshoot, %OS (pos); settling time, Ts; and peak
time, Tp from the following diagram using MATLAB software
Theory:
Software: MATLAB
Program:
clc
clear
p1=[1 3+7*i]; % Define polynomial containing
% first pole.
p2=[1 3-7*i]; % Define polynomial containing
% second pole.
deng=conv(p1,p2); % Multiply the two polynomials to
% find the 2nd order polynomial,
% as^2+bs+c.
omegan=sqrt(deng(3)/deng(1)) % Calculate the natural frequency,
% sqrt(c/a).
zeta=(deng(2)/deng(1))/(2*omegan)
% Calculate damping ratio,
% ((b/a)/2*wn).
Ts=4/(zeta*omegan) % Calculate settling time,
% (4/z*wn).
Tp=pi/(omegan*sqrt(1 -zeta^2)) % Calculate peak time,
% pi/wn*sqrt(1 -z^2).
pos=100*exp(-zeta* pi/sqrt(1 -zeta^2))
% Calculate percent overshoot
% (100*e^(-z*pi/sqrt(1-z^2)).
pause
Result:
omegan =
7.6158
zeta =
0.3939
Ts =
1.3333
Tp =
0.4488
pos =
26.0176
Discussion:
Experiment No: 07
Experiment Name: For the system shown in the following figure, find the peak time, percent
overshoot, and settling time by using MATLAB software.
Theory: We know for the unity feedback system closed loop transfer function
T(s)=G(s)/1+G(s)
Peak time
Percent overshoot
Settling time
Software: MATLAB
Program:
clc
clear
numg=[25]; % Define numerator of G(s).
deng=poly([0 -5]); % Define denominator of G(s).
' G(s)' % Display label.
G=tf(numg,deng) % Create and display G(s).
' T(s)' % Display label.
T=feedback(G,1) % Find T(s).
[numt,dent]=tfdata(T,'v'); % Extract numerator & denominator
% of T(s).
wn=sqrt(dent(3)) % Find natural frequency.
z=dent(2)/(2*wn) % Find damping ratio.
Ts=4/(z*wn) % Find settling time.
Tp=pi/(wn*sqrt(1-z^2)) % Find peak time.
pos=exp(-z*pi/sqrt(1-z^2))*100 % Find percent overshoot.
step(T) % Generate step response.
pause
Result:
T=
25
--------------
s^2 + 5 s + 25
Continuous-time transfer function.
wn = 5
z =0.5000
Ts =1.6000
Tp = 0.7255
pos =16.3034
Step Response
1.4
1.2
System: T
Peak amplitude: 1.16
1 Overshoot (%): 16.3
At time (seconds): 0.718
0.8
Amplitude
0.6
0.4
0.2
0
0 0.5 1 1.5 2 2.5
Time (seconds)
Discussion:
Experiment No: 08
Experiment Name: Find the range of gain, K, for the system of Figure 6.10 that will cause
the system to be stable, unstable, and marginally stable. Assume K > 0 by using MATLAB
software.
Theory: Routh-Hurwitz Criterion
Software: MATLAB
Program:
clc
clear
K=[1:1:2000];
for n=1:length(K);
dent=[1 18 77 K(n)];
poles=roots(dent);
r=real(poles);
if max(r) >=0,
poles
K=K(n)
break
end
end
pause
Result:
poles =
-18.0025 + 0.0000i
0.0012 + 8.7775i
0.0012 - 8.7775i
K=
1387
Discussion:
Experiment No: 09
Experiment Name: For the system of the following figure, evaluate the static error constants
and find the expected error for the standard step, ramp, and parabolic inputs.
Theory:
Software:
Program:
clc
clear
numg=500*poly([-2 -5 -6]);
deng=poly([0 -8 -10 -12]);
G=tf(numg,deng);
' Check Stability'
T=feedback(G,1);
poles=pole(T)
' Step Input'
Kp=dcgain(G)
ess=1/(1+Kp)
' Ramp Input'
numsg=conv([1 0],numg);
densg=poly([0 -8 -10 -12]);
sG=tf(numsg,densg);
sG=minreal(sG);%if any s or higher order of s is common between numerator
and denominator
Kv=dcgain(sG)
ess=1/Kv
' Parabolic Input'
nums2g=conv([1 0 0],numg);
dens2g=poly([0 -8 -10-12]);
s2G=tf(nums2g,dens2g);
s2G=minreal(s2G);
Ka=dcgain(s2G)
ess=1/Ka
pause
Result:
Check Stability
poles =
-516.9544
-5.7623
-5.4278
-1.8554
ans =
Step Input
Kp =
Inf
ess =
ans =
Ramp Input
Kv =
31.2500
ess =
0.0320
ans =
Parabolic Input
Ka =
ess =
Inf
Discussion:
Exp No: 10
Experiment Name: Sketch root locus by using MATLAB script and find out the limit of the
Gain (K) for remaining stable. Assume the following system
Theory:
Program:
clc
clear
numgh=[1 -4 20];
dengh=poly([-4 -6]);
'G(s)H(s)'
GH=tf(numgh,dengh)
rlocus(GH)
Output:
Result: From the figure it is cleared that for the stability gain K≤2.5
Discussion:
Exp No: 11
Experiment Name: Sketch the root locus for the system shown in the following figure
find the following:
a. The exact point and gain where the locus crosses the 0.45 damping ratio line
b. The exact point and gain where the locus crosses the jω-axis
c. The breakaway point on the real axis
d. The range of K within which the system is stable
Theory:
Program:
clc
clear
numgh=[1 -4 20];
dengh=poly([-2 -4]);
'G(s)H(s)'
GH=tf(numgh,dengh)
rlocus(GH)
z=0.2:0.05:0.5;
wn=0:1:10;
sgrid(z,wn)
title('Root Locus')
pause
rlocus(GH)
axis([-3 1 -4 4])
title('Close Up')
z=0.45;
wn=0;
sgrid(z,wn)
for k=1:3
[K,p]=rlocfind(GH)
end
pause
selected_point =
-1.5241 + 3.0190i
K=
0.4153
p=
-1.5327 + 3.0286i
-1.5327 - 3.0286i
Selected point 0.00 + 3.8991i
K = 1.5017
p=
0.00 + 3.8991i
0.00 - 3.8991i
selected_point = -2.8792 - 0.000i
K = 0.0248
p = -2.8792 + 0.000i
-2.8792 - 0.000i
Result:
Discussion:
Exp No: 12
Experiment Name: Study of an Ideal PI Controller/Compensator and compare the following
system before and after compensate.
Fig. : Closed-loop system for Example 9.1: a. before compensation; b. after ideal integral compensation
Theory:
Program:
clc
clear
num=[0 1];
numt=conv([1 .1],[0 0 1]);
den=conv([1 1],[1 2]);
dent=conv(den,[1 10]);
rlocus(num,dent); % Before Compensation
v=[-15 1 -6 6];
axis(v);
axis('square');
sgrid(0.174,[])
[K,r]=rlocfind(num,dent)
pause
dentt=conv(dent,[1 0]);
rlocus(numt,dentt); % After Compensation
v=[-15 1 -6 6];
axis(v);
axis('square');
sgrid(0.174,[])
[K,r]=rlocfind(numt,dentt)
pause
Output:
Brfore Compensation:
selected_point = -0.6891 + 3.8943i
K = 162.0131
r=
-11.5940 + 0.0000i
-0.7030 + 3.8993i
-0.7030 - 3.8993i
After Compensation:
selected_point =-0.6673 + 3.8455i
K = 159.2149
r=
-11.5626 + 0.0000i
-0.6736 + 3.8477i
-0.6736 - 3.8477i
-0.0902 + 0.0000i
Discussion:
Step Input Response (Exp 12) Observation Using MATLAB Simulink:
Discussion:
Exp No 13: Lag Compensator :
Exp No 14: PD Controller (Ideal)
Exp No 15: PD Controller (Lead)
Exp No 16: PID Controller
Exp No 17:
Nyquest
Exp No: 18
Bode Plot
Exp No 19:
PLC