一、问题一
复现二输入单输出模糊控制系统，改动其中一到两个环节（隶属度设置、规则等），对比修改前后控制效果。
定义模糊
%Fuzzy Control for water tank
clear all;
close all;
a=newfis('fuzz_tank');%Fuzzy Inference System structure
a=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 base
2 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==1
showrule(a) %Show fuzzy rule base
ruleview('tank'); %Dynamic Simulation
end
disp('-------------------------------------------------------');
disp(' fuzzy controller table:e=[-3,+3],u=[-4,+4] ');
disp('-------------------------------------------------------');
for i=1:1:7
e(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+1
f(i)=100/N*(i-1);
end
u=trimf(x,[f(1),f(1),f(2)]);
figure(1);
plot(x,u);
for j=2:N
u=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 base
1 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 base
a1=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 Simulation
x=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=[];
otherwise
error(['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:5
h(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:5
sys(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:5
h(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)>delta
sats=sign(s);
else
sats=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个点，每段长度(步长)为ts
G=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)<TE
thr(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:Size
for j=1:D
Path(i,j)=rand*(thd-th0)+th0;
end
end
%**********差分进化计算***************%
for N=1:Nmax
%(2)变异
for i=1:Size
r1=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，且不等于i
r1=ceil(Size*rand);
r2=ceil(Size*rand);
r3=ceil(Size*rand);
end
for j=1:D
mutate_Path(i,j)=Path(r1,j)+F.*(Path(r2,j)-Path(r3,j));%选择前半部分产生变异个体
end
%(3)交叉
for j=1:D
if rand<=CR
cross_Path(i,j)=mutate_Path(i,j);
else
cross_Path(i,j)=Path(i,j);
end
end
%(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)); %计算插值后的轨迹与参考轨迹的距离值
end
new_object = chap15_8obj(cross_Path_spline,distance_cross(i,:),0); %计算交叉后的能量消耗最低及路径逼近最佳值的和
formal_object = chap15_8obj(Path_spline,distance_Path(i,:),0); %计算插值后的能量消耗最低及路径逼近最佳值的和
if new_object<=formal_object
Fitness(i)=new_object;
Path(i,:)=cross_Path(i,:);
else
Fitness(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)<=TE
thop(k)=path(k); %要逼近的最优轨迹
dthop(k)=(thop(k)-thop_1)/ts;
ddthop(k)=(dthop(k)-dthop_1)/ts;
else
thop(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:G
energy_all=energy_all+energy(k);
end
dis=sum(distance);%参考轨迹的逼近误差
%********计算目标********%
Object=w*energy_all+(1-w)*dis; %used for main.m
if flag==1
t(1)=0;
th0=0;
for k=1:1:G %>TE 不包含原点
t(k)=k*ts;
if t(k)<TE
thr(k)=(thd-th0)*(t(k)/TE-1/(2*pi)*sin(2*pi*t(k)/TE))+th0; %不含原点的参考轨迹
else
thr(k)=thd;
end
end
figure(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 function
for k=1:1:2000
times(k)=k;
%Step 4: Calculate du/dt
dU=DeltaU(V,DistanceCity,A,D);
%Step 5: Calculate u(t)
U=U+dU*Step;
%Step 6: Calculate output of NN
V=1./(1+exp(-Mu*U)); % S function
%Step 7: Calculate energy function
E=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 length
disp('Iteration times');k
disp(' the optimization route is');V1
disp('Final optimization engergy function:');Final_E
disp('Initial length');Initial_Length
disp('Final optimization length');Final_Length
PlotR(V1,cities);
else
disp('the optimization route is');V1
disp('the route is invalid');
end
figure(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:cols
V1(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:c
L0=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:cs
L=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
