sieng_parity_a.m

%----System Matrices for linearized SI engine model----%
A1=[-1.6688 4.1250;-.2926 -15.8177];
B1=[0 410.3077;55.6064 0];
C1=[1 0; 0 1;0 -.6655];
D1=[0 0 ;0 0;10.3995 0];
d=D1(:,1);
H=[0 ;0 ];

%-----Disturbance Matrix---%
E=[-23.3822; 0];

%-----Conversion to Discrete time model----%
[A B C D]=c2dM(A1,B1,C1,D1,.01);
p=[-2 -18];
K = place(A1,B1,p);

%----History of states, inputs and outputs----%
s=4;
t=1000;
x=zeros(2,1001);
y=zeros(3,1001);
u=zeros(2,1001);
x(:,1)=[1;1];
u(:,1)=-K*x(:,1);
y(:,1)=C*x(:,1)+D*u(:,1);
for i=1:t;
    x(:,i+1)=A*x(:,i)+B*u(:,i);
    u(:,i+1)=-K*x(:,i+1);
    y(:,i+1)=C*x(:,i+1)+D*u(:,i+1);
end
Y=zeros(15,501);
for i=1:501
    Y(:,i)=[y(:,i);y(:,i+1);y(:,i+2);y(:,i+3);y(:,i+4)];
end
U=zeros(10,501);
for i=1:501
    U(:,i)=[u(:,i);u(:,i+1);u(:,i+2);u(:,i+3);u(:,i+4)];
end

%----Calculation of different Fault distribution matices----%
Q1=[D;C*B;C*A*B;C*A^2*B;C*A^3*B];
Q2=[zeros(size(D));D;C*B;C*A*B;C*A^2*B];
Q3=[zeros(size(D));zeros(size(D));D;C*B;C*A*B];
Q4=[zeros(size(D));zeros(size(D));zeros(size(D));D;C*B];
Q5=[zeros(size(D));zeros(size(D));zeros(size(D));zeros(size(D));D];
Q=[Q1 Q2 Q3 Q4 Q5];

%-----Actuator Fault Distribution matrix---%
Qa1=[D;C*B;C*A*B;C*A^2*B;C*A^3*B];
Qa2=[zeros(size(D));D;C*B;C*A*B;C*A^2*B];
Qa3=[zeros(size(D));zeros(size(D));D;C*B;C*A*B];
Qa4=[zeros(size(D));zeros(size(D));zeros(size(D));D;C*B];
Qa5=[zeros(size(D));zeros(size(D));zeros(size(D));zeros(size(D));D];
Qa=[Qa1 Qa2 Qa3 Qa4 Qa5];

%-----Componenet Fault Distribution matrix---%
Tc1=[d;C*H;C*A*H;C*A^2*H;C*A^3*H];
Tc2=[zeros(size(d));d;C*H;C*A*H;C*A^2*H];
Tc3=[zeros(size(d));zeros(size(d));d;C*H;C*A*H];
Tc4=[zeros(size(d));zeros(size(d));zeros(size(d));d;C*H];
Tc5=[zeros(size(d));zeros(size(d));zeros(size(d));zeros(size(d));d];
Tc=[Tc1 Tc2 Tc3 Tc4 Tc5];
 
%-----Diturbance Distribution matrix---%
Td1=[d;C*E;C*A*E;C*A^2*E;C*A^3*E];
Td2=[zeros(size(d));d;C*E;C*A*E;C*A^2*E];
Td3=[zeros(size(d));zeros(size(d));d;C*E;C*A*E];
Td4=[zeros(size(d));zeros(size(d));zeros(size(d));d;C*E];
Td5=[zeros(size(d));zeros(size(d));zeros(size(d));zeros(size(d));d];
Td=[Td1 Td2 Td3 Td4 Td5];

%-----Calculation of weight matrix 'w'---%
R=[C;C*A;C*A^2;C*A^3;C*A^4];
 w=null(R');
 w1=w(:,3);
 w1([1,4,7,10,13],1)=0;
 w2=w(:,2);
 w2([2,5,8,11,14],1)=0;
 w3=w(:,3);

%------Fault Matrices----%
fault=zeros(1,501);
Fa=zeros(10,501);
Fc=zeros(5,501);
Dt=zeros(5,501);
S=zeros(15,501);
%------Fault Creation in Throttle acuator----%
fault(:,350:501)=5;
Fa(1,:)=fault;

%-----Residual Calculation---%
r_x=w1'*((Y-Q*U)-Qa*Fa-Tc*Fc-Td*Dt-S);
r_phi=w2'*((Y-Q*U)-Qa*Fa-Tc*Fc-Td*Dt-S);
r_a=w3'*((Y-Q*U)-Tc*Fc-Td*Dt-S);

%-----Plot of Residuals---%
figure;
subplot(2,2,1);plot(rest,r_x,rest,2.5e-3*l,'r--',rest,-2.5e-3*l,'r--'),grid, ylabel('res_s1'),title('Residual for Speed sensor');
subplot(2,2,2);plot(rest,r_s2,rest,2.5e-3*l,'r--',rest,-6e-3*l,'r--'),grid, ylabel('res_s2'),title('Residual for pman sensor');
subplot(2,2,3);plot(rest,r_a),grid, ylabel('res_a'),xlabel('time'),title('Residual for Throttle actuator fault');
subplot(2,2,4);plot(rest,fault),grid, ylabel('fault'),xlabel('time'),title('Fault in Throttle Acutator');