用遗传算法优化模糊控制规则，有部分代码，不懂计算适应度那块调用子函数，子函数要怎么写？以及M程序如何和simulink相互调用？

clear all;

close all;

b=newfis('yichuan');

b=addvar(b,'input','e',[-3,3]); %Parametere

b=addmf(b,'input',1,'NB','zmf',[-3,-1]);

b=addmf(b,'input',1,'NM','trimf',[-3,-2,0]);

b=addmf(b,'input',1,'NS','trimf',[-3,-1,1]);

b=addmf(b,'input',1,'Z','trimf',[-2,0,2]);

b=addmf(b,'input',1,'PS','trimf',[-1,1,3]);

b=addmf(b,'input',1,'PM','trimf',[0,2,3]);

b=addmf(b,'input',1,'PB','smf',[1,3]);

b=addvar(b,'input','ec',[-3,3]); %Parameterec

b=addmf(b,'input',2,'NB','zmf',[-3,-1]);

b=addmf(b,'input',2,'NM','trimf',[-3,-2,0]);

b=addmf(b,'input',2,'NS','trimf',[-3,-1,1]);

b=addmf(b,'input',2,'Z','trimf',[-2,0,2]);

b=addmf(b,'input',2,'PS','trimf',[-1,1,3]);

b=addmf(b,'input',2,'PM','trimf',[0,2,3]);

b=addmf(b,'input',2,'PB','smf',[1,3]);

b=addvar(b,'output','kp',[-0.3,0.3]);%Parameter kp

b=addmf(b,'output',1,'NB','zmf',[-0.3,-0.1]);

b=addmf(b,'output',1,'NM','trimf',[-0.3,-0.2,0]);

b=addmf(b,'output',1,'NS','trimf',[-0.3,-0.1,0.1]);

b=addmf(b,'output',1,'Z','trimf',[-0.2,0,0.2]);

b=addmf(b,'output',1,'PS','trimf',[-0.1,0.1,0.3]);

b=addmf(b,'output',1,'PM','trimf',[0,0.2,0.3]);

b=addmf(b,'output',1,'PB','smf',[0.1,0.3]);

b=addvar(b,'output','ki',[-0.06,0.06]);%Parameter ki

b=addmf(b,'output',2,'NB','zmf',[-0.06,-0.02]);

b=addmf(b,'output',2,'NM','trimf',[-0.06,-0.04,0]);

b=addmf(b,'output',2,'NS','trimf',[-0.06,-0.02,0.02]);

b=addmf(b,'output',2,'Z','trimf',[-0.04,0,0.04]);

b=addmf(b,'output',2,'PS','trimf',[-0.02,0.02,0.06]);

b=addmf(b,'output',2,'PM','trimf',[0,0.04,0.06]);

b=addmf(b,'output',2,'PB','smf',[0.02,0.06]);

b=addvar(b,'output','kd',[-3,3]);%Parameter kd

b=addmf(b,'output',3,'NB','zmf',[-3,-1]);

b=addmf(b,'output',3,'NM','trimf',[-3,-2,0]);

b=addmf(b,'output',3,'NS','trimf',[-3,-1,1]);

b=addmf(b,'output',3,'Z','trimf',[-2,0,2]);

b=addmf(b,'output',3,'PS','trimf',[-1,1,3]);

b=addmf(b,'output',3,'PM','trimf',[0,2,3]);

b=addmf(b,'output',3,'PB','smf',[1,3]);

b=setfis(b,'DefuzzMethod','centroid');

writefis(b,'yichuan');

b=readfis('yichuan');

shuru0=[1 1;1 2;1 3;1 4;1 5;1 6;1 7;

21;2 2;2 3;2 4;2 5;2 6;2 7;

31;3 2;3 3;3 4;3 5;3 6;3 7;

4 1;4 2;4 3;44;4 5;4 6;4 7;

5 1;5 2;5 3;54;5 5;5 6;5 7;

6 1;6 2;6 3;64;6 5;6 6;6 7;

7 1;7 2;7 3;74;7 5;7 6;7 7];%输入规则

lianjie0=[1 1;1 1;1 1;1 1;1 1;1 1;1 1;

11;1 1;1 1;1 1;1 1;1 1;1 1;

11;1 1;1 1;1 1;1 1;1 1;1 1;

1 1;1 1;1 1;1 1;11;1 1;1 1;

1 1;1 1;1 1;11;1 1;1 1;1 1;

1 1;1 1;1 1;11;1 1;1 1;1 1;

1 1;1 1;1 1;11;1 1;1 1;1 1];%模糊规则中的权重和连接信息

shuchu0=ones(49,3);%生成输出规则，其中Outi是1*128的矩阵，重新生成16*8的输出矩阵

rulelist1=[shuru0 shuchu0 lianjie0];%将个体所携带的规则信息导入到模糊控制规则

b=readfis('yichuan');

b=addrule(a,rulelist1);

writefis(b,'yichuan');

%%%%%%%%%  initial code   %%%%%%%%%%

Size=30;%种群个体数

CodeL=128;%基因长度

MinX=zeros(1,CodeL);

MaxX=ones(1,CodeL);%输出规则的最大最小值

for sc=1:1:128%sc代表128个基因的序数，生成初始群体用

Out(:,sc)=MinX(1)+round((MaxX(1)-MinX(1))*rand(Size,1))+1;%Out是一个30*128的初始群体

end%每个个体里面对应的每个基因生成的基因范围都是1~4

G=100;%总代数

BsJ=0;%ITAE最优指标初始化

%%%%%%%%%%     Start Running    %%%%%%%%%%%%

for kg=1:1:G%kg是代数

time(kg)=kg;

%****** Step 1 : Evaluate BestJ ******

for i=1:1:Size

Outi=Out(i,:);

[Outi,BsJ]=zixitong(Outi,BsJ);%调用子函数

BsJi(i)=BsJ;    %最优代价值

end

[OderJi,IndexJi]=sort(BsJi);  %排序

BestJ(kg)=OderJi(1);   %取最小值

BJ=BestJ(kg);

Ji=BsJi+1e-10;    %Avoiding deviding zero，1的-10次方很小 可以忽略不计 避免分母为零

fi=1./Ji;  %适应函数值

% Cm=max(Ji);

% fi=Cm-Ji;

[Oderfi,Indexfi]=sort(fi);   %Arranging fi small to bigger 适应函数值排序

Bestfi=Oderfi(Size);          %LetBestfi=max(fi)取最大值

BestS=Out(Indexfi(Size),:);  %LetBestS=E(m), m is the Indexfi belong to max(fi)

%****** Step 2 : Select and Reproduct Operation******

fi_sum=sum(fi);

fi_Size=(Oderfi/fi_sum)*Size;

fi_S=floor(fi_Size);              % Selecting Bigger fi value取较大的适应值确定其位置

r=Size-sum(fi_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;     % Adding rest to equal Size

end

k=1;

for i=Size:-1:1       % Select the Sizeth and Reproducefirstly

for j=1:1:fi_S(i)

TempE(k,:)=Out(Indexfi(i),:);      % Select and Reproduce

k=k+1;                            % k is used toreproduce

end

end

%************ Step 3 : Crossover Operation************

Pc=0.90;

for i=1:2:(Size-1)

temp=rand;

if Pc>temp                      %Crossover Condition

alfa=rand;

TempE(i,:)=round(alfa*Out(i+1,:)+(1-alfa)*Out(i,:));

TempE(i+1,:)=round(alfa*Out(i,:)+(1-alfa)*Out(i+1,:));%因为控制规则矩阵元素必须是整数

end

end

TempE(Size,:)=BestS;

Out=TempE;

%************ Step 4: Mutation Operation**************变异的结果也可能不为整数！

Pm=0.10-[1:1:Size]*(0.01)/Size;       %Bigger fi,smaller Pm

Pm_rand=rand(Size,CodeL);

Mean=(MaxX+1+ MinX+1)/2;

Dif=(MaxX-MinX);

for i=1:1:Size

for j=1:1:CodeL

if Pm(i)>Pm_rand(i,j)        %Mutation Condition

TempE(i,j)=round(Mean(j)+Dif(j)*(rand-0.5));

end

end

end

%Guarantee TempE(Size,:) belong to the best individual

TempE(Size,:)=BestS;

Out=TempE;

end

Bestfi

BestS

Best_J=BestJ(G)

figure(1);

plot(time,BestJ,'r','linewidth',2);

xlabel('Times');ylabel('Best J');

figure(2);

plot(timef,yd,'r',timef,y,'b:','linewidth',2);

xlabel('Time(s)');ylabel('yd,y');

legend('Ideal position signal','Positionsignal tracking');