1、T-S模糊系統結論參數的MATLAB仿真研究    摘  要  介紹了應用最小二乘法對T-S結論參數進行粗略辨識，確定參數的大致范圍，再應用遺傳算法對前提參數和結論參數同時優化的參數辨識方法。對非線性函數進行逼近實驗，給出了用MATLAB編程進行仿真的具體實現方法，結果證明該方法的可行性和有效性。     關鍵詞  最小二乘法；參數辨識；遺傳算法  0  引言    對T-S模糊系統參數辨識的過程大致分為結構辨識和參數辨識，而參數辨識則是整個系統

2、辨識的關鍵所在。遺傳算法以其在解空間內進行高效啟發式搜索，尋優速度快，不易陷入局部最優解等優點成為近來應用較多的優化方法。將遺傳算法用于解決T-S模糊模型的參數辨識問題，在應用最小二乘法進行粗略辨識的前提下，用遺傳算法對結論參數進行尋優，用MATLAB進行仿真，取得了較好的效果。 1  用最小二乘法對T-S模糊模型參數的初步辨識    T-S模糊模型辨識的過程一般分為以下幾個階段：前提結構辨識；前提參數辨識；結論結構辨識；結論參數辨識，直到模型滿足要求為止。結構辨識的方法，在此不再詳細說明，只對結論參數辨識問題展開討論。   

3、0; 在確定了前提結構和結論結構之后，對模糊模型的結論參數進行粗略的辨識，以確定遺傳算法尋優的范圍。在眾多的參數辨識方法中，最小二乘法是最基本的一種，Gauss于1795年就用最小二乘法，由觀測結果估算了行星的運行軌道。此后，這種方法被廣泛應用，并根據實際問題提出了許多改進的最小二乘法，如正交最小二乘法，廣義最小二乘法，增廣最小二乘法等。這里所用的是線性最小二乘法，將前提結構劃分的各個范圍中的輸入輸出數據擬合為一次多項式函數。從而得出粗略的結論參數。以此來大致確定遺傳算法要優化的結論參數范圍。前提參數的大致范圍可根據所選的隸屬函數來確定。     為簡單起見，考慮

4、一維的單輸入非線性系統。對下列函數進行逼近1： 設定輸入范圍為-1，1，將它模糊分割為五個區，隸屬度函數采用廣義的鐘形函數，這里只有一個輸入變量，輸出為y = a x + b的線性方程，待優化的結論參數有2×5 = 10個，用MATLAB編程來初步得出待優化的結論參數，主要代碼如下：     data_n=100;newdata_n=1001;x=linspace(-1，-0.6，data_n);     y=0.7*sin(pi*x)+0.3*sin(3*pi*x)+0.1*sin(5*pi*x);   

5、;  polyfit(x，y，1)     由上述方法仿真可得出的10個參數，確定參數范圍，如表1中所示：表1 用最小二乘法估計的參數及優化所選的參數范圍    參數 a1 b1 a2 b2 a3 b3 a4 b4 a5 b5 估計參數 -0.6499 -1.0459 0.9010 0.0518 8.3065 0.0000 0.9010 -0.0518 -0.6499 1.0459 參數范圍 -5，5 -5，5 -5，5 -5，5 -2，15 -5.5， 5.5 -5，5 -5， 5 -5.5， 5.5 -1，5

6、2 用遺傳算法來優化T-S模型的結論參數    由于同時優化的參數的數量較多，故采用實數編碼的方式對參數進行編碼。編碼過程是通過計算機產生所要優化的參數范圍內的隨機數，對每一個參數進行編碼后，連接在一起形成一條染色體，然后就可以對它進行遺傳操作。用MATLAB編程，確定尋優范圍的主要代碼如下： MinX(1)=-5.0*ones(1); MaxX(1)=5.0*ones(1); MinX(2)=-5.0*ones(1); MaxX(2)=5.5*ones(1); MinX(10)=-1*ones(1); MaxX(10)=5*ones(1); Kpar(：，1)=

7、MinX(1)+(MaxX(1)-MinX(1)*rand(Size，1); Kpar(：，2)=MinX(2)+(MaxX(2)-MinX(2)*rand(Size，1); Kpar(：，9)=MinX(9)+(MaxX(9)-MinX(9)*rand(Size，1); Kpar(：，10)=MinX(10)+(MaxX(10)-MinX(10)*rand (Size，1);     在產生大種群后，對個體進行初步篩選，去掉一些適應度差的個體，剩下的作為初始種群進行遺傳操作，這種方法可以使種群在保持多樣性的同時，節省計算時間。進化代數的確定，一般是根據問題所要求

8、的精度來確定，精度低，進化的代數就可以少一些，反之則多一些。     適應度函數是衡量個體優劣的指標，為了達到尋優的目標，適應度函數一般是通過目標函數變換而來的，這里對T-S模型的參數進行辨識，采用的目標函數為均方誤差：        （1）     其中：t(k)為由辨識的模糊模型計算出的第k個采樣時刻的輸出值，(k)為第k 個采樣時刻的實際輸出值，Q為總的采樣次數。     在進行選擇操作時，是按照適應度越大，被選擇的概率越大，所以

9、，這里選用的適應度函數為：           （2）     其中是一個較小的正實數，目的是為了避免除數為零的情況發生3。其主要代碼如下： for i=1：1：Size  %以下為初始篩選         cansu=Kpar(i，：);         code=cansu;    &

10、#160;    rmse=computeitae(code);%計算均方誤差;         BsJ=rmse;         BsJi(i)=BsJ;%將均誤差賦給BsJi     end         fi=1./BsJi;% 求適應度的值     Oderfi，Indexfi=sort

11、(fi)%對適應度值由小到大排列        sortrsa=Kpar(Indexfi(Size)，：);     for s=2：1：selectsize          shuijirsa=Kpar(Indexfi(Size-s+1)，：);         sortrsa=sortrsa  shuijirsa;   &

12、#160;      end          sortrsa          Kpar=sortrsa          Size=selectsize     遺傳操作一般包括選擇，交叉和變異。選擇方法采用蒙特卡羅法，按比例的適應度分配。若某個個體i ，其適應度為f i ，

13、 則其被選擇的概率表示為：                （3）     由于在進行遺傳操作前，已經對個體進行了初步的篩選，所以為了避免對種群中優良個體的破壞，這里采用單點交叉的方法，隨機選擇交叉點之后，將兩個個體的交叉點后面的基因進行交換。變異采用實值變異的方法，隨機選擇染色體的某個基因，由于用實數編碼的方法，每個基因就是一個要辨識的參數值，所以可以用一個函數實現在參數范圍內適當改變該參數值的大小，從而達到保持種群

14、多樣性的目的。為了使尋優不過早的收斂到次優解，隨著進化代數的增加，需要適當增大變異率，其實現方法只需用一個函數來表示變異率： pm=0.1 + 1：1：G×0.1/G       （4）     式中：pm表示變異率，G代表進化代數。1：1：G表示一個數組，變化的范圍為1，G，步長為1。為了避免破壞優良個體，變異率不宜取的過大。通過實驗得知，在第一代時，可取變異率 pm=0.1+1×0.1/G，第二代時，變異率pm=0.1+2×0.1/G，以此類推。這樣隨著進化代數的增加，p

15、m也隨著增加。     另外，為了防止遺傳操作對最優個體的破壞，采取保留最優個體的方法。將每一代產生的最優個體放在該種群的最后，再繼續進行下一代的操作。整個算法的流程圖如圖所示。具體實現的遺傳操作的主要代碼如下： G=300;%進化的代數 BsJ=0; for kg=1：1：G   time(kg)=kg;       %*step 1：計算誤差*        for i=1：1：Size     

16、60;   cansu1=Kpar(i，：);         selectcode=cansu1;         error=computeitae(selectcode);%計算均方誤差;         BsJ1=error;         BsJi1(i)=BsJ1;  

17、60;  end        OderJi，IndexJi=sort(BsJi1);    BestJ(kg)=OderJi(1);    BJ=BestJ(kg);     Ji=BsJi1+1e-10;%避免除零        fi=1./Ji;% 求適應度的值     Oderfi，Indexfi=sort(fi)%對適應度值由小到大排列  

18、0;  Bestfi=Oderfi(Size);     BestS=Kpar(Indexfi(Size)，：)%保存最大適應度值對應的染色體(參數)         %*step2：選擇和復制* 圖 算法流程圖 fi_sum=sum(fi);     fi_Size=(Oderfi/fi_sum)*Size;     fi_S=floor(fi_Size);     r=Size-sum(fi_

19、S);     Rest=fi_Size-fi_S;     RestValue，Index=sort(Rest);        for i=Size：-1：Size-r+1         fi_S(Index(i)=fi_S(Index(i)+1;     end        k=1;   

20、  for i=Size：-1：1         for j=1：1：fi_S(i)             TempE(k，：)=Kpar(Indexfi(i)，：);             k=k+1;         e

21、nd     end         % step 3交叉率不能太大        Pc=0.7;     for i=1：2：(Size-1)         temp=rand;         if Pc>temp     &#

22、160;       alfa=rand; TempE(i，：)=alfa*Kpar(i+1，：)+(1-alfa)*Kpar(i，：); TempE(i+1，：)=alfa*Kpar(i，：)+(1-alfa)*Kpar(i+1，：);         end     end     TempE(Size，：)=BestS;     Kpar=TempE; 

23、0;      %*Step 4：變異*變異率不能太大     Pm=0.15+1：1：Size*(0.01)/Size;     Pm_rand=rand(Size，CodeL);     Mean=(MaxX+MinX)/2;     Dif=(MaxX-MinX)/2;        for i=1：1：Size      

24、;   for j=1：1：CodeL             if Pm(i)>Pm_rand(i，j)                TempE(i，j)=Mean(j)+Dif(j)*(rand- 0.4);             End         end     end     TempE(Size，：)=BestS;%保留最優個體     Kpar=TempE;     end BestS Best_J=BestJ(G) figure(1); plot(time，BestJ); xlabel('T