本文章以最简单的二阶系统为例,介绍其simulink仿真实现和m代码实现
案例中的二阶系统如下所示
经典ADRC的基本结构如下:
本案例中的simulink仿真整体结构(为便于理解,结构图与上述ADRC整体结构类似)
仿真参数初始化所需m文件程序(文章尾部附有本仿真模型及m,文件支持matlabb以上)
%-----------ADRC参数初始化------------%%参数初始化%跟踪微分器r=100;%表示跟踪快慢h0=5*h;%h0代表信号的平滑程度(滤波效果)v1_last=0;v2_last=0;v0_last=0;%扩张状态观测器beta01=10;beta02=200;beta03=30;alpha1=0.5;%文献里给定值alpha2=0.25;%给定值delta=0.0025;b=1;z1_last=0;z2_last=0;z3_last=0;%非线性误差反馈nlsef_alpha1=0.7;nlsef_alpha2=1;%被控对象初始化temp_y=[0.5;0];u_last=0;
TD微分跟踪器结构如下所示
图中fst函数利用matlab_function搭建,内部代码如下所示
%fst函数function fn=fst(x1,x2,r,h)d=h*r;d0=h*d;y=x1+h*x2;a0=sqrt(d^2+8*r*abs(y));if abs(y)<=d0a=x2+y/h;elsea=x2+0.5*(a0-d)*sgn(y);endfn=-r*sat(a,d);end%符号函数function y=sgn(x)if x>0y=1;elseif x<0y=-1;elsey=0endend%sat函数function y=sat(a,d)if abs(a)<=dy=a/d;elsey=sgn(a);endend%符号函数
ESO扩张状态观测器结构如下所示
图中fal函数利用matlab_function搭建,内部代码如下所示
function y=fal(e,alpha,delta)if abs(e)>deltay=sign(e)*abs(e)^alpha;elsey=e/(delta^(1-alpha));endend
BLSEF非线性误差反馈控制律结构如下所示
图中fal函数利用matlab_function搭建,内部代码同上
典型被控系统结构如下所示
阶跃信号下输入/输出/跟踪信号仿真结果
系统输出的微分观测效果
系统非线性项的观测效果
系统仿真的m文件实现(想要深入理解实现过程的朋友可以尝试复现该代码)
m代码(ADRC控制器)
clc;clear all;close all;%设定运行时间time=10;%设定仿真步长h=0.01;%时间定义t=0.01:h:time;%跟踪信号v0=zeros(1,time/h);for i=time/h/2+1:time/h;v0(i)=1;endrand_noise=0.05*randn(1,time/h);%跟踪信号中加入随机噪声vn=v0+rand_noise;%-----------ADRC------------%%参数初始化%跟踪微分器r=100;%表示跟踪快慢h0=5*h;%h0代表信号的平滑程度(滤波效果)v1_last=0;v2_last=0;v0_last=0;%扩张状态观测器beta01=10;beta02=200;beta03=30;alpha1=0.5;%文献里给定值alpha2=0.25;%给定值delta=0.0025;b=1;z1_last=0;z2_last=0;z3_last=0;%非线性误差反馈nlsef_alpha1=0.7;nlsef_alpha2=1;%被控对象初始化temp_y=[0.5;0];u_last=0;%----ADRC正式开始------%for k=1:time/h%第一轮迭代处理%两个参数分别为控制量和当前时间parameter1=u_last;parameter2=k*h;tSpan=[0 h];%利用龙格库塔法求解微分方程[~,total_y]=ode45('PlantModel',tSpan,temp_y,[],parameter1,parameter2);%total_state里面的元素都是龙格塔库一点点计算的结果,直接使用最后一列,即计算结果即可temp_y=total_y(length(total_y),:);%寻访最后一行,全部列的元素%记录下输出和输出的微分y(k)=temp_y(1);dy(k)=temp_y(2);%---跟踪微分器TD----%v1(k)=v1_last+h*v2_last;v2(k)=v2_last+h*fst(v1_last-vn(k),v2_last,r,h0);x3(k)=-v1_last^2;v1_last=v1(k);v2_last=v2(k);v0_last=vn(k);%----扩张状态观测器--%e=z1_last-y(k);z1(k)=z1_last+h*(z2_last-beta01*e);z2(k)=z2_last+h*(z3_last-beta02*(fal(e,alpha1,delta))+b*u_last);z3(k)=z3_last-h*beta03*(fal(e,alpha2,delta));z1_last=z1(k);z2_last=z2(k);z3_last=z3(k);%---非线性误差反馈----%e1(k)=v1(k)-z1(k);e2(k)=v2(k)-z2(k);u0(k)=beta01*fal(e1(k),nlsef_alpha1,delta)+beta02*fal(e2(k),nlsef_alpha2,delta);u(k)=u0(k)-z3(k)/b;u_last=u(k);endfigure(1);plot(t,u,'r');figure(2);subplot(311);plot(t,z1,'r',t,y,'k',t,vn,'b','linewidth',2);xlabel('time(s)');ylabel('z1,y');legend('目标输出信号','估计输出信号','实际输出信号');subplot(312);plot(t,z2,'r',t,dy,'k','linewidth',2);xlabel('time(s)'),ylabel('z2,dy');legend('估计输出微分信号','实际输出微分信号');subplot(313);plot(t,z3,'r',t,x3,'k','linewidth',2);xlabel('time(s)'),ylabel('z3,x3');legend('估计扰动','实际扰动');%---------函数部分---------%%sat函数function y=sat(a,d)if abs(a)<=dy=a/d;elsey=sgn(a);endend%符号函数function y=sgn(x)if x>0y=1;elseif x<0y=-1;elsey=0endend%fst函数function fn=fst(x1,x2,r,h)d=h*r;d0=h*d;y=x1+h*x2;a0=sqrt(d^2+8*r*abs(y));if abs(y)<=d0a=x2+y/h;elsea=x2+0.5*(a0-d)*sgn(y);endfn=-r*sat(a,d);end%fal函数function y=fal(e,alpha,delta)if abs(e)>deltay=abs(e)^alpha*sign(e);elsey=e/(delta^(1-alpha));endend
系统模型
%系统方程function dy=PlantModel(t,y,flag,p1,p2)u=p1;time=p2;dy=zeros(2,1);dy(1)=y(2);dy(2)=-y(1)^2+u;end
m文件实现的仿真结果
仿真源文件网盘链接
链接:/s/1bZ-teW6aMX2XjmLkBW1ozg
提取码:4m5p