伺服系统刚性模型的建立

一.系统工作原理

为了实现对运动控制系统精准的位置控制，首先要对伺服进给系统进行准确建模和模型辨识。人们对于运动控制系统的研究中已经提出了多种多样的系统建模和辨识方法。

图1 伺服电机滚珠丝杠传动系统刚性模型

下面对整个系统的工作原理进行解释：

1.$u^{{}^{\prime }}(V)$是输入电机伺服单元的控制电压信号，由控制器生成；

2.$K_a(A/V)$是伺服单元的放大增益，生成电机的电流指令$i$，电流$i$输入给电机本体；

3.电机模型，电流$i$通过电机模型转换为电机轴上的扭矩$T_m$。电机的电流动态可以用通过L和R来描述。$K_t(Nm/A)$是电机的转矩常量，定义了电流到扭矩的转换比例。

4.除了电机输出的转矩$T_m(Nm)$外，同时伺服轴还收到外部扰动$T_d(Nm)$的作用，外部扰动通常由导轨、丝杠的摩擦力和切削力等，它可以直接作用于电机产生的扭矩。

5.机械系统的模型参数主要包括各轴对一个的等效惯量$J(kg\cdot m^2)$和粘滞阻尼系数$B(kg\cdot m^2/s)$，电机的机械动态（扭矩到角速度的转换）由传递函数$\frac{1}{Js+B}$描述，它将电机扭矩$T_m$转换电机角速度$w$

6.$T_m$和$T_d$的差值作用于机械系统，转化为电机轴的角速度$w(rad/s)$，并通过导程为$R_g(mm/rad)$滚珠丝杠的速度$\dot{x}(mm/s)$,其积分就是工作台或者刀具的位置$x(mm)$。

7.输出积分，积分环节$\frac{1}{s}$将速度$\dot{x}(mm/s)$转换为位移$x(mm)$，即系统的最终输出。

二.复数域系统模型表示

伺服进给系统连续函数模型可以表示为：
$$x(s)=\frac{r_g}{s}\frac{1}{Js+B}(K_t{K}_a{u}^{{}^{\prime }}(s)-T_d(s))\text{\hspace{1em}(1)}$$
其中$T_d(s)$为外部扰动，由两部分组成如下式，$T_f(s)$表示摩擦力，$T_c(s)$表示切削力，在机床空切时，切削力为。
$$T_d(s)=T_f(s)+T_c(s)\text{\hspace{1em}(2)}$$
定义增益$K_w=\frac{K_t{K}_a}{J}$和极点$p_w=\frac{-B}{J}$以及外部扰动$T_d(s)$可以假设为电压值，用$d(s)=\frac{T_d(s)}{K_t{K}_a}$表示。

将外部扰动$T_d(s)$假设为电压值。我们可以将扭矩扰动$d(s)$通过$K_t$ 转换成电流扰动，然后用$K_a$ 将它转换为等效的电压扰动，以便在控制系统设计中进行处理。这样，我们可以在控制器设计中考虑到扰动，而不需要直接处理机械扭矩扰动。通过这种方式，控制器可以设计出来抵抗这些电压扰动，从而间接抵抗机械扰动。

根据上述定义的关系式，式（1）中的电机角速度$w$和伺服轴的位置$x$可以表示为：
$$\begin{gathered} w(s)=\frac{K_w}{s-p_w}(u^{{}^{\prime }}(s)-d(s)) \\ x(s)=\frac{r_g}{s}\frac{K_w}{s-p_w}(u^{{}^{\prime }}(s)-d(s))\text{\hspace{1em}(3)} \end{gathered}$$
上述式子$x(s)$分子和分母同时除以$J$,然后把定义的增益$K_w$、极点$p_w$以及将扰动转换为电压的$d(s)$带进去。

三.状态空间法表示

在经典控制理论中，对一个线性定长系统，可用常微分方程或者传递函数加以描述，可将某个单变量作为输出，直接和输入联系起来。实际上系统除了输出量这个变量之外，还包含其他相互独立的变量，而微分方程或者传递函数对这些内部的中间变量是不便描述的，因而不能包含系统的所有信息。显然，从能否完全揭示系统的全部运动状态来说，用微分方程或传递函数来描述一个线性定长系统有其不足之处。

在用状态空间法分析系统时，系统的动态特性是用由状态变量构成的一阶微分方程组来描述的。他能反应系统的全部独立变量的变化，从而能同时确定系统的全部内部运动状态，而且还可以方便地处理初始条件。这样，在设计控制系统时，不再只局限于输入量、输出量、误差量，为提高系统性能提供了有力的工具。价值可利用计算机进行分析设计及实时控制，因而可以应用与非线性系统、时变系统、多输入——多输出系统以及随机过程等。

用状态空间表示滚柱丝杠进给系统，如下式：
$$\left[\begin{array}{c}\dot{x}(t)\\ \dot{w}(t)\end{array}\right]=A_c\left[\begin{array}{c}x(t)\\ w(t)\end{array}\right]+\left[\begin{array}{c}{B}_c-B_c\end{array}\right]\left[\begin{array}{c}{u}^{{}^{\prime }}(t)\\ d(t)\end{array}\right],其中A_c=\left[\begin{array}{cc}0& r_g\\ 0& p_w\end{array}\right],B_c=\left[\begin{array}{c}0\\ K_w\end{array}\right]\text{\hspace{1em}(4)}$$
我们对这个式子进行正反推导：

1.反向推导

式（4）即：
$$\begin{gathered} \left[\begin{array}{c}\dot{x}(t)\\ \dot{w}(t)\end{array}\right]=\left[\begin{array}{cc}0& r_g\\ 0& p_w\end{array}\right]\left[\begin{array}{c}x(t)\\ w(t)\end{array}\right]+\left[\begin{array}{cc}0& 0\\ K_w& -K_w\end{array}\right]\left[\begin{array}{c}{u}^{{}^{\prime }}(t)\\ d(t)\end{array}\right] \\ =\left[\begin{array}{cc}0& r_g\\ 0& p_w\end{array}\right]\left[\begin{array}{c}x(t)\\ w(t)\end{array}\right]+\left[\begin{array}{c}0\\ K_w\end{array}\right]u^{{}^{\prime }}(t)+\left[\begin{array}{c}0\\ -K_w\end{array}\right]d(t) \\ \text{\hspace{1em}(5)} \end{gathered}$$

即：

$$\begin{gathered} \dot{x}(t)=r_{g}w(t) \\ \dot{w}(t)=p_{w}w(t)+K_w{u}^{{}^{\prime }}(t)-K_{w}d(t) \end{gathered}$$
将式（3）拉普拉斯反转换成常微分方程：
$$\begin{gathered} w(s)=\frac{K_w}{s-p_w}(u^{{}^{\prime }}(s)-d(s)) \\ w(s)(s-p_w)=K_w{u}^{{}^{\prime }}(s)-K_{w}d(s) \\ w(s)s-w(s)p_w=K_w{u}^{{}^{\prime }}(s)-K_{w}d(s) \\ 进行拉普拉斯反变换结果为：\dot{w}(t)-w(t)p_w=K_w{u}^{{}^{\prime }}(t)-K_{w}d(t) \\ \dot{w}(t)=w(t)p_w+K_w{u}^{{}^{\prime }}(t)-K_{w}d(t)\text{\hspace{1em}(5)} \end{gathered}$$
对于$x(t)$，因为有一个额外的$\frac{1}{s}$,它是$w(t)$的积分，所以我们有：
$$x(t)=\int r_{g}w(t)dt\text{\hspace{1em}(6)}$$
这里$x(t)$不能直接进行拉普拉斯反变换，$x(t)$中包含两个分母项$\frac{1}{s}$和$\frac{1}{s-p_w}$，这意味着它可能对应一个复合函数，直接应用拉普拉斯反变换需要复杂的部分分式分解。

尽管可以对 $x(t)$ 直接进行反拉普拉斯变换，但实际操作中通常会根据输入信号的类型和系统分析的需求，采用更分解、具体的方式来进行，下式是直接反变换的错误结果。
$$\begin{gathered} x(s)=\frac{r_g}{s}\frac{K_w}{s-p_w}(u^{{}^{\prime }}(s)-d(s)) \\ x(s)(s^2-s{p}_w)=r_g{K}_w(u^{{}^{\prime }}(s)-d(s)) \\ 进行拉普拉斯反变换错误的结果为：\ddot{x}(t)-\dot{x}(t)p_w=r_g{K}_w(u^{{}^{\prime }}(t)-d(t))\text{\hspace{1em}(7)} \end{gathered}$$

1.正向推导

现在根据式（3）写出它的状态空间表达式

状态空间模型的标准形式为：
$$\begin{gathered} \dot{X}(t)=AX(t)+Bu(t) \\ y(t)=CX(t)+Du(t)\text{\hspace{1em}(8)} \end{gathered}$$
其中$X(t)$是状态向量，$u(t)$是输入向量，$y(t)$是输出向量，以及A，B，C和D是状态空间模型的矩阵。

选择状态向量$X(t)$为：
$$X(t)=\left[\begin{array}{c}x(t)\\ w(t)\end{array}\right]\text{\hspace{1em}(9)}$$
然后，选择输入向量$u(t)$为：
$$u(t)=\left[\begin{array}{c}{u}^{{}^{\prime }}(t)\\ d(t)\end{array}\right]\text{\hspace{1em}(10)}$$
现在可以写出状态空间模型的矩阵A，B，C和D

（1）矩阵A由系统的自然动态决定：
$$A=\left[\begin{array}{cc}0& r_g\\ 0& p_w\end{array}\right]\text{\hspace{1em}(11)}$$
（2）矩阵B描述了输入对状态的影响：
$$B=\left[\begin{array}{cc}0& 0\\ K_w& -K_w\end{array}\right]\text{\hspace{1em}(12)}$$
（3）矩阵C用于从状态向量获得输出$y(t)$:
$$C=\left[\begin{array}{cc}{r}_g& 0\end{array}\right]\text{\hspace{1em}(13)}$$
（4）矩阵D通常是零，除非有直接从输入到输出的路径：
$$D=\left[\begin{array}{cc}0& 0\end{array}\right]\text{\hspace{1em}(14)}$$
然后我们使用拉普拉斯反变换来获得$w(t)$和$x(t)$的时间域方程，对于$w(t)$，我们有:
$$\dot{w}(t)=w(t)p_w+K_w{u}^{{}^{\prime }}(t)-K_{w}d(t)\text{\hspace{1em}(15)}$$
对于$x(t)$，因为有一个额外的$\frac{1}{s}$,它是$w(t)$的积分，所以我们有：
$$x(t)=\int r_{g}w(t)dt\text{\hspace{1em}(16)}$$

$$\begin{gathered} \dot{X}(t)=\left[\begin{array}{c}\dot{x}(t)\\ \dot{w}(t)\end{array}\right]=AX(t)+Bu(t)=\left[\begin{array}{cc}0& r_g\\ 0& p_w\end{array}\right]\left[\begin{array}{c}x(t)\\ w(t)\end{array}\right]+\left[\begin{array}{cc}0& 0\\ K_w& -K_w\end{array}\right]\left[\begin{array}{c}{u}^{{}^{\prime }}(t)\\ d(t)\end{array}\right] \\ y(t)=CX(t)+Du(t)=\left[\begin{array}{cc}1& 0\end{array}\right]\left[\begin{array}{c}x(t)\\ w(t)\end{array}\right]\text{\hspace{1em}(17)} \end{gathered}$$

四.离散域表示

如果控制系统中的所有信号都是时间变量的连续函数，换句话说，这些信号在全部时间上都是已知的，则这样的系统称为连续时间系统，简称连续系统；如果控制系统中有一处或者几处信号是一串脉冲或者数码，换句话说，这些信号仅定义在离散时间上，则这样的系统称为离散时间系统，简称离散系统。

通常，把系统中的离散信号是脉冲和序列形式的离散系统，称为采样控制系统或脉冲控制系统；而把数字序列形式的离散系统，称为数字控制系统或计算机控制系统。

把来自控制器的数字信号转换为连续信号，用于这一转换过程的装置叫做保持器，从数学上说，保持器的任务是解决各采样点之间的插值问题。

1.保持器的数学描述

在采样时刻上，连续信号的函数值与脉冲序列的脉冲强度相等，在nT时刻有：
$$e(t){|}_{t=nT}=e(nT)=e^{\ast }(nT)\text{\hspace{1em}(18)}$$
在(n+1)T时刻，有：
$$e(t){|}_{t=(n+1)T}=e((n+1)T)=e^{\ast }((n+1)T)\text{\hspace{1em}(19)}$$
当脉冲序列$e^{\ast }(t)$向连续信号$e(t)$的转换过程中，在nT与(n+1)T时刻之间，即当$0<△t<T$时，连续信号$e(nT+△t)$究竟多大，他与$e(nT)$的关系如何？这就是保持器要解决的问题。

保持器是具有外推功能的元件，保持器的外推作用，表现为现在时刻的输出信号取决于过去时刻离散信号的外推。通常采用下面多项式外推公式描述保持器：
$$e(nT+△t)=a_0+a_1△t+a_2(△t{)}^2+\dots \dots +a_m(△t{)}^m\text{\hspace{1em}(20)}$$
式中，$△t$是以nT时刻为原点的坐标，式（20）表示：现在时刻的输出$e(nT+△t)$值，取决于△t=0，-T，-2T，……，-mT各过去时刻的离散信号$e^{\ast }(nT)$，$e^{\ast }((n-1)T)$，$e^{\ast }((n-2)T)$，……，$e^{\ast }((n-m)T)$的（m+1）个值。外推公式中（m+1）个待定系数$a_i$由过去各采样时刻（m+1）个离散信号值确定，所以$a_i$有唯一解。这样的保持器称为m阶保持器，若取m=0，则称零阶保持器。

2.零阶保持器

零阶保持器的外推公式为：
$$e(nT+△t)=a_0=e(nT)\text{\hspace{1em}(21)}$$
上式说明，零阶保持器是一种按常值外推的保持器，它把前一采样时刻nT的采样值一直保持到下一采样时刻之前，然后使采样信号变成阶梯信号。如果把阶梯信号的中点连起来，可以得到与连续信号形状一致但是在时间上落后T/2的响应。

式（21）还表明，零阶保持器过程是由于理想脉冲$e(nT)\delta (t-nT)$的作用结果。如果给零阶保持器输入一个理想单位脉冲$\delta (t)$，则其脉冲过度函数$g_h(t)$是幅值为1，持续时间为T的矩形脉冲，并可分解为两个单位阶跃函数的和：
$$g_h(t)=1(t)-(t-T)\text{\hspace{1em}(22)}$$
对脉冲过渡函数$g_h(t)$取拉氏变换，可得零阶保持器的传递函数：
$$G_h(t)=\frac{1}{s}-\frac{e^{-Ts}}{s}=\frac{1-e^{-Ts}}{s}\text{\hspace{1em}(23)}$$
伺服轴的输入信号$u^{{}^{\prime }}$经过D/A转换，式（3）的速度表达式加入零阶保持，然后根据z变换表，进行z变换
$$\begin{gathered} w(k)=Z[\frac{K_w}{s-p_w}(u^{{}^{\prime }}(s)-d(s))\cdot \frac{1-e^{-Ts}}{s}] \\ =Z[\frac{-p_w}{(s-p_w)s}\cdot \frac{K_w(1-e^{-Ts})}{-p_w}\cdot (u^{{}^{\prime }}(s)-d(s))] \\ =\frac{(1-e^{p_w{T}_s})z}{(z-1)(z-e^{p_w{T}_s})}\cdot \frac{K_w(1-z^{-1})}{-p_w}\cdot (u^{{}^{\prime }}(k)-d(k)) \\ =\frac{(1-e^{p_w{T}_s})z}{(z-1)(z-e^{p_w{T}_s})}\cdot \frac{K_w(z-1)}{-z{p}_w}\cdot (u^{{}^{\prime }}(k)-d(k)) \\ =\frac{1}{(z-e^{p_w{T}_s})}\cdot \frac{K_w(1-e^{p_w{T}_s})}{-p_w}\cdot (u^{{}^{\prime }}(k)-d(k))\text{\hspace{1em}(24)} \end{gathered}$$
之后离散域表达为：
$$w(k)=\frac{K_{wd}}{z-p_w}(u^{{}^{\prime }}(k)-d(k)),其中K_{wd}=\frac{K_w}{-p_w}(1-e^{p_{wd}{T}_s}),p_{wd}=e^{p_w{T}_s}\text{\hspace{1em}(25)}$$

连续时间状态空间表达式离散化以后如下式：
$$\left[\begin{array}{c}x(k+1)\\ w(k+1)\end{array}\right]=A_d\left[\begin{array}{c}x(k)\\ w(k)\end{array}\right]+\left[\begin{array}{c}{B}_d-B_d\end{array}\right]\left[\begin{array}{c}{u}^{{}^{\prime }}(k)\\ d(k)\end{array}\right],其中A_d=e^{A_c{T}_s},B_d={\int }_0^T{e}^{A_c\lambda }d\lambda \cdot B_c\text{\hspace{1em}(26)}$$

五.考虑噪声

虽然简单的刚性模型可以表示机床伺服驱动的过程，实际情况中存在多种多样的非线性误差，例如 D/A 转换的量化误差，过程噪声和测量噪声等。这些非线性误差最终都会影响跟踪精度，为了更好地控制效果，应该辨识出准确的模型参数等效惯量$J$和粘滞阻尼系数$B$ ，以及摩擦力的扰动。

加入电压输入量化噪声、位置测量噪声、转速测量噪声和摩擦力扰动噪声。

$$\begin{gathered} \left[\begin{array}{c}x(k+1)\\ w(k+1)\\ d(k+1)\end{array}\right]=A\left[\begin{array}{c}x(k)\\ w(k)\\ d(k)\end{array}\right]+Bu(k)+W\left[\begin{array}{c}\tilde{u}(k)\\ w_d(k)\end{array}\right] \\ 其中A=\left[\begin{array}{cc}{A}_d& -B_d\\ 0& 01\end{array}\right]，B=\left[\begin{array}{c}{B}_d\\ 0\end{array}\right],W=\left[\begin{array}{cc}{B}_d& 0\\ & 0\\ 0& 1\end{array}\right] \end{gathered}$$
$$\left[\begin{array}{c}{x}_m(k)\\ w_m(k)\end{array}\right]=C\left[\begin{array}{c}x(k)\\ w(k)\\ d(k)\end{array}\right]+V\left[\begin{array}{c}\tilde{x}(k)\\ \tilde{w}(k)\end{array}\right],其中C=\left[\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\end{array}\right],V=\left[\begin{array}{cc}1& 0\\ 0& 1\end{array}\right]$$
$d(k)$是由卡尔曼滤波器估计的值，既用于摩擦力辨识，也用于控制回路里面的反馈。
\end{bmatrix}
$$

$$\left[\begin{array}{c}{x}_m(k)\\ w_m(k)\end{array}\right]=C\left[\begin{array}{c}x(k)\\ w(k)\\ d(k)\end{array}\right]+V\left[\begin{array}{c}\tilde{x}(k)\\ \tilde{w}(k)\end{array}\right],其中C=\left[\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\end{array}\right],V=\left[\begin{array}{cc}1& 0\\ 0& 1\end{array}\right]$$

$d(k)$是由卡尔曼滤波器估计的值，既用于摩擦力辨识，也用于控制回路里面的反馈。

参考文献
[1]史振辉. 基于振动抑制的高带宽多轴轮廓跟踪控制研究[D].华中科技大学,2018.