% This program is used to explain how PCA works and the idea behind the% method. This example is given in the ISA book (Chapter 8, simple example)% The process is composed of a valve a thermometer and a fermenter. The% valve is used to control the flow rate of cooling water. The thermometer% is used to measure the temperature of the fermenter.% The steady state of the system is: valve=10% max. Temperature T=25 D.% YANG ZHANG Email: yz@che.utexas.edu% 2006-April, UTexas-Austinclear all% Give valve set points and fermenter temperature.% Valve=10+rand(1,10)-0.5; % ten value of Valve around mean.% Temp=25-(Valve-10)+(rand(1,10)-0.5)*0.5; % 10 value of Temp: T=2.5*valve + noiseTemp = [2.5000e+001 2.5000e+001 2.5000e+001 2.5300e+001 2.5200e+001 2.5200e+001 2.4800e+001 2.4700e+001 2.5100e+001 2.4700e+001];Valve = [1.0100e+001 9.9000e+000 1.0000e+001 9.8000e+000 9.9000e+000 9.8000e+000 1.0100e+001 1.0300e+001 1.0000e+001 1.0100e+001];% Define databasex=[Valve;Temp]';figure(1)plot(x(:,1),x(:,2),'ko')xlabel('Valve Position (%)')ylabel('Fermenter Temperature (C)')% Calculate the covariance of dataset x (relationship between Valve and% Temperature)covariance=cov(x);% Singular Value Decomposition x=T*E*P% E is the eignvalue of x; U is orthnormal; V=U'.%[U E V]=svd(covariance);%U=U'; V=V';[U E] = eig(covariance);E = diag([E(2,2),E(1,1)]); % arrange the Eignvalue in a decending orderU = [U(:,2),U(:,1)]; % arrange the Eignvector to make it compatable with Eignvalue.% Mean center the raw dataValvebar=x(:,1)-mean(x(:,1));Tempbar=x(:,2)-mean(x(:,2));xbar=[Valvebar,Tempbar];% Principal axis rotation of the covariance matrix.z=U'*xbar';z=z';%covz is the covariance matrix of principal components which is equal to E%matrix. The variance of transformed variable (z1 & z2) will have variance%0.2283 and 0.0123 respectively.covz=U'*covariance*U;% The first column of U factor is - and + with nearly equal value which% means the first principal component is related to variability which both% measurements have difference. The second column is both positive which% means the 2nd PC is concentrate on the varaiability of the common between% the two.% scaling of PCS.for i=1:2Vs(:,i)=sqrt(E(i,i))*U(:,i);Ws(:,i)=U(:,i)/sqrt(E(i,i));end% after rescale: Vs'*Vs=E; V'*COV*Vs=E^2; Ws'*Ws=inv(E); Ws'*COV*Ws=Iy=Ws'*xbar';%define T2T2=diag(y'*y);% find the control limit of T2p=2; %PC numbern=10; % sample numberT2Upper=p*(n-1)/(n-p)*finv(0.95,p,n-p);% draw the T2 control chartsfigure(2)plot(T2,'ko')hold onplot(T2Upper*ones(1,13),'k')xlabel('Sample Number')ylabel('T2 Score')hold off% Draw the control ellipse% s1^2*s2^2/(s1^2*s2^2-s12^2)*[(x1-x1mean)^2/s1^2+(x2-x2mean)^2/s2^2-2*s12(% x1-x1mean)(x2-x2mean)/s1^2/s2^2]=T2^2(upper control limit)% S=covariance;% a=S(1)*S(4)/(S(1)*S(4)-S(2)^2);% b=xbar(:,1).^2/S(1)+xbar(:,2).^2/S(4)-2*S(2)*xbar(:,1).*xbar(:,2)./S(1)/S(4);% error detection% two new points is added to be detected which are [10.1,25.2];[9.9;24];% plot two values with the original 10.Valve2=[10.2,9.9];Temp2=[25.1,25.5];figure(1)hold onplot(Valve2,Temp2,'k*')hold off% plot the control charts for valve and temperature separately.figure(3)hold onplot(1:10,Valve,'ko')plot(11:12,Valve2,'k*')plot(0:13, mean(Valve)*ones(1,14),'k')plot(0:13, (mean(Valve)-std(Valve)*2)*ones(1,14),'k')plot(0:13, (mean(Valve)+std(Valve)*2)*ones(1,14),'k')xlabel('Sample Number')ylabel('Valve Position (%)')hold offfigure(4)hold onplot(1:10,Temp,'ko')plot(11:12,Temp2,'k*')plot(0:13, mean(Temp)*ones(1,14),'k')plot(0:13, (mean(Temp)+std(Temp)*2)*ones(1,14),'k')plot(0:13, (mean(Temp)-std(Temp)*2)*ones(1,14),'k')xlabel('Sample Number')ylabel('Temperature (C)')hold off% transform new observations to PC axisy2=Ws'*([Valve2-mean(Valve);Temp2-mean(Temp)]);% plot the control charts for principal componentsfigure(5)hold onplot(1:10,y(1,:),'ko')plot(11:12,y2(1,:),'k*')plot(0:13, mean(y(1,:))*ones(1,14),'k')plot(0:13, -2*ones(1,14),'k')plot(0:13, 2*ones(1,14),'k')xlabel ('Sample Number')ylabel ('Principal Component 1')hold offfigure(6)hold onplot(1:10,y(2,:),'ko')plot(11:12,y2(2,:),'k*')plot(0:13,mean(y(2,:))*ones(1,14),'k')plot(0:13, -2*ones(1,14),'k')plot(0:13, 2*ones(1,14),'k')xlabel ('Sample Number')ylabel ('Principal Component 2')hold off% plot the T2 chart for the two new pointsfigure(2)T22=diag(y2'*y2);hold onplot([11,12],T22,'k*')hold off% draw the control ellipsera = sqrt(T2Upper*E(1,1));rb = sqrt(T2Upper*E(2,2));V=[sqrt(E(1,1))*U(:,1),sqrt(E(2,2))*-U(:,2)];k1 = V(2,1)/V(1,1);k2 = V(2,2)/V(1,2);ang = atan(k1); %atan(k1);x0 = mean(x(:,1));y0 = mean(x(:,2));% i = 9.5:0.01:10.5;% s = 9.72:0.01:10.3;% j = k1*i + 42.6303;% l = k2*s + 19.3279;figure (7)h=ellipse(ra,rb,ang,x0,y0,'k',1000);hold onplot(x(:,1),x(:,2),'ko')plot(Valve2, Temp2,'k*')% plot(i,j)% plot(s,l)xlabel ('Valve Position')ylabel ('Temperature')hold off
