operator_est.m

一个用MATLAB编写的优化控制工具箱

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Construction of linear and fuzzy controllers from operator data for a chemical 
% plant.
%
% By: Kevin Passino
% Version: 8/20/99
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear

load operator.dat		% For information on the data set see this file.

yy=operator(:,6);  % The settings made by the operator
y=yy';

xx=operator(:,1:5); 
x=xx';
                % The data input from the plant that the operator uses to make decisions, and in particular
		%  Inputs x holds u1(k),u2(k),u3(k),u4(k),u5(k)
		
time=0:length(y)-1;  % Set up the time vector for the data

M=length(y);

%%%%%%%%%%%%%%%
% Plot the data (this is what is used for training):

figure(1)
clf
subplot(321)
plot(time,x(1,:),'k-')  
ylabel('u_1(k)')
axis([min(time) max(time) min(x(1,:)) max(x(1,:))])
zoom
subplot(322)
plot(time,x(2,:),'k-')  
ylabel('u_2(k)')
axis([min(time) max(time) min(x(2,:)) max(x(2,:))])
zoom
subplot(323)
plot(time,x(3,:),'k-')  
ylabel('u_3(k)')
axis([min(time) max(time) min(x(3,:)) max(x(3,:))])
zoom
subplot(324)
plot(time,x(4,:),'k-')  
ylabel('u_4(k)')
axis([min(time) max(time) min(x(4,:)) max(x(4,:))])
zoom
subplot(325)
plot(time,x(5,:),'k-')  
ylabel('u_5(k)')
xlabel('Time, k')
axis([min(time) max(time) min(x(5,:)) max(x(5,:))])
zoom
subplot(326)
plot(time,y,'k-')
ylabel('y(k)')
xlabel('Time, k')
axis([min(time) max(time) min(y) max(y)])


% Next, we form the test data.  Note that we only have M=70 training data points.
% We do not have an ability to generate more data from the physical system; hence,
% we will generate MGamma=69 points which are "in between" the given data.  To 
% do this we simply place points at the average of the two values they lie between:

for i=1:M-1
	xt(:,i)=0.5*(x(:,i)+x(:,i+1));
	yt(i)=0.5*(y(i)+y(i+1));
end

MGamma=length(yt);
timet=0:length(yt)-1;  

%%%%%%%%%%%%%%%
% Plot the data (this is what is used for testing):

figure(2)
clf
subplot(321)
plot(timet,xt(1,:),'k-')  
ylabel('u_1(k)')
axis([min(timet) max(timet) min(xt(1,:)) max(xt(1,:))])
zoom
subplot(322)
plot(timet,xt(2,:),'k-')  
ylabel('u_2(k)')
axis([min(timet) max(timet) min(xt(2,:)) max(xt(2,:))])
zoom
subplot(323)
plot(timet,xt(3,:),'k-')  
ylabel('u_3(k)')
axis([min(timet) max(timet) min(xt(3,:)) max(xt(3,:))])
zoom
subplot(324)
plot(timet,xt(4,:),'k-')  
ylabel('u_4(k)')
axis([min(timet) max(timet) min(xt(4,:)) max(xt(4,:))])
zoom
subplot(325)
plot(timet,xt(5,:),'k-')  
ylabel('u_5(k)')
xlabel('Time, k')
axis([min(time) max(time) min(xt(5,:)) max(xt(5,:))])
zoom
subplot(326)
plot(timet,yt,'k-')
ylabel('y(k)')
xlabel('Time, k')
axis([min(timet) max(timet) min(yt) max(yt)])


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Perform correlation analysis to determine which inputs to use
% to the estimator
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% The i,j element of Rho is the correlation coefficient between the ith and jth
% variable (of course it is symetric with 1s all down the diagonal since every
% variable is perfectly correlated with itself).  Hence, the mth row shows the 
% how every variable is correlated with the mth variable.

Rho=corrcoef([yy xx])

% To determine the effect of each of the inputs of the estimator (i.e., the 
% elements of x(k)) on the output we need to look at the 1,j elements.  To
% do this we plot the first row as a histogram:

figure(3)
clf
bar(Rho(1,2:6),'w')
ylabel('Correlation coefficient values')
title('Correlation coefficient between inputs and y(k)')
xlabel('i indices for u_i(k), i=1,2,3,4,5')
grid
axis([0.5 5.5 -1 1])

% To view this use the following 3d plot (this helps to visualize the cross-correlations)

figure(4)
clf
surf(Rho);
colormap(jet)
% Use next line for generating plots to put in black and white documents.
%colormap(white);
xlabel('i');
ylabel('j');
zlabel('Correlation coefficient');
title('Matrix of correlation coefficient values');
rotate3d  % For rotating for better viewing

% From this you can see that the operator does not seem to use u4 and u5 too strongly 
% in deciding what set point to pick.  We may use this as a basis for deciding that 
% we should not use these variables as inputs to the controller. Also, note that
% the input u2 does not have that big of an impact - so if we want, we could not
% use it as an input either.


%%%%%%%%%%%%%%%%%%%%%%
% Controller estimator development:
%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Batch linear least squares training of an affine function 
% using all the inputs
%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Y=y(1,1:M)';

% Initialize Phi

phi(:,1)=[x(:,1)', 1]';
Phi=phi(:,1)';

% Form the rest of Phi

for i=2:M
	phi(:,i)=[x(:,i)', 1]';
	Phi=[Phi ; phi(:,i)'];
end

theta=Phi\Y

% Next, compute the approximator values and the error vector
% at the training data

for i=1:M,
	phi(:,i)=[x(:,i)', 1]';
	Fltrain(i)=theta'*phi(:,i);
	etrain(i)=y(i)-Fltrain(i);
end

% Print out the training error

training_error_bls=(1/M)*etrain*etrain'

% Next, compute the approximator values and the error vector
% at the testing data

for i=1:MGamma,
	phit(:,i)=[xt(:,i)', 1]';
	Fltest(i)=theta'*phit(:,i);
	etest(i)=yt(i)-Fltest(i);
end

% Print out the testing error

testing_error_bls=(1/MGamma)*etest*etest'

% Next, plot the data and the approximator to compare

figure(5)
clf
subplot(211)
plot(time,y,'k-',timet,Fltest,'k--')
ylabel('y and its estimate')
title('Estimator performance, y (solid line), estimate of y (dashed line)')
grid
subplot(212)
plot(timet,etest,'k-')
xlabel('Time, k')
title('Estimation error')
grid

%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Batch linear least squares training of an affine function 
% using only inputs u1, u2, u3
%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Y=y(1,1:M)';

x3=x(1:3,:); % Picks the first three inputs
x3t=xt(1:3,:);

% Initialize Phi

phi3(:,1)=[x3(:,1)', 1]';
Phi3=phi3(:,1)';

% Form the rest of Phi

for i=2:M
	phi3(:,i)=[x3(:,i)', 1]';
	Phi3=[Phi3 ; phi3(:,i)'];
end

theta3=Phi3\Y

% Next, compute the approximator values and the error vector
% at the training data

for i=1:M,
	phi3(:,i)=[x3(:,i)', 1]';
	Fltrain3(i)=theta3'*phi3(:,i);
	etrain3(i)=y(i)-Fltrain3(i);
end

% Print out the training error

training_error_bls_3inputs=(1/M)*etrain3*etrain3'

% Next, compute the approximator values and the error vector
% at the test data

for i=1:MGamma,
	phi3t(:,i)=[x3t(:,i)', 1]';
	Fltest3(i)=theta3'*phi3t(:,i);
	etest3(i)=yt(i)-Fltest3(i);
end

% Print out the testing error

testing_error_bls_3inputs=(1/MGamma)*etest3*etest3'

% Next, plot the data and the approximator to compare

figure(6)
clf
subplot(211)
plot(time,y,'k-',timet,Fltest3,'k--')
ylabel('y and its estimate')
title('Estimator performance, y (solid line), estimate of y (dashed line) (3 inputs)')
grid
subplot(212)
plot(timet,etest3,'k-')
xlabel('Time, k')
ylabel('Estimation errors')
grid

%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Batch linear least squares training of an affine function 
% using only inputs u1 and u3
%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Y=y(1,1:M)';

x2=x([1,3],:); % Picks the u1 and u3 inputs
x2t=xt([1,3],:);

% Initialize Phi

phi2(:,1)=[x2(:,1)', 1]';
Phi2=phi2(:,1)';

% Form the rest of Phi

for i=2:M
	phi2(:,i)=[x2(:,i)', 1]';
	Phi2=[Phi2 ; phi2(:,i)'];
end

theta2=Phi2\Y

% Next, compute the approximator values and the error vector
% at the training data

for i=1:M,
	phi2(:,i)=[x2(:,i)', 1]';
	Fltrain2(i)=theta2'*phi2(:,i);
	etrain2(i)=y(i)-Fltrain2(i);
end

% Print out the training error

training_error_bls_2inputs=(1/M)*etrain2*etrain2'

% Next, compute the approximator values and the error vector
% at the test data

for i=1:MGamma,
	phi2t(:,i)=[x2t(:,i)', 1]';
	Fltest2(i)=theta2'*phi2t(:,i);
	etest2(i)=yt(i)-Fltest2(i);
end

% Print out the training error

testing_error_bls_2inputs=(1/MGamma)*etest2*etest2'

% Next, plot the data and the approximator to compare

figure(7)
clf
subplot(211)
plot(time,y,'k-',timet,Fltest2,'k--')
ylabel('y and its estimate')
title('Estimator performance, y (solid line), estimate of y (dashed line) (2 inputs)')
grid
subplot(212)
plot(timet,etest2,'k-')
xlabel('Time, k')
ylabel('Estimation errors')
grid

%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Batch linear least squares training of an affine function 
% using only u1 as an input
%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Y=y(1,1:M)';

x1=x(1,:); % Picks the u1 input
x1t=xt(1,:);

% Initialize Phi

phi1(:,1)=[x1(:,1)', 1]';
Phi1=phi1(:,1)';

% Form the rest of Phi

for i=2:M
	phi1(:,i)=[x1(:,i)', 1]';
	Phi1=[Phi1 ; phi1(:,i)'];
end

theta1=Phi1\Y

% Next, compute the approximator values and the error vector
% at the training data

for i=1:M,
	phi1(:,i)=[x1(:,i)', 1]';
	Fltrain1(i)=theta1'*phi1(:,i);
	etrain1(i)=y(i)-Fltrain1(i);
end

% Print out the training error

training_error_bls_1input=(1/M)*etrain1*etrain1'

% Next, compute the approximator values and the error vector
% at the training data

for i=1:MGamma,
	phi1t(:,i)=[x1t(:,i)', 1]';
	Fltest1(i)=theta1'*phi1t(:,i);
	etest1(i)=yt(i)-Fltest1(i);
end

% Print out the training error

testing_error_bls_1input=(1/MGamma)*etest1*etest1'

% Next, plot the data and the approximator to compare

figure(8)
clf
subplot(211)
plot(time,y,'k-',timet,Fltest1,'k--')
ylabel('y and its estimate')
title('Estimator performance, y (solid line), estimate of y (dashed line) (1 input)')
grid
subplot(212)
plot(timet,etest1,'k-')
xlabel('Time, k')
ylabel('Estimation errors')
grid

%%%%%%%%%%%%%%%%%%%%%
% Train a fuzzy system - 1 input to premise (1, 3, or all for consequent functions)
%%%%%%%%%%%%%%%%%%%%%

% Next, find Phi, which involves processing x through phi

R=9;  % The number of rules.
n=1; % One input (to membership functions)

% Check ranges of data for u1and grid on that region

c1(1,:)=4.5:.25:6.5;
sigma1(1,:)=0.25*ones(1,R);

% Initialize Phif

	for j=1:R
		  mu1(j,1)=exp(-0.5*((x1(1,1)-c1(1,j))/sigma1(1,j))^2);
	end
	
	denominator1(1)=sum(mu1(:,1)); 
	
	for j=1:R
		xi1(j,1)=mu1(j,1)/denominator1(1);
	end

% To use only one input to the consequent functions:

%Phif1=[xi1(:,1)', x1(1,1)*xi1(:,1)'];

% To use all the inputs for the consequents, except u4 and u5:

Phif1=[xi1(:,1)', x(1,1)*xi1(:,1)', x(2,1)*xi1(:,1)', x(3,1)*xi1(:,1)'];

% To use all the inputs for the consequents:

%Phif1=[xi1(:,1)', x(1,1)*xi1(:,1)', x(2,1)*xi1(:,1)', x(3,1)*xi1(:,1)', x(4,1)*xi1(:,1)', x(5,1)*xi1(:,1)'];

% Form the rest of Phif

for i=2:M
	
	for j=1:R
		  mu1(j,i)=exp(-0.5*((x1(1,i)-c1(1,j))/sigma1(1,j))^2);
	end
	
	denominator1(i)=sum(mu1(:,i)); 
	
	for j=1:R
		xi1(j,i)=mu1(j,i)/denominator1(i);
	end

% To use only one input to the consequent functions: