一、问题一
复现二输入单输出模糊控制系统，改动其中一到两个环节（隶属度设置、规则等），对比修改前后控制效果。
定义模糊
%Fuzzy Control for water tank
clear all;
close all;a=newfis('fuzz_tank');%Fuzzy Inference System structurea=addvar(a,'input','e',[-3,3]); %Parameter e
a=addmf(a,'input',1,'NB','zmf',[-3,-1]);
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,'PB','smf',[1,3]);a=addvar(a,'output','u',[-4,4]); %Parameter u
a=addmf(a,'output',1,'NB','zmf',[-4,-1]);
a=addmf(a,'output',1,'NS','trimf',[-4,-2,1]);
a=addmf(a,'output',1,'Z','trimf',[-2,0,2]);
a=addmf(a,'output',1,'PS','trimf',[-1,2,4]);
a=addmf(a,'output',1,'PB','smf',[1,4]);rulelist=[1 1 1 1; %Edit rule base2 2 1 1;3 3 1 1;4 4 1 1;5 5 1 1];a=addrule(a,rulelist);a1=setfis(a,'DefuzzMethod','mom'); %Defuzzy，解模糊的方法等（centroid，bisector，middle of maximum，largest of maximum，smallest of maximum）
writefis(a1,'tank'); %Save to fuzzy file "tank.fis"
a2=readfis('tank');figure(1);
plotfis(a2);
figure(2);
plotmf(a,'input',1);
figure(3);
plotmf(a,'output',1);flag=1;
if flag==1showrule(a) %Show fuzzy rule baseruleview('tank'); %Dynamic Simulation
end
disp('-------------------------------------------------------');
disp(' fuzzy controller table:e=[-3,+3],u=[-4,+4] ');
disp('-------------------------------------------------------');for i=1:1:7e(i)=i-4;Ulist(i)=evalfis([e(i)],a2);
end
Ulist=round(Ulist)e=-3; % Error
u=evalfis([e],a2) %Using fuzzy inference进行求解
%Define N+1 triangle membership function
clear all;
close all;
N=2;x=0:0.1:100;
for i=1:N+1f(i)=100/N*(i-1);
endu=trimf(x,[f(1),f(1),f(2)]);
figure(1);
plot(x,u);for j=2:Nu=trimf(x,[f(j-1),f(j),f(j+1)]);hold on;plot(x,u);
end
u=trimf(x,[f(N),f(N+1),f(N+1)]);
hold on;
plot(x,u);
xlabel('x');
ylabel('Degree of membership');%Define N+1 triangle membership function
clear all;
close all;
z=0:0.1:60;u=trimf(z,[0,0,10]);
figure(1);
plot(z,u);u=trimf(z,[0,10,25]);
hold on;
plot(z,u);u=trimf(z,[10,25,40]);
hold on;
plot(z,u);u=trimf(z,[25,40,60]);
hold on;
plot(z,u);u=trimf(z,[40,60,60]);
hold on;
plot(z,u);xlabel('z');
ylabel('Degree of membership');修改隶属关系
%Fuzzy Control for washer
clear all;
close all;a=newfis('fuzz_wash');a=addvar(a,'input','x',[0,100]); %Fuzzy Stain
a=addmf(a,'input',1,'SD','trimf',[0,0,50]);
a=addmf(a,'input',1,'MD','trimf',[0,50,100]);
a=addmf(a,'input',1,'LD','trimf',[50,100,100]);a=addvar(a,'input','y',[0,100]); %Fuzzy Axunge
a=addmf(a,'input',2,'NG','trimf',[0,0,50]);
a=addmf(a,'input',2,'MG','trimf',[0,50,100]);
a=addmf(a,'input',2,'LG','trimf',[50,100,100]);a=addvar(a,'output','z',[0,60]); %Fuzzy Time
a=addmf(a,'output',1,'VS','trimf',[0,0,10]);
a=addmf(a,'output',1,'S','trimf',[0,10,25]);
a=addmf(a,'output',1,'M','trimf',[10,25,40]);
a=addmf(a,'output',1,'L','trimf',[25,40,60]);
a=addmf(a,'output',1,'VL','trimf',[40,60,60]);rulelist=[1 1 1 1 1; %Edit rule base1 2 3 1 1;1 3 4 1 1;2 1 2 1 1;2 2 3 1 1;2 3 4 1 1;3 1 3 1 1;3 2 4 1 1;3 3 5 1 1];a=addrule(a,rulelist);
showrule(a) %Show fuzzy rule basea1=setfis(a,'DefuzzMethod','mom'); %Defuzzy
% （1）centroid：面积重心法；
% （2）bisector：面积等分法；
% （3）mom：最大隶属度平均法；
% （4）som最大隶属度取小法；
% （5）lom：大隶属度取大法；% a1=setfis(a,'DefuzzMethod','centroid'); %Defuzzy
writefis(a1,'wash'); %Save to fuzzy file "wash.fis"
a2=readfis('wash');figure(1);
plotfis(a2);
figure(2);
plotmf(a,'input',1);
figure(3);
plotmf(a,'input',2);
figure(4);
plotmf(a,'output',1);ruleview('wash'); %Dynamic Simulationx=60;
y=70;
z=evalfis([x,y],a2) %Using fuzzy inference二、问题二
1.改变一个控制参数大小，其他参数取默认值，看控制效果
2.修改delta值看效果
3.在simulink模型上添加不同干扰、噪声，查看默认控制参数下的控制效果
作业框架（套路）：改了哪些参数、怎么改的、效果如何、为什么
控制函数
function [sys,x0,str,ts] = spacemodel(t,x,u,flag)
switch flag,
case 0,[sys,x0,str,ts]=mdlInitializeSizes;
case 1,sys=mdlDerivatives(t,x,u);
case 3,sys=mdlOutputs(t,x,u);
case {2,4,9}sys=[];
otherwiseerror(['Unhandled flag = ',num2str(flag)]);
end
function [sys,x0,str,ts]=mdlInitializeSizes
global cij bj c
sizes = simsizes;
sizes.NumContStates = 5;
sizes.NumDiscStates = 0;
sizes.NumOutputs = 2;
sizes.NumInputs = 4;
sizes.DirFeedthrough = 1;
sizes.NumSampleTimes = 0;
sys = simsizes(sizes);
x0 = 0*ones(1,5);
str = [];
ts = [];
cij=[-1 -0.5 0 0.5 1;-1 -0.5 0 0.5 1];
bj=1.0;
c=15;
function sys=mdlDerivatives(t,x,u)
global cij bj c
x1d=u(1);dx1d=cos(t);
x1=u(2);x2=u(3);
e=x1d-x1;
de=dx1d-x2;
s=c*e+de;xi=[x1;x2];
h=zeros(5,1);
for j=1:1:5h(j)=exp(-norm(xi-cij(:,j))^2/(2*bj^2));
end
gama=0.015;
W=[x(1) x(2) x(3) x(4) x(5)]';
for i=1:1:5sys(i)=-1/gama*s*h(i);
end
function sys=mdlOutputs(t,x,u)
global cij bj c
x1d=u(1);
dx1d=cos(t);ddx1d=-sin(t);
x1=u(2);
x2=u(3);
e=x1d-x1;
de=dx1d-x2;s=c*e+de;
W=[x(1) x(2) x(3) x(4) x(5)]';
xi=[x1;x2];
h=zeros(5,1);
for j=1:1:5h(j)=exp(-norm(xi-cij(:,j))^2/(2*bj^2));
end
fn=W'*h;
b=133;
xite=1.10;
%ut=1/b*(-fn+ddx1d+c*de+xite*sign(s));delta=0.05;
kk=1/delta;
if abs(s)>deltasats=sign(s);
elsesats=kk*s;
end
% ut=1/b*(-fn+ddx1d+c*de+xite*sats);
ut=1/b*(-fn+ddx1d+c*de+xite*sats+1*s);sys(1)=ut;
sys(2)=fn;绘图函数
close all;
%
% figure(1);
% subplot(211);
% plot(t,y(:,1),'k',t,y(:,2),'r:','linewidth',2);
% xlabel('time(s)');ylabel('Position tracking');
% legend('ideal position signal','practical signal');
% subplot(212);
% plot(t,cos(t),'k',t,y(:,3),'r:','linewidth',2);
% xlabel('time(s)');ylabel('Speed tracking');
% legend('ideal speed signal','practical signal');
%
% figure(2);
% plot(t,u(:,1),'k','linewidth',2);
% xlabel('time(s)');ylabel('Control input');
%
% figure(3);
% plot(t,fx(:,1),'k',t,fx(:,2),'r:','linewidth',2);
% xlabel('time(s)');ylabel('fx and estiamted fx');
% legend('fx','estiamted fx');subplot(2,3,1);
plot(t,y(:,1),'k',t,y(:,2),'r:','linewidth',2);
xlabel('time(s)');ylabel('Position tracking');
legend('ideal position signal','practical signal');
daspect([10,4,1])
subplot(2,3,4);
plot(t,cos(t),'k',t,y(:,3),'r:','linewidth',2);
xlabel('time(s)');ylabel('Speed tracking');
legend('ideal speed signal','practical signal');
daspect([10,4,1])subplot(2,3,[2 5]);
plot(t,u(:,1),'k','linewidth',2);
xlabel('time(s)');ylabel('Control input');
daspect([10,0.2,1])subplot(2,3,[3 6]);
plot(t,fx(:,1),'k',t,fx(:,2),'r:','linewidth',2);
xlabel('time(s)');ylabel('fx and estiamted fx');
legend('fx','estiamted fx');
daspect([10,40,1])三、问题三
比较遗传算法和差分进化解决同一优化问题的速度，哪个快？注明你电脑的CPU型号和主频，内存大小。轨迹规划，修改变异因子F（也称放大因子）、交叉因子（CR）、群体规模（M），查看轨迹生成效果和运行时间。
优化主程序：
clear all;
close all;
global TE G ts
Size=50; %样本个数
D=4; %每个样本有４个固定点,即分成4段
F=0.5; %变异因子
CR=0.9; %交叉因子Nmax=30; %DE优化次数TE=1; %参考轨迹参数TE
thd=0.50;
aim=[TE;thd];%摆线路径终点start=[0;0];%路径起点
tmax=3*TE; %仿真时间n=500;
ts=TE/(2*n); %将TE分为1000个点，每段长度(步长)为tsG=tmax/ts; %仿真时间为G=3000
%***************摆线参考轨迹*************%
th0=0;
for k=1:1:G
t(k)=k*ts; %t(1)=0.001;t(2)=0.002;.....
if t(k)<TEthr(k)=(thd-th0)*(t(k)/TE-1/(2*pi)*sin(2*pi*t(k)/TE))+th0; %不含原点的参考轨迹(1)
else thr(k)=thd;
end
end
%(1)初始化路径
for i=1:Sizefor j=1:DPath(i,j)=rand*(thd-th0)+th0;end
end%**********差分进化计算***************%
for N=1:Nmax
%(2)变异for i=1:Sizer1=ceil(Size*rand);r2=ceil(Size*rand);r3=ceil(Size*rand);while(r1==r2||r1==r3||r2==r3||r1==i||r2==i||r3==i)%选取不同的r1,r2,r3，且不等于ir1=ceil(Size*rand);r2=ceil(Size*rand);r3=ceil(Size*rand);endfor j=1:Dmutate_Path(i,j)=Path(r1,j)+F.*(Path(r2,j)-Path(r3,j));%选择前半部分产生变异个体end
%(3)交叉for j=1:Dif rand<=CRcross_Path(i,j)=mutate_Path(i,j);elsecross_Path(i,j)=Path(i,j);endend
%(4)选择算法 %先进行三次样条插值，此为D=4时的特殊情况%XX(1)=0;XX(2)=200*ts;XX(3)=400*ts;XX(4)=600*ts;XX(5)=800*ts;XX(6)=1000*ts;YY(1)=th0;YY(2)=cross_Path(i,1);YY(3)=cross_Path(i,2);YY(4)=cross_Path(i,3);YY(5)=cross_Path(i,4);YY(6)=thd;dY=[0 0];cross_Path_spline=spline(XX,YY,linspace(0,1,1000));%输出插值拟合后的曲线，注意步长nt的一致,此时输出1000个点YY(2)=Path(i,1);YY(3)=Path(i,2);YY(4)=Path(i,3);YY(5)=Path(i,4);Path_spline=spline(XX,YY,linspace(0,1,1000));%*** 计算指标并比较***%for k=1:1000 distance_cross(i,k)=abs(cross_Path_spline(k)-thr(k)); %计算交叉后的轨迹与参考轨迹的距离值distance_Path(i,k)=abs(Path_spline(k)-thr(k)); %计算插值后的轨迹与参考轨迹的距离值endnew_object = chap15_8obj(cross_Path_spline,distance_cross(i,:),0); %计算交叉后的能量消耗最低及路径逼近最佳值的和formal_object = chap15_8obj(Path_spline,distance_Path(i,:),0); %计算插值后的能量消耗最低及路径逼近最佳值的和if new_object<=formal_objectFitness(i)=new_object;Path(i,:)=cross_Path(i,:);elseFitness(i)=formal_object;Path(i,:)=Path(i,:);end
%结束end[iteraion_fitness(N),flag]=min(Fitness);%记下第NC次迭代的最小数值及其维数lujing(N,:)=Path(flag,:) %第NC次迭代的最佳路径fprintf('N=%d Jmin=%g\n',N,iteraion_fitness(N));
end
[Best_fitness,flag1]=min(iteraion_fitness);
Best_solution=lujing(flag1,:);
YY(2)=Best_solution(1);YY(3)=Best_solution(2);YY(4)=Best_solution(3);YY(5)=Best_solution(4);Finally_spline=spline(XX,YY,linspace(0,1,1000));
chap15_8obj(Finally_spline,distance_Path(Size,:),1);figure(3);
plot((0:0.001:tmax),[0,thr(1:1:3000)],'k','linewidth',2);
xlabel('Time (s)');ylabel('Ideal Path');
hold on;
plot((0:0.2:1), YY,'ko','linewidth',2);
hold on;
plot((0:0.001:1),[0,Finally_spline],'k-.','linewidth',2);
xlabel('Time (s)');ylabel('Optimized Path');
legend('Ideal Path','Interpolation points','Optimized Path');figure(4);
plot((1:Nmax),iteraion_fitness,'k','linewidth',2);
xlabel('Time (s)');ylabel('Fitness Change');目标函数程序
%***********计算控制输入能量消耗最低及路径逼近最佳值之和的子函数*************%
function Object=object(path,distance,flag) %path,distance是2000维
global TE G ts
w=0.60;
th_1=0;tol_1=0;e_1=0;
tmax=3*TE; %目标函数积分上限为3TE
thd=0.5;
thop_1=0;dthop_1=0;
x1_1=0;x2_1=0;
for k=1:1:G %Begin th(k)从2开始和thop（1）对应t(k)=k*ts;if t(k)<=TEthop(k)=path(k); %要逼近的最优轨迹dthop(k)=(thop(k)-thop_1)/ts;ddthop(k)=(dthop(k)-dthop_1)/ts;elsethop(k)=thd;dthop(k)=0;ddthop(k)=0; end %离散模型
I=1/133;b=25/133;
d(k)=1*sin(k*ts);x2(k)=x2_1+ts*1/I*(tol_1-b*x2_1+d(k));
x1(k)=x1_1+ts*x2(k);th(k)=x1(k);
dth(k)=x2(k); e(k)=thop(k)-th(k);
de(k)=(e(k)-e_1)/ts; kp=300;kd=0.30;tol(k)=kp*e(k)+kd*de(k);energy(k)=abs(tol(k)*dth(k));tol_1=tol(k);x1_1=x1(k);x2_1=x2(k);e_1=e(k);thop_1=thop(k);dthop_1=dthop(k);
end
%************计算总能量******************%
energy_all=0;
for k=1:1:Genergy_all=energy_all+energy(k);
end
dis=sum(distance);%参考轨迹的逼近误差
%********计算目标********%
Object=w*energy_all+(1-w)*dis; %used for main.m
if flag==1t(1)=0;th0=0;for k=1:1:G %>TE 不包含原点t(k)=k*ts;if t(k)<TEthr(k)=(thd-th0)*(t(k)/TE-1/(2*pi)*sin(2*pi*t(k)/TE))+th0; %不含原点的参考轨迹else thr(k)=thd;endendfigure(1);plot(t,thr,'r.-',t,thop,'k',t,th,'k-.','linewidth',2);legend('Ideal trajectory','Optimal trajectory', 'Trajectory tracking');xlabel('Time (s)');ylabel('Joint angle tracking');figure(2);plot(t,tol,'k','linewidth',2);xlabel('Time (s)');ylabel('Control input,tol');
end
end四、问题四
修改权值等参数查看算法收敛性和求解结果，修改城市数量，查看算法增长时间
% TSP Solving by Hopfield Neural Network
function TSP_hopfield()
clear all;
close all;%Step 1: Initialization
A=1.5;
D=1;
Mu=50;
Step=0.01;%Step 2: %Calculate initial route length
N=8;
cityfile = fopen('city8.txt','rt' );
cities = fscanf(cityfile, '%f %f',[ 2,inf] )
fclose(cityfile);
Initial_Length=Initial_RouteLength(cities); DistanceCity=dist(cities',cities);
%Step 3: Initialization NN
U=rands(N,N);
V=1./(1+exp(-Mu*U)); % S functionfor k=1:1:2000
times(k)=k;
%Step 4: Calculate du/dtdU=DeltaU(V,DistanceCity,A,D);
%Step 5: Calculate u(t)U=U+dU*Step;
%Step 6: Calculate output of NNV=1./(1+exp(-Mu*U)); % S function
%Step 7: Calculate energy functionE=Energy(V,DistanceCity,A,D);Ep(k)=E;
%Step 8: Check validity of the route[V1,CheckR]=RouteCheck(V);
end%Step 9: Results
if(CheckR==0)Final_E=Energy(V1,DistanceCity,A,D);Final_Length=Final_RouteLength(V1,cities); %Give final lengthdisp('Iteration times');kdisp(' the optimization route is');V1disp('Final optimization engergy function:');Final_Edisp('Initial length');Initial_Lengthdisp('Final optimization length');Final_LengthPlotR(V1,cities);
elsedisp('the optimization route is');V1disp('the route is invalid');
endfigure(2);
plot(times,Ep,'r','linewidth',2);
title('Energy Function Change');
xlabel('k');ylabel('E');% Calculate energy function
function E=Energy(V,d,A,D)
[n,n]=size(V);
t1=sumsqr(sum(V,2)-1);
t2=sumsqr(sum(V,1)-1);
PermitV=V(:,2:n);
PermitV=[PermitV,V(:,1)];
temp=d*PermitV;
t3=sum(sum(V.*temp));
E=0.5*(A*t1+A*t2+D*t3);%%%%%%% Calculate du/dt
function du=DeltaU(V,d,A,D)
[n,n]=size(V);
t1=repmat(sum(V,2)-1,1,n);
t2=repmat(sum(V,1)-1,n,1);
PermitV=V(:,2:n);
PermitV=[PermitV, V(:,1)];
t3=d*PermitV;
du=-1*(A*t1+A*t2+D*t3);%Check the validity of route
function [V1,CheckR]=RouteCheck(V)
[rows,cols]=size(V);
V1=zeros(rows,cols);
[XC,Order]=max(V);
for j=1:colsV1(Order(j),j)=1;
end
C=sum(V1);
R=sum(V1');
CheckR=sumsqr(C-R);% Calculate Initial Route Length
function L0=Initial_RouteLength(cities)
[r,c]=size(cities);
L0=0;
for i=2:cL0=L0+dist(cities(:,i-1)',cities(:,i));
end
% Calculate Final Route Length
function L=Final_RouteLength(V,cities)
[xxx,order]=max(V);
New=cities(:,order);
New=[New New(:,1)];
[rows,cs]=size(New);L=0;
for i=2:csL=L+dist(New(:,i-1)',New(:,i));
end% Give Path optimization plot
function PlotR(V,cities)
figure;cities=[cities cities(:,1)];[xxx,order]=max(V);
New=cities(:,order);
New=[New New(:,1)];subplot(1,2,1);
plot( cities(1,1), cities(2,1),'r*' ); %First city
hold on;
plot( cities(1,2), cities(2,2),'+' ); %Second city
hold on;
plot( cities(1,:), cities(2,:),'o-' ), xlabel('X axis'), ylabel('Y axis'), title('Original Route');
axis([0,1,0,1]);subplot(1,2,2);
plot( New(1,1), New(2,1),'r*' ); %First city
hold on;
plot( New(1,2), New(2,2),'+' ); %Second city
hold on;
plot(New(1,:),New(2,:),'o-');
title('TSP solution');
xlabel('X axis');ylabel('Y axis');
axis([0,1,0,1]);
axis on
