% Generate a fuzzy model for the pressure system

clear FM

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% define constants
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
FM.c = 3;          	   % number of clusters
FM.seed = 0;            % seed
FM.m = 2;               % fuzziness
FM.ante = 2;            % antecedent:   1 - product-space MFS
                        %               2 - projected MFS
FM.Ny = {{[1]};{[1]}};              % denominator order
FM.Nu = {{[1]};{[1]}};              % numerator orders
%FM.Ny = {{[1]};{[1]}};              % denominator order
%FM.Nu = {{[1]};{[1]}};              % numerator orders
%FM.Nd = 1;              % transport delays
                        % (set to 1 for y(k+1) = f(u(k),....))

FM.Ts = 5;              % sample time [s]

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% load identification data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% DATA(:,1) - time                 DATA(:,2) - output pressure
% DATA(:,3) - measured gas flow    DATA(:,4) - disturbance
% DATA(:,5) - valve position       DATA(:,6) - setpoint for gas flow
% DATA(:,7) - n.a.

load trsart3; Dat.U{1} = DATA(:,5); Dat.Y{1} = DATA(:,2);
load trsart4; Dat.U{2} = DATA(:,5); Dat.Y{2} = DATA(:,2);

Dat.InputName{1}  = 'valve position';
Dat.OutputName{1} = 'pressure';

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% load validation data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
load tssart2; ue = DATA(:,5); ye = DATA(:,2);
load tssart3; ue = [ue;DATA(:,5)]; ye = [ye;DATA(:,2)];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% make fuzzy model by means of fuzzy clustering
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[FM,Part] = fmclust(Dat,FM);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% simulate the fuzzy model for validation data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(1); clf
[ym,VAF,dof,yl,ylm] = fmsim(ue,ye,FM,[],[],2); VAF
title('Process output (blue) and model output (magenta)');
ylabel('Output');
disp('Hit any key to plot the local models ...'); pause

figure(1); clf
subplot(211); plot([ylm{1}]); set(gca,'ylim',[1 2])
title('Individual local models');
xlabel('Time'); ylabel('Output');
subplot(212); plot(dof{1})
title('Degrees of fulfillment');
xlabel('Time'); ylabel('Membership grade');

disp('Hit any key to plot the membership functions ...'); pause
plotmfs(FM)

if exist('fuzblock','file') == 4,
    disp('Hit any key to open a Simulink model ...'); pause
    t = 5*(0:length(u)-1)';
    FIS = fm2fis(FM);
    %pressure
end;