MFAC 算法基本原理是在每个工作点处,建立非线性系统等价的动态线性数据模型,利用受控系统的I/O数据在线估计系统的伪偏导数,然后设计加权一步向前的控制器,进而实现非线性系统数据驱动的无模型自适应控制。MFAC 有三种不同动态线性化方法的算法设计,即基于紧格式动态线性化的无模型自适应控制(Compact Form Dynamic Linearization based MFAC,CFDL-MFAC),基于偏格式动态线性化的无模型自适应控制(Partial Form Dynamic Linearization based MFAC,PFDL-MFAC)以及基于全格式动态线性化的无模型自适应控制(Full Form Dynamic Linearization based MFAC,FFDL-MFAC)。
本文控制方法使用紧格式动态线性化的无模型自适应控制。
参考文献:《无人艇重定义无模型自适应艏向控制方法与试验》
控制律:
控制原理:
MATLAB 实现:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 使用CFDL-MFAC(紧格式动态线性化方法)进行无人艇航向控制
%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
close all;
clear;
%% 算法参数设定
%仿真时间步长
ts=0.5;
%非线性KT方程参数k,t,α
Kit=0.701;
Tit=0.332;
afa=0.001;
% Kit=0.18;
% Tit=4.65;
% afa=1.07;
% Kit=2.36;
% Tit=5.489;
% afa=0.000094;
%风浪流干扰参数
Wn=0.958;
yip=0.9;
%风流浪干扰
%浪干扰
langrr=wgn(1,800,(0.1/0.1),'linear');
%4级海况
w0=0.808; %波浪频率
tw=6; %波浪周期
kw=0.255; %增益常数
citam=0.527; %波浪强度系数
yipp=0.3; %阻尼系数
%MFAC算法参数设置
nata=0.5; %λ
miu=1; %μ
yita=0.1; %η
ru=1; %ρ
Kz=8; %K1
yipsl=0.0000001; %ε
PI=3.1415926;
%% 参数初始化
num=[kw 0];
den=[1 (2*yipp*w0) (w0*w0)];
hs=tf(num,den);
tt=0.1:0.1:80;
[langr,xr]=lsim(hs,langrr,tt);
Y0 = [0;0];%初值
t = 0.1:ts:80;
t0=0;
Y = cell(1,length(t));
z = zeros(2,length(t));
Y0_m=[0;0]; %参考模型
Y_m = cell(1,length(t));
z_m = zeros(2,length(t));
kp=1;
kd=1;
% kp=0.98;
% kd=0.95;
ki=0.0;
error_1=0;error_2=0;
u_1=0;u_2=0;
y_1=0;y_2=0;
r_1=0;r_2=0;
sf_1=0.1;
%% 循环
for k=1:200
% d=(langr(k,1)+0.5)*PI/180;
d=0;
time(k)=k*ts;
M=0;
if M==1
rin(k)=30;
fai_r=30;
if time(k)>=20&&time(k)<40
rin(k)=0;
fai_r=0;
else
if time(k)>=40&&time(k)<60
rin(k)=30;
fai_r=30;
else
if time(k)>=60&&time(k)<80
rin(k)=0;
fai_r=0;
end
end
end
else
rin(k)=PI/3;
fai_r=PI/3;
end
Y_m{k}=rkfa_m( ts,t0,Y0_m,Wn,yip,fai_r); %龙格库塔法解航向一阶非线性方程
z_m(1,k) = Y_m{k}(1); %参考fai
z_m(2,k) = Y_m{k}(2);
yrout(k)=z_m(1,k);
% yr(1,k)=z_m(2,k); %参考r
% yr(2,k)=(yr(1,k)-dyr_1)/ts; %参考r'
% u_m(k)=(Tit*yr(2,k)+yr(1,k)+afa*yr(1,k)*yr(1,k)*yr(1,k))/Kit;
Y{k}=rkfa_GR( ts,t0,Y0,Kit,Tit,u_1,afa,d); %龙格库塔法解航向一阶非线性方程
z(1,k) = Y{k}(1); %fai
z(2,k) = Y{k}(2); %r
yout(k)=z(1,k);
r(k)=z(2,k);
error(k)=rin(k)-z(1,k);
M2=1; %0:PD控制 1:MFAC控制
if M2==0
%PD控制部分
z_error(k)=yrout(k)-z(1,k);
m=1;
if m==1
%位置式pid
u(k)=kp*error(k)+kd*(error(k)-error_1)/ts;
else
%增量式pid
xc(1)=error(k)-error_1; %Calculating P
xc(2)=error(k)-2*error_1+error_2; %Calculating D
xc(3)=error(k); %Calculating I
du(k)=kp*xc(1)+kd*xc(2)/ts+ki*xc(3)*ts;
if du(k)>PI/6
du(k)=PI/6;
end
if du(k)<-PI/6
du(k)=-PI/6;
end
u(k)=u_1+du(k);
end
else
%MFAC方法
sf(k)=sf_1+((yita*(u_1-u_2))*((yout(k)-y_1+(Kz*r(k)-(Kz*r_1)))-(sf_1*(u_1-u_2))))/(miu+(u_1-u_2)*(u_1-u_2));
if ((abs(sf(k))<yipsl)||(abs(u_1-u_2)<yipsl))
sf(k)=sf_1;
end
u(k)=u_1+((ru*sf(k))*(rin(k)-(yout(k)+Kz*r(k))))/(nata+(sf(k)*sf(k)));
sf_2=sf_1;sf_1=sf(k);
end
if u(k)>PI/6
u(k)=PI/6;
end
if u(k)<-PI/6
u(k)=-PI/6;
end
%舵速控制
% if (((u(k)-u_1)/ts)>10)
% u(k)=10*ts+u_1;
% else
% if(((u(k)-u_1)/ts)<-10)
% u(k)=-10*ts+u_1;
% end
% end
error_2=error_1;
error_1=error(k);
u_2=u_1;u_1=u(k);
y_2=y_1;y_1=yout(k);
r_2=r_1;r_1=r(k);
Y0=Y{k};
z_1=z(2,k);
Y0_m=Y_m{k};
end
%% 结果可视化
kk=1:199;
figure;
plot(kk*ts,yout(1:199)*180/PI,'r-',kk*ts,rin(1:199)*180/PI,'--');
title('船舶实际航向与设定航向');
xlabel('时间 \fontname{Times New Roman}t/s');
legend('实际航向','设定航向')
ylabel('角度 \fontname{Times New Roman}/° ');
figure;
plot(kk*ts,u(1:199)*180/PI,'r-');
title('舵角');
xlabel('时间 \fontname{Times New Roman}t/s');
ylabel('角度 \fontname{Times New Roman}/° ');
figure
plot(kk*ts,error(1:199)*180/PI,'r');
title('航向误差');
xlabel('时间 \fontname{Times New Roman}t/s');
ylabel('角度 \fontname{Times New Roman}/° ');
% figure
% plot(kk*ts,rin(1:199),'--',kk*ts,yrout(1:199),'b',kk*ts,yout(1:199),'r-',kk*ts,u(1:199),'g');
% title('船舶航向角与舵角');
% axis([0 80 -20 40])
% xlabel('时间 \fontname{Times New Roman}t/s');
% legend('设定航向角','参考航向角','实际航向角','实际舵角')
% ylabel('角度 \fontname{Times New Roman}/° ');
function Y = rkfa_GR( h,t0,Y0,K,T,u,alfa,d )
k1 = fai_GR(t0,Y0,K,T,u,alfa,d);
k2 = fai_GR(t0+h/2,Y0+h/2*k1,K,T,u,alfa,d);
k3 = fai_GR(t0+h/2,Y0+h/2*k2,K,T,u,alfa,d);
k4 = fai_GR(t0+h,Y0+h*k3,K,T,u,alfa,d);
Y = Y0+h*(k1+2*k2+2*k3+k4)/6;
end
function Y = fai_GR( ~,Y0,K,T,u,alfa,d )
% 此处Y0为一个列向量,因为时间t未显含在一阶方程组中
% 所以ode函数的第一个参数为空,要根据具体情况而定。
Y = [Y0(2);
(K*(u+d)-Y0(2)-alfa*(Y0(2))^3)/T;];
end
仿真结果分析:
PID:
MFAC:
PID参数手动调优的效果比MFAC的最优效果要好。
当加上干扰后:
PID:
MFAC:
加上干扰后,PID在稳定后依然有误差,且无法真正稳定,MFAC则不需要调整参数也能保证优秀的控制效果。
但是MFAC在调整时间上始终比PID控制慢,需要进一步研究解决。