模糊控制——WM演算法

在現代智慧控制演算法中，模糊控制是在實際控制系統設計中使用比較成熟的一種方法。模糊控制可以使用在一些無法建立系統模型的場合，根據專家經驗確定模糊規則，實現對系統的控制。

這裡主要介紹一種基本的模糊控制演算法，WM演算法。該演算法的思想是根據取樣的資料對（一組輸入、輸出資料），確定出模糊規則，通常是一條資料對就可以確定一條規則。

首先我們需要確定系統的輸入輸出數量，假設系統為單輸入單輸出。對輸入變數x，輸出變數y分別劃分模糊集合，可以使用正態分佈隸屬度函式，或者三角隸屬度函式來劃分。這叫做變數的模糊化。如x的論域為[0,2],劃分13個模糊集合，分別為A1,A2,A3,...,A13,如下圖：

  對於輸出y,論域為[-1.5,1.5],劃分13個模糊集合，為B1,B2,B3,...,B13，如下圖所示：

現在有0-2論域上均勻分佈的樣本點共21個，利用它們來確定模糊規則。

需要分別計算每一個數據點在模糊集合上的隸屬度，選取最高的隸屬度值作為確定一條模糊規則的依據。如樣本點（0.2,1），需要計算0.2在輸入隸屬度函式中的隸屬度值，需要計算13個值，找出其中最大的值如A5，則輸入為A5；再計算輸出1在13個模糊集合中的隸屬度函式值，找出最大的那個，如B2，則輸出為B2；由此可以確定一條模糊規則：IF x=A5 THEN y=B2;由此可以確定21個規則；但是這些規則有大量的重複和衝突的規則，需要計算它們置信度

conf=u(Ax)*u(By);

由此公式可以求出矛盾規則的置信度，把置信度低的規則去掉；按照WM演算法的提出則王立新的做法，還應該乘上一個專家經驗係數，也就是專家認為這條規則的可信度大不大。上面的公式改寫為：

conf=u*u(Ax)*u(By);

由此可以建立模糊規則庫；

上面表中的第一行代表輸入x的隸屬集合的下標，第二行代表輸出y的隸屬集合的下標。

利用模糊規則庫，計算輸出y;根據去模糊化公式：

即可計算輸出。

選用函式

y=0.9*sin(PI*x)+0.3*cos(3*PI*x);

以下是在matlab中的模擬程式碼：

（1）計算輸出變數y的隸屬度函式

function u = u_y_B(y,a,left,right,step)
%u_y_B
% 計算輸出變數y的隸屬度函式值
%y:輸出變數的值
%a：區間中點的值
%left:表示論域區間的左端點
%right：表示論域區間的右端點
%step:三角形底邊長的一半
%論域為[-1.5,1.5]，共有5個模糊區間：B1,B2,B3,B4,B5
b=a-step;
c=a+step;
len=length(y);
u=zeros(1,len);
for i=1:len
    if a==left+step
        if y(i)>=b&&y(i)<=a
        u(i)=1;
        end
        if y(i)>a&&y(i)<=c
         u(i)=(c-y(i))/(c-a);
        end
        if y(i)<b||y(i)>c
          u(i)=0;
        end
    elseif a==right-step
            if y(i)>=b&&y(i)<=a
            u(i)=(y(i)-b)/(a-b);
            end
            if y(i)>a&&y(i)<=c
                u(i)=1;
            end
            if y(i)<b||y(i)>c
                u(i)=0;
            end  
    else
            if y(i)>=b&&y(i)<=a
            u(i)=(y(i)-b)/(a-b);
            end
            if y(i)>a&&y(i)<=c
                u(i)=(c-y(i))/(c-a);
            end
            if y(i)<b||y(i)>c
                u(i)=0;
            end
    end
end
end

（2）計算輸入隸屬度（輸入模糊化）

function u=u_x_input(x)
%計算任意一個輸入在輸入模糊區間A1，A2，A3，A4，A5，A6，A7上各自的隸屬度值
% x：輸入值
% u：輸出一個有隸屬度組成的陣列
a=0;
u=zeros(1,10);
for i=1:10
    a=a+0.2;
   u(i)=u_x_A(x,a); 
end
end

function u= u_x_A(x,a)
%u_x_A
% 計算輸入變數x的隸屬度函式值
%x:輸入變數的值
%a：區間中點的值
%論域為[0,2]，共有10個模糊區間A1,A2,A3,A4,A5,A6,A7，A8，A9，A10
len=length(x);
u=zeros(1,len);
b=a-0.2;
c=a+0.2;
for i=1:len
    if a==0.2
        if x(i)>=b&&x(i)<=a
        u(i)=1;
        end
        if x(i)>a&&x(i)<=c
         u(i)=(c-x(i))/(c-a);
        end
        if x(i)<b||x(i)>c
          u(i)=0;
        end
    elseif a==1.8
            if x(i)>=b&&x(i)<=a
            u(i)=(x(i)-b)/(a-b);
            end
            if x(i)>a&&x(i)<=c
                u(i)=1;
            end
            if x(i)<b||x(i)>c
                u(i)=0;
            end  
    else
            if x(i)>=b&&x(i)<=a
            u(i)=(x(i)-b)/(a-b);
            end
            if x(i)>a&&x(i)<=c
                u(i)=(c-x(i))/(c-a);
            end
            if x(i)<b||x(i)>c
                u(i)=0;
            end
    end
end
end

clc;
clear;
PI=3.1415926;
t=0:0.01:2;
y=0.9*sin(PI*t)+0.3*cos(3*PI*t);
plot(t,y);%原始影象
xlabel('輸出值x');
ylabel('輸入值y');
grid on;
title('y=0.9*sin(PI*t)+0.3*cos(3*PI*t)');

%獲取取樣點,取樣21組資料
sample_x=0:0.1:2;
sample_y=0.9*sin(PI*sample_x)+0.3*cos(3*PI*sample_x);
sample_num=length(sample_x);%取樣個數

%論域x劃分set_X個模糊區間，使用正態（高斯）形隸屬函式,論域[0,2]
set_X=13;
xmin=0;
xmax=2;
x_step=(xmax-xmin)/(set_X-1);%x模糊集合的步長
av_x=xmin:x_step:xmax;       %計算高斯分佈均值
sigma_x=sqrt(-x_step^2/(8*log(0.5)));%計算高斯分佈方差
%sigma_x=0.09;
x=xmin:0.01:xmax;
figure(2)
for i=1:set_X
    plot(gaussmf(x,[sigma_x,av_x(i)]));%繪製x的模糊函式曲線
    hold on;
end
  legend('A1','A2','A3','A4','A5','A6','A7','A8','A9','A10','A11','A12','A13');
  xlabel('輸入值x');
  ylabel('隸屬度值u(x)');
  set(gca,'XTick',0:50:250);
  set(gca,'XTickLabel',{'0','0.5','1.0','1.5','2','2.5'});
  title('輸入變數x的模糊區間劃分以及隸屬度');
  
%論域y劃分set_Y模糊區間，使用三角隸屬函式,論域[-1.5,1.5]
figure(3)
set_Y=13;
ymin=-1.5;%論域下限
ymax=1.5; %論域上限
y_step=(ymax-ymin)/(set_Y+1);%三角形兩個尖點之間的步長
a=ymin;%儲存論域下限，方便後面的隸屬度計算
y=ymin:0.01:ymax;%獲取一組y的數值
for i=1:set_Y
    a=a+y_step;
    plot(u_y_B(y,a,ymin,ymax,y_step));%繪製y的模糊函式曲線
    hold on;
end
 legend('B1','B2','B3','B4','B5','B6','B7','B8','B9','B10','B11','B12','B13');
  xlabel('輸出值y');
  ylabel('隸屬度值u(y)');
  set(gca,'XTick',0:50:350);
  set(gca,'XTickLabel',{'-1.5','-1.0','-0.5','0','0.5','1.0','1.5','2.0'});
  title('輸出變數y的模糊區間劃分以及隸屬度');

  %WM演算法
 uxA=zeros(sample_num,set_X);%儲存每條樣本資料x的隸屬度函式值
 uyB=zeros(sample_num,set_Y);%儲存每條樣本資料y的隸屬度函式值
 for i=1:set_X
     uxA(:,i)=gaussmf(sample_x,[sigma_x,av_x(i)]);%sample_num個樣本x的在第i個模糊區間中的隸屬度值
 end
  a=ymin;
 for j=1:set_Y
     a=a+y_step;
     uyB(:,j)=u_y_B(sample_y,a,ymin,ymax,y_step);%sample_num個樣本y的在第j個模糊區間中的隸屬度值
 end
 WM_rule=zeros(3,sample_num);%儲存每個樣本資料所在的模糊集合的下標
 [~,WM_rule(1,:)]=max(uxA,[],2);%計算每個樣本x所在的模糊集合下標
 [~,WM_rule(2,:)]=max(uyB,[],2);%計算每個樣本y所在的模糊集合下標
 for i=1:sample_num
      WM_rule(3,i)=uxA(i,WM_rule(1,i))*uyB(i,WM_rule(2,i));  %計算每條規則的支援度 
 end
 %對於屬於x屬於同一個模糊區間，輸出卻不同的規則，去除信任度低的規則
 for i=2:sample_num
     if(WM_rule(1,i-1)==WM_rule(1,i))
         if(WM_rule(3,i-1)<=WM_rule(3,i))
             WM_rule(:,i-1)=0;
         else
             WM_rule(:,i)=0;
         end
     end
 end
 
 WM_rule(:,all(WM_rule==0,1))=[];%去除多於的規則（把每列全為0的刪除）
 %WM演算法的模糊規則庫已經建立完成

p_value=zeros(1,set_Y);%用於儲存y模糊函式尖點所對應的橫座標的值
a=ymin;
for i=1:set_Y
    a=a+y_step;
    p_value(i)=a;%儲存y模糊函式尖點所對應的橫座標的值
end

%測試規則
x=0;
y_x=zeros(1,201);
WM_y_x=zeros(1,201);
for i=1:201
    x=x+0.01;
    ux=zeros(1,set_X);
    for m=1:set_X
      ux(m)=gaussmf(x,[sigma_x,av_x(m)]);
    end
   num=0;
   num1=0;
   den=0;
   for j=1:set_X
       num=num+p_value(B_index(j))*ux(j);
       num1=num1+p_value(WM_rule(2,j))*ux(j);
       den=den+ux(j);
   end
   y_x(i)=num/den;
   WM_y_x(i)=num1/den;
end
figure(4);
x=0:0.01:2;
plot(x,WM_y_x,'-.b');%畫出WM演算法的輸出曲線
hold on;
y=0.9*sin(PI*x)+0.3*cos(3*PI*x);
plot(x,y,'-g');  %畫出原始函式的曲線
 xlabel('輸入值x');
 ylabel('輸出值y');
 title('使用WM、DM模糊控制演算法的結果');
 legend('WM演算法輸出曲線','原始輸出曲線');
grid on;

模擬結果如下：

由上圖可知：WM演算法能夠較好的擬合原始曲線。