% 本程序用來實現CSTR模型的模糊PID控制仿真
% 
% 作者：

clear

% 創建一個模糊推理系統
a = newfis('fuzzcontrol');

% 添加輸入變量e和隸屬度函數
a = addvar(a,'input','e',[-3,3]);

a = addmf(a,'input',1,'NB','zmf',[-3,-1]);
a = addmf(a,'input',1,'NM','trimf',[-3,-2,0]);
a = addmf(a,'input',1,'NS','trimf',[-3,-1,1]);
a = addmf(a,'input',1,'Z','trimf',[-2,0,2]);
a = addmf(a,'input',1,'PS','trimf',[-1,1,3]);
a = addmf(a,'input',1,'PM','trimf',[0,2,3]);
a = addmf(a,'input',1,'PB','smf',[1,3]);

% figure(1);
% plotmf(a,'input',1);

% 添加輸入變量ec和隸屬度函數
a = addvar(a,'input','ec',[-3,3]);

a = addmf(a,'input',2,'NB','zmf',[-3,-1]);
a = addmf(a,'input',2,'NM','trimf',[-3,-2,0]);
a = addmf(a,'input',2,'NS','trimf',[-3,-1,1]);
a = addmf(a,'input',2,'Z','trimf',[-2,0,2]);
a = addmf(a,'input',2,'PS','trimf',[-1,1,3]);
a = addmf(a,'input',2,'PM','trimf',[0,2,3]);
a = addmf(a,'input',2,'PB','smf',[1,3]);

% figure(2);
% plotmf(a,'input',2);

% 添加輸出變量U和隸屬度函數

a = addvar(a,'output','U',[-3,3]);

a = addmf(a,'output',1,'NB','trimf',[-4,-3,-2]);
a = addmf(a,'output',1,'NM','trimf',[-3,-2,-1]);
a = addmf(a,'output',1,'NS','trimf',[-2,-1,0]);
a = addmf(a,'output',1,'Z','trimf',[-1,0,1]);
a = addmf(a,'output',1,'PS','trimf',[0,1,2]);
a = addmf(a,'output',1,'PM','trimf',[1,2,3]);
a = addmf(a,'output',1,'PB','trimf',[2,3,4]); 

% figure(3);
% plotmf(a,'output',1);

% 推理規則表
rulelist = [ 
     1     1     1     1     1;
     1     2     1     1     1;
     1     3     2     1     1;
     1     4     2     1     1;
     1     5     3     1     1;
     1     6     3     1     1;
     1     7     4     1     1;
     2     1     1     1     1;
     2     2     2     1     1;
     2     3     2     1     1;
     2     4     3     1     1;
     2     5     3     1     1;
     2     6     4     1     1;
     2     7     5     1     1;
     3     1     2     1     1;
     3     2     2     1     1;
     3     3     3     1     1;
     3     4     3     1     1;
     3     5     4     1     1;
     3     6     5     1     1;
     3     7     5     1     1;
     4     1     2     1     1;
     4     2     3     1     1;
     4     3     3     1     1;
     4     4     4     1     1;
     4     5     5     1     1;
     4     6     5     1     1;
     4     7     6     1     1;
     5     1     3     1     1;
     5     2     3     1     1;
     5     3     4     1     1;
     5     4     5     1     1;
     5     5     5     1     1;
     5     6     6     1     1;
     5     7     6     1     1;
     6     1     3     1     1;
     6     2     4     1     1;
     6     3     5     1     1;
     6     4     5     1     1;
     6     5     6     1     1;
     6     6     6     1     1;
     6     7     7     1     1;
     7     1     4     1     1;
     7     2     5     1     1;
     7     3     5     1     1;
     7     4     6     1     1;
     7     5     6     1     1;
     7     6     7     1     1;
     7     7     7     1     1;
];

a =  addrule(a,rulelist);
 
%  解模糊方法
a = setfis(a,'DefuzzMethod','centroid');

% 模糊推理系統創建完畢
 
% *****************************************************************

% 仿真參數
N = 200;
h = 0.1;

% 控制變量初值
Q(1) = 100.0;
Qc(1) = 15.0;

% 狀態變量初值
Ca(1) = 0.037;
T(1) = 402.35;
Tc(1) = 345.44;
H(1) = 6.0; % 單位是分米

% 被控變量設定值
Cas = 0.06;
Ts = 405.0;
Tcs = 350.0;

% PID調節參數、模糊控制比例因子
% 控制Ca
% Kp = 20.0;  
% Ti = 0.05;
% Td = 0.0; 
% Ke = 120; 
% Kc = 10;  
% Ku = 0.1; 

% 控制T
% Kp = 1.5;  
% Ti = 1.0;
% Td = 0.0;  
% Ke = 0.3; 
% Kc = 0.1;  
% Ku = 0.2; 

% 控制Tc
Kp = 0.5;  
Ti = 0.2;
Td = 0.0;
Ke = 0.5; 
Kc = 0.05;  
Ku = 0.1; 

q0 = Kp*(1+h/Ti+Td/h);
q1 = -1*Kp*(1+2*Td/h);
q2 = Kp*Td/h;

% 誤差變量
ek = 0;
ek_1 = 0;
ek_2 = 0;
  
% 控制仿真
for i = 1:N 
%     % 控制生成物濃度Ca
%     ek_2 = ek_1;
%     ek_1 = ek;
%     ek = Cas-Ca(i);
%     e(i) = ek;
%     ec(i) = (ek-ek_1)/h;
%     
%     if abs(ek)>0.005
%         dQ =(evalfis([ek*Ke,(ek-ek_1)/h*Kc],a))*Ku;
%         Q(i+1) = Q(i)+dQ;
%     else 
%         dQ = q0*ek+q1*ek_1+q2*ek_2;
%         Q(i+1) = Q(i)+dQ;
%     end 
%     
%     [Ca(i+1), T(i+1), Tc(i+1), H(i+1)] = CSTR_RK4(Q(i+1), Qc, Ca(i), T(i), Tc(i), H(i));
    
%     % 控制反應溫度 T
%     ek_2 = ek_1;
%     ek_1 = ek;
%     ek = Ts-T(i);
%     e(i) = ek;
%     ec(i) = (ek-ek_1)/h;
%     
%     if abs(ek)>0.5
%         dQc =(evalfis([ek*Ke,(ek-ek_1)/h*Kc],a))*Ku;
%         Qc(i+1) = Qc(i)-dQc;
%     else 
%         dQc = q0*ek+q1*ek_1+q2*ek_2;
%         Qc(i+1) = Qc(i)-dQc;
%     end 
%     
%     [Ca(i+1), T(i+1), Tc(i+1), H(i+1)] = CSTR_RK4(Q, Qc(i+1), Ca(i), T(i), Tc(i), H(i));

% 控制冷卻液溫度 Tc
    
    ek_2 = ek_1;
    ek_1 = ek;
    ek = Tcs-Tc(i);
    e(i) = ek;
    ec(i) = (ek-ek_1)/h;
    
    if abs(ek)>0.5
        dQc =(evalfis([ek*Ke,(ek-ek_1)/h*Kc],a))*Ku;
        Qc(i+1) = Qc(i)-dQc;
    else 
        dQc = q0*ek+q1*ek_1+q2*ek_2;
        Qc(i+1) = Qc(i)-dQc;
    end 
    
    [Ca(i+1), T(i+1), Tc(i+1), H(i+1)] = CSTR_RK4(Q, Qc(i+1), Ca(i), T(i), Tc(i), H(i));
end
% **********************************************

% figure(1)
% plot(Ca,'b.-')
% grid on

% figure(2)
% plot(T,'b.-')
% grid on

figure(3)
plot(Tc,'b.-')
grid on

% figure(2)
% plot(e)
% grid on
% 
% figure(4)
% plot(ec)
% grid on

% figure(4)
% plot(Qc)
% grid on

% figure(5)
% plot(Q)
% grid on