1

EXP 1: Design and simulation of mathematical model for

a linear mechanical system

FORMULA:

The corresponding differential equations are as follows:

 2

𝑑2 𝑑 𝑑
𝑀2 2 (𝑦2) + 𝑓2 ( (𝑦2) − (𝑦1)) + 𝐾2(𝑦2 − 𝑦1) = 𝐹(𝑡)
𝑑𝑡 𝑑𝑡 𝑑𝑡

𝑑2 𝑑 𝑑 𝑑
𝑀1 (𝑦1) − 𝑓 2 ( (𝑦2) − (𝑦1)) − 𝐾 (𝑦
2 2 − 𝑦1 ) + 𝑓 1 ( (𝑦1)) + 𝐾1𝑦1 = 0
𝑑𝑡2 𝑑𝑡 𝑑𝑡 𝑑𝑡

Consider the parameters as given below:

M1 = M2 = 5 kg
K1 = 2 N/m
K2 = 6 N/m
f1 = 2N/m/s
f2 = 3N/m/s
The transfer function in s domain for the corresponding differential equation is:
𝑦2(𝑠) 5𝑠2 + 5𝑠 + 11
=
(𝑠) 25𝑠4 + 40𝑠3 + 91𝑠2 + 27𝑠 + 30

1. Simulink model:
 3

1. Simulink:
 4

MEASUREMENTS VALUES
60.484%
Overshoot
0.806
Preshoot
5.563
Undershoot
1.8585 seconds
Rise Time
1.848 seconds
Fall Time
45.5 seconds
Settling Time
0.347
Final Value
0.7317
Steady state Error

INFERENCE:

 The differential equations for a linear mechanical system were established and it
was realized using Simulink.
 The corresponding transfer function graph was plotted for step input, in Simulink
and MATLAB, and the properties were analyzed.
 The damping effect of the spring system was established, and the corresponding
data was tabulated.
 5

EXP 2: Design and simulation of mathematical model for

a non-linear mechanical (Ball and Beam) system

FORMULA/DERIVATION:
The Lagrange method is an energy-based approach for deriving the equations of
motion of a dynamical system. This is convenient since it does not require the use
of vectors. We will now derive the equations of motions for the Ball and Beam
using this method. Consider the Ball and Beam system:

The ball rolls on the beam without slipping under the action of the force of
gravity. The beam is tilted from an external torque to control the position of the
ball on the beam. We first define a set of generalized coordinates which fully
describe the system. Here, p(t) is the position of the ball and 𝜃(𝑡) is the angle of
the beam. The lagrangian is described as:
 6

L=K-U

where K is the kinetic energy and U is the potential energy of the system.

The kinetic energy of the beam is given by:

1
𝐾1 = 𝐽𝜃′2
2

Where J is the moment of inertia of the beam. The kinetic energy of the ball is:

1 ′2
1 2
𝐾2 = 𝐽𝑏𝜃𝑏 + 𝑚𝑣𝑏
2 2

where θ is the angular velocity of the ball and VB is the linear velocity of the ball.
The quantity θb can be expressed in terms of the generalized coordinates as

where r is the radius of the ball. We can also express vb in terms of the
generalized coordinates.
𝑣𝑏2 = 𝑥2 + 𝑦2
 7

𝑥 = 𝑝 𝑐𝑜𝑠 𝜃

ⅆ𝑥 ⅆ𝑝
= ( ) cos 𝜃 − 𝑝𝜃 sin 𝜃
ⅆ𝑡 ⅆ𝑡

ⅆ𝑥 2 ⅆ𝑝 2 ⅆ𝑝
( ) = ( ) cos2 𝜃 − 2𝑝 ( ) 𝜃 cos 𝜃 sin 𝜃 + 𝑝2𝜃2 sin2 𝜃
ⅆ𝑡 ⅆ𝑡 ⅆ𝑡

𝑦 = 𝑝 𝑠𝑖𝑛 𝜃

ⅆ𝑦 ⅆ𝑝
( )=( ) 𝑠𝑖𝑛 𝜃 + 𝑝𝜃 𝑐𝑜𝑠 𝜃
ⅆ𝑡 ⅆ𝑡

𝑑𝑦 2 𝑑𝑝 2 𝑑𝑝
( ) = ( ) sin2 𝜃 + 2𝑝 ( ) cos 𝜃 sin 𝜃 + 𝑝2𝜃2 cos2 𝜃
𝑑𝑡 𝑑𝑡 𝑑𝑡

𝑑𝑝 2
𝑣2
𝑏 = ( ) + 𝑝2 𝜃2
𝑑𝑡
Substituting for kinetic energy,

1
𝑘 = 1 ( 𝐽𝐵 2 2 2
+ 𝑚𝑝 𝐴
2 2
+ 𝑚) 𝑝 2
2 𝑟

The potential energy is given by:

𝑈 = 𝑚𝑔𝑝𝑠𝑖𝑛𝜃
 8

The resultant Lagrangian for this system is,

1 𝐽𝑏 1
2 + 𝑚𝑝2 2
𝐿= ( + 𝑚) 𝑝 𝐴 − 𝑚𝑔𝑝𝑠𝑖𝑛𝜃
2 𝑟2 2

The first Lagrange is given by:

𝑑 𝜕𝐿 𝜕𝐿
( ) − =𝑈
𝑑𝑡 𝜕𝑝′ 𝜕𝑝

𝜕𝐿 𝐽𝐵
= ( + 𝑚) 𝑝′′
𝜕𝑡 𝑟2

𝜕𝐿
= 𝑚𝑝𝜃 − 𝑚𝑔𝑠𝑖𝑛𝜃
𝜕𝑝

Thus, the equation of the ball beam system is given by:

𝑱𝑩
( + 𝒎) 𝒑′′ + 𝒎𝒈𝒔𝒊𝒏𝜽 − 𝒎𝒑𝜽𝟐 = 𝟎
𝒓𝟐

BLOCK DIAGRAM/CODE:
 9

Here the lever is connected to a servo gear. Thus, as the servo angle changes by
theta, the lever changes by an angle alpha.
A controller will be designed for this system so that the ball's position can be
manipulated.
𝐽𝑏
( + 𝑚) 𝑟′′ + 𝑚𝑔𝑠𝑖𝑛𝛼 − 𝑚𝑟𝛼′2 = 0
𝑅2

𝑑
𝛼= 𝜃
𝐿
Linearization of this equation about the beam angle, alpha = 0, gives us the
following linear approximation of the system:

𝐽 + 𝑚) 𝑟′′ = − 𝑚𝑔𝑑 𝜃
( 2
𝑅 𝐿

 r - Distance of the ball from the center.

 α – angle made by the beam with respect to the horizontal
 R – radius of the ball
 JB – Moment of inertia of the ball

1. MATLAB initialization parameters and Code:

Parameter initialization in MATLAB:

a) m=0.1 kg.
b) R=0.01 m (radius of the ball).
c) g=-9.8 m/s2
d) L=1 m (length of beam).
e) d=0.01 m (level arm offset).
f) J=4 * 10-6 kg/m2(Moment of inertia of the ball).
 10

Code:
>>m = 0.1;
>>R = 0.01;
>>g = -9.8;
>>L = 1;
>>d = 0.01;
>>J = 4e-6;

>>K = (m*g*d)/(L*(J/R^2+m)); %simplifies input

>>num = [-K];
>>den = [1 0 0];
>>step (0.25*num, den);

2. Simulink model:
 11

OBSERVATIONS:

INFERENCE:

 We can observe that as θ increases, beam becomes steeper and as a result, the
system becomes highly unstable;
 In addition to the existing model, we can design an external controller which
helps in maintaining the θ such that the stability of the system is maintained.
 A lead compensator can be designed which prevents the ball from
rolling off the beam and hence making the system stable.
 12

EXP 3: Design and simulation of mathematical model for

Electrical system
1. Series RLC circuit

FORMULA/DERIVATION:

Here, let:
 L = inductance (H)
 R = resistance ( R)
 C = capacitance(F)
 I(t) = current through the circuit(clockwise)(A)
 q = charge
 Vs = Unit step input.
 13

The current is moving in a clockwise manner. Using Kirchhoff’s voltage law:

𝐿𝑑(𝑡) 1
𝑉𝑠 (𝑡) − − (𝑡)𝑅 − ∫ 𝑖(𝑡)𝑑𝑡 = 0
𝑑𝑡 𝐶
𝑑𝑞(𝑡)
Here we know that (𝑡) = ,
𝑑𝑡

𝐿𝑑 𝑑𝑞(𝑡) 𝑅𝑑(𝑡) 1
𝑉(𝑡)𝑠 = ( )+ + ∫ 𝑑𝑞
𝑑𝑡 𝑑𝑡 𝑑𝑡 𝐶

𝐿𝑑2(𝑡) 𝑅𝑑𝑞 𝑞
𝑉𝑠 (𝑡) = + +
𝑑𝑡2 𝑑𝑡 𝐶
𝑑2(𝑡) 𝑑𝑞(𝑡) 1
𝐿 = 𝑅( ) + 𝑞(𝑡) = 𝑉𝑠 (𝑡)
𝑑𝑡2 𝑡 𝐶

Taking Laplace transform, assuming Vs = u(t) [unit step response];

1
𝐿𝑠2(𝑠) + 𝑅𝑠𝑄(𝑠) + 𝑄(𝑠) = 𝑉 (𝑠)
𝑠
𝐶

𝑄(𝑠) = (1/𝐿)
𝑉 (𝑠) 𝑅𝑠 1
𝑠 𝑠2 + 𝐿 + 𝐿𝐶
 14

BLOCK DIAGRAM/CODE:

1. MATLAB initialization parameters and Code:

Parameter initialization in MATLAB:

Parameters:
 R=3 Ohm
 L=1 H
 C= 0.05 F
Code:
>>R=3;
>>L=1;
>>C=0.05;

>>num=[1./L];
>>den=[1 (R/L) 1./(L*C) ];
>>G=tf(num,den);
>>step(G);

>>[w,zeta,p]=damp(G); %damp() gives the frequency domain

parameters
>>tau=1./w; %time constant

2. Simulink model:
 15

OBSERVATIONS:

1. MATLAB code:

2. Simulink:
 16

MEASUREMENTS VALUES
30.921%
Overshoot
0.658%
Pre shoot
9.302%
Undershoot
317.030 ms
Rise Time

INFERENCE:

 We know that current lags in the inductor, leads in capacitor and is in phase in
resistor.
 Here, the graph is divided into transient and steady state regions, where the
transient state is influenced by initial conditions(R,L,C), while the steady state
component remains unchanged.
 Hence from the graphs, we can observe the underdamping of the wave,
before reaching a steady state value. By changing the values of R, L and C, we
can change the damping ratio and in turn the response as well.
 17

EXP 3: Design and simulation of mathematical model for

Electrical system

2. Parallel RLC circuit

BLOCK DIAGRAM/CODE:

1. MATLAB initialization parameters and Code:

Parameter initialization in MATLAB:

Parameters:
 R=2 Ohm
 L=1 H
 C= 0.5 F
Code:
>>R=2;
>>L=1
>>C=0.5
>>num=[0 (R*L) 0 ];
>>den=[(R*L*C) L R];
>>G=tf(num,den);
>>step(G);
>>[w,zeta,p]=damp(G); %damp() gives the frequency domain
parameters
>>tau=1. /w; %time constant

2. Simulink model:
 18

OBSERVATIONS:

1. MATLAB code:

2. Simulink:
 19

MEASUREMENTS VALUES
6.956%
Overshoot
0.658%
Pre shoot
30.921%
Undershoot
545.708 ms
Rise Time

INFERENCE:

 We know that current lags in the inductor, leads in capacitor and is in phase in
resistor.
 Here, we can see that a sinusoidal, gradually declining response signal is initially
observed, after which it attains steady state.
 From the graph, we can observe the underdamped response of the parallel RLC
circuit. By changing the values of R, L and C, we can change the damping ratio
and in turn the response as well.
 20

EXP 4: System Identification through Process Reaction

Curves

FORMULA/THEORY:

To generate a process reaction curve, the process can reach steady state or as
close to steady state as possible. Then, in open loop, so that there is no control
action, a small step disturbance is introduced, and the reaction of the process
variable is b recorded.

For an oscillatory reaction curve, the dimensionless equation for model 1 will
be adopted. To distinguish between process models that do not have a process
zero (where 𝑎̅ = 0) and those which do have a zero (where 𝑎̅≠ 0) the following
lemma is needed.
21
22
 23

BLOCK DIAGRAM:

1. Simulink model:

OBSERVATIONS:

1. Simulink:
Process A:

Process B:
 24

[Link] Parameters Process A Process B

tm,0 2.897s 3.297s
1.
tp,1 10.33s 8.531s
2.
yp,1 1.206 1.385
3.
tp,2 24.73s 22.73s
4.
yp,2 1.005 1.01
5.
ym,1 0.97097 0.95816
6.
p 14.4 14.3
7.

CALCULATIONS:

1. Process A:

1 1.005 − 1
𝜒= 𝑙𝑛 ( )
10.33 − 24.73 1.206 − 1

𝝌 = 𝟎. 𝟐𝟓𝟖

ln(1.206 − 1)2
𝜍= √
𝜋2 + (𝑙𝑛1.206 − 1)2

𝝇 = 𝟎. 𝟒𝟒𝟗

14.4
𝑟=
√4𝜋2 + 14.42 ∗ 0.2582

𝑟 = 𝟏. 𝟗𝟕𝟐𝟕
 25

1.9727
𝜃 = 10.33 + ln(1.206 − 1)
0.449

𝜽 = 𝟑. 𝟑𝟗

2. Process B:

1 1.01 − 1
𝜒= 𝑙𝑛 ( )
8.331 − 22.73 1.385 − 1

𝝌 = 𝟎. 𝟐𝟓𝟕

14.199 ∗ 0.257
𝜍=
√4𝜋2 + 14.1992 ∗ 0.2572

𝝇 = 𝟎. 𝟓

14.199
𝑟=
√4𝜋2 + 14.1992 ∗ 0.2572

𝑟 = 𝟏. 𝟗𝟓𝟒

1.385 − 1 2
𝑎′ = 0.5 + √0.52 + (1 − ( ) )
𝑒−0.5∗8.531

𝒂′ = 𝟏. 𝟓𝟑𝟗

14.199 1.539√1 − 0.52

𝜃 = 8.537 − (π − tan−1 ( ))
6.28 1 − 1.539 ∗ 0.5
𝜽 = 𝟒. 𝟔𝟓𝟕
 26

[Link] Parameters Process A Process B

1. ς 0.449 0.5

2. χ 0.258 0.257

3. τ 1.9727 1.954

4. a’ 0 1.539

5. θ 3.39 4.657

Thus, the equation of process A is:

𝑒−3.39𝑠
3.88𝑠2 + 1.771𝑠 + 1
The equation of process B is:
(1 + 1.539𝑠)𝑒−4.6575
3.818𝑠2 + 1.954𝑠 + 1

INFERENCE:
 From the graphs, we can analyze the process curve and obtain the necessary
parameters.
 The damping ratio for both the processes reinforce the fact that the curves are
underdamped.
 The delay time θ is higher for process B and it is indicative for the fact that
process reaction curve for process B takes longer than that of A.
 27

EXP 5: Stability Analysis using Pole Zero maps and

Routh Horwitz Criterion in simulation platform

THEORY:

The poles and zeros of the equation can be plotted on the s-plane, which is
known as the pole zero map.

Thus, if all the roots are placed in the left side of s-plane, then the system is
stable. Even if some roots are placed on the right side of the s-plane, then the
system becomes unstable.

Here, we take 3 different equations and obtain the characteristic equations. We

will model a MATLAB program and plug in the values.

CALCULATION:
𝟐𝟓
1. (𝒔) =
𝒔𝟐+𝟔𝒔

The characteristic equation is: 𝑠2 + 6𝑠 + 25 = 0

s2 1 25
6
s1
25
s0
Since there are no sign changes in the first column of the Routh array, the
system is ‘Stable’. Both roots lie on the left half of s-plane.

Number of Zeros: 0
Number of Roots: 2 [0, -6]
 28

𝟖
2. (𝒔) =
𝐬𝟐−𝟒𝐬+𝟖

The characteristic equation is: 𝑠2 − 4𝑠 + 16 = 0

s2 1 16
-4
s1
16
s0
Since there are two sign changes in the first column of the Routh array, the
system is ‘Unstable’ and there are two roots present in the right half of the
s-plane.

Solving the denominator equation gives the poles as follows:

(4 ± √16 − 4(8))
𝑠=
2
 29

𝑠 = 2 ± 2𝑗

Number of zeros: 0
Number of Roots: 2 (2 ± 2𝑗)

𝟏𝟔
3. (𝒔) =
𝒔𝟒+𝟖𝒔𝟐+𝟏𝟔

The characteristic equation is: 𝑠4 + 0𝑠3 + 8𝑠2 + 0𝑠 + 32 = 0

s4 1 8 32
0 0 0
s3
Since there are a row of zeros, we will consider an auxiliary equation,
differentiate wrt s and use those coefficient values.

𝑑𝐴
= 4𝑠3 + 16𝑠 = 0
𝑑𝑠

s4 1 8 32
1 4
s3
1 8
s2
-4
s1
8
s0
Since there are two sign changes in the first column of the Routh array, the
system is ‘Unstable’ and there are two roots present in the right half of the
s-plane, while the others lie on the left half of s plane.
Solving the denominator equation gives:

𝑠 = 2𝑗, 2𝑗, −2𝑗, −2𝑗

 30

Number of zeros: 0
Number of Roots: 4 (2𝑗, 2𝑗, −2𝑗, −2𝑗)

MATLAB CODE:

𝟐𝟓
a. (𝒔) = 𝒔𝟐+𝟔𝒔

clear;
close all;
clc;
format compact;

% Taking coefficients vector and organizing the first two rows

coeffVector = input (‘input vector of your system coefficients: \n i.e. [an an-
1 an-2 ... a0] = ‘);
ceoffLength = length(coeffVector);
rhTableColumn = round(ceoffLength/2);

% Initialize Routh-Hurwitz table with empty zero array

rhTable = zeros(ceoffLength,rhTableColumn);

% Compute first row of the table

rhTable (1,: ) = coeffVector(1,1:2:ceoffLength);

% Check if length of coefficients vector is even or odd

if (rem(ceoffLength,2) ~= 0)
% if odd, second row of table will be
rhTable(2,1:rhTableColumn – 1) = coeffVector(1,2:2:ceoffLength);
else
% if even, second row of table will be
rhTable (2,: ) = coeffVector(1,2:2:ceoffLength);
end
 31

%% Calculate Routh-Hurwitz table’s rows

% Set epss as a small value

epss = 0.01;

% Calculate other elements of the table

for I = 3:ceoffLength

% special case: row of all zeros

if rhTable(i-1, : ) == 0
order = (ceoffLength – i);
cnt1 = 0;
cnt2 = 1;
for j = 1: rhTableColumn – 1
rhTable(i-1,j) = (order – cnt1) * rhTable(i-2,cnt2);
cnt2 = cnt2 + 1;
cnt1 = cnt1 + 2;
end
end

for j = 1:rhTableColumn – 1
% first element of upper row
firstElemUpperRow = rhTable(i-1,1);

% compute each element of the table

rhTable(I,j) = ((rhTable(i-1,1) * rhTable(i-2,j+1)) - ....
(rhTable(i-2,1) * rhTable(i-1,j+1))) / firstElemUpperRow;
end

% special case: zero in the first column

if rhTable(I,1) == 0
rhTable(I,1) = epss;
end
end
%% Compute number of right hand side poles(unstable poles)

% Initialize unstable poles with zero

 32

unstablePoles = 0;

% Check change in signs

for I = 1:ceoffLength – 1
if sign(rhTable(I,1)) * sign(rhTable(i+1,1)) == -1
unstablePoles = unstablePoles + 1;
end
end

% Print calculated data on screen

fprintf (‘\n Routh-Hurwitz Table:\n’)

rhTable

% Print the stability result on screen

if unstablePoles == 0
fprintf (‘~~~~~> it is a stable system! <~~~~~\n’)
else
fprintf (‘~~~~~> it is an unstable system! <~~~~~\n’)
end
fprintf (‘\n Number of right hand side poles =%2.0f\n’,unstablePoles)
reply = input (‘Do you want roots of system be shown? Y/N ‘, ‘s’);
if reply == ‘y’ || reply == ‘Y’
sysRoots = roots(coeffVector);
fprintf (‘\n Given polynomial coefficients roots :\n’)
sysRoots
end

H1=tf([0 0 25],[1 6 0])

pzmap(H1)
grid on
 33

𝟖
b. (𝒔) = 𝒔𝟐−𝟒𝒔+𝟖

clear ;
close all;
clc;
format compact;

% Taking coefficients vector and organizing the first two rows

coeffVector = input('input vector of your system coefficients: \n i.e. [an an-

1 an-2 ... a0] = ');
ceoffLength = length(coeffVector);
rhTableColumn = round(ceoffLength/2);

% Initialize Routh-Hurwitz table with empty zero array

rhTable = zeros(ceoffLength,rhTableColumn);

% Compute first row of the table

rhTable(1,:) = coeffVector(1,1:2:ceoffLength);

% Check if length of coefficients vector is even or odd

if (rem(ceoffLength,2) ~= 0)
% if odd, second row of table will be
rhTable(2,1:rhTableColumn - 1) = coeffVector(1,2:2:ceoffLength);
else
% if even, second row of table will be
rhTable(2,:) = coeffVector(1,2:2:ceoffLength);
end
%% Calculate Routh-Hurwitz table's rows

% Set epss as a small value

epss = 0.01;

% Calculate other elements of the table

 34

for i = 3:ceoffLength

% special case: row of all zeros

if rhTable(i-1,:) == 0
order = (ceoffLength - i);
cnt1 = 0;
cnt2 = 1;
for j = 1:rhTableColumn - 1
rhTable(i-1,j) = (order - cnt1) * rhTable(i-2,cnt2);
cnt2 = cnt2 + 1;
cnt1 = cnt1 + 2;
end
end

for j = 1:rhTableColumn - 1
% first element of upper row
firstElemUpperRow = rhTable(i-1,1);

% compute each element of the table

rhTable(i,j) = ((rhTable(i-1,1) * rhTable(i-2,j+1)) - ....
(rhTable(i-2,1) * rhTable(i-1,j+1))) / firstElemUpperRow;
end

% special case: zero in the first column

if rhTable(i,1) == 0
rhTable(i,1) = epss;
end
end
%% Compute number of right hand side poles(unstable poles)

% Initialize unstable poles with zero

unstablePoles = 0;

% Check change in signs

for i = 1:ceoffLength - 1
if sign(rhTable(i,1)) * sign(rhTable(i+1,1)) == -1
 35

unstablePoles = unstablePoles + 1;
end
end

% Print calculated data on screen

fprintf('\n Routh-Hurwitz Table:\n')

rhTable

% Print the stability result on screen

if unstablePoles == 0
fprintf('~~~~~> it is a stable system! <~~~~~\n')
else
fprintf('~~~~~> it is an unstable system! <~~~~~\n')
end
fprintf('\n Number of right hand side poles =%2.0f\n',unstablePoles)
reply = input('Do you want roots of system be shown? Y/N ', 's');
if reply == 'y' || reply == 'Y'
sysRoots = roots(coeffVector);
fprintf('\n Given polynomial coefficients roots :\n')
sysRoots
end

H1=tf([0 0 8],[1 -4 8])

pzmap(H1)
grid on

𝟏𝟔
c. (𝒔) = 𝒔𝟒+𝟖𝒔𝟐+𝟏𝟔

clear ;
close all;
clc;
format compact;

% Taking coefficients vector and organizing the first two rows

 36

coeffVector = input('input vector of your system coefficients: \n i.e. [an an-

1 an-2 ... a0] = ');
ceoffLength = length(coeffVector);
rhTableColumn = round(ceoffLength/2);

% Initialize Routh-Hurwitz table with empty zero array

rhTable = zeros(ceoffLength,rhTableColumn);

% Compute first row of the table

rhTable(1,:) = coeffVector(1,1:2:ceoffLength);

% Check if length of coefficients vector is even or odd

if (rem(ceoffLength,2) ~= 0)
% if odd, second row of table will be
rhTable(2,1:rhTableColumn - 1) = coeffVector(1,2:2:ceoffLength);
else
% if even, second row of table will be
rhTable(2,:) = coeffVector(1,2:2:ceoffLength);
end
%% Calculate Routh-Hurwitz table's rows

% Set epss as a small value

epss = 0.01;

% Calculate other elements of the table

for i = 3:ceoffLength

% special case: row of all zeros

if rhTable(i-1,:) == 0
order = (ceoffLength - i);
cnt1 = 0;
cnt2 = 1;
for j = 1:rhTableColumn - 1
rhTable(i-1,j) = (order - cnt1) * rhTable(i-2,cnt2);
cnt2 = cnt2 + 1;
 37

cnt1 = cnt1 + 2;
end
end

for j = 1:rhTableColumn - 1
% first element of upper row
firstElemUpperRow = rhTable(i-1,1);

% compute each element of the table

rhTable(i,j) = ((rhTable(i-1,1) * rhTable(i-2,j+1)) - ....
(rhTable(i-2,1) * rhTable(i-1,j+1))) / firstElemUpperRow;
end

% special case: zero in the first column

if rhTable(i,1) == 0
rhTable(i,1) = epss;
end
end
%% Compute number of right hand side poles(unstable poles)
% Initialize unstable poles with zero

unstablePoles = 0;

% Check change in signs

for i = 1:ceoffLength - 1
if sign(rhTable(i,1)) * sign(rhTable(i+1,1)) == -1
unstablePoles = unstablePoles + 1;
end
end

% Print calculated data on screen

fprintf('\n Routh-Hurwitz Table:\n')

rhTable

% Print the stability result on screen

if unstablePoles == 0
fprintf('~~~~~> it is a stable system! <~~~~~\n')
 38

else
fprintf('~~~~~> it is an unstable system! <~~~~~\n')
end
fprintf('\n Number of right hand side poles =%2.0f\n',unstablePoles)
reply = input('Do you want roots of system be shown? Y/N ', 's');
if reply == 'y' || reply == 'Y'
sysRoots = roots(coeffVector);
fprintf('\n Given polynomial coefficients roots :\n')
sysRoots
end

H1=tf([16],[1 0 8 0 16])
pzmap(H1)
grid on

OBSERVATIONS:

1. Code:

𝟐𝟓
a. (𝒔) = 𝒔𝟐+𝟔𝒔

input vector of your system coefficients:

i.e. [an an-1 an-2 ... a0] = [1 6 25]

Routh-Hurwitz Table:
rhTable =
1 25
6 0
25 0
~~~~~> it is a stable system! <~~~~~

Number of right hand side poles = 0

Do you want roots of system be shown? Y/N y

Given polynomial coefficients roots :

sysRoots =
-3.0000 + 4.0000i
 39

-3.0000 - 4.0000i

H1 =

25
--- - -
s^2 + 6 s

Continuous-time transfer function.

𝟖
b. (𝒔) = 𝒔𝟐−𝟒𝒔+𝟖

input vector of your system coefficients:

i.e. [an an-1 an-2 ... a0] = [1 -4 16]

Routh-Hurwitz Table:
rhTable =
1 16
-4 0
16 0
~~~~~> it is an unstable system! <~~~~~
 40

Number of right-hand side poles = 2

Do you want roots of system be shown? Y/N Y

Given polynomial coefficients roots:

sysRoots =
2.0000 + 3.4641i
2.0000 - 3.4641i

H1 =

8
--- - -
s^2 - 4 s + 8

Continuous-time transfer function.

𝟏𝟔
c. (𝐬) = 𝐬𝟒+𝟖𝐬𝟐+𝟏𝟔

input vector of your system coefficients:

i.e. [an an-1 an-2 ... a0] = [1 0 8 0 32]
 41

Routh-Hurwitz Table:
rhTable =
1 8 32
2 0 0
8 32 0
-8 0 0
32 0 0
~~~~~> it is an unstable system! <~~~~~

Number of right hand side poles = 2

Do you want roots of system be shown? Y/N y

Given polynomial coefficients roots :

sysRoots =
-0.9102 + 2.1974i
-0.9102 - 2.1974i
0.9102 + 2.1974i
0.9102 - 2.1974i

H1 =

16
--- - --
s^4 + 8 s^2 + 16

Continuous-time transfer function.

 42

INFERENCE:

 For a stable system, the roots of the characteristic equation must have negative
real parts.
 If any root of the characteristic equation has a positive real part or if there is a
repeated root on the imaginary axis, then the system is unstable.
 If all the roots of the characteristic equation have negative real parts except for
the presence of one or more non repeated roots on the imaginary axis, then the
system is limitedly or marginally stable.
 The pole zero map also reinforces the above conditions, in a graphical manner.
 43

EXP 6: Root Locus Based Analysis in Simulation Platform

THEORY:
The root locus technique in control system was first introduced in the
year 1948 by Evans. Any physical system is represented by a transfer
function in the form of

We can find poles and zeros from G(s). The location of poles and zeros are crucial
keeping view stability, relative stability, transient response, and error analysis.
When the system is put to service stray inductance and capacitance get into the
system, thus changes the location of poles and zeros.

In root locus technique in control system, we will evaluate the position of the
roots, their locus of movement and associated information. This information will
be used to comment upon the system performance.

Some of the advantages of root locus technique are:

1. Root locus technique in control system is easy to implement as

compared to other methods.
2. With the help of root locus, we can easily predict the performance of
the whole system.
3. Root locus provides the better way to indicate the parameters.
 44

Now there are various terms related to root locus technique that we will use
frequently:
1. Characteristic Equation Related to Root Locus
Technique:

1 + G(s)H(s) = 0 is known as characteristic equation.

Now on differentiating the characteristic equation
and on equating dk/ds equals to zero, we can get
break away points.

2. Break away Points:

Suppose two root loci which start from pole and moves
in opposite direction collide with each other such that
after collision they start moving in different directions in
the symmetrical way. Or the breakaway points at which
multiple roots of the characteristic equation 1 + G(s)H(s)
= 0 occur. The value of K is maximum at the points
where the branches of root loci break away. Break away
points may be real, imaginary, or complex.

3. Break in Point:

Condition of break in to be there on the plot is

written below: Root locus must be present between
two adjacent zeros on the real axis.

4. Centre of Gravity:

It is also known centroid and is defined as the point on the

plot from where all the asymptotes start. Mathematically, it
is calculated by the difference of summation of poles and
zeros in the transfer function when divided by the
difference of total number of poles and total number of
zeros. Centre of gravity is always real and it is denoted by
σA.

Where, N is number of poles and M is number of zeros.

 45

5. Asymptotes of Root Loci:

Asymptote originates from the center of gravity or

centroid and goes to infinity at definite some angle.
Asymptotes provide direction to the root locus when
they depart break away points.

6. Angle of Asymptotes:

Asymptotes makes some angle with the real axis and

this angle can be calculated from the given formula,

Where, p = 0, 1, 2 ……. (N-M-1) N is the total number of poles M is

the total number of zeros.

7. Angle of Arrival or Departure:

We calculate angle of departure when there exist

complex poles in the system. Angle of departure can
be calculated as 180-{(sum of angles to a complex pole
from the other poles)-(sum of angle to a complex pole
from the zeros)}.

8. Intersection of Root Locus with the Imaginary Axis:

To find out the point of intersection root locus with

imaginary axis, we must use Routh Hurwitz criterion.
First, we find the auxiliary equation then the
corresponding value of K will give the value of the
point of intersection.

9. Gain Margin:

We define gain margin by which the design value of the

gain factor can be multiplied before the system becomes
unstable. Mathematically it is given by the formula
 46

10. Phase Margin:

Phase margin can be calculated from the given formula:

11. Symmetry of Root Locus:

Root locus is symmetric about the x axis or the real axis

Now there are two ways of determining the value of K, each

way is described below.

1. Magnitude Criteria:

At any points on the root locus, we can apply

magnitude criteria as,

Using this formula, we can calculate the value of K at

any desired point.

2. Using Root Locus Plot:

The value of K at any s on the root locus is given by:

Root Locus Plot:

This is also known as root locus technique in control system and is
used for determining the stability of the given system. Now to
 47

determine the stability of the system using the root locus technique
we find the range of values of K for which the complete performance
of the system will be satisfactory, and the operation is stable.

Now there are some results that one should remember to plot the
root locus. These results are written below:

1. Region where root locus exists:

After plotting all the poles and zeros on the plane, we can
easily find out the region of existence of the root locus by
using one simple rule which is written below, only that
segment will be considered in making root locus if the
total number of poles and zeros at the right hand side of
the segment is odd.

2. How to calculate the number of separate root loci:

Several separate root loci are equal to the total number

of roots if number of roots are greater than the number
of poles otherwise number of separate root loci is equal
to the total number of poles if number of roots are
greater than the number of zeros.

Procedure to Plot Root Locus:

Keeping all these points in mind we are able to draw the root locus plot for any
kind of system. Now let us discuss the procedure of making a root locus:

 Find out all the roots and poles from the open loop transfer function and
then plot them on the complex plane.

 All the root loci start from the poles where k = 0 and terminates at the
zeros where K tends to infinity. The number of branches terminating at
infinity equals to the difference between the number of poles & number
of zeros of G(s)H(s).
 48

 Find the region of existence of the root loci from the method described
above after finding the values of M and N.

 Calculate break away points and break in points if any.

 Plot the asymptotes and centroid point on the complex plane for the root
loci by calculating the slope of the asymptotes.

 Now calculate angle of departure and the intersection of root loci with
imaginary axis.

 Now determine the value of K by using any one method that I have
described above. By following above procedure, you can easily draw the
root locus plot for any open loop transfer function.

 Calculate the gain margin.

 Calculate the phase margin.

 You can easily comment on the stability of the system by using Routh
Array.

MATLAB CODE:

1. Code:
1
1. 𝑠3+6𝑠2+45𝑠
 49

0.075𝑠2+𝑠+1
2.
𝑠3+3𝑠2+5𝑠

OBSERVATIONS:

1. Code:

1
1.
𝑠3+6𝑠2+45𝑠
 50

0.075𝑠2+𝑠+1
2.
𝑠3+3𝑠2+5𝑠

INFERENCE:
 A slight change in the pole-zero configuration may cause significant
changes in the root- locus configurations.
 The cancelled pole of G(s)H(s) is a closed-loop pole of the system and this must
be added to the closed-loop poles obtained from the root-locus plot of G(s)H(s).
 General effects of the addition of zeros:
o pull root locus to the left
o makes system more stable
o speed up the settling of the response.

 General effects of the addition of poles:

 pull root locus to the right
 lower system’s relative stability
 slow down the settling of the response
 51

EXP 7: Frequency Response Analysis of a 3rd Order

System Using Bode Plot

THEORY:

The stability of linear systems can be assessed by checking the poles ofthe
transfer function of the system, i.e., the roots of the characteristic equation in
the s-plane. A linear system is stable if all the poles lie in theleft-half s-plane.
However, merely knowing whether a system is stable or unstable is not
sufficient. We need to find out how stable the system is.

A Bode plot is a graph commonly used in control system engineering to

determine the stability of a control system. A Bode plot maps the frequency
response of the system through two graphs – the Bode magnitude plot
(expressingthe magnitude in decibels) and the Bode phase plot (expressing the
phase shift in degrees). We can measure the relative stability of a systemby
examining the gain margin and the phase margin obtained from
the Bode plot of the frequency response.

Phase margin and gain margin are the measures of closed-loop control
systems stability. They are illustrated in the figure below.
 52

1. Gain margin:

Gain margin is defined as the amount of change in open-loop gain needed

to make a closed-loop system unstable.

 The gain margin is the difference between 0 dB and thegain at the

phase cross-over frequency that gives a phase of −180°. If the gain
|GH(jω) | at the frequency of ∠GH(jω)=−180° is greater than 0 dB as
shown in the left of Fig. 2.25, meaning a positive gain margin, then
theclosed-loop system is stable.

 We can usually read the gain margin directly from the Bodeplot (as
shown in the diagram above).

 This is done by calculating the vertical distance between the magnitude

curve (on the Bode magnitude plot) and the x-axisat the frequency
where the Bode phase plot = 180°. This point is known as the phase
crossover frequency.

 It is important to realize that the Gain and the Gain Margin arenot the
same things. In fact, the Gain Margin is the negative of the gain (in
decibels, dB). This will make sense when we look at the Gain margin
formula.

 The formula for Gain Margin (GM) can be expressed as:

Where G is the gain. This is the magnitude (in dB) as read fromthe vertical
axis of the magnitude plot at the phase crossover frequency.

2. Phase margin:

 Phase margin is defined as the amount of change in open-loop phase

needed to make a closed-loop systemunstable.

 The phase margin is the difference in phase between −180° and the
phase at the gain cross-over frequencythat gives a gain of 0 dB.
 53

 We can usually read the phase margin directly from the Bode plot (as
shown in the diagram above). This is doneby calculating the vertical
distance between the phase curve (on the Bode phase plot) and the x-
axis at the frequency where the Bode magnitude plot = 0 dB. This point
is known as the gain crossover frequency.

 It is important to realize that the phase lag and the Phase Margin
are not the same things. This will makesense when we look at the
phase margin formula.

 The formula for Phase Margin (PM) can be expressed as:

 Where is the phase lag (a number less than 0). This is thephase as read
from the vertical axis of the phase plot at thegain crossover frequency.

In general, as the gain of a system increases, the system becomes lessstable.

The gain margin and the phase margin indicate how much the gain increases
until the system becomes unstable.
Let us look at the relationship between the gain and phase margins and the
stability of a system. For example, assume that the phase margin of a system
is zero degree, i.e., ∠GH(jω)=180° at a frequency ω of |G(jω)|=1dB. In this
case, as shown in Fig. 2.26, the phase difference between Y(s) and R(s) is 180°,
and this closed-loopsystem thus becomes a positive feedback system. This
system will continue to amplify its output and become unstable. Example of
an unstable system is illustrated below.

If a closed-loop system is stable, both the gain margin and the phase margin
need to be positive. In general, the phase margin of 30–60 degrees and the
gain margin of 2–10 dB is desirable in theclosed-loop system design. A
system with a large gain margin and phase margin is stable but has a sluggish
response, while the one with a small gain margin and phase margin has a less
sluggish response but is oscillatory.

For the graphs below the stabilities are analyzed

 54

A) Bode plot A

B) Bode plot B
 55

TABULATION:

[Link] Parameters Bode Plot A Bode Plot B

Gain Margin 48Db infinite
1.
𝜔𝑐𝑔 6.7 rad/s -
2.
Phase Margin 89.8° 94.5°
3.
𝜔𝑐𝑝 0.022rad/s 0.203rad/s
4.

INFERENCE:

 For both the system the phase starts at -90 and hence there is a pole at the
origin for both the systems.

 From the tabular column Gain margin for A is 48dB and that for B it is infinity.

Generally, for stability Gain margin is expected to be more.

 Since for the magnitude Bode plot B the Gain margin is infinity it is inferred
that it is a very stable system, but the response is sluggish with no
oscillations. It is a very stable system because infinite gain can be increased
or decreased, and the stability of the system won’t be affected

 For the bode plot B 48db can be increased or decreased till the system
becomes unstable.

 Greater the phase margin greater the stability of the system. The phase
margin of B is > than the phase margin of A.

 The greater stability of B over A is because B can tolerate a maximum

phase change of 94 degrees without making the system unstable
whereas A can tolerate a maximum of 84degrees without making the
system unstable.
 56

EXP 8: Determination of Transfer Function of a Physical

System using Frequency Response and Bode’s
Asymptotes

FORMULA/THEORY:
It is sometimes useful, while working on experimental data obtained
from physicalsystems, to reverse engineer the frequency response
analysis, i.e., obtaining the transfer function of a system from a given
Bode plot. This is possible because, the Bode plot completely gives the
frequency domain parameters. From the magnitude and phase plots, it
is possible to determine the no. of poles andzeros, the relative order of
the system, the break frequencies of the poles and zeros and the open-
loop system gain.

Method for the estimation of transfer function from the Bode plot:
1. Relative Order of the system:

It is known that each pole in a system contributes −90° to the phase,

whereas eachpole in a system contributes +90° to the phase. Thus, by
looking at the phase plot, one can determine the relative order of the
system,i.e., how many extra poles are there than zeros in a system or
vice versa.
2. Number of zeros and poles:

We utilize the magnitude plot to determine the number of poles and

zeros in a system. The slope changes in the magnitude plot will help
determine the number of poles and zeros. For a pole, the slope change
is −20 𝑑𝐵/𝑑𝑒𝑐 whereas for a zero, it is +20 𝑑𝐵/𝑑𝑒𝑐.
 57

3. Location of poles and zeros:

To obtain the locations of poles and zeros, it is needed to find the break
frequencies of the poles and zeros. Before that, the starting slope of the
magnitude plot will tell us the type of the system.

𝑆𝑡𝑎𝑟𝑡𝑖𝑛𝑔 𝑠𝑙𝑜𝑝𝑒 = 0 => 𝑇𝑦𝑝𝑒 − 0 𝑠𝑦𝑠𝑡𝑒𝑚

𝑆𝑡𝑎𝑟𝑡𝑖𝑛𝑔 𝑠𝑙𝑜𝑝𝑒 = −20 𝑑𝐵⁄𝑑𝑒𝑐 => 𝑇𝑦𝑝𝑒 − 1 𝑠𝑦𝑠𝑡𝑒𝑚
𝑆𝑡𝑎𝑟𝑡𝑖𝑛𝑔 𝑠𝑙𝑜𝑝𝑒 = −40 𝑑𝐵⁄𝑑𝑒𝑐 => 𝑇𝑦𝑝𝑒 − 2 𝑠𝑦𝑠𝑡𝑒𝑚

To find the break frequencies of the poles/zeros, it is necessary to see

the bode asymptotic plot. The frequencies where the slope changes in
the magnitude plot will give the break frequencies of the poles/zeros. At
this point, we can estimate the transfer function as shown below:

Where:
K – open loop system gain
𝑤𝑐𝑧𝑛 − 𝑏𝑟𝑒𝑎𝑘 𝑓𝑟𝑒𝑞𝑢𝑒𝑛𝑐𝑖𝑒𝑠 𝑜𝑓 𝑧𝑒𝑟𝑜𝑠
𝑤𝑐𝑝𝑛 − 𝑏𝑟𝑒𝑎𝑘 𝑓𝑟𝑒𝑞𝑢𝑒𝑛𝑐𝑖𝑒𝑠 𝑜𝑓 𝑝𝑜𝑙𝑒𝑠

4. Finding Gain:
To find the open loop system gain 𝐾, it is necessary to look at the start of the
Bode magnitude plot. If the magnitude at the starting frequency 𝜔𝑠 is 𝑀 (𝑖𝑛 𝑑𝐵),
then:

From this, system gain can be measured, and the final transfer function can be
found out.
 58

CALCULATION:

To determine the Transfer function of the above plot:

1. Relative order of the system:

 The phase plot ends at −180.

 Thus, the system has 2 more poles than zeros.

2. Number of zeros and poles:

 From the magnitude plot, there are no slope changes from

negative to positive (+20 dB/dec). Thus, the system has no zeros.

Number of poles=2
Number of zeros=0
 59

3. Location of poles and zeros:

 From the magnitude plot, it can be seen that the slope with
which it starts is -20 dB/dec. Thus, the system is a type-1 system,
i.e., it has a pole at 𝑠 = 0. This can also be said from the fact that
the phase plot begins at −90°.

 Secondly, we need to look for change in slope and the corresponding

frequency. From the magnitude plot, it is clear that at 𝜔𝑐𝑝1 = 0.5 𝑟𝑎𝑑/𝑠,
there is a change in slope from −20 𝑑𝐵/𝑑𝑒𝑐 𝑡𝑜−40 𝑑𝐵/𝑑𝑒𝑐. Thus, 𝜔𝑐𝑝1is a
break frequency.

4. Finding Gain:

 From the magnitude plot, 𝜔𝑠 = 0.1 𝑟𝑎𝑑/𝑠 and 𝑀 = 34 𝑑𝐵. Thus:

 Thus, the transfer function is given by:

 60

MATLAB CODE:

>> sys = tf([5], [2 1 0])

>> figure (1)
>>bode(sys)
>>title ('Bode Plot of G(s)')
>>figure (2)
>>asymp(sys)
>>title ('Asymptotic Bode Plot of G(s)')

OBSERVATIONS:
MATLAB:
 61

INFERENCE:
 A pole contributes −90° to the phase plot because:

 The time constant form of the transfer function is:

 his method of finding the transfer function is immensely helpful when working
with experimental data. From the data, the Bode plot can be obtained and thus,
the transfer function of the system can be easily obtained.
 62

EXP 9: Design of Lag, Lead compensators and Evaluation

of Closed Loop Performance

FORMULA/THEORY:

1. Lead compensator:

The lead compensator is an electrical network which produces a sinusoidal output

having phase lead when a sinusoidal input is applied. The lead compensator circuit
in the ‘s’ domain is shown in the following figure. Here, the capacitor is parallel to
the resistor R1 and the output is measured across resistor R2.

The transfer function of this lead compensator is –

 63

We know that the phase of the output sinusoidal signal is equal to the sum of the
phase angles of input sinusoidal signal and the transfer function. So, to produce
the phase lead at the output of this compensator, the phase angle of the transfer
function should be positive. This will happen when 0<β <1. Therefore, zero will be
nearer to origin in pole-zero configuration of the lead compensator.

2. Lag Compensator:

The Lag Compensator is an electrical network which produces a sinusoidal output

having the phase lag when a sinusoidal input is applied. The lag compensator
circuit in the ‘s’ domain is shown in the following figure.

Here, the capacitor is in series with the resistor R2 and the output is measured
across this combination.
 64

From the transfer function, we can conclude that the lag compensator has one
1 1
pole at 𝑠 = − and one zero at 𝑠 = − . This means, the pole will be nearer to
𝑐𝛼 𝛼
origin in the pole-zero configuration of the lag compensator.
Substitute, 𝑠 = 𝑗𝜔 in the transfer function.

We know that the phase of the output sinusoidal signal is equal to the sum of the
phase angles of input sinusoidal signal and the transfer function. So, to produce
the phase lag at the output of this compensator, the phase angle of the transfer
function should be negative. This will happen when 𝛼 > 1.

PROCEDURE:

Q1. The open loop transfer function of certain unity feedback control system
is given byG(s) = K/s(s+4) (s+80). It is desired to have the phase margin to be at least 33°
and the velocity error constant Kv= 30 sec-1. Design a phase lag series compensator.
 65

ANS:

1. Determining the value of K.

2. The bode plot of the uncompensated system.

3. The phase margin is calculated.
4. Suitable values for the phase margin of the compensated system.
5. New gain crossover frequency is determined.

6. The value 𝛽 of the compensator is determined.

7. The transfer function of lag compensator is determined.

The zero of the compensator is arbitrarily placed at one tenth of the new
gain crossover frequency. Where Zc and Pc is the zero and pole of the
compensator.

The TF is given as follows

8. The open loop response is generated by placing the plant in series with the lag
compensator.

The actual phase margin is calculated at the new gain crossover frequency.

Q2. Design a lead compensator for a unity feedback system with open loop transfer
function, G(s) = K/s(s+8) to satisfy the following specifications (i) percentage peak
overshoot = 9.5% (ii) natural frequency response, ωn = 12 rad/sec (iii) Velocity error
constant Kv ≥ 10.
 66

ANS:
1. Determination of the dominant pole with the given parameters

2. The pole zero plot is drawn.

3. Angle contributed by the lead network

4. The poles and zeroes of the compensator are found

5. Determination of the TF of the compensator

6. Determining the open loop transfer function of the lead compensated system

The lead compensator is connected in series with the plant the value of K is
determined by the above formula. And the vector lengths from all poles are
represented as l1,l2,l4 and for zeroes it is l3.
7. The check for error requirement is done.
 67

MATLAB CODE:

Q1.
>> Kv = 30
>> K = Kv * 4 * 80
>> G=tf(K, conv([1 4 0],[1 80]))
>> [Gm Pm] = margin(G)
>> margin(G)
>> Pd = 33
>> Pd = Pd + 5
>> Phgcm = -142
>> wgcm= 4.58 % from the first bode plot use the mouse pointer and locate wgcm
corresponding to phgcm.
>> [Beta,p]=bode(G,wgcm)
>> T=10/wgcm
>> Zc=1/T
>> Pc=1/(Beta*T)
>> Gc=tf([1 Zc],[1 Pc])
>> sys=Gc*G/Beta
>> figure(1)
>> margin(G)
>> hold on
>> margin(sys)
>> legend('Without Compensator', 'With Lag compensator')
>> figure(2)
>> step(feedback(G,1))
>> hold on
>> step(feedback(sys,1))
>> legend('Without Compensator', 'With Lag compensator')

Q2.
>> Mp=9.5
>> e=log(Mp)*sqrt(1/(pi^2+(log(Mp)^2)));
>> w=12;

%open loop poles

>> s1 = 0
>> s2 = -8
 68

%Dominant pole

>> sd1=-e*w+j*w*sqrt(1-e^2)
>> sd2=-e*w-j*w*sqrt(1-e^2)
>> theta1=angle(sd1 - s1)
>> theta2=angle(sd1 - s2)
>> phi = theta1 + theta2 - 3.14
>> phi1 = (phi)/2

%finding the pole and the zero

>> m2 = tan(theta1/2 - phi1)

>> x = [linspace(-100,100)]
>> syms k
>> eqn = m2*(k + 7.2) + 9.6 == 0;
>> Pc = double(solve(eqn)); % X values where y=0
>> m3 = tan((theta1/2) + phi1)
>> syms k
>> eqn = m3*(k + 7.2) + 9.6 == 0;
>> Zc = double(solve(eqn)); % X values where y=0
>> T=-1/Zc;
>> a=-1/(Pc*T);
>> s=tf('s');
>> G=(s+1/T)/(s+1/(a*T));
>> l1 = abs(sd1 - s1)
>> l2 = abs(sd1 - s2)
>> l3 = abs(sd1 - Zc)
>> l4 = abs(sd1 - Pc)

>> K=(l1*l2*l4)/l3;
>> G=K*G/(s*(s+8))
>> sys=K/(s*(s+8))
>> figure(1)
>> margin(sys)
>> hold on
>> margin(G)
>> legend('Uncompensated','Compensated')
>> figure(2)
>> step(feedback(sys,1))
>> hold on
 69

>> step(feedback(G,1))
>> legend('Without Compensator', 'With Lead compensator')

OBSERVATION:

Q1.
 70

Q2.
 71

TABULATION:

Parameters Q1 Q2
Uncompensated Compensated Uncompensated Compensated
0.1064 0.2525 0.1110 0.1432
Rise time(s)
2.9940 2.3096 0.9089 0.4917
Settling
Time(s)
68.2050 40.4974 33.7293 12.6234
Overshoot(%)
0.2895 0.649 0.2763 0.3137
Peak time(s)
1.6821 1.4050 1.3373 1.1262
Peak Value
 72

INFERENCE:

 Lag compensators reduce the gain crossover frequency of the system. Thus,
the system response becomes slower.

 Lag compensators reduce the high frequency gain of the system. Hence, the
high frequency noise attenuation becomes better.

 The time constant of the lag compensator should not be unnecessarily large.
Otherwise, the system's response will be too slow.

 From the observation table for the first question Lag compensators allow a
high gain at low frequencies and thus improve the steady state error.

 From the observation table the presence of a lead compensator increases the
gain of the system at high frequencies (the amount of this gain is equal to a).

 This can increase the crossover frequency, which will help to decrease the rise
time and settling time of the system (but may amplify high frequency noise).

 In a lead compensated design additional positive phase increases the phase

margin and thus increases the stability of the system.

 And it is seen that the gain margin of lead compensated system is infinity.
A gain margin of infinity means that no matter how much you increase
the gain, the system will always be stable.

 Due to the presence of phase lead compensation maximum overshoot of the

system decreases.

 The general comparison is summarized below:

o If the design specifications are given in terms of steady state error and
phase margin, either the lead or the lag compensator may be able to
meet the requirement.
o However, if specification on rise time or gain crossover frequency is
also to be met, then only one of the two compensators would be able
to do the job. Or we need to combine lag and lead compensator.
 73

EXP 10: Design of PID controllers and Evaluation of

Closed Loop performance

FORMULA/THEORY:
A controller is a mechanism that seeks to minimize the difference
between the actual value of a system (i.e., the process variable) and
the desired value of the system (i.e., the setpoint). Controllers are a
fundamental part of control engineering and used in all complex
control systems. Some of the uses are:

o Controllers improve the steady-state accuracy by decreasing

the steady state error. As the steady-state accuracy improves,
the stability also improves.
o Controllers also help in reducing the unwanted offsets
produced by the system, while reducing the noise signals
produced by the system.
o Controllers can help to speed up the slow response of an
overdamped system.

Controllers are of different types, but the most common and popular one is the
Proportional-Integral-Derivative (PID) controller. The general form of a PID
controller in s-domain is as follows:
 74

𝑲𝒊
𝑪(𝒔) = 𝑲𝒑 + 𝒔𝑲𝒅 +
𝒔
where Kp, Kd and Ks are the gain constants of the different elements.

1. Proportional Control:

𝐾
Consider the below block diagram, where is the system that needs to be
1+𝑠𝑐
controlled and τ is the system time constant. The closed loop transfer
function of the above system is:

𝐾𝐾𝑝 𝐴𝐾𝐾𝑃
𝐶(𝑠) = ∗
1 + 𝐾𝐾𝑝 1 + 𝑠𝑟

And the system response c(t) for a unit step input is:

𝐴𝐾𝐾𝑃 𝑠𝑡
𝑐(𝑡) = (1 − 𝑒 − 𝑐 )
1+
𝐾𝐾𝑃

The system response is given below:

By using proportional control, the time response improves by a factor of

1
1+𝐾𝐾𝑃
which is good for the system, but we get a steady state offset between
the desired response and the output response. This can be reduced by
increasing the gain of the controller, but it introduces oscillations.

2. Integral Control:

Closed loop block diagram of a system with only integral controller action is
given by:
 75

1
where 𝐾𝑖 = . The closed loop transfer function of this system is given by:
𝑐𝑖

(𝑠) = 𝐾
𝐾 + 𝑟 𝑠 + 𝑟𝑟 𝑠2
𝑖 𝑖

With integral control, order of system increases by one. For a step input,
using the above transfer function, the steady state error turns out to be zero.
Thus, a major advantage of the I-controller is that it eliminates steady state
error, but the output response is generally slow and oscillatory.

3. Proportional - Integral (PI) control:

Closed loop block diagram of a system with only a PI controller action is given
by:

1
where 𝐾𝑖 = 𝑐 . From the discussions of P and I controllers, a PI controller
𝑖
provides dual advantages of fast response due to P-action and the zero
steady state error due to I-action.
 76

By varying 𝐾𝑝, it is possible to vary the damping ratio of the system

4. Proportional Derivative (PD) Control:

For higher order processes, a PD controller improves the stability of the

system. It also improves transient response of the system.

5. Proportional-Integral-Derivative (PID) Control:

It is clear from above discussions that a suitable combination of proportional,

integral, and derivative actions can provide all the desired performances of a
closed loop system. The transfer function of a P-I-D controller is given by:

𝑲𝒊
𝑪(𝒔) = 𝑲𝒑 + 𝒔𝑲𝒅 +
𝒔

The order of the controller is low, but this controller has universal
applicability; it can be used in any type of SISO system. However, for proper
 77

use, a controller must be tuned for a particular process; i.e., selection of P, I,

D parameters are especially important and process dependent.

Unless the parameters are carefully chosen, a controller may cause instability
to the closed loop system. It is not always necessary that all the features of
proportional, derivative, and integral actions should be incorporated in the
controller. In fact, in most of the cases, a simple P-I structure will suffice.

Guideline for selection of controller mode:

1. Proportional Controller:

It is simple regulating type; tuning is easy. But it normally introduces steady

state error. It is recommended for process transfer functions having a pole at
origin, or for transfer functions having a single dominating pole; for example,
with:

(𝑠) = 𝐾
(1 + 𝑟 𝑠)(1 + 𝑟 𝑠)(1 + 𝑟 𝑠)
1 2 3

2. Integral Control:
It does not exhibit steady state error but is relatively slow responding.
It is particularly effective for: fast process, with high noise level,
process dominated by dead time, high order system with all-time
constants of the same magnitude.

3. Proportional plus Integral (P-I) Control:

It does not because offset associated with proportional control. It

also yields much faster response than integral action alone. It is
widely used for process industries for controlling variables like level,
flow, pressure, etc., those do not have large time constants.

4. Proportional plus Derivative (P-D) Control:

It is effective for systems having large number of time constants. It

results in a more rapid response and less offset than is possible by
pure proportional control. But one must be careful while using
derivative action in control of fast processes, or if the measurement
 78

is noisy (ex: flow measurement).

5. Proportional plus Integral plus Derivative (P-I-D) Control:

It finds universal application. But proper tuning of the controller is

difficult. It is particularly useful for controlling slow variables, like pH,
temperature, etc. in process industries.

DESIGNING A PI/PID CONTROLLER ACCORDING TO GIVEN

SPECIFICATIONS:
 The most common way of tuning a P/PI/PID controller is using Ziegler-
Nichol’s technique, but that method cannot be used if design
specifications are given.
 Two types of specifications can be given:
 In frequency domain, like phase margin
 In time domain, like damping ratio, steady state error
 We employ two different methods of design depending on the type of
specification given.
For frequency domain specifications:
When the phase margin and gain crossover frequency are the design
parameters, the following steps are followed:
Step 1:
Determine magnitude and phase of uncompensated open loop sinusoidal
transfer function G(jw)
𝐿𝑒𝑡 ∅d 𝑡ℎ𝑒 𝑑𝑒𝑠𝑖𝑟𝑒𝑑 𝑝ℎ𝑎𝑠𝑒 𝑚𝑎𝑟𝑔𝑖𝑛 𝑎𝑡 𝑔𝑎𝑖𝑛 𝑐𝑟𝑜𝑠𝑠𝑜𝑣𝑒𝑟 𝑓𝑟𝑒𝑞𝑢𝑒𝑛𝑐𝑦 𝜔1.
𝐿𝑒𝑡 𝐾1 = |(𝑗𝑤)| 𝑎𝑡 𝜔 = 𝜔1 𝑎𝑛𝑑 ∅1 = ∠𝐺(𝑗𝑤) 𝑎𝑡 𝜔 = 𝜔1
Step 2:

Determine the phase margin of uncompensated system and the angle to be

contributed by the controller to achieve the desired phase margin.
𝐿𝑒𝑡 ∅u − 𝑃ℎ𝑎𝑠𝑒 𝑚𝑎𝑟𝑔𝑖𝑛 𝑜𝑓 𝑢𝑛𝑐𝑜𝑚𝑝𝑒𝑛𝑠𝑎𝑡𝑒𝑑 𝑠𝑦𝑠𝑡𝑒𝑚 𝑎𝑛𝑑 𝜃 − 𝑃ℎ𝑎𝑠𝑒 𝑎𝑛𝑔𝑙𝑒 𝑜𝑓
𝑐𝑜𝑛𝑡𝑟𝑜𝑙𝑙𝑒𝑟 𝑎𝑡 𝜔 = 𝜔1
Now, ∅u = 180° + ∅1 𝑎𝑛𝑑 𝜃 = ∅d − ∅u.
 79

Step 3:
Determine the Transfer function of the controller.

For a PI controller:

Step 4:
Determine open loop transfer function of compensated system.

For time domain specifications:

When the damping ratio 𝜁, error constant and natural frequency 𝜔d are the
design parameters, the following steps are followed:

Step 1:

Determine the dominant pole 𝑠d and calculate its magnitude and phase:

By considering the dominant pole

Step 2:

Determine magnitude and phase of G(s) at 𝑠=𝑠d.

 80

Step 3:
Determination of transfer function of PI/PID controller

For PI controller:

For PID controller:

𝐺 (𝑠) = 𝐾 + 𝑠𝐾 + 𝐾𝑖
𝑐 𝑝 𝑑
𝑠

Determine 𝐾i from the specified error constant, such that the compensated
system meets the error requirement.
Calculate the parameter 𝐾D and 𝐾P using the following equations

Step 4:
Determine open loop transfer function of compensated system
Questions Given:

1. Consider a unity feedback system with open loop transfer function G(s) =
100/(s+1)(s+2)(s+5). Design a PI controller, so that the phase margin of the
system is 60° at a frequency of 0.5 rad/sec.
 81

2. Consider a unity feedback system with open loop transfer function G(s) =
4/(s+1)(s+5). Design a PI controller so that the closed loop has a damping ratio
of 0.9 and natural frequency of oscillation as 2.5 rad/sec.

3. Consider a unity feedback system with open loop transfer function G(s) = 75/
(s+1)(s+3)(s+8). Design a PID controller to satisfy the following specifications
(a) steady state error for unit ramp input should be less than 0.08, (b)
Damping ratio = 0.8 (c) Natural frequency of oscillation = 2.5 rad/sec.

MATLAB CODE:
For 1st question:

>> PM=60/57; %60 degrees converted to radians

>> w=0.5;
>> s=tf('s');
>> G=10/((s+1)*(1+0.5*s)*(1+0.2*s))
>> G1=10/((j*w+1)*(1+0.5*j*w)*(1+0.2*j*w))
>> M1=abs(G1)
>> theta1=(angle(G1))

%finding phase margin and theta2

>> PM1=3.14+(theta1)
>> theta2=PM-PM1

%determining constants

>> Ki=(-w*sin(theta2))/M1
>> Kp=cos(theta2)/M1
>> contr=pid(Kp,Ki)
>> OL=contr*G
>> figure(1)
>> step(feedback(G,1))
>> hold on
>> step(feedback(OL,1))
>> legend('Without controller', 'With controller')
>> figure(2)
>> subplot 211
>> rlocus(OL)
>> title('Root Locus of G(s)')
 82

>> subplot 212

>> rlocus(G)
>> title('Root locus of compensated system')

For 2nd question:

>> e=0.9;
>> w=2.5;
>> s=tf('s');
>> G=4/((s+1) *(s+5));

%Dominant pole

>> sd1=-e*w + j*w*sqrt(1-e^2)

>> sd2=-e*w-j*w*sqrt(1-e^2)
>> M1=abs(sd1)
>> theta1=angle(sd1)

%Magnitude and Phase of G(s) at s=sd

>> G1=4/((sd1+1) *(sd1+5));

>> M2=abs(G1)
>> theta2=angle(G1)

%Determining the constants

>> Ki=(-M1*sin(theta2)/(M2*sin(theta1)))
>> Kp=(-sin(theta2+theta1)/(M2*sin(theta1)))-(2*Ki*cos(theta1)/M1)

%Finding the TF

>> cont=pid(Kp,Ki)
>> OL1=G*cont
>> figure(1)
>> step(feedback(G, 1))
>> hold on
>> step(feedback(OL, 1))
>> legend('No controller', 'With controller')
>> figure(2)
>> subplot 211
>> rlocus(OL)
>> title('Root Locus of G(s)')
>> subplot 212
>> rlocus(G)
 83

>> title('Root locus of compensated system')

For 3rd question:

>> e=0.8;
>> ess=0.08;
>> w=2.5;
>> s=tf('s');
>> G=75/((s+1)*(s+3)*(s+8));

%Dominant pole

>> sd1=-e*w+j*w*sqrt(1-e^2)
>> sd2=-e*w-j*w*sqrt(1-e^2)
>> M1=abs(sd1)
>> theta1=angle(sd1)

%Magnitude and Phase of G(s) at s=sd

>> G1=75/((sd1+1)*(sd1+3)*(sd1+8));
>> M2=abs(G1)
>> theta2=angle(G1)

%Determining the constants

>> Kv=1/ess
>> Ki=Kv/3.125
>> Kd=(sin(theta2)/(M1*M2*sin(theta1)))+(Ki/M1^2)
>> Kp=(-sin(theta2+theta1)/(M2*sin(theta1)))-(2*Ki*cos(theta1)/M1)

%Finding the TF

>> cont=pid(Kp,Ki,Kd)
>> OL=G*cont
>> figure(1)
>> step(feedback(G, 1))
>> hold on
>> step(feedback(OL, 1))
>> legend('No controller', 'With controller')
>> figure(2)
>> subplot 211
>> rlocus(OL)
>> title('Root Locus of G(s)')
 84

>> subplot 212

>> rlocus(G)
>> title('Root locus of compensated system')

OBSERVATIONS:

First question:
 85

Second question:
 86

Third question:

TABULATIONS:

s no Settling Time(s) Steady state error

Without With Without With

Controller Controller Controller Controller

1. 21.8 8.04 0.091 0

2. 1.94 2.73 0.55 0

3. 2.7 1.28 0.242 0

 87

INFERENCE:

 From the observations, a PI/PID controller always makes the steady state
error vanish.
 From the plots obtained, the system’s closed loop step response reaches
steady state more quickly, except in the case of controller 2, where the design
specification led to the settling time being higher than the uncompensated
system.

 From the Bode plot of controller 1, we can see that the design specification of
phase margin has been met, thus, the controller is functioning properly.
 From the root locus plots of the compensated and uncompensated systems,
we can see that the addition of the poles/zeros by the controller pulls the root
locus towards the left so that it can accommodate the necessary gains.
 The step responses inform us about the closed loop performance of the
system. The primary goal of any system is to obtain the desired response, i.e.,
the input. More the system is close to the input, more closed loop stable it is.
 88

EXP 11: Discretization of a Continuous System and the

Effect of Sampling

Let a unit impulse train be defined as follows:

If a continuous-time signal x(t) is sampled in a periodic manner, mathematically,

the sampled signal can be represented as follows:

where:
𝑥∗(𝑡) − 𝑠𝑎𝑚𝑝𝑙𝑒𝑑 𝑠𝑖𝑔𝑛𝑎𝑙
𝑇 − 𝑠𝑎𝑚𝑝𝑙𝑖𝑛𝑔 𝑡𝑖𝑚𝑒
 89

The disadvantage with using this method to sample an actual continuous system is
not always ideal. This is because, impulse trains are not usually the type of input
that act as a discrete input when you have a discrete system that interacts with a
continuous system.
For example, consider a discrete system that constantly inputs a motor with the
speed for a defined sample time. With impulse sampling, the motor will throttle
only for the period of the impulse and then slow down, only to speed up in the next
impulse. But this is not what is ideal.
Ideally, we would want the motor to hold its speed constant until the next sample
arrives from the discrete system. This is achieved through data hold.

Data Hold:
Data hold is a process of generating a continuous-time signal h(t) from a
discrete-time sequence x(nT).
The signal h(t) during the time interval 𝑛𝑇 ≤ 𝑡 ≤ (𝑘 + 1)𝑇 (𝑇 − 𝑠𝑎𝑚𝑝𝑙𝑒 𝑡𝑖𝑚𝑒),
maybe approximated as a polynomial as follows:
ℎ(𝑛𝑇 + 𝑟) = 𝑎𝑛 𝑟𝑛 + 𝑎𝑛−1 𝑟𝑛 −1 + ⋯ + 𝑎1𝑟 + 𝑥(𝑛𝑇)
If the data hold circuit is an n-th order polynomial extrapolator, it is called
an n-th order hold. It uses the past (n+1) discrete data to generate the
polynomial.

Sampling:
Sampling is the reduction of a continuous-time signal to a discrete-time signal. A
sample is a value or set of values at a point in time and/or space. A sampler is a
subsystem or operation that extracts samples from a continuous signal.
A theoretical ideal sampler produces samples equivalent to the instantaneous
value of the continuous signal at the desired points. While sampling, we come
across several problems, one of which is called aliasing.

Aliasing:
Aliasing refers to the distortion or artifact that results when a signal
reconstructed from samples is different from the original continuous
signal. In addition to this, by altering the sampling rate, we also decimate
or extrapolate the discretized signal.
 90

Decimation:
Consider a signal x(n) compressed by a factor M. This yields: 𝑥(𝑛) = 𝑥(𝑀𝑛)

This means that, only every M-th sample is considered. So, we lose
information while we decimate a signal.

Interpolation:
Consider a signal x(n) expanded by a factor L. These yields: 𝑥(𝑛) = 𝑥(𝑛/𝐿)
Doing this, every 2nd sample will become zero. So, we expand the sampled
data. The newly inserted sampling points are calculated using a proper
interpolating function.

QUESTIONS:
Q1. Determine the z-domain transfer function for the following s-domain transfer
function.

Discretize the given continuous system and study the effect of sampling.

PROCEDURE:
Conversion to s-domain:
Step 1: Take Inverse Laplace transform of H(s)

Step 2: Sample the continuous response by setting t = nT where T – sampling

time:
 91

Step 3: Finding the z-transform

OBSERVATION:

Q1.

𝑃𝑙𝑜𝑡 1: 𝐹𝑖𝑥𝑒𝑑 𝑎, 𝑣𝑎𝑟𝑦𝑖𝑛𝑔 𝑠𝑎𝑚𝑝𝑙𝑖𝑛𝑔

𝑡𝑖𝑚𝑒

𝑃𝑙𝑜𝑡 2: 𝑎 = 2, 𝑇𝑠 = 1
 92

𝑃𝑙𝑜𝑡 3: 𝑎 = 4, 𝑇𝑠 = 1

𝑃𝑙𝑜𝑡 4: 𝑎 = 6, 𝑇𝑠 = 1
 93

INFERENCE:

 If the signal is slowly varying, then fewer samples per second will be required
than if the waveform is rapidly varying.
 The major drawback of random undersampling is that this method can
discard potentially useful data that could be important for the induction
process.
 It can be seen that undersampling causes poor resolution and the step
response is not smooth as well.
 As the sampling time increases the resolution gets clearer. This is because of
Oversampling. It can be inferred that oversampling helps in improving the
resolution and also data is not discarded
 From the above it can be inferred that:
1. Oversampling methods duplicate or create new synthetic examples in the
minority class.
2. Undersampling methods delete or merge examples in the majority class.

 As the ‘a’ in the equation increases, the resolution increases in the step response
of the z domain
 94

EXP 12: Test of Observability and Controllability in

Continuous and Discrete Domain in Simulation Platform

FORMULA/THEORY:

It is impractical to know all the states of a system. So, sometimes it is limited to

the critical states. Thus, if a state is not observable then the controller will not
be able to determine its behavior from the system output and hence not be
able to use that state to stabilize the system.

Controllability and observability work together, as stated in a further theorem.

Controllability:

Consider a dynamic system as follows:

The state equation (or the pair (A, B)) is said to be completely state
controllable or simply state controllable if for any initial state x(0) and any
final state x(N), there exists an input sequence u(k), k = 0, 1, 2, · · · , N, which
transfers x(0) to x(N) for some finite N. Otherwise the state equation is state
uncontrollable.

1. The state equation (1) or the pair (A, B) is state controllable if and
only if the n × np state controllability matrix has rank n, i.e., full row
rank.
The state controllability matrix is given as follows:
𝑈𝑐 = [𝐵 𝐴𝐵 𝐴2𝐵 … … … 𝐴𝑛−1 𝐵]
The full rank condition can also be stated as follows: If the
determinant of the the controllability matrix is zero, then the matrix
is full rank.
 95

2. When the matrix ‘A’ has distinct eigenvalues and in Jordan/Diagonal

canonical form, the state model is controllable if and only if all the rows
of matrix ‘B’ are non-zero.

Observability:
The state model (or the pair (A, C) ) is said to be observable if any initial state x(0)
can be uniquely determined from the knowledge of output y(k) and input
sequence u(k), for k = 0, 1, 2, · · · , N, where N is some finite time. Otherwise, the
state model is unobservable.

1. The state model (1) or the pair (A, C) is observable if the np × n

observability matrix has rank n, i.e., full column rank.

2. If the state model is in observable canonical form, then the system is

observable.

Theorem of Duality:

The pair (A, B) is controllable if and only if the pair ( 𝑇,𝐵𝑇 ) is observable.

QUESTIONS:
Q1. Test the controllability and observability of the given system in both continuous
and discrete domain.
Given characteristic equation: 5𝑧2 − 2𝑧 + 2 = 0

Assume: It is a unity feedback system.

 96

PROCEDURE:
Finding the transfer function:
We know that any transfer function is of the form:

Here, as it is unity feedback, H(z) = 1.

The characteristic equation of the above transfer function is:

Thus,
The transfer function in z-domain is:
 97

MATLAB CODE:

Q1.
disp('To check controllability and observability in z-domain')
disp(' --- ---------------------------------------------------------------- ')
disp('Given characteristic equation: 5z^2 - 2z + 2 = 0')

H=tf([5,-2,1],[5,-2,2],0.1)
[b1 a1] = tfdata(H, 'v');
[A1,B1,C1,D1]=tf2ss(a1,b1);
g1=rank(ctrb(A1,B1));
g2=rank(obsv(A1,C1));
n = size(A1, 1);

if (g1==n || g2==n)
disp('Given system is both fully controllable and observable')
elseif (g1==n)
disp('Given system is fully controllable only')
elseif (g2==n)
disp('Given system is fully observable only')
else
disp('Given system is neither fully controllable nor observable')
end

fprintf('\n')
disp('To check controllability and obsevability in s-domain')
 98

disp(' --- -------------------------------------------------------------- ')

sys=d2c(H) %converting given z domain tf to s domain

[a b] = tfdata(H, 'v');
[A,B,C,D]=tf2ss(b,a);
f1=rank(ctrb(A,B));
f2=rank(obsv(A,C));
n = size(A, 1);

if (f1==n || f2==n)
disp('Given system is both fully controllable and observable')
elseif (f1==n)
disp('Given system is fully controllable only')
elseif (f2==n)
disp('Given system is fully observable only')
else
disp('Given system is neither fully controllable nor observable')
end

OBSERVATION:

Q1.
 99

INFERENCE:

 If a continuous time system is controllable or observable, then its

discrete time version, with any sampling period, is controllable or
observable.
 If a continuous time system is not controllable or observable, then its
discrete time version, with any sampling period, is not controllable or
observable.
 Loss of controllability and/or observability occurs only in presence of
oscillatory modes of the system.
 Controllability is a black and white issue. That is a system can either be
completely controllable or not controllable.
 It is not always possible to observe all the states of a system in a real-life
scenario, in that case, only the critical cases are considered for observability.
 100

EXP 13: Design of State Feedback Controller and

Evaluation of Closed Loop Performance

FORMULA/THEORY:
The state of a system is a minimal set of variables such that the knowledge of
these variables at a particular time t0, with the knowledge of the input(s) to the
system for t>t0, can completely determine the behavior of the system. The state
space approach is a method for modelling, analyzing, and designing a wide range
of both linear and non-linear systems. State space models use state variables to
describe a system by a set of first order differential equations.
A system is represented in state space by the following equations:

Where:
𝑥(𝑡): 𝑠𝑡𝑎𝑡𝑒 𝑣𝑒𝑐𝑡𝑜𝑟
𝑢(𝑡): 𝑖𝑛𝑝𝑢𝑡/𝑐𝑜𝑛𝑡𝑟𝑜𝑙 𝑣𝑒𝑐𝑡𝑜𝑟
𝑥(𝑡): 𝑜𝑢𝑡𝑝𝑢𝑡 𝑣𝑒𝑐𝑡𝑜𝑟
A: 𝑆𝑦𝑠𝑡𝑒𝑚 𝑚𝑎𝑡𝑟𝑖𝑥 𝑡ℎ𝑎𝑡 𝑟𝑒𝑙𝑎𝑡𝑒𝑠 ℎ𝑜𝑤 𝑡ℎ𝑒 𝑐𝑢𝑟𝑟𝑒𝑛𝑡 𝑠𝑡𝑎𝑡𝑒 𝑎𝑓𝑓𝑒𝑐𝑡𝑠 𝑡ℎ𝑒 𝑠𝑡𝑎𝑡𝑒
𝑐ℎ𝑎𝑛𝑔𝑒
 101

B: 𝐼𝑛𝑝𝑢𝑡 𝑜𝑟 𝑐𝑜𝑛𝑡𝑟𝑜𝑙 𝑚𝑎𝑡𝑟𝑖𝑥 − 𝑑𝑒𝑡𝑒𝑟𝑚𝑖𝑛𝑒𝑠 ℎ𝑜𝑤 𝑖𝑛𝑝𝑢𝑡 𝑎𝑓𝑓𝑒𝑐𝑡𝑠 𝑠𝑦𝑠𝑡𝑒𝑚 𝑐ℎ𝑎𝑛𝑔𝑒

C: 𝐷𝑒𝑡𝑒𝑟𝑚𝑖𝑛𝑒𝑠 𝑡ℎ𝑒 𝑟𝑒𝑙𝑎𝑡𝑖𝑜𝑛𝑠ℎ𝑖𝑝 𝑏𝑒𝑡𝑤𝑒𝑒𝑛 𝑡ℎ𝑒 𝑠𝑦𝑠𝑡𝑒𝑚 𝑠𝑡𝑎𝑡𝑒𝑠 𝑎𝑛𝑑 𝑡ℎ𝑒 𝑠𝑦𝑠𝑡𝑒𝑚
𝑜𝑢𝑡𝑝𝑢𝑡
D: 𝐴𝑙𝑙𝑜𝑤𝑠 𝑖𝑛𝑝𝑢𝑡 𝑡𝑜 𝑑𝑖𝑟𝑒𝑐𝑡𝑙𝑦 𝑎𝑓𝑓𝑒𝑐𝑡 𝑡ℎ𝑒 𝑠𝑦𝑠𝑡𝑒𝑚 𝑜𝑢𝑡𝑝𝑢𝑡

Full state feedback:

Full state feedback, or pole placement, is a method employed in feedback

control system theory to place the closed-loop poles of a plant in pre-
determined locations in the s-plane. This state feedback is done so that the
dynamic response of the system becomes more stable. As the pole locations can
be chosen by us, it is possible to get highly stable responses.
For this method to work, the system needs to be completely state
controllable.

Complete state controllability:

Complete state controllability describes the ability of an external input (the

vector of control variables) to move the internal state of a system from any initial
state to any final state in a finite time interval.

From the theory of state space model, the matrix ‘A’ captures the dynamics of
the system. Thus, a controller that we design will have to modify that matrix if we
want to change the dynamic response of the system. One important result here
is:

Thus, if we can change the eigen values w.r.t our need, we can stabilize the
system.
 102

Assume D=0.

Now, as already stated, the x(t) vector contains all the possible states of the
system. So, from the block diagram all the states of the system are fed back to
the reference input via an appropriate gain matrix ‘K’. Hence, this method of
control is known as full-state feedback control.

Here:

K: Gain matrix obtained according to the position of the poles

R: reference input

U: input

In some cases, there can be an additional block before the summing point called
‘Kr’, which adjusts for the steady-state error of the system. Together with the
blocks K, Kr and the summing point, we have a full-state feedback controller that
can be designed appropriately with reference to the poles given.

State feedback design (via pole placement):

Consider the system:

𝑥̇(𝑡) = 𝐴𝑥(𝑡) + 𝐵𝑢(𝑡)
𝑦(𝑡) = 𝐶𝑥(𝑡)

We define a new input (𝑡) = −𝐾𝑥(𝑡) + 𝑟(𝑡), where K is a row vector given by:

𝐾 = (𝑘0 𝑘1 … … 𝑘𝑛−1)

This ‘K’ is the feedback gain matrix that we need to find, and r(t) is the
 103

new input signal. Thus,

(𝑡) = −𝑘0𝑥1 − 𝑘1𝑥2 … … + 𝑟(𝑡)

This is called state feedback. As it can be seen, we feed back to the input u(t) the
states x1, x2, …., xn multiplied by the elements of the matrix K, which are called
gains.

By replacing u(t) in our original equation, we get:

𝑥̇(𝑡) = (𝐴 − 𝐵𝐾)𝑥(𝑡) + 𝐵𝑟(𝑡)

𝑦(𝑡) = 𝐶𝑥(𝑡)

As we have seen, the system stability now depends on the eigenvalues of

the matrix (A - BK). Now, to design a state feedback controller, we will be
given the position where the poles need to be. This means, we know the
characteristic equation of the matrix (A - BK).

Using these eigenvalues, we can find the appropriate gains, i.e., the
matrix K, to have desired placement of eigenvalues of the matrix (A - BK).
In MATLAB, this is achieved via the place( ) function.
 104

Question:
Q1. A single input system is described by the following equations:

Design a state feedback controller which will give closed loop poles at -1±j2, -6.
From the question, we can see that the given system will have 3 state variables. The
gain matrix ‘K’ can be found easily, as the closed loop poles, i.e., eigenvalues of the
matrix (A – BK) is given. An additional block ‘Kr’ can be added in cascade with the
reference input ‘r’ so that the steady state error of the total controlled system goes to
zero.

MATLAB CODE/BLOCK DIAGRAM:

1. Code:

Q1.
>> A=[-1 0 0;1 -2 0;2 1 -3];
>> B=[10;1;0];
>> C=[1 0 1];
>> D=0;

%converts the matrices to state space form

>> sys=ss(A,B,C,D);

>> 'The eigenvectors/poles of the given system is:'

>> disp(eig(sys));

>> %pole placement

>> dp1=complex(-1,2);
 105

>> dp2=complex(-1,-2);
>> dp3=-6;
>> dp=[dp1,dp2,dp3];
>> %Calculating Gain K

>> gainK=place(A,B,dp);

>> %new sys

>> sys_new=ss(A-B*gainK,B,C,D);
>> Kr = 1/dcgain(sys_new);
>> sys_new1=ss(A-B*gainK,B*Kr,C,D);

>> 'The eigenvectors/poles of the new state feedback-controlled system is:'

>> disp(eig(sys_new));
>> step(sys)
>> hold on
>> step(sys_new1)
>> hold on
>> plot([Link], [Link], 'r')
>> legend('Without feedback','With state feedback');

Explanation:

 To convert given matrices to state space form, we use ss( ) function.

 The function eig( ) displays the eigenvalues of the system.
 The place( ) function computes the requires gain matrix K, in
accordance to the needed design parameter, i.e., the required pole
placement, which is given as input.
 The sys_new element is the new state space description of the
state-feedback controlled system.
 The constant Kr is computed to offset for the steady state error.
 Also, we plot the step response of the new system by transferring
datapoints fromSimulink via the block To Workspace.
 106

2. Simulink Diagram:

OBSERVATION:

Q1.
 107

TABULATION:

[Link] Type of Rise time(s) Settling Steady state

system Time(s) error
Without state 2.42 4.21 -17.5
1. feedback
With state 0.372 3.35 0
feedback

INFERENCE:

 The state feedback will work only for fully controllable systems. This is
because all the states of the system are fed back to the reference input,
and if the system is not controllable, then this feedback won’t properly
function.

 From the step response plot for state feedback, the step response for the
uncontrolled system has a higher rise time and has a negative steady-state
error, i.e., the output is higher than the value of the step input. Whereas,
with the state feedback, the rise time & settling time gets reduced, the
steady state error also goes to zero (This is due to the corrective gain term
Kr).

 Only one integrator is used in the Simulink model because all the gain
terms are modelled as Matrix K*u, so that the output can be
demultiplexed to obtain the different states.

 The method of pole placement is basically an upgraded version of the root locus
compensator design. In root locus, we have one gain that we can alter (i.e., ‘K’),
and we can move the positions of the poles only along the root locus. But here,
we have a gain matrix K and we have the entire s-plane to place our poles.
 108

EXP 14: State Observer and Evaluation of Closed Loop

Performance

FORMULA/THEORY:
The state of a system is a minimal set of variables such that the knowledge of
these variables at a particular time t0, with the knowledge of the input(s) to the
system for t>t0, can completely determine the behavior of the system. The state
space approach is a method for modelling, analyzing, and designing a wide range
of both linear and non-linear systems. State space models use state variables to
describe a system by a set of first order differential equations.
A system is represented in state space by the following equations:

Where:
𝑥(𝑡): 𝑠𝑡𝑎𝑡𝑒 𝑣𝑒𝑐𝑡𝑜𝑟
𝑢(𝑡): 𝑖𝑛𝑝𝑢𝑡/𝑐𝑜𝑛𝑡𝑟𝑜𝑙 𝑣𝑒𝑐𝑡𝑜𝑟
𝑥(𝑡): 𝑜𝑢𝑡𝑝𝑢𝑡 𝑣𝑒𝑐𝑡𝑜𝑟
A: 𝑆𝑦𝑠𝑡𝑒𝑚 𝑚𝑎𝑡𝑟𝑖𝑥 𝑡ℎ𝑎𝑡 𝑟𝑒𝑙𝑎𝑡𝑒𝑠 ℎ𝑜𝑤 𝑡ℎ𝑒 𝑐𝑢𝑟𝑟𝑒𝑛𝑡 𝑠𝑡𝑎𝑡𝑒 𝑎𝑓𝑓𝑒𝑐𝑡𝑠 𝑡ℎ𝑒 𝑠𝑡𝑎𝑡𝑒
𝑐ℎ𝑎𝑛𝑔𝑒
 109

B: 𝐼𝑛𝑝𝑢𝑡 𝑜𝑟 𝑐𝑜𝑛𝑡𝑟𝑜𝑙 𝑚𝑎𝑡𝑟𝑖𝑥 − 𝑑𝑒𝑡𝑒𝑟𝑚𝑖𝑛𝑒𝑠 ℎ𝑜𝑤 𝑖𝑛𝑝𝑢𝑡 𝑎𝑓𝑓𝑒𝑐𝑡𝑠 𝑠𝑦𝑠𝑡𝑒𝑚 𝑐ℎ𝑎𝑛𝑔𝑒

C: 𝐷𝑒𝑡𝑒𝑟𝑚𝑖𝑛𝑒𝑠 𝑡ℎ𝑒 𝑟𝑒𝑙𝑎𝑡𝑖𝑜𝑛𝑠ℎ𝑖𝑝 𝑏𝑒𝑡𝑤𝑒𝑒𝑛 𝑡ℎ𝑒 𝑠𝑦𝑠𝑡𝑒𝑚 𝑠𝑡𝑎𝑡𝑒𝑠 𝑎𝑛𝑑 𝑡ℎ𝑒 𝑠𝑦𝑠𝑡𝑒𝑚
𝑜𝑢𝑡𝑝𝑢𝑡
D: 𝐴𝑙𝑙𝑜𝑤𝑠 𝑖𝑛𝑝𝑢𝑡 𝑡𝑜 𝑑𝑖𝑟𝑒𝑐𝑡𝑙𝑦 𝑎𝑓𝑓𝑒𝑐𝑡 𝑡ℎ𝑒 𝑠𝑦𝑠𝑡𝑒𝑚 𝑜𝑢𝑡𝑝𝑢𝑡

State Observer:

An observer is a dynamic system that is used to estimate the state of a system or

someof the states of a system. A full-state observer is used to estimate the
states of the system. The observer can be designed as either a continuous-time
system or a discrete-timesystem.
Basically, in full-state feedback, we assumed that all the states of the
system are available for feedback. But this is not always the case. In
such situations, it becomes necessaryto estimate the unavailable states.
Estimation of unmeasurable state variables is commonly known as
observation and adevice/program that does this is known as a state
observer.

A full-order state observer is one that observes all the states of the system.

For this observer to work, the system needs to be completely state observable.
Observability is a measure of how well all the internal states of a system can be
inferred from knowledge of its external outputs.

The observer is a subsystem to reconstruct the state vectors of

the plant (givensystem). The mathematical model of the observer
is basically the same as that of the plant,except that we include
an additional term that includes the estimation error to
compensate for inaccuracies in matrices A and B and the lack of
the initial error.
The estimation error or observation error is the difference between the
measuredoutput and the estimated output.
 110

State observer design:

As we have already seen, an observer is designed to estimate the
unmeasurable states of a given system. The observer is an equivalent
system with the following state space description, whose sole purpose
is to estimate all the states of the system.

To ensure that the estimation of the state variables is accurate, the

observer is also fed with the error signal between the outputs, i.e., (𝑡)
− 𝑦̂(𝑡), which is then multiplied with the gain matrix, L. Thus, the new
system takes the form:

The use of output feedback can be understood by considering the

approximation error. Let (𝑡) denote the error of the approximation of
the state variables.
 111

Substituting for the above equation from the previous equations, we get:

Converting to time domain:

This above equation means that, we choose the feedback gain matrix ‘L’
such that our system equation (A - LC) becomes stable, so that e(t) goes
to zero and the approximation is accurate.

In case of designing an observer, the eigenvalues of the observer will

have already been given. In that situation, with the matrices ‘A’ and
‘C’ known, and the characteristic equation of the observer system
known (as eigenvalues are known), we can find the feedback gain
matrix ‘L’ using linear algebra.

In MATLAB, this is done using the function acker( ) .

The quicker the approximation, better the observer system is. So,
based on the value of the gain matrix L, we will get quicker or slower
approximations. The process described above is same as the process
adopted for the state feedback controller design. So, it can be said
that both problems are mathematically the same. This property is
called duality.

Question:
Q1.
 112

 From the question, the given system will have 2 state variables.

 We will assume an initial condition of [1;1] for the system

because the output of the system is zero.

 We design the state observer exactly as the given

system without the initial conditions.

 The errors are calculated by feeding back the difference of the

outputs of the plant and observer via a gain matrix L.

 This matrix L is found in a similar way compared to the state feedback

controller.

MATLAB CODE/BLOCK DIAGRAM:

1. Code:

Q1.
>> % Original Plant

>> a=[-1 1;1 -2];

>> b=[0;0];
>> c=[1 0];
>> d=0;

>> sys=ss(a,b,c,d);
>> eig(sys);
>> rank(obsv(sys));

>> % Observer pole placement at -5 and -5

>> L = acker(a',c',[-5;-5])
 113

Explanation:

 The observability of the system is checked first. Only if the system is

observable, it is possible to design a state observer.
 The function acker( ) is very similar to the place( ) function. Here, as we are
designing an observer, we consider transposes of the original matrices, and
finally take the transpose of the total acker( ) function output. This is our
required gain matrix ‘L’.

2. Simulink Diagram:
 114

OBSERVATION FOR SIMULINK:

Q1.

𝑷𝒍𝒐𝒕 𝟏: 𝒙𝟏(𝒕),𝒙̂𝟏(𝒕) 𝒗𝒔 𝒕𝒊𝒎𝒆

𝑷𝒍𝒐𝒕 𝟐:(𝒙𝟏(𝒕)−𝒙^𝟏(𝒕)) 𝒗𝒔 𝒕𝒊𝒎𝒆

 115

𝑷𝒍𝒐𝒕 𝟑: 𝒙𝟐(𝒕),𝒙^𝟐(𝒕) 𝒗𝒔 𝒕𝒊𝒎𝒆

𝑷𝒍𝒐𝒕 𝟒:(𝒙𝟐(𝒕)−𝒙^𝟐(𝒕)) 𝒗𝒔 𝒕𝒊𝒎𝒆

 116

𝑷𝒍𝒐𝒕 𝟓: 𝒚(𝒕),𝒚^(𝒕) 𝒗𝒔 𝒕𝒊𝒎𝒆

INFERENCE:

 From the Simulink plots of the state observer, we can see that the
observer estimatesthe states of the given system perfectly, and that is
exactly the job of an observer system; to estimate the unmeasurable
states of a given system.

 Plots 2 and 4 show the error signal between the original plant and the
observer system. It decays to zero within 2 seconds. This means that it is a
particularly good approximation, and the placement of poles/eigenvalues
is good.

 Plots 1 and 3 show the actual states of the system, and the states which
the observer estimates. As we have taken the plant with initial conditions,
the states of the plant decay to zero from that initial condition. Here too,
we can see that the observer startsto follow the state of the system very
quicky (2-3 seconds). So, the observer is tracking the states of the system
perfectly.

 Plot 5 shows the output of the plant and output of the observer. Similar to
 117

the statesplot, the observer quickly starts following the output of the
system y(t).

 If the system isn’t completely observable, it is not possible to design a full

state observer, because some states can’t be observed, that is, it is not
possible to say howthe system states vary with the input.

 From the Simulink model, the difference between the outputs of theplant
and observer system is fed back to the observer via a gain matrix ‘L’. This
process continuously refines the error until it becomes zero, i.e., the
observer properlyestimates the states/output of the given system/plant. 