Table of Contents

................................................................................................................................. 1
Part 1 ........................................................................................................................ 1
Part 3 ........................................................................................................................ 8
Part 4 ...................................................................................................................... 10
Part 5 ...................................................................................................................... 11
Part 6 ...................................................................................................................... 12

%Initialization
clear all; close all; clc;
dim = 6;
n = dim;
m = 2;
p = n/m;
r = 3;
I = 1; L = 6; GI = 1;
M = I*L/18*[4 1 0;1 4 1;0 1 2];
K = 3*GI/L*[2 -1 0;-1 2 -1;0 -1 1];
C = 0.4*K;
B1 = L/6*[1 1 0;0 1 1;0 0 1];
A = [zeros(3,3) eye(3);-inv(M)*K -inv(M)*C];
B = [zeros(3,3);B1];
C1 = -inv(M)*K; C2 = -inv(M)*C; %C3 = inv(M)*B1;
C = [C1(2:3,1:3) C2(2:3,1:3)]; D = B1(2:3,1:3);

Time sampling and measurements.

Np = 10240; Ne = 100;
fs = 1;
Tf = 1024*50;
dt = 1/fs;
t = 0: dt : Tf;
Np1 = length(t);

Part 1
discrete time models

Ad = expm(A*dt);
eig(Ad);
Bd = inv(A)*(Ad - eye(dim))*B;
combination = [1 0 0;0 1 0;0 0 1];
G_Matrix = zeros(2,3,Np);
ur1 = zeros(1, Np1);
ur1(1,1:Ne) = randn(1,Ne).*hanning(Ne).';
ur2 = zeros(1, Np1);
ur2(1,1:Ne) = randn(1,Ne).*hanning(Ne).';
ur3 = zeros(1, Np1);
ur3(1,1:Ne) = randn(1,Ne).*hanning(Ne).';
for i = 1:3

1
index = combination(i,:);
ur = [index(1)*ur1;index(2)*ur2;index(3)*ur3];
%generation of measurements
xk = zeros(dim,1);
for j = 1 : Np1
y(1:2,j) = C*xk + D*ur(:, j);
xk = Ad*xk + Bd*ur(:, j);
end
figure
subplot(311)
plot(t, ur(i,1:Np1));
xlim([0,Tf])
ylabel(['input(u',num2str(i),')']);
subplot(312)
plot(t, y(1,:))
ylabel('output(y1)');
xlim([0,Tf])
subplot(313)
plot(t, y(2,:))
ylabel('output(y2)');
xlim([0,Tf])
% Frequency response Function
Us = fft(ur(i,1:end), Np)/Np;
Ys1 = fft(y(1,1:end), Np)/Np;
Ys2 = fft(y(2,1:end), Np)/Np;
w = 0.5*fs*linspace(0,1,Np/2+1);
figure
subplot(311)
semilogy(w, abs(Us(1:Np/2+1)));
ylabel(['fft(input', num2str(i),')']);
subplot(312)
semilogy(w, abs(Ys1(1:Np/2+1)));
ylabel('fft(output 1)');
subplot(313)
semilogy(w, abs(Ys2(1:Np/2+1)));
ylabel('fft(output 2)');
YU1 = Ys1.*conj(Us);
YU2 = Ys2.*conj(Us);
UU = Us.*conj(Us);
Gs1 = YU1./UU;
Gs2 = YU2./UU;
Gs_phi1 = angle(Gs1);
Gsa1 = abs(Gs1);
Gs_phi2 = angle(Gs2);
Gsa2 = abs(Gs2);
figure
subplot(221)
semilogy(w, Gsa1(1:Np/2+1));
ylabel(['Frequency Response(y_1 vs u_',num2str(i),')']);
subplot(222)
semilogy(w, Gsa2(1:Np/2+1));
ylabel(['Frequency Response(y_2 vs u_',num2str(i),')']);
subplot(223)
plot(w, Gs_phi1(1:Np/2+1));
ylabel(['Phase(y_1 vs u_',num2str(i),')']);
subplot(224)

2
 plot(w, Gs_phi2(1:Np/2+1));
ylabel(['Phase(y_2 vs u_',num2str(i),')']);
axis('auto')
G_Matrix(1,i,:) = Gs1;
G_Matrix(2,i,:) = Gs2;
end

3
4
5
6
7
Part 3
Generation of Phi Matrix

Gz = cell(Np,1);
Phi1 = cell(2*p+1,Np);
for i = 1:Np
Phi1{p+1,i} = eye(r);
end
for i = 1:p
for j = 1:Np
z = exp(sqrt(-1)*2*pi*(j-1)/Np);
Gz{j} = reshape(G_Matrix(:,:,j),m,r)*z^-i;
%Gz{j} = (C*inv(z.*eye(n) - A)*B + D)*z^-i;
Phi1{i,j} = Gz{j};
Phi1{i+p+1,j} = eye(r)*z^-i;
end
end
% Real part of Phi Matrix
for i = 1:p
for j = 1:Np
Phi_R(m*(i-1)+1:m*i,r*(j-1)+1:r*j) = real(Phi1{i,j});
end
end
row = m*p;
for i = p+1:2*p+1
for j = 1:Np
Phi_R(row+1:row+r,r*(j-1)+1:r*j) = real(Phi1{i,j});
end
row = row + r;
end
% Imaginary part of Phi Matrix
for i = 1:p
for j = 1:Np
Phi_I(m*(i-1)+1:m*i,r*(j-1)+1:r*j) = imag(Phi1{i,j});
end
end
row = m*p;
for i = p+1:2*p+1
for j = 1:Np
Phi_I(row+1:row+r,r*(j-1)+1:r*j) = imag(Phi1{i,j});
end
row = row + r;
end
% Real part of Psi Matrix
for i = 1:Np
Psi_R(1:m,r*(i-1)+1:r*i) = reshape(real(G_Matrix(:,:,i)),m,r);
end
% Imaginary part of Phi Matrix
for i = 1:Np
Psi_I(1:m,r*(i-1)+1:r*i) = reshape(imag(G_Matrix(:,:,i)),m,r);
end

Calculation of Theta Matrix(Coefficients)

Theta = [Psi_R Psi_I]/([Phi_R Phi_I]);

8
for i = 1:p
Q{i} = -Theta(1:m,m*(i-1)+1:m*i);
end
R0 = Theta(1:m,m*p+1:m*p+r);
col = m*p+r;
for i = 1:p
R{i} = Theta(1:m,col+1:col+r);
col = col + 1;
end
Q{1}
Q{2}
Q{3}
R0
R{1}
R{2}
R{3}

ans =
0.763352861591927 14.395711989723697
-1.428488653716280 -3.895932709614300

ans =
-16.925545562969486 -16.300314477740056
3.752569388875708 3.870085180206786

ans =
6.441753869252735 10.664605741673174
-1.261825917586711 -2.047346904217889

R0 =
0.000000000000002 1.000000000000012 1.000000000000006
-0.000000000000001 -0.000000000000000 1.000000000000008

ans =
0.314998527415210 0.470890489449626 14.866602479173133
-0.039733936074097 -0.838225534959949 -5.282018240020592

ans =
0.470890489449626 14.866602479173133 -0.299885755312902
-0.838225534959949 -5.282018240020592 0.039176408484353

ans =
14.866602479173133 -0.299885755312902 -8.799359478370604
-5.282018240020592 0.039176408484353 2.268385701736137

9
Part 4
State Space Realizations

D_hat = R0;
C_hat = zeros(m,m*p);
C_hat(1:m,m*p-m+1:m*p) = eye(m);
B_hat = zeros(m*p,r);
for i = 1:p
B_hat(m*(i-1)+1:m*i,1:r) = R{p-i+1} - Q{p-i+1}*R0;
end
A_hat = zeros(m*p,m*p);
A_hat(1:m,m*p-m+1:m*p) = -Q{p};
for i = 2:p
A_hat(m*(i-1)+1:m*(i-1)+m,m*p-m+1:m*p) = -Q{p-i+1};
A_hat(m*(i-1)+1:m*(i-1)+m,(i-1)*m-m+1:(i-1)*m) = eye(m);
end
fprintf(1,'\n\t\t\tState Space Matrices\n\t\t**************************************
A_hat
B_hat
C_hat
D_hat

State Space Matrices

****************************************
A_hat =
Columns 1 through 4
0 0 0 0
0 0 0 0
1.000000000000000 0 0 0
0 1.000000000000000 0 0
0 0 1.000000000000000 0
0 0 0 1.000000000000000
Columns 5 through 6
-6.441753869252735 -10.664605741673174
1.261825917586711 2.047346904217889
16.925545562969486 16.300314477740056
-3.752569388875708 -3.870085180206786
-0.763352861591927 -14.395711989723697
1.428488653716280 3.895932709614300

B_hat =
14.866602479173125 -6.741639624565713 -25.905719089296642
-5.282018240020590 1.301002326071079 5.577558523540761
0.470890489449654 31.792148042142820 32.925974285396883
-0.838225534959955 -9.034587628896343 -7.583478160598196
0.314998527415216 -0.292462372142308 -0.292462372142616
-0.039733936074096 0.590263118756347 0.042403123310030

C_hat =
0 0 0 0 1 0
0 0 0 0 0 1

10
 D_hat =
0.000000000000002 1.000000000000012 1.000000000000006
-0.000000000000001 -0.000000000000000 1.000000000000008

Part 5
Error Analysis in EigenValues

fprintf(1,'\n\t\t\tEigenValues of the Original Matrix(Discrete Time)\n\t\t*********

eig(Ad)
fprintf(1,'\n\t\t\tEigenValues of the Estimated Matrix(Discrete Time)\n\t\t********
eig(A_hat)
fprintf(1,'\n\t\t\tError in Real part of EigenValues\n\t\t*************************
diffEig = real(eig(Ad))-real(eig(A_hat))
fprintf(1,'\n\tNorm of error in real Part of Eigenvalues = %g',norm(diffEig,'fro'))
fprintf(1,'\n\n\tRelative Norm of error in real Part of Eigenvalues = %g',norm(diff
fprintf(1,'\n\n\t\t\tError in Imaginary part of EigenValues\n\t\t******************
diffEig = imag(eig(Ad))-imag(eig(A_hat))
fprintf(1,'\n\tNorm of error in imaginary Part of Eigenvalues = %g',norm(diffEig,'f
fprintf(1,'\n\n\tRelative Norm of error in imaginary Part of Eigenvalues = %g\n',no

EigenValues of the Original Matrix(Discrete Time)

********************************************************
ans =
0.048337788201464 + 0.608466040156158i
0.048337788201464 - 0.608466040156158i
0.566151729396097 + 0.648298110427239i
0.566151729396097 - 0.648298110427239i
0.951800406413586 + 0.257717279711803i
0.951800406413586 - 0.257717279711803i

EigenValues of the Estimated Matrix(Discrete Time)

********************************************************
ans =
0.048337788201437 + 0.608466040156100i
0.048337788201437 - 0.608466040156100i
0.566151729396612 + 0.648298110427191i
0.566151729396612 - 0.648298110427191i
0.951800406413140 + 0.257717279711516i
0.951800406413140 - 0.257717279711516i

Error in Real part of EigenValues

*************************************************
diffEig =
1.0e-012 *

11
 0.027512714328992
0.027512714328992
-0.514699394216223
-0.514699394216223
0.445754544387000
0.445754544387000

Norm of error in real Part of Eigenvalues = 9.63711e-013

Relative Norm of error in real Part of Eigenvalues = 6.14743e-013
Error in Imaginary part of EigenValues
*************************************************
diffEig =
1.0e-012 *
0.058619775700208
-0.058619775700208
0.047517545453957
-0.047517545453957
0.287159185319297
-0.287159185319297

Norm of error in imaginary Part of Eigenvalues = 4.19892e-013

Relative Norm of error in imaginary Part of Eigenvalues = 3.20736e-013

Part 6
Error Analysis in Direct Transmission Term

D_hat
D
fprintf(1,'\n\t\t\tError in Direct Transmission\n\t\t******************************
diffD = D_hat - D
fprintf(1,'\n\tNorm of error in Direct Transmission = %g',norm(diffD,'fro'));
fprintf(1,'\n\n\tRelative Norm of error in Direct Transmission = %g\n',norm(diffD,'

D_hat =
0.000000000000002 1.000000000000012 1.000000000000006
-0.000000000000001 -0.000000000000000 1.000000000000008

D =
0 1 1
0 0 1

Error in Direct Transmission

*************************************************
diffD =

12
 1.0e-013 *
0.021901643927287 0.119904086659517 0.059952043329758
-0.005356080577116 -0.001492109116083 0.084376949871512

Norm of error in Direct Transmission = 1.60004e-014

Relative Norm of error in Direct Transmission = 9.23784e-015

Published with MATLAB® 7.9

13