operator_est.m

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

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

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

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

% To use all inputs to the consequents:

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

end

% Find the least squares estimate

thetaf1=Phif1\Y


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

for i=1:M
	for j=1:R
		  mu1tr(j,i)=exp(-0.5*((x1(1,i)-c1(1,j))/sigma1(1,j))^2);
	end
	
	denominator1tr(i)=sum(mu1tr(:,i)); 
	
	for j=1:R
		xi1tr(j,i)=mu1tr(j,i)/denominator1tr(i);
	end
	
% To use only one input to the consequent functions:

%phif1tr(:,i)=[xi1tr(:,i)', x1(1,i)*xi1tr(:,i)']';

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

	phif1tr(:,i)=[xi1tr(:,i)', x(1,i)*xi1tr(:,i)', x(2,i)*xi1tr(:,i)', x(3,i)*xi1tr(:,i)']';

% To use all inputs to the consequents:

%	phif1tr(:,i)=[xi1tr(:,i)', x(1,i)*xi1tr(:,i)', x(2,i)*xi1tr(:,i)', x(3,i)*xi1tr(:,i)', x(4,i)*xi1tr(:,i)', x(5,i)*xi1tr(:,i)']';

	Fftrain1(i)=thetaf1'*phif1tr(:,i);
	etrainf1(i)=y(i)-Fftrain1(i);
end

% Print out the training error

training_error_TS_1input=(1/M)*etrainf1*etrainf1'

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

for i=1:MGamma
	for j=1:R
		  mu1t(j,i)=exp(-0.5*((x1t(1,i)-c1(1,j))/sigma1(1,j))^2);
	end
	
	denominator1t(i)=sum(mu1t(:,i)); 
	
	for j=1:R
		xi1t(j,i)=mu1t(j,i)/denominator1t(i);
	end
	
% To use only one input to the consequent functions:

%phif1t(:,i)=[xi1t(:,i)', x1t(1,i)*xi1t(:,i)']';

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

	phif1t(:,i)=[xi1t(:,i)', xt(1,i)*xi1t(:,i)', xt(2,i)*xi1t(:,i)', xt(3,i)*xi1t(:,i)']';

% To use all inputs to the consequents:

%	phif1t(:,i)=[xi1t(:,i)', xt(1,i)*xi1t(:,i)', xt(2,i)*xi1t(:,i)', xt(3,i)*xi1t(:,i)', xt(4,i)*xi1t(:,i)', xt(5,i)*xi1t(:,i)']';

	Fftest1(i)=thetaf1'*phif1t(:,i);
	etestf1(i)=yt(i)-Fftest1(i);
end

% Print out the testing error

testing_error_TS_1input=(1/MGamma)*etestf1*etestf1'

% Next, plot the data and the approximator to compare

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

% Next, test to see if there is correlation between any input and the error
% and if there is then that input should be added as an input to the system.

Rhof1=corrcoef([etrainf1' 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(10)
clf
bar(Rhof1(1,2:6),'w')
ylabel('Correlation coefficient values')
title('Correlation coefficient between input and error (for training data)')
xlabel('i indices for u_i(k), i=1,2,3,4,5')
grid
axis([0.5 5.5 -1 1])


%%%%%%%%%%%%%%%%%%%%%
% Train a fuzzy system - 2 inputs (2 or all for consequent functions)
%%%%%%%%%%%%%%%%%%%%%

% In this case we use two inputs (u1, u3) to the membership functions, and
% to the consequent functions (or all the inputs to the consequent functions).

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

R=9;  % The number of rules.
n=2; % Two inputs (to membership functions)

% Check ranges of data for u1 and u3 and grid on that region

c(1,:)=[4.5,4.5,4.5,5.5,5.5,5.5,6.5,6.5,6.5];
c(2,:)=[1000,3500,6000,1000,3500,6000,1000,3500,6000];
sigma(1,:)=1*ones(1,R);
sigma(2,:)=2500*ones(1,R);

% Initialize Phif

	for j=1:R
		mu(j,1)=1; % Initialize
		for ii=1:n
		  mu(j,1)=mu(j,1)*exp(-0.5*((x2(ii,1)-c(ii,j))/sigma(ii,j))^2);
		end
	end
	
	denominator(1)=sum(mu(:,1)); 
	
	for j=1:R
		xi(j,1)=mu(j,1)/denominator(1);
	end

% To use only two inputs to the consequent functions:

Phif=[xi(:,1)', x2(1,1)*xi(:,1)', x2(2,1)*xi(:,1)'];

% To use all the inputs for the consequents:

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

% Form the rest of Phif

for i=2:M
	
	for j=1:R
		mu(j,i)=1; % Initialize
		for ii=1:n
		  mu(j,i)=mu(j,i)*exp(-0.5*((x2(ii,i)-c(ii,j))/sigma(ii,j))^2);
		end
	end
	
	denominator(i)=sum(mu(:,i)); 
	
	for j=1:R
		xi(j,i)=mu(j,i)/denominator(i);
	end

% To use only two inputs to the consequent functions:

Phif=[Phif; xi(:,i)', x2(1,i)*xi(:,i)', x2(2,i)*xi(:,i)'];

% To use all inputs to the consequents:

%	Phif=[Phif; xi(:,i)', x(1,i)*xi(:,i)', x(2,i)*xi(:,i)', x(3,i)*xi(:,i)', x(4,i)*xi(:,i)', x(5,i)*xi(:,i)'];

end

% Find the least squares estimate

thetaf=Phif\Y


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

for i=1:M
	for j=1:R
		mu(j,i)=1; % Initialize
		for ii=1:n
		  mu(j,i)=mu(j,i)*exp(-0.5*((x2(ii,i)-c(ii,j))/sigma(ii,j))^2);
		end
	end
	
	denominator(i)=sum(mu(:,i)); 
	
	for j=1:R
		xi(j,i)=mu(j,i)/denominator(i);
	end
	
% To use only two inputs to the consequent functions:

phif(:,i)=[xi(:,i)', x2(1,i)*xi(:,i)', x2(2,i)*xi(:,i)']';

% To use all inputs to the consequents:

%	phif(:,i)=[xi(:,i)', x(1,i)*xi(:,i)', x(2,i)*xi(:,i)', x(3,i)*xi(:,i)', x(4,i)*xi(:,i)', x(5,i)*xi(:,i)']';

	Fftrain(i)=thetaf'*phif(:,i);
	etrainf(i)=y(i)-Fftrain(i);
end

% Print out the training error

training_error_TS_2inputs=(1/M)*etrainf*etrainf'

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

for i=1:MGamma
	for j=1:R
		mut(j,i)=1; % Initialize
		for ii=1:n
		  mut(j,i)=mut(j,i)*exp(-0.5*((x2t(ii,i)-c(ii,j))/sigma(ii,j))^2);
		end
	end
	
	denominatort(i)=sum(mut(:,i)); 
	
	for j=1:R
		xit(j,i)=mut(j,i)/denominatort(i);
	end
	
% To use only two inputs to the consequent functions:

phift(:,i)=[xit(:,i)', x2t(1,i)*xit(:,i)', x2t(2,i)*xit(:,i)']';

% To use all inputs to the consequents:

%	phift(:,i)=[xit(:,i)', xt(1,i)*xit(:,i)', xt(2,i)*xit(:,i)', xt(3,i)*xit(:,i)', xt(4,i)*xit(:,i)', xt(5,i)*xit(:,i)']';

	Fftest(i)=thetaf'*phift(:,i);
	etestf(i)=yt(i)-Fftest(i);
end

% Print out the testing error

testing_error_TS_2inputs=(1/MGamma)*etestf*etestf'


% Next, plot the data and the approximator to compare

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

% Next, test to see if there is correlation between any input and the error
% and if there is then that input should be added as an input to the system.

Rhof2=corrcoef([etrainf' 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(12)
clf
bar(Rhof2(1,2:6),'w')
ylabel('Correlation coefficient values')
title('Correlation coefficient between input and error (for training data)')
xlabel('i indices for u_i(k), i=1,2,3,4,5')
grid
axis([0.5 5.5 -1 1])


%%%%%%%%%%%%%%%%%%%%%
% Train a fuzzy system - 2 inputs - but u1 and u4 (2 or all for consequent functions)
%%%%%%%%%%%%%%%%%%%%%

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

R=9;  % The number of rules.
n=2; % Two inputs (to membership functions)

% Check ranges of data for u1 and u4 and grid on that region

c(1,:)=[4.5,4.5,4.5,5.5,5.5,5.5,6.5,6.5,6.5];
c(2,:)=[-.3,-.1,.1,-.3,-.1,.1,-.3,-.1,.1];
sigma(1,:)=1*ones(1,R);
sigma(2,:)=0.2*ones(1,R);

x4=x([1,4],:); % Picks the u1 and u4 inputs
x4t=xt([1,4],:);

% Initialize Phif

	for j=1:R
		mu(j,1)=1; % Initialize
		for ii=1:n
		  mu(j,1)=mu(j,1)*exp(-0.5*((x4(ii,1)-c(ii,j))/sigma(ii,j))^2);
		end
	end
	
	denominator(1)=sum(mu(:,1)); 
	
	for j=1:R
		xi(j,1)=mu(j,1)/denominator(1);
	end

% To use only two inputs to the consequent functions:

Phif4=[xi(:,1)', x4(1,1)*xi(:,1)', x4(2,1)*xi(:,1)'];

% To use all the inputs for the consequents:

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

% Form the rest of Phif

for i=2:M
	
	for j=1:R
		mu(j,i)=1; % Initialize
		for ii=1:n
		  mu(j,i)=mu(j,i)*exp(-0.5*((x4(ii,i)-c(ii,j))/sigma(ii,j))^2);
		end
	end
	
	denominator(i)=sum(mu(:,i)); 
	
	for j=1:R
		xi(j,i)=mu(j,i)/denominator(i);
	end

% To use only two inputs to the consequent functions:

Phif4=[Phif4; xi(:,i)', x4(1,i)*xi(:,i)', x4(2,i)*xi(:,i)'];

% To use all inputs to the consequents:

%	Phif4=[Phif4; xi(:,i)', x(1,i)*xi(:,i)', x(2,i)*xi(:,i)', x(3,i)*xi(:,i)', x(4,i)*xi(:,i)', x(5,i)*xi(:,i)'];

end

% Find the least squares estimate

thetaf4=Phif4\Y


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

for i=1:M
	for j=1:R
		mu(j,i)=1; % Initialize
		for ii=1:n
		  mu(j,i)=mu(j,i)*exp(-0.5*((x4(ii,i)-c(ii,j))/sigma(ii,j))^2);
		end
	end
	
	denominator(i)=sum(mu(:,i)); 
	
	for j=1:R
		xi(j,i)=mu(j,i)/denominator(i);
	end
	
% To use only two inputs to the consequent functions:

phif4(:,i)=[xi(:,i)', x4(1,i)*xi(:,i)', x4(2,i)*xi(:,i)']';

% To use all inputs to the consequents:

%	phif4(:,i)=[xi(:,i)', x(1,i)*xi(:,i)', x(2,i)*xi(:,i)', x(3,i)*xi(:,i)', x(4,i)*xi(:,i)', x(5,i)*xi(:,i)']';

	Fftrain4(i)=thetaf4'*phif4(:,i);
	etrainf4(i)=y(i)-Fftrain4(i);
end

% Print out the training error

training_error_TS_2inputs4=(1/M)*etrainf4*etrainf4'

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

for i=1:MGamma
	for j=1:R
		mut(j,i)=1; % Initialize
		for ii=1:n
		  mut(j,i)=mut(j,i)*exp(-0.5*((x4t(ii,i)-c(ii,j))/sigma(ii,j))^2);
		end
	end
	
	denominatort(i)=sum(mut(:,i)); 
	
	for j=1:R
		xit(j,i)=mut(j,i)/denominatort(i);
	end
	
% To use only two inputs to the consequent functions:

phif4t(:,i)=[xit(:,i)', x4t(1,i)*xit(:,i)', x4t(2,i)*xit(:,i)']';

% To use all inputs to the consequents:

%	phif4t(:,i)=[xit(:,i)', xt(1,i)*xit(:,i)', xt(2,i)*xit(:,i)', xt(3,i)*xit(:,i)', xt(4,i)*xit(:,i)', xt(5,i)*xit(:,i)']';

	Fftest4(i)=thetaf4'*phif4t(:,i);
	etestf4(i)=yt(i)-Fftest4(i);
end

% Print out the testing error

testing_error_TS_2inputs4=(1/MGamma)*etestf4*etestf4'


% Next, plot the data and the approximator to compare

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

% Next, test to see if there is correlation between any input and the error
% and if there is then that input should be added as an input to the system.

Rhof24=corrcoef([etrainf4' 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(14)
clf
bar(Rhof24(1,2:6),'w')
ylabel('Correlation coefficient values')
title('Correlation coefficient between input and error (for training data)')
xlabel('i indices for u_i(k), i=1,2,3,4,5')
grid
axis([0.5 5.5 -1 1])



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% End of Program
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%