用s函数实现Lugre摩擦模型

1. Lugre摩擦模型

2. Lugre摩擦模型的状态空间方程

3. Lugre摩擦模型的simulink建模

4. Lugre摩擦模型的s函数建模

5. simulink模型和s函数模型对比

6. Lugre摩擦模型的s函数源代码

1. Lugre摩擦模型

摩擦模型主要分为动态和静态摩擦模型。Lugre摩擦模型能够很好的反应出摩擦的静态和动态效应。Lugre摩擦模型：

$\dot{z} = \dot{\theta } - \sigma _{0}\frac{| \dot{\theta }|}{g( \dot{\theta })}z$ ①

$F = \sigma _{0}z+\sigma _{1}\dot{z}+\sigma _{2}\dot{\theta }$ ②

其中 ${g( \dot{\theta })}$ 为非线性函数：

${g( \dot{\theta })} = F_c+(F_s-F_c)e^{-(\frac{\dot{\theta }}{w_s})^2}$ ③

表1 参数意义及取值
参数	意义	本文取值
$\dot{\theta }$	接触面相对角速度	——
z	鬃毛的平均形变	——
$F_c$	库伦摩擦力	2.4440
$F_s$	最大静摩擦力	0.5991
$w_s$	Stribeck速度	0.0103
$\sigma_0$	鬃毛刚性系数	0.4766
	滑动阻尼系数	0.2701
$\sigma_2$	黏性摩擦因数	0.0049

这些参数根据不同的系统取不同的值，且会在外界温度、润滑条件、材料磨损及接触面作用力等因素的影响下发生改变（取值方法：查论文，一般通过参数辨识来计算取值；以上取值仅为示例）

2. Lugre摩擦模型的状态空间方程

$\dot{z} = \dot{\theta } - \sigma _{0}\frac{| \dot{\theta }|}{g( \dot{\theta })}z$ ④

$z = z + \dot{z}T$ ⑤

其中，T为离散采样时间， $T=\Delta t$

$\because \dot{z} = \frac{dz}{dt}\approx \frac{\Delta z}{\Delta t}$

$\because z\approx z+\Delta z$

$\therefore z = z+\dot{z}\Delta t = z+\dot{z}T$

取状态变量：

$x = \begin{bmatrix} z\\ \dot{z} \end{bmatrix}$ ⑥

则根据表达式④、⑤有状态方程：

$\begin{bmatrix} z\\ \dot{z} \end{bmatrix} = \begin{bmatrix} 1 & T\\ -\sigma _{0}\frac{| \dot{\theta }|}{g( \dot{\theta })} & 0 \end{bmatrix}\begin{bmatrix} z\\ \dot{z} \end{bmatrix}+\begin{bmatrix} 0\\ 1 \end{bmatrix}\dot{\theta }$ ⑦

根据表达式②有输出方程：

$F =\begin{bmatrix} \sigma_0&\sigma_1 \end{bmatrix} \begin{bmatrix} z\\\dot{z} \end{bmatrix}+\sigma_2\dot{\theta }$ ⑧

3. Lugre摩擦模型的simulink建模

根据表达式 ①②③搭建simulink仿真模型如图1所示：

参数初始化：（按照表1中参数初始化，仿真前先运行init.m文件，将参数加载到工作空间）

输入信号：

输出波形：

4. Lugre摩擦模型的s函数建模

根据表达式⑦⑧和表1中的参数编写s函数，取离散采样时间为0.001s：

$\begin{bmatrix} z\\ \dot{z} \end{bmatrix} = \begin{bmatrix} 1 & 0.001\\ -\frac{0.4766| \dot{\theta }|}{2.4440+(0.5991-2.4440)e^{-(\frac{\dot{\theta }}{0.0103})^2}} & 0 \end{bmatrix}\begin{bmatrix} z\\ \dot{z} \end{bmatrix}+\begin{bmatrix} 0\\ 1 \end{bmatrix}\dot{\theta }$ ⑨

$F =\begin{bmatrix} 0.4766&0.2701 \end{bmatrix} \begin{bmatrix} z\\\dot{z} \end{bmatrix}+0.0049\dot{\theta }$ ⑩

离散变量、输入、输出设置：

sizes.NumDiscStates  = 2;   %离散变量个数
sizes.NumOutputs     = 1;   %输出个数
sizes.NumInputs      = 1;   %输入个数

初始条件设置：

x0  = [0;0];

采样时间设置：

ts  = [0.001 0];

表达式⑨：

function sys=mdlUpdate(t,x,u)
w = u(1);
gw = 2.4440 + (0.5991 - 2.4440)*exp(-(w/0.0103)^2);
sys = [1 0.001;-0.4766*((abs(w))/gw) 0]*x+[0;1]*w;

表达式⑩：

function sys=mdlOutputs(t,x,u)
w = u(1);
sys = [0.4766 0.2701]*x+0.0049*w;

输入波形：

输出波形：

5. simulink模型和s函数模型对比

图8中上半部分为simulink模型，下半部分s-function模型

根据图9可知s函数模型和simulink模型输出波形基本一致。s函数离散变量的初始状态会影响0时刻的仿真值，修改s函数离散变量初始值会减小0时刻的误差。

x0  = [0;-10];

误差在0.001级，可以忽略。

6. Lugre摩擦模型的s函数源代码

function [sys,x0,str,ts,simStateCompliance] = Lugre(t,x,u,flag)switch flagcase 0   %初始化[sys,x0,str,ts,simStateCompliance]=mdlInitializeSizes;case 1   %计算连续状态的微分sys=mdlDerivatives(t,x,u);case 2   %计算下一个离散状态  sys=mdlUpdate(t,x,u);case 3   %计算输出sys=mdlOutputs(t,x,u);case 4   %计算下一次采样的时间sys=mdlGetTimeOfNextVarHit(t,x,u);case 9   %系统结束sys=mdlTerminate(t,x,u);otherwise %其他情况DAStudio.error('Simulink:blocks:unhandledFlag', num2str(flag));end%case 0:初始化
function [sys,x0,str,ts,simStateCompliance]=mdlInitializeSizessizes = simsizes;           %用于设置模块参数的结构体sizes.NumContStates  = 0;   %连续变量个数
sizes.NumDiscStates  = 2;   %离散变量个数
sizes.NumOutputs     = 1;   %输出个数
sizes.NumInputs      = 1;   %输入个数
sizes.DirFeedthrough = 1;   %输入是否直接影响输出
sizes.NumSampleTimes = 1;   %至少一个采样时间sys = simsizes(sizes);%初始化初始条件
x0  = [0;0];%保留项
str = [];% 初始化采样时间
ts  = [0.001 0]; %[0 0] 连续采样%[0 1] 小步长的连续采样%[PERIOD OFFSET] 离散采用时间 [采样时间 步长]%[-2 0] 变步长的采样时间 FLAG=4时获取下一次采样时间%指定simStateCompliance块值
simStateCompliance = 'UnknownSimState'; %'UnknownSimState'  默认设置%'DefaultSimState'  与内置块的模拟状态相同%'HasNoSimState'    没有模拟状态%'DisallowSimState' 保存或恢复模型模拟状态时出错%case 1:计算连续状态的微分
%sys=AX+BU
function sys=mdlDerivatives(t,x,u)
sys = [];%case 2:计算下一个离散状态
%sys=AX+BU
function sys=mdlUpdate(t,x,u)
w = u(1);
gw = 2.4440 + (0.5991 - 2.4440)*exp(-(w/0.0103)^2);
sys = [1 0.001;-0.4766*((abs(w))/gw) 0]*x+[0;1]*w;%case 3:计算输出
%sys = CX+DU
function sys=mdlOutputs(t,x,u)
w = u(1);
sys = [0.4766 0.2701]*x+0.0049*w;%case 4:计算下一次采样的时间[-2 0]时使用该函数
function sys=mdlGetTimeOfNextVarHit(t,x,u)
sampleTime = 0.01;
sys = t + sampleTime;%case 9:系统结束
function sys=mdlTerminate(t,x,u)
sys = [];