CN101261147A - 基于模型化测量和动态重量变化率预测的动态称重方法 - Google Patents

Abstract

本发明公开了一种基于模型化测量和动态重量变化率预测的动态称重方法。涉及运用嵌入式(DSP)计算机这类哈佛总线的新型微处理机，构成测量和控制系统，采用模型化测量和动态重量变化率预测相结合的动态信号处理方式，以解决动态称重中既高精度又高速的问题。为低成本地实现固体粉状或颗粒状产品物流高速定量称重包装。

Description

基于模型化测量和动态重量变化率预测的动态称重方法

技术领域

本发明涉及称重技术，尤其是涉及一种基于模型化测量和动态重量变化率预测的动态称重方法。

背景技术

在食品加工、轻化工、粮食加工、饲料加工、制药、冶金等行业中有大量的固体粉状或颗粒状产品是以定量包装的形式与消费者见面的，如糖、盐、味精、洗衣粉、药品、化工原料等。目前固体粉状或颗粒状产品在最后的定量称重包装阶段，有相当一部分企业还是采用人工称重、人工包装，或用速度较慢的自动定量称重包装设备实现。这对于企业大规模生产、提高劳动生产率来说已成一种瓶颈或制约，出现这类情况主要是因为现已有的定量称重包装设备要么称重、包装过程的速度慢。本发明就是为低成本地实现高速定量称重包装。

发明内容

本发明的目的在于提供一种基于模型化测量和动态重量变化率预测的动态称重方法，为低成本地实现固体粉状或颗粒状产品物流高速定量称重包装。

本发明解决其技术问题所采用的技术方案是：

一种基于模型化测量和动态重量变化率预测的动态称重方法，该方法的步骤如下：

1)模型化测量：

设入测量系统的输入即被测量为u_i(k)，测量系统的输出即测量输出值为u_o(k)；u_i(k)是固体粉状或颗粒状物料落入称重料斗的变化质量，u_o(k)是测量系统的测量输出值，当测量系统是理想系统时，u_i(k)、u_o(k)的动态关系描述为

u_o(k)＝K·u_i(k)......................(1)

式(1)中K是常数。而实际的测量系统是一个线性定常离散动态系统，其差分方程的描述为

a₀u_o(k)+a₁u_o(k-1)+…+a_nu_o(k-n)＝b₀·u_i(k)+b₁·u_i(k-1)+…+b_m·u_i(k-m)

(m≤n)...........................(2)

为了获得式(2)差分方程，采用阶跃输入方法，观察并获取输出信号，将获取的输出信号和阶跃输入信号代入式(2)，采用最小二乘法估算出a_o、a₁、…、a_n和b_o、b₁、…、b_n，将式(2)写成算子形式，并令n＝m、a_o＝1，则：

A(q^-1)u_o(k)＝B(q^-1)·u_i(k)....................(5)

式(5)中

A(q^-1)＝1+a₁q^-1+a₂q^-2+…+a_nq^-n

B(q^-1)＝b₀+b₁q^-1+b₂q^-2+…+b_nq^-n

考虑到模型的误差和测量噪声，式(5)方程改写为

A(q^-1)u_o(k)＝B(q^-1)·u_i(k)+e(k)....................(6)

式(6)中e(k)称为模型残差或方程误差，将式(6)改写成

u_o(k)＝-a₁u_o(k-1)-…-a_nu_o(k-n)+b₀·u_i(k)+b₁·u_i(k-1)+…+b_n·u_i(k-n)+e(k)

...........................(7)

定义

$\stackrel{\→}{x}\left(k\right)={[{-u}_o(k-1),\·\·\·,{-u}_o(k-n),u_i\left(k\right),u_i(k-1),\·\·\·,u_i(k-n)]}^T$

$\stackrel{\→}{\mathrm{\θ}}={[a_1,\·\·\·,a_n,b_0,b_1,\·\·\·,b_n]}^T$

式(7)可写为

$u_o\left(k\right)={\stackrel{\→}{x}}^T\left(k\right)\·\stackrel{\→}{\mathrm{\θ}}+e\left(k\right)\·\·\·\left(8\right)$

将测得的(N+n)个数据代入式(8)建立N个方程组，则写成矩阵方程的形式

${\stackrel{\→}{u}}_o=\stackrel{\‾}{X}\·\stackrel{\→}{\mathrm{\θ}}+\stackrel{\→}{e}\·\·\·\left(9\right)$

其中

${\stackrel{\→}{u}}_o=\left[\begin{array}{c}{u}_o(n+1)\\ u_o(n+2)\\ \·\\ \·\\ \·\\ u_o(n+N)\end{array}\right],$ $\stackrel{\→}{e}=\left[\begin{array}{c}e(n+1)\\ e(n+2)\\ \·\\ \·\\ \·\\ e(n+N)\end{array}\right],$

$\stackrel{\‾}{X}=\left[\begin{array}{c}{\stackrel{\→}{x}}^T(n+1)\\ {\stackrel{\→}{x}}^T(n+2)\\ \·\\ \·\\ \·\\ \stackrel{\→}{x}_T(n+N)\end{array}\right]=\left[\begin{array}{cccccc}-u_o\left(n\right)& \·\·\·& -u_o\left(1\right)& u_i(n+1)& \·\·\·& u_i\left(1\right)\\ -u_o(n+1)& \·\·\·& -u_o\left(2\right)& u_i(n+2)& \·\·\·& u_i\left(2\right)\\ \·& & \·& \·& & \·\\ \·& & \·& \·& & \·\\ \·& & \·& \·& & \·\\ -u_o(n+N-1)& \·\·\·& -u_o\left(N\right)& u_i(n+N)& \·\·\·& u_i\left(N\right)\end{array}\right]$

根据最小二乘法，选择

使得误差准则函数J为最小：

$J=\underset{k=n+1}{\overset{N+n}{\mathrm{\Σ}}}{e}^2\left(k\right)={\stackrel{\→}{e}}^T\stackrel{\→}{e}={(\stackrel{\→}{u}-\stackrel{\‾}{X}\·\stackrel{\→}{\mathrm{\θ}})}^T(\stackrel{\→}{u}-\stackrel{\‾}{X}\·\stackrel{\→}{\mathrm{\θ}})$

令 $\begin{array}{c}{\frac{\∂J}{\∂\stackrel{\→}{\mathrm{\θ}}}|}_{\stackrel{\→}{\mathrm{\θ}}=\widehat{\stackrel{\→}{\mathrm{\θ}}}}=0\\ {\frac{\∂J}{\∂\stackrel{\→}{\mathrm{\θ}}}|}_{\stackrel{\→}{\mathrm{\θ}}=\widehat{\stackrel{\→}{\mathrm{\θ}}}}=\left[\frac{\∂J}{\∂\stackrel{\→}{\mathrm{\θ}}}{(\stackrel{\→}{u}-\stackrel{\‾}{X}\·\stackrel{\→}{\mathrm{\θ}})}^T\right](\stackrel{\→}{u}-\stackrel{\‾}{X}\·\stackrel{\→}{\mathrm{\θ}})+{(\stackrel{\→}{u}-\stackrel{\‾}{X}\·\stackrel{\→}{\mathrm{\θ}})}^T\left[\frac{\∂J}{\∂\stackrel{\→}{\mathrm{\θ}}}(\stackrel{\→}{u}-\stackrel{\‾}{X}\·\stackrel{\→}{\mathrm{\θ}})\right]\end{array}$

则 $0=-2{\stackrel{\‾}{X}}^T\·\stackrel{\→}{u}+2{\stackrel{\‾}{X}}^T\·\stackrel{\‾}{X}\·\stackrel{\→}{\mathrm{\θ}}$

求得

的估计值

$\widehat{\stackrel{\→}{\mathrm{\θ}}}={({\stackrel{\‾}{X}}^T\·\stackrel{\‾}{X})}^{-1}{\stackrel{\‾}{X}}^T\·\stackrel{\→}{u}={[{\hat{a}}_1,\·\·\·,{\hat{a}}_n,{\hat{b}}_0,{\hat{b}}_1,\·\·\·,{\hat{b}}_n]}^T\·\·\·\left(10\right)$

将式(10)代入式(7)，令e(k)＝0，并改写为：

$u_i\left(k\right)=\frac{1}{{\hat{b}}_0}[u_o\left(k\right)+{\hat{a}}_1{u}_o(k-1)+\·\·\·+{\hat{a}}_n{u}_o(k-n)-{\hat{b}}_1\·u_i(k-1)-\·\·\·-{\hat{b}}_n\·u_i(k-n)]\·\·\·\left(11\right)$

根据式(11)，可由u_o(k)、u_o(k-1)、…、u_i(k-1)、…，得出u_i(k)。

2)动态重量变化率预测：

实际物料落入称重料斗的重量信号是u_i(kΔt)，u_i(kΔt)的变化是不均匀的，Δt是采样间隔时间，kΔt＝k_jΔt处是送料装置关断点，但由于物料送料机构的惯性和空中物料的存在，实际关断点的重量与稳态重量存在Δu_i的偏差；最大重量大于稳态重量，是因为物料自由落体的瞬时冲击力造成，从重量变化看，物料重量信号u_i(kΔt)的变化均匀性用u_i(kΔt)的变化率du_i/dt衡量，du_i/dt为大于零的常数时，说明物料重量是匀速上升，而实际表明du_i/dt不是常数，即送料速度不均匀是一定存在的；且du_i/dt大，则将造成实际关断点的重量与稳态重量偏差值Δu_i大；du_i/dt小，则将造成实际关断点的重量与稳态重量偏差值Δu_i小；因此，物料动态重量的变化率du_i/dt与偏差值Δu_i存在对应关系，根据实测实验表明，偏差值Δu_i与变化率du_i/dt为线性关系：

${\mathrm{\Δu}}_i=K\frac{{\mathrm{du}}_i}{\mathrm{dt}}\·\·\·\left(12\right)$

式(12)中K是经实验确定的常数，其量纲是S(秒)；

对于实际离散测量系统，采用

$\frac{{\mathrm{du}}_i}{\mathrm{dt}}=\frac{u_i\left(k_j\mathrm{\Δt}\right)-u_i\left(k_h\mathrm{\Δt}\right)}{(k_j-k_h)\mathrm{\Δt}}\·\·\·\left(13\right)$

计算du_i/dt，式(13)中k_j＞k_h，k_j是到实际关断点时的第k_j次采样，若k_j-k_h＝m，则k_h是在k_j前早m次的采样时刻。将式(13)代入式(12)可估算出将出现的偏差值Δu_i，或者在处k_j可更精确地估算出稳态的重量值，从而提高动态定量称重的精度。

本发明具有的有益效果是：

本发明涉及运用嵌入式(DSP)计算机这类哈佛总线的新型微处理机，构成测量和控制系统，采用模型化测量和动态重量变化率预测相结合的动态信号处理方式，以解决动态称重中既高精度又高速的问题。

附图说明

图1是测量系统。

图2是测量系统阶跃响应曲线。

图3是料斗中物料重量信号变化曲线。

图4是基于模型化测量和DSP技术的高速动态粉状物流定量称重系统示意框图。

具体实施方式

下面结合附图和实施例对本发明作进一步的说明。

1.1模型化测量(根据测量系统模型估算动态输入值)

由于本测量系统(包括传感器)是以嵌入式(DSP)计算机为核心的数字化测量系统，因此测量系统是一个离散系统，设测量系统的输入即被测量为ui(k)，测量系统的输出即测量输出值为u_o(k)。u_i(k)、u_o(k)以及与测量系统G(z)的关系如图1所示，图1中u_i(k)是固体粉状或颗粒状物料落入称重料斗的变化质量，u_o(k)是测量系统的测量输出值，当测量系统是理想系统(仅为比例系统)时，u_i(k)、u_o(k)的动态关系可以描述为

u_o(k)＝K·u_i(k)......................(1)

式(1)中K是常数。而实际的测量系统是一个线性定常离散动态系统，其差分方程的可描述为

a₀u_o(k)+a₁u_o(k-1)+…+a_nu_o(k-n)＝b₀·u_i(k)+b₁·u_i(k-1)+…+b_m·u_i(k-m)

(m≤n)...........................(2)

实际测量时直接获得的是u_o(k)，而不是u_i(k)，这对于u_i(k)发生瞬间变化时u_o(k)不能及时反应，从而使测量系统在u_i(k)发生瞬变过程时不能通过观察u_o(k)获得u_i(k)的变化。

如某测量系统的阶跃响应如图2所示，在0-1000ms间测量系统输出不能紧跟阶跃输入，此时的测量输出值不能反应输入的变化。

对式(2)的差分方程进行Z变换可得线性定常离散系统的脉冲传递函数描述为

$G\left(z\right)=\frac{b_0+b_1{z}^{-1}+\·\·\·+b_m{z}^{-m}}{a_0+a_1{z}^{-1}+\·\·\·+a_n{z}^{-n}},(m\≤n)\·\·\·\left(3\right)$

测量系统的脉冲传递函数G(z)又是测量系统脉冲响应序列g(k)的Z变换

G(z)＝Z[g(k)]............................(4)

当G(z)已知时，即a_o、a₁、…、a_n和b_o、b₁、…、b_m已知，u_i(k)、u_o(k)的动态对应关系就明确了，也就有可能通过测量输出uo(k)推算出u_i(k)。然而，通常测量系统的a_o、a₁、…、a_n和b_o、b₁、…、b_m是不知道的。

为了获得式(2)差分方程，可采用阶跃输入方法，观察并获取输出信号，将获取的输出信号和阶跃输入信号代入式(2)，采用最小二乘法可估算出a_o、a₁、…、a_n和b_o、b₁、…、b_n。将式(2)写成算子形式，并令n＝m、a_o＝1，则：

A(q^-1)u_o(k)＝B(q^-1)·u_i(k)....................(5)

式(5)中

A(q^-1)＝1+a₁q^-1+a₂q^-2+…+a_nq^-n

B(q^-1)＝b₀+b₁q^-1+b₂q^-2+…+b_nq^-n

考虑到模型的误差和测量噪声，式(5)方程改写为

A(q^-1)u_o(k)＝B(q^-1)·u_i(k)+e(k)....................(6)

式(6)中e(k)称为模型残差或方程误差。

将式(6)改写成

u_o(k)＝-a₁u_o(k-1)-…-a_nu_o(k-n)+b₀·u_i(k)+b₁·u_i(k-1)+…+b_n·u_i(k-n)+e(k)

...........................(7)

定义

$\stackrel{\→}{x}\left(k\right)={[{-u}_o(k-1),\·\·\·,{-u}_o(k-n),u_i\left(k\right),u_i(k-1),\·\·\·,u_i(k-n)]}^T$

$\stackrel{\→}{\mathrm{\θ}}={[a_1,\·\·\·,a_n,b_0,b_1,\·\·\·,b_n]}^T$

式(7)可写为

$u_o\left(k\right)={\stackrel{\→}{x}}^T\left(k\right)\·\stackrel{\→}{\mathrm{\θ}}+e\left(k\right)\·\·\·\left(8\right)$

将测得的(N+n)个数据代入式(8)，可以建立N个方程组，则可写成矩阵方程的形式

${\stackrel{\→}{u}}_o=\stackrel{\‾}{X}\·\stackrel{\→}{\mathrm{\θ}}+\stackrel{\→}{e}\·\·\·\left(9\right)$

其中

${\stackrel{\→}{u}}_o=\left[\begin{array}{c}{u}_o(n+1)\\ u_o(n+2)\\ \·\\ \·\\ \·\\ u_o(n+N)\end{array}\right],$ $\stackrel{\→}{e}=\left[\begin{array}{c}e(n+1)\\ e(n+2)\\ \·\\ \·\\ \·\\ e(n+N)\end{array}\right],$

$\stackrel{\‾}{X}=\left[\begin{array}{c}{\stackrel{\→}{x}}^T(n+1)\\ {\stackrel{\→}{x}}^T(n+2)\\ \·\\ \·\\ \·\\ \stackrel{\→}{x}_T(n+N)\end{array}\right]=\left[\begin{array}{cccccc}-u_o\left(n\right)& \·\·\·& -u_o\left(1\right)& u_i(n+1)& \·\·\·& u_i\left(1\right)\\ -u_o(n+1)& \·\·\·& -u_o\left(2\right)& u_i(n+2)& \·\·\·& u_i\left(2\right)\\ \·& & \·& \·& & \·\\ \·& & \·& \·& & \·\\ \·& & \·& \·& & \·\\ -u_o(n+N-1)& \·\·\·& -u_o\left(N\right)& u_i(n+N)& \·\·\·& u_i\left(N\right)\end{array}\right]$

根据最小二乘法，选择

使得误差准则函数J为最小：

$J=\underset{k=n+1}{\overset{N+n}{\mathrm{\Σ}}}{e}^2\left(k\right)={\stackrel{\→}{e}}^T\stackrel{\→}{e}={(\stackrel{\→}{u}-\stackrel{\‾}{X}\·\stackrel{\→}{\mathrm{\θ}})}^T(\stackrel{\→}{u}-\stackrel{\‾}{X}\·\stackrel{\→}{\mathrm{\θ}})$

令 $\begin{array}{c}{\frac{\∂J}{\∂\stackrel{\→}{\mathrm{\θ}}}|}_{\stackrel{\→}{\mathrm{\θ}}=\widehat{\stackrel{\→}{\mathrm{\θ}}}}=0\\ {\frac{\∂J}{\∂\stackrel{\→}{\mathrm{\θ}}}|}_{\stackrel{\→}{\mathrm{\θ}}=\widehat{\stackrel{\→}{\mathrm{\θ}}}}=\left[\frac{\∂J}{\∂\stackrel{\→}{\mathrm{\θ}}}{(\stackrel{\→}{u}-\stackrel{\‾}{X}\·\stackrel{\→}{\mathrm{\θ}})}^T\right](\stackrel{\→}{u}-\stackrel{\‾}{X}\·\stackrel{\→}{\mathrm{\θ}})+{(\stackrel{\→}{u}-\stackrel{\‾}{X}\·\stackrel{\→}{\mathrm{\θ}})}^T\left[\frac{\∂J}{\∂\stackrel{\→}{\mathrm{\θ}}}(\stackrel{\→}{u}-\stackrel{\‾}{X}\·\stackrel{\→}{\mathrm{\θ}})\right]\end{array}$

则 $0=-2{\stackrel{\‾}{X}}^T\·\stackrel{\→}{u}+2{\stackrel{\‾}{X}}^T\·\stackrel{\‾}{X}\·\stackrel{\→}{\mathrm{\θ}}$

求得

的估计值

$\widehat{\stackrel{\→}{\mathrm{\θ}}}={({\stackrel{\‾}{X}}^T\·\stackrel{\‾}{X})}^{-1}{\stackrel{\‾}{X}}^T\·\stackrel{\→}{u}={[{\hat{a}}_1,\·\·\·,{\hat{a}}_n,{\hat{b}}_0,{\hat{b}}_1,\·\·\·,{\hat{b}}_n]}^T\·\·\·\left(10\right)$

将式(10)代入式(7)，令e(k)＝0，并改写为：

$u_i\left(k\right)=\frac{1}{{\hat{b}}_0}[u_o\left(k\right)+{\hat{a}}_1{u}_o(k-1)+\·\·\·+{\hat{a}}_n{u}_o(k-n)-{\hat{b}}_1\·u_i(k-1)-\·\·\·-{\hat{b}}_n\·u_i(k-n)]\·\·\·\left(11\right)$

根据式(11)，可由u_o(k)、u_o(k-1)、…、u_i(k-1)、…，得出u_i(k)。

综上所述，本发明中的模型化测量方法，就是先用已知的阶跃输入信号测出测量系统的输出响应值，再用最小二乘法辨识测量系统的差分方程(如式(2)所示)，在实际测量时用辨识估计值形成式(11)，从而实现对u_i(k)的及时测量。

1.2根据动态重量变化率预测稳态重量

实际粉状或颗粒状物料落入料斗的过程，使称重传感器获取的、与物料重量成正比的称重信号值不断增大，若以当时测得的信号作为称重定量的控制，则实际入袋包装的物料将多于欲实现的定量值，目前一般采用估计一个提前量，提前将送料装置电源关断，停止送料，从而使定量称重保证一定的精度。这种提前停止送料方式在动态物流速度较慢时，且送料速度均匀，是可以保证精度的，但高速送料和送料速度不均匀时，定量称重精度就不能保证。

本发明提出一种动态重量变化率预测方法可以实现在高速送料和送料速度不均匀时的高精度定量称重。实际物料落入称重料斗的重量信号u_i(kΔt)的曲线如图3所示，从图3可见物料重量信号u_i(kΔt)的变化是不均匀的，Δt是采样间隔时间。图中kΔt＝k_jΔt处是送料装置关断点，但由于物料送料机构的惯性和空中物料的存在，实际关断点的重量与稳态重量存在Δu_i的偏差；曲线的最大重量大于稳态重量，是因为物料自由落体的瞬时冲击力造成。从重量变化曲线看，物料重量信号u_i(kΔt)的变化均匀性可用u_i(kΔt)的变化率du_i/dt衡量，du_i/dt为大于零的常数时，说明物料重量是匀速上升，而实际表明du_i/dt不可能是常数，即送料速度不均匀是一定存在的；且du_i/dt大，则将造成实际关断点的重量与稳态重量偏差值Δu_i大；du_i/dt小，则将造成实际关断点的重量与稳态重量偏差值Δu_i小。因此，物料动态重量的变化率du_i/dt与偏差值Δu_i存在对应关系，根据的实测实验表明，偏差值Δu_i与变化率du_i/dt为线性关系：

${\mathrm{\Δu}}_i=K\frac{{\mathrm{du}}_i}{\mathrm{dt}}\·\·\·\left(12\right)$

式(12)中K是经实验确定的常数，其量纲是S(秒)。

对于实际离散测量系统，采用

$\frac{{\mathrm{du}}_i}{\mathrm{dt}}=\frac{u_i\left(k_j\mathrm{\Δt}\right)-u_i\left(k_h\mathrm{\Δt}\right)}{(k_j-k_h)\mathrm{\Δt}}\·\·\·\left(13\right)$

计算du_i/dt，式(13)中k_j＞k_h，k_j是到实际关断点时的第k_j次采样，若k_j-k_h＝m，则k_h是在k_j前早m次的采样时刻。将式(13)代入式(12)可估算出将出现的偏差值Δu_i，或者在处k_j可更精确地估算出稳态的重量值，从而提高动态定量称重的精度。

2、测量装置：

运用嵌入式(DSP)计算机这类哈佛总线的新型微处理机为核心，构成测量和控制系统。图4是本发明的测量系统装置的构成示意图。

本系统装置工作时，粉状或颗粒状物流从散装到定量包装过程的机械装置由上料提升机、振动喂料机、称重料斗、下料包装机构成；系统的信号获取、信号处理、模型化测量和重量预测、料斗门控制、上料控制、振动控制、包装机控制，由传感器、模拟信号处理电路、嵌入式(DSP)计算机为核心的数据处理和控制系统、各接口电路和控制电路构成的电子电气系统实现。系统工作时，由“嵌入式(DSP)计算机为核心的数据处理和控制系统”发出控制信号，启动上料提升机、振动喂料机和下料包装机，待称重包装的粉状或颗粒状物料经皮带式上料提升机送到振动喂料机的入口，在振动喂料机一定振动幅度的控制下，物料按一定的速度向称重料斗送料，料斗下部的传感器将物料的重力和冲击力的合力信号转变为电压信号，该电压信号与物料重力和冲击力相关，电压信号经模拟信号处理电路处理后送到含有A/D电路的“嵌入式(DSP)计算机为核心的数据处理和控制系统”，“嵌入式(DSP)计算机为核心的数据处理和控制系统”接到信号后先进行滤波处理，然后进行基于测量系统模型的模型化测量处理和基于重力变化率的预测，若预测的值与设定值相等时，“嵌入式(DSP)计算机为核心的数据处理和控制系统”发出使振动喂料机停止运行的信号，并发出称重料斗门开启的信号，物料从称重料斗进入包装机并由包装机自动完成装袋封口的工作，至此，物料定量称重和包装完成。可见，该称重系统可用于食品加工、轻化工、粮食加工、饲料加工、制药、冶金等行业中的粉状和颗粒状产品自动定量称重和包装的全过程。

Claims (1)

1、一种基于模型化测量和动态重量变化率预测的动态称重方法，其特征在于该方法的步骤如下：

1)模型化测量：

设入测量系统的输入即被测量为u_i(k)，测量系统的输出即测量输出值为u_o(k)；u_i(k)是固体粉状或颗粒状物料落入称重料斗的变化质量，u_o(k)是测量系统的测量输出值，当测量系统是理想系统时，u_i(k)、u_o(k)的动态关系描述为

u_o(k)＝K·u_i(k)......................(1)

式(1)中K是常数。而实际的测量系统是一个线性定常离散动态系统，其差分方程的描述为

a₀u_o(k)+a₁u_o(k-1)+…+a_nu_o(k-n)＝b₀·u_i(k)+b₁·u_i(k-1)+…+b_m·u_i(k-m)

(m≤n)...........................(2)

为了获得式(2)差分方程，采用阶跃输入方法，观察并获取输出信号，将获取的输出信号和阶跃输入信号代入式(2)，采用最小二乘法估算出a_o、a₁、…、a_n和b_o、b₁、…、b_n，将式(2)写成算子形式，并令n＝m、a_o＝1，则：

A(q^-1)u_o(k)＝B(q^-1)·u_i(k)....................(5)

式(5)中

A(q^-1)＝1+a₁q^-1+a₂q^-2+…+a_nq^-n

B(q^-1)＝b₀+b₁q^-1+b₂q^-2+…+b_nq^-n

考虑到模型的误差和测量噪声，式(5)方程改写为

A(q^-1)u_o(k)＝B(q^-1)·u_i(k)+e(k)....................(6)

式(6)中e(k)称为模型残差或方程误差，将式(6)改写成

u_o(k)＝-a₁u_o(k-1)-…-a_nu_o(k-n)+b₀·u_i(k)+b₁·u_i(k-1)+…+b_n·u_i(k-n)+e(k)

...........................(7)

定义

$\stackrel{\→}{x}\left(k\right)={[{-u}_o(k-1),\·\·\·,{-u}_o(k-n),u_i\left(k\right),u_i(k-1),\·\·\·,u_i(k-n)]}^T$

$\stackrel{\→}{\mathrm{\θ}}={[a_1,\·\·\·,a_n,b_0,b_1,\·\·\·,b_n]}^T$

式(7)可写为

$u_o\left(k\right)={\stackrel{\→}{x}}^T\left(k\right)\·\stackrel{\→}{\mathrm{\θ}}+e\left(k\right)\·\·\·\left(8\right)$

将测得的(N+n)个数据代入式(8)，建立N个方程组，则写成矩阵方程的形式

${\stackrel{\→}{u}}_o=\stackrel{\‾}{X}\·\stackrel{\→}{\mathrm{\θ}}+\stackrel{\→}{e}\·\·\·\left(9\right)$

其中

${\stackrel{\→}{u}}_o=\left[\begin{array}{c}{u}_o(n+1)\\ u_o(n+2)\\ \·\\ \·\\ \·\\ u_o(n+N)\end{array}\right],$ $\stackrel{\→}{e}=\left[\begin{array}{c}e(n+1)\\ e(n+2)\\ \·\\ \·\\ \·\\ e(n+N)\end{array}\right],$

$\stackrel{\‾}{X}=\left[\begin{array}{c}{\stackrel{\→}{x}}^T(n+1)\\ \stackrel{\→}{x}_T(n+2)\\ \·\\ \·\\ \·\\ {\stackrel{\→}{x}}^T(n+N)\end{array}\right]=\left[\begin{array}{cccccc}-u_o\left(n\right)& \·\·\·& -u_o\left(1\right)& u_i(n+1)& \·\·\·& u_i\left(1\right)\\ -u_o(n+1)& \·\·\·& -u_o\left(2\right)& u_i(n+2)& \·\·\·& u_i\left(2\right)\\ \·& & \·& \·& & \·\\ \·& & \·& \·& & \·\\ \·& & \·& \·& & \·\\ -u_o(n+N-1)& \·\·\·& -u_o\left(N\right)& u_i(n+N)& \·\·\·& u_i\left(N\right)\end{array}\right]$

根据最小二乘法，选择

使得误差准则函数J为最小：

$J=\underset{k=n+1}{\overset{N+n}{\mathrm{\Σ}}}{e}^2\left(k\right)={\stackrel{\→}{e}}^T\stackrel{\→}{e}={(\stackrel{\→}{u}-\stackrel{\‾}{X}\·\stackrel{\→}{\mathrm{\θ}})}^T(\stackrel{\→}{u}-\stackrel{\‾}{X}\·\stackrel{\→}{\mathrm{\θ}})$

令 $\begin{array}{c}{\frac{\∂J}{\∂\stackrel{\→}{\mathrm{\θ}}}|}_{\stackrel{\→}{\mathrm{\θ}}=\widehat{\stackrel{\→}{\mathrm{\θ}}}}=0\\ {\frac{\∂J}{\∂\stackrel{\→}{\mathrm{\θ}}}|}_{\stackrel{\→}{\mathrm{\θ}}=\widehat{\stackrel{\→}{\mathrm{\θ}}}}=\left[\frac{\∂J}{\∂\stackrel{\→}{\mathrm{\θ}}}{(\stackrel{\→}{u}-\stackrel{\‾}{X}\·\stackrel{\→}{\mathrm{\θ}})}^T\right](\stackrel{\→}{u}-\stackrel{\‾}{X}\·\stackrel{\→}{\mathrm{\θ}})+{(\stackrel{\→}{u}-\stackrel{\‾}{X}\·\stackrel{\→}{\mathrm{\θ}})}^T\left[\frac{\∂J}{\∂\stackrel{\→}{\mathrm{\θ}}}(\stackrel{\→}{u}-\stackrel{\‾}{X}\·\stackrel{\→}{\mathrm{\θ}})\right]\end{array}$

则 $0=-2{\stackrel{\‾}{X}}^T\·\stackrel{\→}{u}+2{\stackrel{\‾}{X}}^T\·\stackrel{\‾}{X}\·\stackrel{\→}{\mathrm{\θ}}$

求得

的估计值

$\widehat{\stackrel{\→}{\mathrm{\θ}}}={({\stackrel{\‾}{X}}^T\·\stackrel{\‾}{X})}^{-1}{\stackrel{\‾}{X}}^T\·\stackrel{\→}{u}={[{\hat{a}}_1,\·\·\·,{\hat{a}}_n,{\hat{b}}_0,{\hat{b}}_1,\·\·\·,{\hat{b}}_n]}^T\·\·\·\left(10\right)$

将式(10)代入式(7)，令e(k)＝0，并改写为：

$u_i\left(k\right)=\frac{1}{{\hat{b}}_0}[u_o\left(k\right)+{\hat{a}}_1{u}_o(k-1)+\·\·\·+{\hat{a}}_n{u}_o(k-n)-{\hat{b}}_1\·u_i(k-1)-\·\·\·-{\hat{b}}_n\·u_i(k-n)]\·\·\·\left(11\right)$

根据式(11)，可由u_o(k)、u_o(k-1)、…、u_i(k-1)、…，得出u_i(k)。

2)动态重量变化率预测：

实际物料落入称重料斗的重量信号是u_i(kΔt)，u_i(kΔt)的变化是不均匀的，Δt是采样间隔时间，kΔt＝k_jΔt处是送料装置关断点，但由于物料送料机构的惯性和空中物料的存在，实际关断点的重量与稳态重量存在Δu_i的偏差；最大重量大于稳态重量，是因为物料自由落体的瞬时冲击力造成，从重量变化看，物料重量信号u_i(kΔt)的变化均匀性用u_i(kΔt)的变化率du_i/dt衡量，du_i/dt为大于零的常数时，说明物料重量是匀速上升，而实际表明du_i/dt不是常数，即送料速度不均匀是一定存在的；且du_i/dt大，则将造成实际关断点的重量与稳态重量偏差值Δu_i大；du_i/dt小，则将造成实际关断点的重量与稳态重量偏差值Δu_i小；因此，物料动态重量的变化率du_i/dt与偏差值Δu_i存在对应关系，根据实测实验表明，偏差值Δu_i与变化率du_i/dt为线性关系：

${\mathrm{\Δu}}_i=K\frac{{\mathrm{du}}_i}{\mathrm{dt}}\·\·\·\left(12\right)$

式(12)中K是经实验确定的常数，其量纲是S(秒)；

对于实际离散测量系统，采用

$\frac{{\mathrm{du}}_i}{\mathrm{dt}}=\frac{u_i\left(k_j\mathrm{\Δt}\right)-u_i\left(k_h\mathrm{\Δt}\right)}{(k_j-k_h)\mathrm{\Δt}}\·\·\·\left(13\right)$

计算du_i/dt，式(13)中k_j＞k_h，k_j是到实际关断点时的第k_j次采样，若k_j-k_h＝m，则k_h是在k_j前早m次的采样时刻。将式(13)代入式(12)可估算出将出现的偏差值Δu_i，或者在处k_j可更精确地估算出稳态的重量值，从而提高动态定量称重的精度。

Images