LABORATORY MANUAL
CONTROL SYSTEM
B.E 4th SEM 2016
Department of Instrumentation and Control Engineering
L D College of Engineering, Ahmedabad
1
LAB 1
%EXERCISE 1:
A=[1 0 1;2 3 4;-1 6 7]
A =
1 0 1
2 3 4
-1 6 7
B=[7 4 2;3 5 6;-1 2 1]
B =
7 4 2
3 5 6
-1 2 1
A+B
ans =
8 4 3
5 8 10
-2 8 8
A*B
ans =
6 6 3
19 31 26
4 40 41
A*A
ans =
0 6 8
4 33 42
4 60 72
2
A'
ans =
1 2 -1
0 3 6
1 4 7
inv(B)
ans =
0.1111 0.0000 -0.2222
0.1429 -0.1429 0.5714
-0.1746 0.2857 -0.3651
B'*A'
ans =
6 19 4
6 31 40
3 26 41
A*A+B*B-A*B
ans =
53 52 45
15 51 58
-2 28 42
det(A)
ans =
12
det(B)
ans =
-63.0000
det(A*B)
ans =
-756.0000
%EXERCISE 2:
A=[4 2 3;-1 1 3;2 5 7]
A =
4 2 3
-1 1 3
2 5 7
B=[1 2 3;8 7 6;5 3 1]
B =
1 2 3
8 7 6
5 3 1
C=eig(A)
C =
-0.8667
3.2339
9.6328
3
D=eig(B)
D =
11.8655
-2.8655
0.0000
E=eig(A*B)
E =
97.3809
-5.3809
0.0000
%EXERCISE 3:
A=[0 1 -3;2 3 -1;4 5 -2]
A =
0 1 -3
2 3 -1
4 5 -2
B=[-5;7;10]
B =
-5
7
10
inv(A)*B
ans =
-1.0000
4.0000
3.0000
%EXERCISE 4:
t=[0:0.1:2*pi];
y1=sin(t);
y2=t;
4
y3=t-(t.^3)/factorial(3)+(t.^5)/factorial(5)+(t.^7)/factorial(7);
plot(y1)
hold all
plot(y2)
hold all
plot(y3)
%EXERCISE 5:
%Ques 1:
t=[0:1:10*pi];
y=t.*cos(t);
plot(y)
5
%QUES 2:
t=[0:0.5:2*pi];
x=exp(t);
y=100+exp(3*t);
plot(x)
hold all
plot(y)
%EXERCISE 6:
%Ques a:
t=[0:0.1:1.0];
x=t;
y=t.^2;
z=t.^3;
plot(x)
hold all
plot(y)
hold all
plot(z)
6
%QUES b:
x=[-5:1:5];
y=[-5:1:5];
z=-7./(1+x.^2+y.^2);
plot(z)
7
LAB 2
%Q1%
N = [1 6 5 4 3];
D = [ 7 6 5 4 7];
H = N./D;
%(a)%
n=polyval(N,-10)
n(-10) = 4463
n(-5)=polyval(N,-5)
n(-1) = -17
n=polyval(N,-3)
n(-3) = -45
n=polyval(N,-1)
n(-1) = -1
%(b)%
d=polyval(D,-10)
d(-10) = 64467
d=polyval(D,-5)
d (-5)= 3737
d=polyval(D,-3)
d(-3) = 445
d=polyval(D,-1)
d(-1) = 9
%(c)%
h=polyval(H,-10)
h(-10) = 519.0000
h=polyval(H,-5)
h(-5) = -15.2857
h=polyval(H,-3)
h(-3) = -9.0000
h=polyval(H,-1)
h(-1) = -0.4286
%Q2%
w= 15%rad/s;
x=[0:0.1:15];
y=exp(-0.7*x).*sin(w*x);
8
plot(x,y)
%Q3%
w= 10%rad/s;
x=[0:0.05:15];
y=exp(-0.6*x).*cos(w*x);
plot(x,y)
%Q4%
9
t = [0 : 6*3.14];
x=sqrt(t).*sin(3*t);
y=sqrt(t).*cos(3*t);
z=0.8*t;
plot(t,x,y,z);
%Q5(a)%
A=[1 2 3 5 ;-2 5 7 -9 ;5 7 2 -5; -1 3 -7 7];
B=[ 21 ; 18 ; 25 ; 30];
X=inv(A)*B
X =
-1.5130
5.9108
0.3762
1.9126
%Q5(b)%
A=[1 2 3 4; 2 -2 -1 -1; 1 -3 4 -4;2 2 -3 4];
B=[8; -3; 8 ;-2];
10
X=inv(A)*B
X =
2.7869
4.4918
2.1311
-2.5410
%Q6%
syms x;
s1=exp(x^8);
y1=diff(s1)
y1 = 8*x^7*exp(x^8)
s2=(3*x^3)*exp(x^5);
y2=diff(s2)
y2 = 9*x^2*exp(x^5) + 15*x^7*exp(x^5)
s3=(5*x^3)-(7*x^2)+(3*x)+6;
y3=diff(s3)
y3 = 15*x^2 - 14*x + 3
%Q7%
syms x y;
y1=int(abs(x),0.2,0.7)
y1 = 9/40
s2=cos(y)+7*(y^2);
y2=int(s2,0.2,pi)
y2 =(7*pi^3)/3 - sin(1/5) - 7/375
y3=int(sqrt(x))
y3 =(2*x^(3/2))/3
y4=int( 7*x^5 - 6*x^4 + 11*x^3 +4*x^2 +8*x + 9)
y4 =(7*x^6)/6 - (6*x^5)/5 + (11*x^4)/4 + (4*x^3)/3 + 4*x^2 + 9*x
LAB 3
%Q1%
syms s t;
Fs= 1 /(s^4 + 5*s^3 + 7*s^2)
Fs =1/(s^4 + 5*s^3 + 7*s^2)
Ft=ilaplace(Fs,s,t)
Ft = t/7 + (5*exp(-(5*t)/2)*(cos((3^(1/2)*t)/2) +
(11*3^(1/2)*sin((3^(1/2)*t)/2))/15))/49 - 5/49
%Q2%
b=[5 3 6];
a=[1 3 7 9 12];
[r,p,k]=residue(b,a)
r =
-0.5357 - 1.0394i
11
-0.5357 + 1.0394i
0.5357 - 0.1856i
0.5357 + 0.1856i
p =
-1.5000 + 1.3229i
-1.5000 - 1.3229i
0.0000 + 1.7321i
0.0000 - 1.7321i
k =
[]
%Q3%
syms s t;
Ft1 = 7*t^3*cos(5*t+ pi*60/180 );
Ft2 = -7*t*exp(-5*t);
Ft3 = -3*cos(5*t);
Ft4 = t*sin(7*t);
Ft5 = 5*exp(-2*t)*cos(5*t);
Ft6 = 3*sin(5*t + 45 );
Ft7 = 5*exp(-3*t)*cos(t-45);
Fs1 = laplace(Ft1);
Fs1 =(7*3^(1/2)*((120*s)/(s^2 + 25)^3 - (240*s^3)/(s^2 + 25)^4))/2 +
21/(s^2 + 25)^2 - (168*s^2)/(s^2 + 25)^3 + (168*s^4)/(s^2 + 25)^4
Fs2 = laplace(Ft2)
Fs2 =-7/(s + 5)^2
Fs3 = laplace(Ft3)
Fs3 =-(3*s)/(s^2 + 25)
Fs4 = laplace(Ft4)
Fs4 =(14*s)/(s^2 + 49)^2
Fs5 = laplace(Ft5)
Fs5 =(5*(s + 2))/((s + 2)^2 + 25)
Fs6 = laplace(Ft6)
Fs6 =(3*(5*cos(45) + s*sin(45)))/(s^2 + 25)
Fs7 = laplace(Ft7)
Fs7 =(5*sin(45))/((s + 3)^2 + 1) + (5*cos(45)*(s + 3))/((s + 3)^2 + 1)
LAB 4
%Q1%
12
%Q1(a)%
Ta = tf(130,[1 15 130])
Ta =
130
----------------
s^2 + 15 s + 130
Continuous-time transfer function.
[W1,z1]=damp(Ta);
Wn=W1(1:1)
Wn = 11.4018
zeta=z1(1:1)
zeta = 0.6578
stepinfo(Ta)
ans =
RiseTime: 0.1758
SettlingTime: 0.5272
SettlingMin: 0.9037
SettlingMax: 1.0643
Overshoot: 6.4307
Undershoot: 0
Peak: 1.0643
PeakTime: 0.3684
%Q1(b)%
Tb = tf(0.045,[1 0.025 0.045])
Tb =
0.045
---------------------
s^2 + 0.025 s + 0.045
Continuous-time transfer function.
[W1,z1]=damp(Tb);
Wn=W1(1:1)
Wn = 0.2121
13
zeta=z1(1:1)
zeta = 0.0589
stepinfo(Tb)
ans =
RiseTime: 5.1423
SettlingTime: 312.2477
SettlingMin: 0.3099
SettlingMax: 1.8307
Overshoot: 83.0724
Undershoot: 0
Peak: 1.8307
PeakTime: 14.8096
%Q1(c)
Tc = tf(10^8,[1 1.325*(10^3) 10^8])
Tc =
1e08
-------------------
s^2 + 1325 s + 1e08
Continuous-time transfer function.
[W1,z1]=damp(Tc);
Wn=W1(1:1)
Wn = 10000
zeta=z1(1:1)
zeta = 0.0663
stepinfo(Tc)
ans =
RiseTime: 1.0972e-04
SettlingTime: 0.0057
SettlingMin: 0.3412
SettlingMax: 1.8117
Overshoot: 81.1710
Undershoot: 0
Peak: 1.8117
PeakTime: 3.1416e-04
%Q2%
G =zpk([],[-5 -7 -9 -11],150)
14
G =
150
------------------------
(s+5) (s+7) (s+9) (s+11)
Continuous-time zero/pole/gain model.
GH = feedback(G,1)
GH =
150
---------------------------------------------
(s^2 + 21.93s + 124.1) (s^2 + 10.07s + 29.13)
Continuous-time zero/pole/gain model.
sys=zpk(GH)
sys =
150
---------------------------------------------
(s^2 + 21.93s + 124.1) (s^2 + 10.07s + 29.13)
Continuous-time zero/pole/gain model.
rlocus(GH)
15
%Q3%
s=0:0.1:10;
num=30.*(s.^2 -5.*s + 3);
den=(s+1).*(s+2).*(s+3).*(s+5);
G=tf(num,den);
Gc=feedback(G,1);
step(G)
%Q4%
16
G = zpk([-1],[0 -1 -5 -6],1)
G =
(s+1)
-------------------
s (s+1) (s+5) (s+6)
Continuous-time zero/pole/gain model.
Gc = feedback(G,1)
Gc =
(s+1)
-------------------------------------
(s+6.142) (s+4.824) (s+1) (s+0.03375)
Continuous-time zero/pole/gain model.
rlocus(Gc)
%Q5%
Gc = zpk([-0.57 -0.57],[0],29.125)
17
Gc =
29.125 (s+0.57)^2
-----------------
s
Continuous-time zero/pole/gain model.
bode(Gc)
%Q6%
t=0:0.2:10;
zeta=[0 0.1 0.2 0.4 0.5 0.6 0.8 1.0];
for n=1:8;
num=[1];
den=[1 2*zeta(n) 1];
[y(1:51,n), x, t]= step(num,den,t);
18
plot(t,y)
end;
%Q7%
For the unity feedback system having forwardpath transfer function
Gs=tf(3,[1 5 9 0])
Gs =
3
-----------------
s^3 + 5 s^2 + 9 s
Continuous-time transfer function.
Gc=feedback(Gs,1)
Gc =
3
---------------------
s^3 + 5 s^2 + 9 s + 3
Continuous-time transfer function.
Gc1=zpk(Gc)
19
Gc1 =
3
---------------------------------
(s+0.4253) (s^2 + 4.575s + 7.055)
Continuous-time zero/pole/gain model.
rlocus(Gc)
20
LAB 5
DC Motor Position: System Modeling
Physical setup
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 equivalent circuit of the armature and the free-body
diagram of the rotor are shown in the following figure.
(J) moment of inertia of the rotor 3.2284E-6 kg.m^2
(b) motor viscous friction constant 3.5077E-6 N.m.s
(Kb) electromotive force constant 0.0274 V/rad/sec
(Kt) motor torque constant 0.0274 N.m/Amp
(R) electric resistance 4 Ohm
(L) electric inductance 2.75E-6H
21
In this example, we assume that the input of the system is the voltage source (V) applied to the motor's armature,
while the output is the position of the shaft (theta). 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.
1. Transfer Function
Applying the Laplace transform, the above modeling equations can be expressed in terms of the Laplace variable s.
(5)
(6)
We arrive at the following open-loop transfer function by eliminating I(s) between the two above equations, where the
rotational speed is considered the output and the armature voltage is considered the input.
(7)
However, during this example we will be looking at the position as the output. We can obtain the position by integrating
the speed, therefore, we just need to divide the above transfer function by s.
(8)
2. State-Space
The differential equations from above can also be expressed in state-space form by choosing the motor position,
motor speed and armature current as the state variables. Again the armature voltage is treated as the input and the
rotational position is chosen as the output.
(9)
(10)
22
Design requirements
We will want to be able to position the motor very precisely, thus the steady-state error of the motor position should
be zero when given a commanded position. We will also want the steady-state error due to a constant disturbance to
be zero as well. The other performance requirement is that the motor reaches its final position very quickly without
excessive overshoot. In this case, we want the system to have a settling time of 40 ms and an overshoot smaller than
16%.
If we simulate the reference input by a unit step input, then the motor position output should have:
Settling time less than 40 milliseconds
Overshoot less than 16%
No steady-state error, even in the presence of a step disturbance input
MATLAB representation
1. Transfer Function
We can represent the above open-loop transfer function of the motor in MATLAB by defining the parameters and
transfer function as follows. Running this code in the command window produces the output shown below.
J = 3.2284E-6;
b = 3.5077E-6;
K = 0.0274;
23
R = 4;
L = 2.75E-6;
s = tf('s');
P_motor = K/(s*((J*s+b)*(L*s+R)+K^2))
P_motor = 0.0274
-------------------------------------
8.878e-12 s^3 + 1.291e-05 s^2 + 0.0007648 s
Continuous-time transfer function.
2. State Space
We can also represent the system using the state-space equations. The following additional MATLAB commands
create a state-space model of the motor and produce the output shown below when run in the MATLAB command
window.
A = [0 1 0
0 -b/J K/J
0 -K/L -R/L];
B = [0 ; 0 ; 1/L];
C = [1 0 0];
D = [0];
24
motor_ss = ss(A,B,C,D)
motor_ss =
a =
x1 x2 x3
x1 0 1 0
x2 0 -1.087 8487
x3 0 -9964 -1.455e+06
b =
u1
x1 0
x2 0
x3 3.636e+05
c =
x1 x2 x3
y1 1 0 0
25
d =
u1
y1 0
Continuous-time state-space model.
The above state-space model can also be generated by converting your existing transfer function model into state-
space form. This is again accomplished with the ss command as shown below.
motor_ss = ss(P_motor);
LAB 6
For unity feedback control system, where system open loop transfer function is given by () =
16
. with the use of simulation and root locus in MATLAB, discuss about the effects of P control
(+1.6)
and PD control on system performance
num = [16];
denum = [1 1.6 0];
Gp = tf(num,denum)
Gp =
16
-----------
s^2 + 1.6 s
Continuous-time transfer function.
H = [1];
M = feedback(Gp,H)
M =
16
----------------
s^2 + 1.6 s + 16
26
Continuous-time transfer function.
%step(M)
hold on
Kp = 2;
Ki = 0;
Kd = 1;
Gc = pid(Kp, Ki, Kd)
Gc =
Kp + Kd * s
with Kp = 2, Kd = 1
Continuous-time PD controller in parallel form.
Mc = feedback(Gc * Gp,H)
Mc =
16 s + 32
-----------------
s^2 + 17.6 s + 32
Continuous-time transfer function.
step(Mc)
grid on
27
