八、基于干扰观测器的单机械臂滑模控制

8.1  单机械臂模型

 通过引入干扰观测器,可精确地估计被控对象的不确定性和外加干扰,从而降低滑模控制中的增益,有效地降低抖振。

        不确定单机械臂的动力学方程为

                                                                                                      \small \left ( I+\Delta I \right )\ddot{\theta }+\left ( d+\Delta d \right )\dot{\theta }+\delta _{0}\theta +mglcos\theta =u-f_{c}\left ( \dot{\theta },u\right )                    (8.25)

其中\small \theta为系统输出转角,\small I=\frac{4}{3}ml^{2}为转动惯量,\small mg为重力,\small u为控制输入,\small f_{c}\left ( \dot{\theta },u \right )为未知的非线性摩擦,质心距连杆的转动中心为\small l,连杆运动的粘性摩擦为\small d\small \Delta I\small \Delta d非别为相应参数的不确定值,\small \delta _{0}为弹性摩擦系数。

        将式(8.25)变为

                                          \small \ddot{\theta }=\frac{1}{I}u-\frac{d}{I}\dot{\theta }-\frac{\Delta I}{I}\ddot{\theta }-\frac{\Delta d}{I}\dot{\theta }-\frac{\delta _{0}}{I}{\theta }-\frac{1}{I}mglcos\theta -\frac{1}{I}f_{c}\left ( \dot{\theta },u \right )

        则不确定单机械臂可采用二阶微分方程来描述:

                                                                       \small \ddot{\theta }=-b\dot{\theta }+au-f                                                 (8.26)

        其中\small b=\frac{d}{I}> 0\small a=\frac{1}{I}> 0\small a\small b为已知值,\small f代表不确定项、重力项和摩擦项的总和, 

                                                   \small f=\frac{\Delta I}{I}\ddot{\theta }+\frac{\Delta d}{I}\dot{\theta }+\frac{\delta _{0}}{I}\theta +\frac{1}{I}mglcos\theta +\frac{1}{I}f_{c}\left ( \hat{\theta },u \right ) 

8.2  单机械臂模型的滑模控制器设计及分析

        滑模面设计为

                                                                           \small s=ce+\dot{e},c> 0                                               (8.27)

其中\small e=r-\theta\small r为位置。

        针对被控对象式(8.26),设计滑模控制律为                                                     

                                                           \small u(t)=\frac{1}{a}\left ( c\dot{e}+\ddot{r}+b\dot{\theta }+\hat{f}+k_{1} sgn\left ( s \right )\right )                              (8.28)

其中\small \hat{f}为通过干扰观测器对\small f项的估计值,\small \hat{f}\small f项的估计误差。

        切换增益系数\small k_{f}设计为

                                                                            \small k_{f}> \left | \tilde{f} \right |                                                  (8.29)

      Lyapunov函数为

                                                                            \small V_{1}=\frac{1}{2}s^{2}

        由于

                                                    \small \dot{s}=c\dot{e}+\ddot{e}=c\dot{e}+\ddot{r}-\ddot{\theta }=c\dot{e}+\ddot{r}+b\dot{\theta }-au+f

        将控制律式(8.28)代入上式,得:

                     \tiny \dot{s}=c\dot{e}+\ddot{r}+b\dot{\theta }-\left ( c\dot{e}+\ddot{r} +b\dot{\theta }+\hat{f}+k_{f}sgn\left ( s \right )\right )+f=-\left ( \hat{f}+k_{1}sgn\left ( s \right ) \right )+f=f-k_{1}sgn\left ( s \right )

                                                                     \small \dot{V}_{1}=s\dot{s}=s\tilde{f}-k_{1}\left | s \right |<0

        可见,为了满足\small \dot{V}_{1}< 0,需要满足\small k_{f}> \left | \tilde{f} \right |。如果对\small f项的估计误差\small \tilde{f}足够小,则切换增益系数\small k_{f}可设计成很小的值,从而有效地降低抖振。

8.3  干扰观测器的设计

         为了观测干扰\small f项,设计观测器为

                                               \small \begin{bmatrix} \dot{\hat{f}}\\ \dot{\hat{x}} \end{bmatrix}=\begin{bmatrix} 0 & 0\\ 1 & 0 \end{bmatrix}\begin{bmatrix} \hat{f}\\ \hat{x} \end{bmatrix}+\begin{bmatrix} 0\\ 1 \end{bmatrix}\left ( \ddot{r}+b\dot{\theta } \right )+\begin{bmatrix} 0\\ -a \end{bmatrix}u+\begin{bmatrix} k_{1}\\ k_{2} \end{bmatrix}\left [ \dot{e}-\hat{x} \right ]              (8.30)

其中\small \hat{f}为对干扰\small f的估计,\small \hat{x}为对\small \dot{e}的估计,\small \tilde{x}=\dot{e}-\hat{x}\small k_{1}\small k_{2}为通过极点配置的增益。

        干扰观测器式(9.30)表示为:

                                                                           \small \dot{\hat{f}}=k_{1}\tilde{x}                                                          (8.31)

                                                              \small \dot{\hat{x}}=\hat{f}-au+k_{2}\tilde{x}+\ddot{r}+b\dot{\theta }                                            (8.32)

8.4  仿真实例

        假设单机械臂的动态方程为

                                                                    \small \ddot{\theta }=-b\dot{\theta }+au-f                                                   (8.33)                            

其中\small a=5,b=15

        被控对象中,取\small f=5+0.15sint。采用控制律为(8.28),取\small c=3.0,干扰观测器取式(8.31)和式(8.32),取\small k_{1}=1500,k_{2}=200。取\small M=1,不采用干扰观测器,为了克服\small f项,需要设计\small k_{f}=6.0,仿真结果如图和如图所示。取\small M=2,采用干扰观测器,取\small k_{f}=0.15,仿真结果如图所示。

                                                       图9.1  控制输入信号(M=2) 

                                     图9.2  位置跟踪及跟踪误差(M=2)

                                                   图9.3  干扰及其观测结果(M=2) 

仿真程序:

simulink主程序:chap9_5sim.mdl

控制器S函数:chap9_5ctrl.m

function [sys,x0,str,ts]=spacemodel(t,x,u,flag)
switch flag,
case 0,
      [sys, x0,str,ts] = mdlInitializeSizes;
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
sizes = simsizes;
sizes.NumOutputs       =1;
sizes.NumInputs        =5;
sizes.DirFeedthrough   =1;
sizes.NumSampleTimes   =1;
sys = simsizes(sizes);
x0 = [];
str = [];
ts = [0 0];
 
function sys = mdlOutputs(t,x,u)
r = u(1);
dr = cos(t);
ddr = -sin(t);
th = u(2);
dth = u(3);
fp = u(5);
 
e = r-th;
de = dr - dth;
c = 3;
s = de+c*e;
 
b = 15;
a = 5;
 
M = 2;
if M ==1             %  Traditional with SMC 
    Kf = 6;
%     Kf = 0.15;
    ut = 1/a*(c*de+ddr+b*dth+Kf*sign(s));
elseif M ==2         % SMC with observer
    Kf = 0.15;
    ut = 1/a*(c*de+ddr+b*dth+1*fp+Kf*sign(s));
end
sys(1) = ut;

干扰观测器S函数:chap9_5obv.m

function [sys,x0,str,ts]=s_function(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
sizes = simsizes;
sizes.NumContStates    =2;
sizes.NumDiscStates    =0;
sizes.NumOutputs       =1;
sizes.NumInputs        =4;
sizes.DirFeedthrough   =0;
sizes.NumSampleTimes   =0;
sys = simsizes(sizes);
x0 = [0;0];
str = [];
ts = [];
 
function sys = mdlDerivatives(t,x,u)
r = sin(t);
dr = cos(t);
ddr = -sin(t);
 
ut = u(1);
dth = u(3);
x2 = dr - dth;
K1 = 1500;
K2 = 200;
a =5;b = 15;
sys(1) = K1*(x2-x(2));
sys(2) = x(1)-b*x(2)-a*ut+K2*(x2-x(2))+ddr+b*dth+b*x(2);
function sys = mdlOutputs(t,x,u)
 
 
sys(1) = x(1);

被控对象S函数:chap9_5plant.m

function [sys,x0,str,ts] = s_function(t,x,u,flag)
switch flag,
case 0,
      [sys,x0,str,ts] = mdlInitializesSizes;
case 1,
      sys = mdlDerivatives(t,x,u);
case 3,
      sys = mdlOutputs(t,x,u);
case{2,4,9}
      sys = [];
otnerwise
      error(['Unhandled flag = ',num2str(flag)]);
end 
 
function [sys,x0,str,ts] = mdlInitializeSizes
sizes = simsizes;
sizes.NumContStates    =2;
sizes.NumDiscStates    =0;
sizes.NumOutputs       =3;
sizes.NumInputs        =1;
sizes.DirFeedthrough   =0;
sizes.NumSampleTimes   =0;
sys = simsizes(sizes);
x0 = [0;0];
str = [];
ts = [];
 
function sys = mdlDerivatives(t,x,u)
ut = u(1);
b = 15;
a = 5;
 
f = 5+0.15*sin(t);
ddth = -b*x(2)+a*ut-f;
 
sys(1)=x(2);
sys(2)=ddth;
function sys = mdlOutputs(t,x,u)
f = 5+0.15*sin(t);
sys(1) = x(1);
sys(2) = x(2);
sys(3) = f;

作图程序:chap9_5plot.m

close all;
 
figure(1);
subplot(211);
plot(t,y(:,1),'r',y(:,2),'b');
xlabel('time(s)');
ylabel('Position tracking');
subplot(212);
plot(t,y(:,1)-y(:,2),'r');
xlabel('time(s)');
ylabel('Position tracking error');
 
figure(2);
plot(t,ut(:,1),'r');
xlabel('time(s)');
ylabel('Control input');
 
figure(3);
plot(t,f(:,3),'r',t,f(:,4),'b');
xlabel('time(s)');
ylabel('f and fp');

 Simulink仿真图和相应S-Function函数的m文件已经打包上传至资源,代码稍有改动,有需要请下载。