M.

Aamir Khan
Roll No – 024

DCS Assignment

Statement
Question
1) Identify the Transfer Function of any established physical system at your workplace/ from your
knowledge or from any research article after going through the following steps.
a) For a chosen system write a detailed problem statement mentioning the use/application of that
physical system very precisely and to the point and find G(s) and C(s)/R(s) using SISO diagram. Write
from your understanding and study and avoid copying the information from sources. [Marks: 06].

Small model of steam turbine plant and its transformation.

It’s a small steam turbine model which have important parts like governer and valve position control
inorder to control the steam injection to turbine, fuel system which circulate the steam in turbine, Ownd
desigend ST dynamices, Load and speed controller (Gain).

This feed back loop system work in this manner when ever steam turbine loaded so it need more power
to drive generator shaft. So the speed regulator sense the decrease speed in term of frequency of
turbine as load increases, so it immediately send the signal to governor to inject more stem by settting
the valve position in this maner that the speed of steam turbine will be maintained as required delta F.

Delta F is taken load frequency it control during load disturbance whenever the variation occure.

1. Governor = Xs+1/Ys+Z.
2. Valve = a/bs+c.
3. Fuel = Kf/1+sTf.
4. Turbine = 1/1+sTcd.
5. Load = Kp/1+sTp1.
6. Speed Regulator = 1/R.
 7. Load frequency = Df

Supposed value are:

1. Kp = 100.
2. Tp = 206.
3. a,b,c = 1, 0.05, 1.
4. Tf = 0.4.
5. kF = 1.
6. Tcd = 0.1.
7. Governor cofficient x,y,z =0, 0.05,1.
8. Speed regulator = 20

= 100 /(0.002S^5+0.1051S^4+1.8505s^3+12.0925S^2+20.6S +1)

= 20

= DF / Dfref = G(s) / 1 +G(s)H(s).

G(s) = 100 /(0.002S^5+0.1051S^4+1.8505s^3+12.0925S^2+20.6S +1).

= T/F = 100 /(0.002S^5+0.1051S^4+1.8505s^3+12.0925S^2+20.6S +1) / 1 +100 *20

/(0.002S^5+0.1051S^4+1.8505s^3+12.0925S^2+20.6S +1) .

= DF/Dref = 100 /(0.002S^5+0.1051S^4+1.8505s^3+12.0925S^2+20.6S +2000).

Question
Clearly mention the input and output signal(s) to be controlled or tracked in your obtained TF. Also
mention why this output signal is necessary to be tracked in your problem statement. [Marks: 02]
c) Develop a SISO diagram in sampled-data domain, z-domain. And convert your transfer function into
G(z), also find C(z)/R(z) where needed. Also find the discrete equation of the system. [Marks: 04]

SISO
When ever the frequency very from its set valve so the system become unstable and in this
manner gain/sensor should send signal to governor to very the voltage inorder to get desire
output by adjusting the valve position.

Matlab code
'For G(S)'
num = [100 0];
den = [0.002 0.1051 1.8505 12.0925 20.6 2000];
G = tf(num,den);
cltf = (G)
step(G,'--',cltf,'-')
cltfz = feedback( cltf,1)
step(cltfz)

‘For Root Locus’

rlocus(cltfz)

‘For Bode’
bode(cltfz)

'For G(P)'
p = [(1+0.1/2*a)/(1-0.1/2*a)]
syms z
x =[(z^16 - 4.517*z^15 + 10.47*z^14 - 13.09*z^13 + 9.62*z^12 - 3.12*z^11 +
0.8553*z^10 - 0.1337*z^9 + 0.01837*z^8 - 0.001209*z^7 + 7.225e-05*z^6 +
4.863e-07*z^5- 1.701e-07*z^4 - 2.788e-10*z^3)/(z^16 - 4.517*z^15 + 10.47*z^14
- 13.09*z^13 + 9.62*z^12 - 3.121*z^11 + 0.8583*z^10- 0.1385*z^9 + 0.01654*z^8
- 0.001115*z^7 - 1.772e-05*z^6 - 4.863e-07*z^5 + 1.701e-07*z^4 + 2.788e-
10*z^3)]
y= [z^(-3) * (8.515e-05*z^5 + 0.009243*z^4 + 0.02988*z^3 + 0.01241*z^2 +
0.0006743*z + 1.068e-06) /(z^5 - 2.258*z^4 + 2.684*z^3 - 0.4858*z^2 + 0.1113*z
- 0.005221)]
C = x*y

'routh_hurwitz'
e=input('enter the coefficients of characteristic equation: ');
disp('------------------------------------------------------------------------
-')
l=length(e);
m=mod(l,2);
if m==0
a=rand(1,(l/2));
b=rand(1,(l/2));
for i=1:(l/2)
a(i)=e((2*i)-1);
b(i)=e(2*i);
end
else
e1=[e 0];
a=rand(1,((l+1)/2));
b=[rand(1,((l-1)/2)),0];
for i=1:((l+1)/2)
a(i)=e1((2*i)-1);
b(i)=e1(2*i);
end
end
%%now we genrate the remaining rows of routh matrix
l1=length(a);
c=zeros(l,l1);
c(1,:)=a;
c(2,:)=b;
for m=3:l
for n=1:l1-1
c(m,n)=-(1/c(m-1,1))*det([c((m-2),1) c((m-2),(n+1));c((m-1),1) c((m-
1),(n+1))]);
end
end
disp('the routh matrix:')
disp(c)
%%now we check the stablity of system
if c(:,1)>0
disp('System is Stable')
else
disp('System is Unstable');
end

Solution
For G(Z)

100
h= exp(-0.25*s) * --------------------------------------------------------------
0.002 s^5 + 0.1051 s^4 + 1.851 s^3 + 12.09 s^2 + 20.6 s + 2000

Continuous-time transfer function.

8.515e-05 z^5 + 0.009243 z^4 + 0.02988 z^3 + 0.01241 z^2 + 0.0006743 z + 1.068e-06

hd = z^(-3) * ---------------------------------------------------------------------------------------------------
z^5 - 2.258 z^4 + 2.684 z^3 - 0.4858 z^2 + 0.1113 z - 0.005221
Sample time: 0.1 seconds
Discrete-time transfer function.

8.515e-05 z^5 + 0.009243 z^4 + 0.02988 z^3 + 0.01241 z^2 + 0.0006743 z+ 1.068e-06
a= z^(-3) * -----------------------------------------------------------------------------------------------------
z^5 - 2.258 z^4 + 2.684 z^3 - 0.4858 z^2 + 0.1113 z - 0.005221

Sample time: 0.1 seconds

Discrete-time transfer function.

Step response
For G(P)

z^16 - 4.517 z^15 + 10.47 z^14 - 13.09 z^13 + 9.62 z^12 - 3.12 z^11 + 0.8553 z^10 - 0.1337 z^9 + 0.01837 z^8 -
0.001209 z^7 + 7.225e-05 z^6 + 4.863e-07 z^5 - 1.701e-07 z^4 - 2.788e-10 z^3
P= ------------------------------------------------------------------------------------------------------------------------------------------------------------------

z^16 - 4.517 z^15 + 10.47 z^14 - 13.09 z^13 + 9.62 z^12 - 3.121 z^11 + 0.8583 z^10 - 0.1385 z^9 + 0.01654 z^8 - 0.001115 3
z^7 - 1.772e-05 z^6 - 4.863e-07 z^5+ 1.701e-07 z^4 + 2.788e-10 z^

Sample time: 0.1 seconds

Discrete-time transfer function.
x=

-(- z^16 + (4517*z^15)/1000 - (1047*z^14)/100 + (1309*z^13)/100 - (481*z^12)/50 +

(78*z^11)/25 - (8553*z^10)/10000 + (1337*z^9)/10000 - (1837*z^8)/100000 +
(696941049534839*z^7)/576460752303423488 -
(1332777259325515*z^6)/18446744073709551616 -
(4592973641239017*z^5)/9444732965739290427392 +
(6426196309889013*z^4)/37778931862957161709568 +
(2696388148068469*z^3)/9671406556917033397649408)/(z^16 - (4517*z^15)/1000 +
(1047*z^14)/100 - (1309*z^13)/100 + (481*z^12)/50 - (3121*z^11)/1000 + (8583*z^10)/10000
- (277*z^9)/2000 + (827*z^8)/50000 - (2571014955273269*z^7)/2305843009213693952 -
(1307505219944533*z^6)/73786976294838206464 -
(4592973641239017*z^5)/9444732965739290427392 +
(6426196309889013*z^4)/37778931862957161709568 +
(2696388148068469*z^3)/9671406556917033397649408)

y=

((6282961031505473*z^5)/73786976294838206464 +
(5328226733540543*z^4)/576460752303423488 + (747*z^3)/25000 + (1241*z^2)/100000 +
(6219319764451175*z)/9223372036854775808 +
5043487403704781/4722366482869645213696)/(z^3*(z^5 - (1129*z^4)/500 + (671*z^3)/250 -
(2429*z^2)/5000 + (1113*z)/10000 - 1504850793888087/288230376151711744))

C=

-(((6282961031505473*z^5)/73786976294838206464 +
(5328226733540543*z^4)/576460752303423488 + (747*z^3)/25000 + (1241*z^2)/100000 +
(6219319764451175*z)/9223372036854775808 +
5043487403704781/4722366482869645213696)*(- z^16 + (4517*z^15)/1000 -
(1047*z^14)/100 + (1309*z^13)/100 - (481*z^12)/50 + (78*z^11)/25 - (8553*z^10)/10000 +
(1337*z^9)/10000 - (1837*z^8)/100000 + (696941049534839*z^7)/576460752303423488 -
(1332777259325515*z^6)/18446744073709551616 -
(4592973641239017*z^5)/9444732965739290427392 +
(6426196309889013*z^4)/37778931862957161709568 +
(2696388148068469*z^3)/9671406556917033397649408))/(z^3*(z^5 - (1129*z^4)/500 +
(671*z^3)/250 - (2429*z^2)/5000 + (1113*z)/10000 -
1504850793888087/288230376151711744)*(z^16 - (4517*z^15)/1000 + (1047*z^14)/100 -
(1309*z^13)/100 + (481*z^12)/50 - (3121*z^11)/1000 + (8583*z^10)/10000 - (277*z^9)/2000 +
(827*z^8)/50000 - (2571014955273269*z^7)/2305843009213693952 -
(1307505219944533*z^6)/73786976294838206464 -
(4592973641239017*z^5)/9444732965739290427392 +
(6426196309889013*z^4)/37778931862957161709568 +
(2696388148068469*z^3)/9671406556917033397649408))
routh_hurwitz
enter the coefficients of characteristic equation:
-(((6282961031505473)/73786976294838206464 +
(5328226733540543)/576460752303423488 + (747)/25000 + (1241)/100000 +
(6219319764451175)/9223372036854775808 +
5043487403704781/4722366482869645213696)*(- 1 + (4517)/1000 - (1047)/100 + (1309)/100
- (481)/50 + (78)/25 - (8553)/10000 + (1337)/10000 - (1837)/100000 +
(696941049534839)/576460752303423488 - (1332777259325515)/18446744073709551616 -
(4592973641239017)/9444732965739290427392 +
(6426196309889013)/37778931862957161709568 +
(2696388148068469)/9671406556917033397649408))/(1*(1 - (1129)/500 + (671)/250 -
(2429)/5000 + (1113)/10000 - 1504850793888087/288230376151711744)*(1 - (4517)/1000 +
(1047)/100 - (1309)/100 + (481)/50 - (3121)/1000 + (8583)/10000 - (277)/2000 + (827)/50000 -
(2571014955273269)/2305843009213693952 - (1307505219944533)/73786976294838206464
- (4592973641239017)/9444732965739290427392 +
(6426196309889013)/37778931862957161709568 +
(2696388148068469)/9671406556917033397649408))
-------------------------------------------------------------------------
the routh matrix:
0.0502
0
System is Unstable

Question

2) By using Matlab/SIMULINK or any digital controller, do the following:

a) Plot or obtain different response by adding a control parameter (K) or to change ζ (zeta) in analog
control and then in digital domain and compare them along with time, sampled-time, s- and z-domain
mapping. [Marks: 07]

'For G(Z)'
h = tf(100,[0.002 0.1051 1.8505 12.0925 20.6 2000],'IODelay',0.25)
hd = c2d(h,0.1)
a=hd

'Step response'
step(h,'--',hd,'-')
'For G(P)'
p = [(1+0.1/2*a)/(1-0.1/2*a)]
syms z
x =[(z^16 - 4.517*z^15 + 10.47*z^14 - 13.09*z^13 + 9.62*z^12 - 3.12*z^11 +
0.8553*z^10 - 0.1337*z^9 + 0.01837*z^8 - 0.001209*z^7 + 7.225e-05*z^6 +
4.863e-07*z^5- 1.701e-07*z^4 - 2.788e-10*z^3)/(z^16 - 4.517*z^15 + 10.47*z^14
- 13.09*z^13 + 9.62*z^12 - 3.121*z^11 + 0.8583*z^10- 0.1385*z^9 + 0.01654*z^8
- 0.001115*z^7 - 1.772e-05*z^6 - 4.863e-07*z^5 + 1.701e-07*z^4 + 2.788e-
10*z^3)]
y= [z^(-3) * (8.515e-05*z^5 + 0.009243*z^4 + 0.02988*z^3 + 0.01241*z^2 +
0.0006743*z + 1.068e-06) /(z^5 - 2.258*z^4 + 2.684*z^3 - 0.4858*z^2 + 0.1113*z
- 0.005221)]
C = x*y

'routh_hurwitz'
e=input('enter the coefficients of characteristic equation: ');
disp('------------------------------------------------------------------------
-')
l=length(e);
m=mod(l,2);
if m==0
a=rand(1,(l/2));
b=rand(1,(l/2));
for i=1:(l/2)
a(i)=e((2*i)-1);
b(i)=e(2*i);
end
else
e1=[e 0];
a=rand(1,((l+1)/2));
b=[rand(1,((l-1)/2)),0];
for i=1:((l+1)/2)
a(i)=e1((2*i)-1);
b(i)=e1(2*i);
end
end
%%now we genrate the remaining rows of routh matrix
l1=length(a);
c=zeros(l,l1);
c(1,:)=a;
c(2,:)=b;
for m=3:l
for n=1:l1-1
c(m,n)=-(1/c(m-1,1))*det([c((m-2),1) c((m-2),(n+1));c((m-1),1) c((m-
1),(n+1))]);
end
end
disp('the routh matrix:')
disp(c)
%%now we check the stablity of system
if c(:,1)>0
disp('System is Stable')
else
disp('System is Unstable');
end
Root Locus & Bode Plot
cltfz = feedback( hd,1)
step(cltfz)
rlocus(cltfz)
bode(cltfz)

Unstable root Locus

Unstable bode plots

Placing PID
Output of the System
Output of the system also track step input .

Root locus

Bode Plot
By changing in PID contoller
Kp=56 ki=15

Kp=5ki=15
Kp=200ki =15

Ki = 2 Kp=20

Ki = 50 Kp=20
Ki = 100 Kp=2

Question

Apply Routh Criteria on your system and find the range of stability for different values of sampling
interval ‘T’ and controller Gains (PID gain values). [Marks: 06]
c) Apply Jury’s Test on your system and find the range of stability for different values of sampling
interval ‘T’ and controller Gains (PID gain values). [Marks: 06]
d) Analyze the steady state error efficiency for your system. [Marks: 04]
e) All your assignments given previously for DCS time-by-time should be submitted. [Marks: 05]
f) Plot root locus (Bonus Marks)
g) Plot bode plot (Bonus Marks)
h) And any other analysis you learnt etc. (Bonus Marks)
'routh_hurwitz'
e=input('enter the coefficients of characteristic equation: ');
disp('------------------------------------------------------------------------
-')
l=length(e);
m=mod(l,2);
if m==0
a=rand(1,(l/2));
b=rand(1,(l/2));
for i=1:(l/2)
a(i)=e((2*i)-1);
b(i)=e(2*i);
end
else
e1=[e 0];
a=rand(1,((l+1)/2));
b=[rand(1,((l-1)/2)),0];
for i=1:((l+1)/2)
a(i)=e1((2*i)-1);
b(i)=e1(2*i);
end
end
%%now we genrate the remaining rows of routh matrix
l1=length(a);
c=zeros(l,l1);
c(1,:)=a;
c(2,:)=b;
for m=3:l
for n=1:l1-1
c(m,n)=-(1/c(m-1,1))*det([c((m-2),1) c((m-2),(n+1));c((m-1),1) c((m-
1),(n+1))]);
end
end
disp('the routh matrix:')
disp(c)
%%now we check the stablity of system
if c(:,1)>0
disp('System is Stable')
else
disp('System is Unstable');
end

Root Locus
cltfz = feedback( hd,1)
step(cltfz)
rlocus(cltfz)

System is stable on Kp=1 and Ki=1.

Kp=20 Ki=15
System is unstable on Kp=20 and Ki=15.

Kp=0.2 ki=1 (on Stable )

Kp=0.2 ki=0.5 (stable )

Kp=50 ki=15 (unstable)

Sampled time
Kp=1 ki=1
Sampled time = 0.001

Kp=1 ki=1
Sampled time = 0.01

Jury’s Test
clc;
P=[1 -1.231 0.2891];
[M,L] = jury(P);
if nargin <2
N=10;
end
P=P(1)/abs(P(1))*P;
P=sym(P);
n=length(P);
M=sym(zeros(2*n,n));
M(1,:)=P(n:-1:1);
M(2,:)=P(1:1:n);
syms epsilon
for i=1:n-1
for j=1:n-i
M(2*i+1,j)=det( [M(2*i-1,1) M(2*i-1,j+1);M(2*i,1) M(2*i,j+1)] )/(-
M(2*i,1));
M(2*i+2,n-i+1-j)=M(2*i+1,j);
end
if isempty(symvar(M(2*i+2,1))) == 1 % There is a special case only if the
first element is not a variable
if isempty(symvar(M(2*i+2,:))) == 1 % Special cases when all values
are numbers
if sum(abs(double(M(2*i+1,:)))) < 10^(-N) % Special case of a row
of zeros
disp(['Row of zeros: ' num2str(2*i+1) ', ' num2str(2*i+2) '.
Order of auxiliar polynomial: ' num2str(n-i) ]);
for j=1:n-i
M(2*i+1,j)=M(2*i-1,j+1)*j; % Derivative of auxiliar
polynomial
M(2*i+2,n-i+1-j)=M(2*i+1,j);
end
sig=M(2*i+1,n-i)/abs(M(2*i+1,n-i)); % Sign change of auxiliar
polynomial
M(2*i+1,:)=sig*M(2*i+1,:);
M(2*i+2,:)=sig*M(2*i+2,:);
elseif abs(double(M(2*i+2,1)))<10^(-N) &&
sum(abs(double(M(2*i,2:n)))) > 10^(-N) % Special case when the first element
of the row iz zero and the others are values
M(2*i+2,1)=epsilon;
M(2*i+1,n-i)=epsilon;
end
elseif abs(double(M(2*i+2,1)))<10^(-N) % Special case when the first
element of the row iz zero and some other elements are variables
M(2*i+2,1)=epsilon;
M(2*i+1,n-i)=epsilon;
end
end
end
M=vpa(M);
M1=limit(M,epsilon,0,'right'); % RD_NINF = -inf, RD_INF = inf
M1=vpa(M1);
for i=1:n-1
L(i,1)=M1(2*i+2,1);