ASSIGNMENT-4
%Assignment submitted by
%Name: Jeswin Tom Joseph
%Roll_no: CH17B052
1. The transfer function is
%initialising the transfer function
Gs= zpk(-0.5,[-1/20 -0.1 -0.2 -1],1/1000,'InputDelay',3)
Gs =
0.001 (s+0.5)
exp(-3*s) * ------------------------------
(s+0.05) (s+0.1) (s+0.2) (s+1)
Continuous-time zero/pole/gain model.
%step response for the given transfer function
[ystep,tk] = step(Gs,(0:0.01:300));
a. For Krishnaswamy and Sundareshan model for FOPTD,
where
Kp = ystep(end);
[~, t1] = min(abs(ystep-.353*ystep(end)));
[~, t2] = min(abs(ystep-0.853*ystep(end)));
D = 1.3*tk(t1) - 0.29*tk(t2);
tau = 0.67*(tk(t2) - tk(t1));
G_ks= tf(ystep(end),[tau 1],'ioDelay',D)
G_ks =
0.5
exp(-15.4*s) * -----------
22.79 s + 1
Continuous-time transfer function.
step(Gs,G_ks)
legend({'Process','Model'},'Location','best')
By calculation simulink, and that of variable ystep,
We have
b. Skogestad's Half rule
The transfer function is
Since we have a positive numerator term, we have to cancel it with a corresponding
denominator term.
Here, nearest neighbour for the term 2s+1 is s+1.
By Rule T1b: Skogestad Half Rule
Hence,
For FOPTD,
For SOPTD
G_sk_F=tf(0.5,[25 1],'ioDelay',13)
G_sk_F =
0.5
exp(-13*s) * --------
25 s + 1
Continuous-time transfer function.
G_sk_S=tf(0.5,[250 32.5 1],'ioDelay',5.5)
G_sk_S =
0.5
exp(-5.5*s) * --------------------
250 s^2 + 32.5 s + 1
Continuous-time transfer function.
step(Gs,G_sk_F,G_sk_S)
legend({'Process','FOPTD Model','SOPTD Model'},'Location','best')
c. For least square approximation of SOPTD,
%Initialising Bode plots and hence finding frequency response
[Gmag,Gphase,wvec]=bode(Gs,0:0.001:300);
mag = lsqcurvefit(@(mpar,wdata) modpred(mpar,wdata),[1 1
1],wvec,squeeze(Gmag));
Local minimum possible.
lsqcurvefit stopped because the final change in the sum of squares relative to
its initial value is less than the value of the function tolerance.
<stopping criteria details>
G_ls = tf(mag(1),[mag(2)*mag(2) 2*mag(2)*mag(3) 1]);
Gain=mag(1);
taup=mag(2);
zeta=mag(3);
Phase1 = lsqcurvefit(@(mpar,wdata)
phasepred(mpar,wdata,Gain,taup,zeta),4,wvec,cos(squeeze(Gmag)));
Solver stopped prematurely.
lsqcurvefit stopped because it exceeded the function evaluation limit,
options.MaxFunctionEvaluations = 1.000000e+02.
G_ls.iodelay=Phase1
G_ls =
0.4953
exp(-4*s) * -----------------------
272.1 s^2 + 31.29 s + 1
Continuous-time transfer function.
step(Gs,G_ls)
legend({'Process','Least Square Model'},'Location','best')
%function for least square fit given below
function yhat = modpred(mpar,wdata)
% Function for least squares fit
Kp = mpar(1);
taup1= mpar(2);
zeta1=mpar(3);
den=sqrt((1-(taup1^2).*(wdata.^2)).^2+4*(zeta1*taup1.*wdata).^2);
yhat=Kp./den;
end
function ph=phasepred(mpar,wdata,~,taup1,zeta1)
D=mpar(1);
ph=cos(-D*wdata-angle(1-(((taup1).*wdata).^2)+1i*2*zeta1*(taup1)));
end
d. Comparison
step(Gs,G_ks,G_sk_F,G_sk_S,G_ls)
legend({'Process','Krishnaswamy & Sundareshan','Skogestad FOPTD','Skogestad
SOPTD','Least Square SOPTD'},'Location','southeast')
bode(Gs,G_ks,G_sk_F,G_sk_S,G_ls)
legend({'Process','Krishnaswamy & Sundareshan','Skogestad FOPTD','Skogestad
SOPTD','Least Square SOPTD'},'Location','southwest')
2.
a. Simulink Block diagram of the process with delay and ZOH controller is given below
b. Codes for creating dataset, and partitioning it
%Creating input datasignal
B=0.2;
Ts=0.2;
u=idinput(2555,'prbs',[0 B],[-1 1]);
%passing the created datasignal to SIMULINK
in_data=[(0:1:length(u)-1)'*Ts (u)];
%SIMULINK is run, and the output dataset generated
%Accessing the first 2555 datapoints
in_dup=out.in.data(1:2555);
out_dup=out.o.data(1:2555);
%creation of iddata object
data=iddata(out_dup,in_dup,1,'outputname','output','inputname','input');
%plotting input-output data profile
plot(data)
%Partioning data into trainingset and testset
train_data=data(1:1500);
test_data=data(1501:end);
%removing mean from the training data set
[ztrain,Tr]=detrend(train_data,0);
ztest=detrend(test_data,Tr);
c. Code for Fitting the given FIR Model to training data set is given below
sysimpest=impulseest(ztrain);
%calculating and plotting the impulse response of the estimated model.
figure
impulse(sysimpest,'sd',6);
%using 6 as sd to get a higher confidence interval(99.999%)
Analysing the IR coefficients of the variable sysimpest, we find that the system is stable with
an input-output delay of 6.
d. Code for calculating the step response and estimating the gain of the data
%calculating and plotting the step response of the estimated model
figure
step(sysimpest);
sysimp=impulse(sysimpest,'sd',6);
systep=step(sysimpest);
%estimating the steady state gain of data
ss_gain=systep(end)
ss_gain = 1.7521
• The steady state gain of the data is 1.7521
• The system is not oscillatory, and is sluggish, hence it is underdamped.
• The system is stable with an input-output delay of 6.
e. Code for estimating parameters of iven model is given below
%using oe routine to find the parameters of given model
model_oe=oe(ztrain,[1,2,6]);
From the analysis of impulse and step responses for the model
we have
m=1; n=2; and d=6;
f. Code for analyzing goodness of estimated model using resid routine and comparison with
test data set
%Analzing of goodness of data
figure
resid(model_oe,ztrain);
From the plots, it is clear that there is no significant correlation between residuals and inputs.
Also, there is no correlation netween samples, as given by auto-correlation function plot.
%Examining errors in parameter estimates
present(model_oe);
model_oe =
Discrete-time OE model: y(t) = [B(z)/F(z)]u(t) + e(t)
B(z) = 0.01621 (+/- 0.0009147) z^-6
F(z) = 1 - 1.749 (+/- 0.01466) z^-1 + 0.7571 (+/- 0.0143) z^-2
Sample time: 1 seconds
Parameterization:
Polynomial orders: nb=1 nf=2 nk=6
Number of free coefficients: 3
Use "polydata", "getpvec", "getcov" for parameters and their uncertainties.
Status:
Termination condition: Near (local) minimum, (norm(g) < tol)..
Number of iterations: 5, Number of function evaluations: 13
Estimated using OE on time domain data "ztrain".
Fit to estimation data: 48.84%
FPE: 0.1059, MSE: 0.1055
More information in model's "Report" property.
g. Code for prediction of estimated model on test data
%Comparison of estimated model with test data
figure
compare(model_oe,ztest);
The OE(1,2,6) is written in the transfer function form as:
where
Hence,
