基于LVW-LM-QN算法的Wiener系统辨识

Identification of Wiener Systems Based on LVW-LM-QN Algorithm

摘要

Wiener系统由线性动态子系统与静态非线性子系统串联组成，广泛应用于石油、化工等过程工业中，获得Wiener系统的模型具有重要意义.本文针对Wiener系统提出一种基于线性变化权重-列文伯格马夸尔特-拟牛顿(linear variable weight-Levenberg Marquardt-quasi Newton，LVW-LM-QN)算法的非线性系统辨识方法.将Wiener系统分成两个子系统分别处理，对于线性动态部分，采用规范变量分析(canonical variate analysis，CVA)算法的子空间识别方法进行参数估计；对于非线性静态部分，采用LVW-LM-QN算法进行辨识处理.最后通过数值例子和双储罐系统液位控制的应用案例来评估该方法，仿真结果验证了所提方法的有效性和精确性.

Abstract

Wiener systems， consisting of a linear dynamic subsystem and a static nonlinear subsystem in series， find extensive application in process industries such as petroleum and chemical engineering. Obtaining the model of Wiener systems holds significant importance. A nonlinear system identification method based on the linear variable weight–Levenberg Marquardt–quasi Newton （LVW-LM-QN） algorithm for Wiener systems was proposed. The Wiener system was divided into two subsystems for separate processing. For the linear dynamic part， the subspace identification method with the canonical variate analysis （CVA） algorithm was used for parameter estimation， whereas for the subsequent nonlinear static part， the LVW-LM-QN algorithm was employed for identification. Finally， the method was evaluated through numerical examples and an application case of liquid level control in a two-tank system， and the effectiveness and accuracy of the proposed method were verified by the simulation results.

Graphical abstract

关键词

Wiener系统 / 子空间方法 / LVW-LM-QN算法 / 神经网络 / CVA算法

Key words

Wiener system / subspace method / LVW-LM-QN algorithm / neural network / CVA algorithm

系统辨识是根据收集到的输入和输出数据集建立数学模型.近几十年，线性系统的研究领域取得了一些重大进展.然而，非线性系统广泛存在于工业过程中，由于非线性系统的复杂性和多样性，建立全局模型有困难，因此非线性系统的辨识仍然具有挑战性.

面向块的模型，例如Wiener^［1］，Hammerstein，Hammerstein-Wiener^［2］等，用简单的块来构建模型，以便找到足够灵活的结构来表示许多实际非线性系统^［3-4］.本文考虑建立Wiener系统，它是面向块的非线性模型的典型块结构，其特点是线性动态部分后面是静态非线性块.

工业过程的许多非线性特征都可以用Wiener系统来描述，例如连续搅拌釜反应器^［5］、热交换系统^［6］等.Wiener系统的优点在于线性块包含了系统动力学的复杂性，而非线性的复杂性仅包含在静态块中，这一优点降低了辨识的难度，因为可以应用线性和非线性系统理论的分而治之策略来分别处理两个子系统^［7］.有许多方法来辨识Wiener系统，区别主要在于动态线性块和静态非线性块的表示方式.用来表示动态线性块的结构包括阶跃响应模型^［8-9］、外源输入自回归（ARX）模型^［8］、状态空间模型^［10］等.对非线性块提出的结构有分段线性函数^［10］、多项式函数^［11］、支持向量机^［12］等.

多项式函数经常被用来表示Wiener系统的非线性部分，但对于具有高度非线性增益的过程，它存在一定限制^［13］.还有许多性能更强大的非线性建模方法，例如在Wiener系统构造中引入B样条网络^［14］和神经网络^［15］.近年来，神经网络已成为分类和回归领域的常用方法^［16］，对于黑盒模型识别具有几个有吸引力的特性：强大的非线性拟合能力、泛化能力良好、可以并行计算、具有鲁棒性.Mohsen等^［17］应用神经网络将非线性行为过程的子空间识别方法扩展到Wiener结构.

本文采用CVA子空间方法辨识的状态空间方程与考虑偏置参数的神经网络相结合构建了Wiener系统，该模型是Mohsen等^［17］提出的Wiener系统的改进版本.并提出了LVW-LM-QN算法来识别Wiener系统的非线性静态部分，引入了线性变化权重，提高了算法效率和性能，数值例子以及应用案例的仿真结果验证了所提方法的有效性和优越性.

1 问题描述

Wiener系统由动态线性块串联静态非线性块组成，如图1所示.其中，u（k），y（k），v（k）分别表示模型的输入、输出和噪声信号；x（k）表示线性状态空间子系统的状态向量； A， B， C 和 D 是适当维度的系统矩阵； （k），s（k）是未知的中间变量；F（\cdot）表示系统的输出非线性函数.

系统的线性子系统以状态空间的形式描述为

输出非线性子系统表示为

s(k)=F( (k)).（2）

系统输出则为

y(k)=s(k)+v(k).（3）

对图1的Wiener系统，有以下假设.

假设1：线性子系统是渐进稳定的，即其特征值严格位于单位圆内.

假设2：线性子系统是最小实现的（即对（ A， C ）是能观的，对（ A， B ）是能控的）.

假设3：系统输入u(k)和系统输出y(k)是可用的持续激励信号^［18］.

假设4：扰动v(k)是均值为零、方差为{\sigma }^{\mathrm{2}}的高斯白噪声^［19］.

用多层感知器神经网络来表示F（\cdot），可以得到非线性动态子系统为神经网络形式的Wiener模型，如图2所示，其中非线性部分的 w₁（k）=［w₁₁（k）…w_1m （k）］^T， w₂（k）=［w₂₁（k）…w_2m （k）］^T表示神经网络的权重，b₁（k），b₂（k）∈\mathbf{R}表示神经网络的偏置，\Phi（\cdot）∈\mathbf{R}^m 表示隐藏层神经元的激活函数.

用该模型来辨识非线性动态系统，模型输出为

\widehat{y}(k)={\mathit{w}}_{\mathrm{2}}{(k)}^{\mathrm{T}}\Phi ({\mathit{w}}_{\mathrm{1}}(k)(k)+b_{\mathrm{1}}(k))+b_{\mathrm{2}}(k).（4）

辨识过程的成本函数定义为

V(k)=\frac{\mathrm{1}}{\mathrm{2}}f{(k)}^{\mathrm{2}}，（5）

f(k)=y(k)-\widehat{y}(k).（6）

f（k）表示系统输出y和Wiener系统输出\widehat{y}之间的误差.

本文的目标是使用所提出的LVW-LM-QN辨识方法最小化式（5），从而得到Wiener模型参数.

2 Wiener系统辨识算法

2.1 Wiener系统线性部分辨识

子空间方法在辨识中对模型结构先验知识需求较少，在数值计算中具有一定的鲁棒性，以及广泛适用于多变量系统，因此子空间方法得到了辨识和控制领域的广泛关注.

考虑到假设1~3，子空间识别方法利用测量的输入输出信号序列，确定式（1）表示的状态空间子系统的系统矩阵 A， B， C 和 D 以及系统阶数.

下列输入输出矩阵方程在子空间识别的过程中发挥了重要作用，它表示未来输出 Y_f块Hankel矩阵与未来输入块Hankel矩阵 U_f和未来状态序列 X_i 线性相关：

{\mathit{Y}}_{\mathrm{f}}={\mathit{\Gamma }}_i{\mathit{X}}_i+{\mathit{H}}_i{\mathit{U}}_{\mathrm{f}}.（7）

其中扩展可观测矩阵{\mathit{\Gamma }}_i：

{\mathit{\Gamma }}_i\stackrel{\mathrm{d}\mathrm{e}\mathrm{f}}{=}\left[\begin{array}{l}\mathit{C}\\ \mathit{C}\mathit{A}\\ \mathit{C}{\mathit{A}}^{\mathrm{2}}\\ ⋮\\ \mathit{C}{\mathit{A}}^{i-\mathrm{1}}\end{array}\right].(8)

下三角托普利茨矩阵 H_i 定义为

H_i\stackrel{\mathrm{d}\mathrm{e}\mathrm{f}}{=}\left[\begin{array}{ccccc}\mathit{D}& \mathrm{0}& \mathrm{0}& \dots & \mathrm{0}\\ \mathit{C}\mathit{B}& \mathit{D}& \mathrm{0}& \dots & \mathrm{0}\\ \mathit{C}\mathit{A}\mathit{B}& \mathit{C}\mathit{B}& \mathit{D}& \dots & \mathrm{0}\\ ⋮& ⋮& ⋮& & ⋮\\ \mathit{C}{\mathit{A}}^{i-\mathrm{2}}\mathit{B}& \mathit{C}{\mathit{A}}^{i-\mathrm{3}}\mathit{B}& \mathit{C}{\mathit{A}}^{i-\mathrm{4}}\mathit{B}& \dots & \mathrm{0}\end{array}\right].(9)

输入输出块Hankel矩阵定义为

{\mathit{U}}_{\mathrm{0}|i-\mathrm{1}}\stackrel{\mathrm{d}\mathrm{e}\mathrm{f}}{=}\left[\begin{array}{cccc}{u}_{\mathrm{0}}& u_{\mathrm{1}}& \dots & u_{j-\mathrm{1}}\\ u_{\mathrm{1}}& u_{\mathrm{2}}& \dots & u_j\\ ⋮& ⋮& & ⋮\\ u_{i-\mathrm{1}}& u_i& \dots & u_{i+j-\mathrm{2}}\end{array}\right]，(10)

{\mathit{Y}}_{\mathrm{0}|i-\mathrm{1}}\stackrel{\mathrm{d}\mathrm{e}\mathrm{f}}{=}\left[\begin{array}{cccc}{y}_{\mathrm{0}}& y_{\mathrm{1}}& \dots & y_{j-\mathrm{1}}\\ y_{\mathrm{1}}& y_{\mathrm{2}}& \dots & y_j\\ ⋮& ⋮& & ⋮\\ y_{i-\mathrm{1}}& y_i& \dots & y_{i+j-\mathrm{2}}\end{array}\right].(11)

简写为

\begin{array}{l}{\mathit{U}}_{\mathrm{p}}\stackrel{\mathrm{d}\mathrm{e}\mathrm{f}}{=}{\mathit{U}}_{\mathrm{0}|i-\mathrm{1}}，{\mathit{U}}_{\mathrm{f}}\stackrel{\mathrm{d}\mathrm{e}\mathrm{f}}{=}{\mathit{U}}_{i|\mathrm{2}i-\mathrm{1}}，\\ {\mathit{Y}}_{\mathrm{p}}\stackrel{\mathrm{d}\mathrm{e}\mathrm{f}}{=}{\mathit{Y}}_{\mathrm{0}|i-\mathrm{1}}，{\mathit{Y}}_{\mathrm{f}}\stackrel{\mathrm{d}\mathrm{e}\mathrm{f}}{=}{\mathit{Y}}_{i|\mathrm{2}i-\mathrm{1}}，\end{array}(12)

其中:下标p和f分别表示过去和未来.包含输入{\mathit{U}}_{\mathrm{p}}和输出{\mathit{Y}}_{\mathrm{p}}的矩阵称为{\mathit{W}}_{\mathrm{p}}:

{\mathit{W}}_{\mathrm{p}}\stackrel{\mathrm{d}\mathrm{e}\mathrm{f}}{=}\left[\begin{array}{l}{\mathit{Y}}_{\mathrm{p}}\\ {\mathit{U}}_{\mathrm{p}}\end{array}\right].(13)

将状态序列 X_i 表示为

{\mathit{X}}_i\stackrel{\mathrm{d}\mathrm{e}\mathrm{f}}{=}(x_i{x}_{i+\mathrm{1}}{x}_{i+\mathrm{2}}\dots x_{i+j-\mathrm{1}}).(14)

子空间识别算法始终包含两个步骤：第1步对先前定义的Hankel矩阵生成的行空间进行加权投影，以得到扩展可观测性矩阵{\mathit{\Gamma }}_i或未知系统状态序列{\mathit{X}}_i的估计{\tilde{\mathit{X}}}_i；第2步从扩展可观测性矩阵{\mathit{\Gamma }}_i或状态序列估计{\tilde{\mathit{X}}}_i中检索出系统矩阵 A， B， C， D .

使用几何工具如下：

运算符{\mathrm{\Pi }}_{\mathit{B}}表示将某一矩阵的行空间投影到矩阵 B\in {\mathbf{R}}^{q\times j}的行空间上：

{\mathrm{\Pi }}_{\mathit{B}}\stackrel{\mathrm{d}\mathrm{e}\mathrm{f}}{=}{\mathit{B}}^{\mathrm{T}}·(\mathit{B}{\mathit{B}}^{\mathrm{T}})⁺·\mathit{B}.(15)

A / B 表示将矩阵 A\in {\mathbf{R}}^{p\times j}的行空间投影到矩阵 B 的行空间：

A / B\stackrel{\mathrm{d}\mathrm{e}\mathrm{f}}{=}A ⋅Π _B = AB^T⋅( BB^T)⁺⋅ B.(16)

而{\mathrm{\Pi }}_{{\mathit{B}}^{\perp }}是将某一矩阵的行空间投影到矩阵 B 的行空间的正交补上的几何运算符：

A / B{}^{\perp }\stackrel{\mathrm{d}\mathrm{e}\mathrm{f}}{=}A ⋅{\mathrm{\Pi }}_{{\mathit{B}}^{\perp }}.(17)

{\mathrm{\Pi }}_{{\mathit{B}}^{\perp }}={\mathit{I}}_j-{\mathrm{\Pi }}_{\mathit{B}}.(18)

投影{\mathrm{\Pi }}_{\mathit{B}}和{\mathrm{\Pi }}_{{\mathit{B}}^{\perp }}的组合将矩阵 A 分解为行空间正交的2个矩阵：

A = A ⋅Π{}_{\mathit{B}}+ A ⋅Π{}_{{\mathit{B}}^{\perp }}= A / B + A / B{}^{\perp }. (19)

子空间识别算法的具体步骤：首先将 Y_f的行空间投影到 U_f行空间的正交补 U{}_{\mathrm{f}}^{\perp }上：

Y_f/ U{}_{\mathrm{f}}^{\perp }= Γ_iX_i / U{}_{\mathrm{f}}^{\perp }+ H_iU_f/ U{}_{\mathrm{f}}^{\perp }.(20)

由式（15）~式（18），有

U_f/ U{}_{\mathrm{f}}^{\perp }= U_f⋅Π{}_{{\mathit{U}}_{\mathrm{f}}^{\perp }}= U_f⋅（I_j-Π{}_{{\mathit{U}}_{\mathrm{f}}}）=

U_f- U_fU{}_{\mathrm{f}}^{\mathrm{T}}⋅（ U_fU{}_{\mathrm{f}}^{\mathrm{T}}）⁺⋅ U_f=0.(21)

于是 Y_f/ U{}_{\mathrm{f}}^{\perp }= Γ_iX_i / U{}_{\mathrm{f}}^{\perp }.

使用矩阵 W₁和 W₂对该投影向左和向右加权：输入 U_f和加权矩阵 W₁和 W₂不能任意选择，应满足以下2个条件：

rank( W₁⋅ Γ_i )=rank Γ_i，(22)

rank（ X_i / U{}_{\mathrm{f}}^{\perp }⋅ W₂）=rank X_i .(23)

保证 Γ_iX_i 在投影到 U{}_{\mathrm{f}}^{\perp }上且由矩阵 W₁， W₂加权后，其秩-n特性不变.如果满足上述两个条件，可得

O_i\stackrel{\mathrm{d}\mathrm{e}\mathrm{f}}{=}W₁⋅ Y_f/ U{}_{\mathrm{f}}^{\perp }⋅ W₂= W₁⋅ Γ_iX_i / U{}_{\mathrm{f}}^{\perp }⋅ W₂.(24)

通过SVD分解：

O_i =（ U₁U₂）\left(\begin{array}{cc}{\mathit{S}}_{\mathrm{1}}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}\end{array}\right)\left(\begin{array}{c}{\mathit{V}}_{\mathrm{1}}^{\mathrm{T}}\\ {\mathit{V}}_{\mathrm{2}}^{\mathrm{T}}\end{array}\right).(25)

得到以下重要属性：

rank O_i =n， Γ_i = W{}_{\mathrm{1}}^{-\mathrm{1}}U₁S{}_{\mathrm{1}}^{\mathrm{1}/\mathrm{2}}， X_i / U{}_{\mathrm{f}}^{\perp }⋅ W₂= S{}_{\mathrm{1}}^{\mathrm{1}/\mathrm{2}}V{}_{\mathrm{1}}^{\mathrm{T}} .(26)

显然，从矩阵 W₁⋅ Y_f / U{}_{\mathrm{f}}^{\perp }⋅ W₂的奇异值分解中可以获得系统阶数n.此外，从对应于非零奇异值的左奇异向量可以找到扩展的可观测性矩阵 Γ_i 的相似变换，而右奇异向量包含了状态序列 X_i 的相关信息.适当选择权重矩阵 W₂，则

{\tilde{\mathit{X}}}_i\stackrel{\mathrm{d}\mathrm{e}\mathrm{f}}{=}X_i / U{}_{\mathrm{f}}^{\perp }⋅ W₂(27)

可被视为状态序列 X_i 的估计.

只要选择合适的权重矩阵 W₁和 W₂，LTI（linear-time-invariant）系统的所有子空间算法都可以在上述框架中解释，包括N4SID，CVA，MOESP，IV-4SID和basic-4SID.

根据系统参数矩阵不同的估计方式，可以将系统参数矩阵的计算方法分为两类：一类是通过提取状态序列估计值{\tilde{\mathit{X}}}_i的列子空间检索系统矩阵，为回归法，如N4SID，CVA；另一类是通过提取扩展可观测矩阵 Γ_i 的行子空间检索系统矩阵，为实现法，如MOESP，IV-4SID，basic-4SID.

本文使用CVA算法来辨识线性动态系统，算法原理如下.

已知CVA算法的权重矩阵取值：

W₁=［（ Y_f/ U{}_{\mathrm{f}}^{\perp }）⋅（ Y_f/ U{}_{\mathrm{f}}^{\perp }）^T］^-1/2，(28)

W₂=（ W_P/ U{}_{\mathrm{f}}^{\perp }）⁺⋅（ W_P/ U{}_{\mathrm{f}}^{\perp }）.(29)

由式（20）~式（27）得到{\tilde{\mathit{X}}}_i，通过类似的推理，容易证明：

O_i+1\stackrel{\mathrm{d}\mathrm{e}\mathrm{f}}{=}{\overline{\mathit{W}}}_{\mathrm{1}}⋅ Y{}_{\mathrm{f}}^{-}/ U{}_{\mathrm{f}}^{-\perp }⋅\overline{\mathit{W}}{}_{\mathrm{2}}={\overline{\mathit{W}}}_{\mathrm{1}}⋅ Γ_i-1⋅{\tilde{\mathit{X}}}_{i+\mathrm{1}}⋅\overline{\mathit{W}}{}_{\mathrm{2}}.(30)

其中：\overline{\mathit{W}}=\mathit{W}o|i； U{}_{\mathrm{f}}^{-}=U_i+1|2i-1； Y{}_{\mathrm{f}}^{-}=Y_i+1|2i-1.

{\mathit{U}}_{i+\mathrm{1}|\mathrm{2}i-\mathrm{1}}\stackrel{\mathrm{d}\mathrm{e}\mathrm{f}}{=}\left(\begin{array}{cccc}{u}_{i+\mathrm{1}}& u_{i+\mathrm{2}}& \cdots & u_{i+j}\\ ⋮& ⋮& & ⋮\\ u_{\mathrm{2}i-\mathrm{1}}& u_{\mathrm{2}i}& \cdots & u_{\mathrm{2}i+j-\mathrm{2}}\end{array}\right)，(31)

{\mathit{Y}}_{i+\mathrm{1}|\mathrm{2}i-\mathrm{1}}\stackrel{\mathrm{d}\mathrm{e}\mathrm{f}}{=}\left(\begin{array}{cccc}{y}_{i+\mathrm{1}}& y_{i+\mathrm{2}}& \cdots & y_{i+j}\\ ⋮& ⋮& & ⋮\\ y_{\mathrm{2}i-\mathrm{1}}& y_{\mathrm{2}i}& \cdots & y_{\mathrm{2}i+j-\mathrm{2}}\end{array}\right).(32)

很容易得出，如果去掉由式（26）计算出的 Γ_i 的最后l（系统输出的数量）行，便得到 Γ_i-1：

Γ_i-1=\underline{{\mathit{\Gamma }}_i}.(33)

此时{\tilde{\mathit{X}}}_{i+\mathrm{1}}为

{\tilde{\mathit{X}}}_{i+\mathrm{1}}=({\overline{\mathit{W}}}_{\mathrm{1}}\cdot {\mathit{\Gamma }}_{i-\mathrm{1}})\times O_{i+\mathrm{1}}{{\overline{\mathit{W}}}_{\mathrm{2}}}^{+}.(34)

已经得到了{\tilde{\mathit{X}}}_i，{\tilde{\mathit{X}}}_{i+\mathrm{1}}，只使用系统的输入输出数据，可以通过最小二乘法求解一组简单的超定方程来得到状态空间矩阵 A， B， C 和 D ：

\left(\begin{array}{c}{\tilde{\mathit{X}}}_{i+\mathrm{1}}\\ {\mathit{Y}}_{i|i}\end{array}\right)=\left(\begin{array}{cc}\mathit{A}& \mathit{B}\\ \mathit{C}& \mathit{D}\end{array}\right)\left(\begin{array}{c}{\tilde{\mathit{X}}}_i\\ {\mathit{U}}_{i|i}\end{array}\right).(35)

在Wiener系统的辨识中，CVA算法利用系统的输入输出数据辨识出模型的系统矩阵 A， B， C 和 D，即模型的线性动态块.

2.2 Wiener系统非线性部分辨识

Gauss-Newton（GN）算法是解决非线性最优问题的常见算法之一^［20］，它不直接计算Hessian矩阵 H，而是通过Jacobi矩阵 J 对 H 进行近似（ H = J^TJ ）：

h_gn=-（ J （k）^TJ （k））^-1g（k），(36)

w_new=w+h_gn.(37)

其中：w为神经网络权重；h_gn为GN算法的迭代步长；k是当前迭代次数；g= J^Tf =\frac{\partial V}{\partial w}为成本函数V对w的梯度.Gauss-Newton算法减少了计算量和所需的存储空间，但如果 J 不是满秩矩阵，Gauss-Newton算法不收敛.

因此引入单位矩阵 H = J^TJ +\muI 来保证近似矩阵的可逆性，就得到了Levenberg Marquardt（LM）算法：

h_{\mathrm{l}\mathrm{m}}=-(\mathit{J}{(k)}^{\mathrm{T}}\mathit{J}{(k)+\mu \mathit{I})}^{-\mathrm{1}}g(k)，(38)

w_new=w+h_lm.(39)

其中：正则化参数\mu \ge \mathrm{0}；h_lm为LM算法的迭代步长.

假设\stackrel{\mathrm{*}}{w}为局部最小值点，如果f（k）在这一点等于0，LM算法在\stackrel{\mathrm{*}}{w}处是二倍收敛^［21］，否则为线性收敛^［21］.

为了得到更好的性能，Loke等^［22］将LM算法与QN算法相结合，提出了一种混合算法（LM-QN），在f（k）于\stackrel{\mathrm{*}}{w}处≠0的条件下仍是超线性收敛^［21］.QN算法基于矩阵 B 对矩阵 H 进行近似：

h_qn=-\mathit{B}{(k)}^{-\mathrm{1}}g(k)，(40)

w_{\mathrm{n}\mathrm{e}\mathrm{w}}=w+h_{\mathrm{q}\mathrm{n}}.(41)

其中： B 为Hessian矩阵 H 的近似矩阵；h_{\mathrm{q}\mathrm{n}}为QN算法的迭代步长.

基本LM-QN算法的学习率是一个常数，通常取1.然而对于一些复杂的问题，由于存在许多局部解，LM-QN算法可能在寻找全局最优解时性能不佳，过早陷入局部最优解.早熟收敛降低了优化效果，为了避免这一缺陷，本文引入线性变化权重LVW（linear variable weight）^［23］策略，应用于学习率\alpha：

\alpha (k)={\alpha }_{\mathrm{m}\mathrm{a}\mathrm{x}}-\frac{{\alpha }_{\mathrm{m}\mathrm{a}\mathrm{x}}-{\alpha }_{\mathrm{m}\mathrm{i}\mathrm{n}}}{k_{\mathrm{m}\mathrm{a}\mathrm{x}}}k，(42)

\mathrm{L}\mathrm{M}:w_{\mathrm{n}\mathrm{e}\mathrm{w}}=w+\alpha (k)h_{\mathrm{l}\mathrm{m}}，(43)

\mathrm{Q}\mathrm{N}:w_{\mathrm{n}\mathrm{e}\mathrm{w}}=w+\alpha (k)h_{\mathrm{q}\mathrm{n}}.(44)

其中：k_{\mathrm{m}\mathrm{a}\mathrm{x}}表示最大迭代次数；{\alpha }_{\mathrm{m}\mathrm{a}\mathrm{x}}和{\alpha }_{\mathrm{m}\mathrm{i}\mathrm{n}}分别表示学习率的最大值和最小值.

LVW-LM-QN辨识算法的迭代过程如下：

LM-QN算法的权重w在第1次迭代时，算法默认设置为LM，根据式（38）计算出迭代步长h_{\mathrm{l}\mathrm{m}}，再通过式（43）得到w_{\mathrm{n}\mathrm{e}\mathrm{w}}，此时根据增益率公式计算出\rho值：

\rho =(V(k)-V(k+\mathrm{1}))/(L(\mathrm{0})-L(h_{\mathrm{l}\mathrm{m}})).(45)

该式表示函数值的实际下降与预测下降之间的比率，L为在当前迭代值w的邻域对V值的近似模型：

V(k+\mathrm{1})\simeq L(h)\equiv V(k)+{\mathit{h}}^{\mathrm{T}}g+\frac{\mathrm{1}}{\mathrm{2}}{h}^{\mathrm{T}}\mathit{B}\mathit{h}.(46)

如果\rho值大于0，说明此次迭代使得V减小，于是w更新为w_new.否则，不会进行w值的更新.再根据式（47）进行判断：

.(47)

如果式（47）连续满足3次，表明获取的权重越来越接近w⁺ （全局最小值即w⁺=argmin _w\left\{F(w)\right\}），其中g（w⁺）=0且V（w⁺）显著非0，这可能导致缓慢的线性收敛^［21］，优化效果差，于是会在下一次迭代时选择算法为QN.

算法为QN时，根据式（40），式（44）计算得到w_new，此时将检验是否满足：

，(48)

(V(k+\mathrm{1})\le (\mathrm{1}+\delta )V(k))，(49)

.(50)

其中\delta为浮点数的相对精度，值为2.220 4e-16.式（48）的满足表明成本函数V减小；式（49）的满足表明此时允许V值略微增加，因为QN算法的主要作用在于使g(k)更接近于0；式（50）的满足说明g(k)的范数值在快速下降.当式（48）和式（49）其中之一满足，同时式（50）也满足时，迭代算法具有超线性收敛的良好性能，于是w更新为w_new，否则不更新.若

{‖g(k+\mathrm{1})‖}_{\mathrm{\infty }}\ge {‖g(k)‖}_{\mathrm{\infty }}，(51)

说明V的梯度的范数下降得不够快，此步迭代QN算法作用不明显，收敛性能不佳，于是将在下一次迭代时再切换算法为LM.

该算法可以为解决许多非线性问题提供一种高效且具有前瞻性的方案，也为寻找最优解提供了新的可能性.

1） 所有非线性优化方法都是迭代的：从起点 w₀（即权重初始值）开始生成一系列向量 w₁， w₂，\dots，逐渐收敛到\stackrel{\mathrm{*}}{w}，即给定函数的局部最小值.当 w₀远离\stackrel{\mathrm{*}}{\mathit{w}}时，希望该方法向\stackrel{\mathrm{*}}{\mathit{w}}稳步迭代，在迭代的全局阶段，仅要求除了最初的迭代步骤之外e_k 没有增加，即，其中e_k={\mathit{w}}_k-\stackrel{\mathrm{*}}{\mathit{w}}.

但在迭代的最后阶段，{\mathit{w}}_k越来越接近\stackrel{\mathrm{*}}{\mathit{w}}，希望得到更快的收敛速度，于是有

线性收敛：

(52)

超线性收敛：

‖e_{k+\mathrm{1}}‖/‖e_k‖\to \mathrm{0}，\mathrm{其中}k\to \mathrm{\infty }.(53)

二倍收敛：

‖e_{k+\mathrm{1}}‖=O\left({‖e_k‖}^{\mathrm{2}}\right)，\mathrm{其中}k\to \mathrm{\infty }.(54)

2） 本文使用拟牛顿法更新矩阵 B，其数值效果好，并且具有全局收敛性和超线性收敛速度：

h（k）=w_new-w，(55)

\mathit{s}(k)=\mathit{J}(k+\mathrm{1})\mathit{J}(k)\mathit{h}(k)+(\mathit{J}(k+\mathrm{1})-\mathit{J}{(k))}^{\mathrm{T}}f(k+\mathrm{1}).(56)

如果 h （k）^Ts （k）>0，则对矩阵 B 作以下更新：

\mathit{v}(k)=\mathit{B}(k)\mathit{h}(k).(57)

\begin{array}{l}\mathit{B}(k+\mathrm{1})=\mathit{B}(k)+\left(\frac{\mathit{s}(k)}{\mathit{h}{(k)}^{\mathrm{T}}\mathit{s}(k)}\right)\mathit{s}{(k)}^{\mathrm{T}}-\\ \left(\frac{\mathit{v}(k)}{\mathit{h}{(k)}^{\mathrm{T}}\mathit{v}(k)}\right)\mathit{v}{(k)}^{\mathrm{T}}.\end{array}(58)

为简便表示，LVW-LM-QN算法的伪代码如表1所示.

所提辨识方法框图如图3所示.辨识过程可概括如下：子空间方法以状态空间形式来估计实际系统输入输出数据的线性动态近似，该初始估计构成了Wiener系统的线性动态子系统（LDS）.前馈神经网络构成非线性静态子系统（NSS），其权重使用LDS的输出数据和系统输入输出数据进行迭代训练.LDS的参数是使用CVA算法确定的，而NSS的参数是使用LVW-LM-QN算法获得的，从而获得Wiener整体模型参数.该方法的主要优点之一是只需要系统的输入输出数据.

3 仿真实验

3.1 数值示例

为了评估所提出的算法，考虑以下带有色噪声干扰的Wiener系统：

x(k+\mathrm{1})=\left[\begin{array}{cc}-\mathrm{1.5}& -\mathrm{3}\\ \mathrm{3}& -\mathrm{1}\end{array}\right]x(k)+\left[\begin{array}{c}\mathrm{1.3}\\ \mathrm{0}\end{array}\right]u(k)，

y（k）=ξ（k）+v（k）.

仿真时，输入u（k）采用单位方差的不相关零均值随机信号序列，v（k）采用方差σ²=0.01的零均值高斯噪声序列.通过使用Matlab软件仿真，收集了2 400对输入输出样本的数据集（见图4）.

使用前600对数据作为训练集，剩余的1 800对数据作为验证集.对于线性动态子系统（LDS）的辨识使用CVA算法^［24］，CVA算法先对系统训练集的输入输出块Hankel矩阵生成的行空间进行加权投影，根据式（20）~式（28）得到未知系统状态序列{\mathit{X}}_i的估计{\tilde{\mathit{X}}}_i，再根据加权投影（式（30）），由式（28）~式（35）得到{\tilde{\mathit{X}}}_{i+\mathrm{1}}，此时求解最小二乘问题即可得到状态空间矩阵.

训练集的输入数据经过CVA辨识出的系统矩阵 A， B， C， D 构成的线性动态块，得到中间变量 (k)，即非线性静态块的输入.NSS块中神经元的激活函数选择为\frac{\mathrm{1}}{\mathrm{1}+{\mathrm{e}}^{-\mathrm{1.5}x}}，突触权重向量（ w₁， w₂）和偏置参数（b₁，b₂）在训练开始时设定为［0，1］之间的随机数.LVW-LM-QN算法使用 (k)和系统输出数据，通过反馈回来的误差f(k)进行迭代训练，目的是最小化成本函数V，最终得到构成NSS的多层感知器神经网络的权重w.只需要训练集的输入输出数据，便得到了Wiener模型在训练集上的线性块和非线性块的参数.

图5为本文所提出的LVW-LM-QN方法、Mohsen等^［17］提出的LM方法以及Gómez等^［25］提出的Subspace-based方法在验证集的输出.如图5所示，本文方法优于LM以及Subspace-based 方法.为了评估所提出方法的准确性，引入了均方根误差（RMSE）和平均绝对百分比误差（MAPE）：

\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}=\sqrt[]{\frac{\mathrm{1}}{N}{\displaystyle \sum _{i=\mathrm{1}}^N}(y_i-{\widehat{y}}_i{)}^{\mathrm{2}}},(59)

\mathrm{M}\mathrm{A}\mathrm{P}\mathrm{E}=\frac{\mathrm{1}}{N}{\displaystyle \sum _{i=\mathrm{1}}^N}\left|\frac{{\widehat{y}}_i-y_i}{y_i}\right|.(60)

其中：N为样本总数；y_i和{\widehat{y}}_i分别是第i个样本的实际值和模型估计值.RMSE和MAPE越小说明辨识效果越好.辨识结果与实际值的均方根误差和平均绝对百分比误差如表2所示.

从表2可以看出，相比其他两种方法，LVW-LM-QN方法的RMSE和MAPE最小，因此，所提方法在数值例子中具有优越的辨识性能.

3.2 双储罐系统的液位控制

双储罐系统中的液位控制^［26］是大多数过程工业的组成部分，例如化工厂、炼油厂、发电厂等.液位控制的主要目的是根据施加到泵上的电压的变化，控制罐中所需的液体液位.在此示例中，采样时间为1 s，使用Matlab软件收集了3 000对输入输出数据样本.图6显示了该过程的输入和输出信号.输入u（k）是施加到泵上的电压，水泵产生流入上部水箱的流量，水经过上部储罐底部的一个小孔进入下储罐，下储罐的液位是双储罐系统的输出y（k）.

使用前500对数据作为训练集，剩余的2 500对数据作为验证集.使用CVA算法辨识线性动态子系统（LDS），先对训练集的输入输出数据构成的块Hankel矩阵生成的行空间进行加权投影，再根据式（20）~式（35）得到未知系统状态序列{\mathit{X}}_i的估计{\tilde{\mathit{X}}}_i和{\tilde{\mathit{X}}}_{i+\mathrm{1}}，此时求解一组超定方程即可得到状态空间矩阵 A， B， C 和 D .

训练集的输入数据经过CVA辨识的系统矩阵 A， B， C， D 构成的线性动态块，得到中间变量\tilde{\omega }（k），即非线性静态块的输入.NSS块中神经元的激活函数选择为\Phi (x)=\frac{\mathrm{1}-{\mathrm{e}}^{-\mathrm{0.5}x}}{\mathrm{1}+{\mathrm{e}}^{-\mathrm{0.5}x}}，突触权重向量（ w₁， w₂）和偏置参数（b₁，b₂）在训练开始时设定为［0，1］之间的随机数.LVW-LM-QN算法使用\tilde{\omega }（k）和系统输出数据以及误差f（k）进行迭代训练，不断减小成本函数V，最终实现V的最小化，同时得到构成NSS的多层感知器神经网络的权重w.Wiener模型在训练集数据上的线性块和非线性块的参数只需要训练集的输入输出数据就可辨识得到.

图7将所提出的LVW-LM-QN方法与LM方法以及Subspace-based方法在验证集上获得的输出进行比较.所提出的方法明显更优.为了进一步验证LVW-LM-QN的有效性，辨识结果与实际值的均方根误差和平均绝对百分比误差如表3所示.

从表3可以看出，所提方法的误差指标明显低于其他两种方法，体现了其在实际应用中的优越辨识性能，这说明LVW-LM-QN方法可以实现Wiener模型的有效辨识.另外，所提算法辨识的Wiener模型结构简单.线性动态子系统（LDS）的最佳近似值为2阶，非线性静态子系统（NSS）具有1个输入、1个输出和3个具有正切双曲激活函数的神经元.

4 结 语

本文研究Wiener系统的辨识，提出了使用面向块的Wiener系统来辨识非线性动态系统的方法，该模型在工程系统中具有广泛的应用前景和很强的实用性.采用子空间方法来辨识模型的线性动态子系统，多层感知器神经网络（MLP）则被用来表示静态非线性子系统.Weiner模型的两个子系统的参数分别使用CVA和LVW-LM-QN方法辨识.所提方法为Wiener系统辨识高度非线性系统提供了有效的辨识算法.仿真实验证明了所提方法具有良好的辨识性能.

参考文献

[1] 陈山， 宋樱， 房胜男， 等. 基于头脑风暴优化算法的Wiener模型参数辨识［J］. 控制与决策， 2017， 32（12）： 2291-2295.

[2] Chen Shan， Song Ying， Fang Sheng-nan， et al. Parameter identification of Wiener systems using brain storm optimization algorithm［J］. Control and Decision， 2017， 32（12）： 2291-2295.

[3] 白晶， 毛志忠， 浦铁成. 多变量Hammerstein-Wiener模型的参数辨识［J］. 东北大学学报（自然科学版）， 2018， 39（1）： 6-10.

[4] Bai Jing， Mao Zhi-zhong， Pu Tie-cheng. Parameter identification of multivariate Hammerstein-Wiener model［J］. Journal of Northeastern University （Natural Science）， 2018， 39（1）： 6-10.

[5] Mao Y W， Ding F. Multi-innovation stochastic gradient identification for Hammerstein controlled autoregressive autoregressive systems based on the filtering technique ［J］. Nonlinear Dynamics， 2015，79（3）： 1745-1755.

[6] Ding F， Liu X M， Liu M M. The recursive least squares identification algorithm for a class of Wiener nonlinear systems ［J］. Journal of the Franklin Institute， 2016， 353（7）： 1518-1526.

[7] Benítez M， Bermúdez A， Rodríguez-Calo J F. Adjoint method for parameter identification problems in models of stirred tank chemical reactors ［J］. Chemical Engineering Research and Design， 2017， 123： 214-229.

[8] Ghani F， Waser R， O’Donovan T S， et al. Non-linear system identification of a latent heat thermal energy storage system［J］. Applied Thermal Engineering， 2018， 134： 585-593.

[9] Hsu Y L， Wang J S. A Wiener-type recurrent neural network and its control strategy for nonlinear dynamic applications ［J］. Journal of Process Control， 2009， 19（6）： 942-953.

[10] Norquay S J， Palazoglu A， Romagnoli J A. Model predictive control based on Wiener models［J］. Chemical Engineering Science， 1998， 53（1）： 75-84.

[11] Norquay S J， Palazoglu A， Romagnoli J A. Application of Wiener model predictive control to a pH neutralization experiment［J］. IEEE Transactions on Control Systems Technology， 1999， 7（4）： 437-445.

[12] Shafiee G， Arefi M M， Jahed-Motlagh M R， et al. Nonlinear predictive control of a polymerization reactor based on piecewise linear Wiener model［J］. Chemical Engineering Journal， 2008， 143： 282-292.

[13] Andrzej J， Józef K. Two-stage instrumental variables identification of polynomial Wiener systems with invertible nonlinearities［J］. International Journal of Applied Mathematics and Computer Science， 2019， 29（3）： 571-580.

[14] Ma J H， Zhao L H， Han Z Z， et al. Identification of Wiener model using least squares support vector machine optimized by adaptive particle swarm optimization［J］. Journal of Control Automation and Electrical Systems， 2015， 26： 609-615.

[15] Qian S Y， Ding Z Y， Li F. Estimation of Wiener model based on neural fuzzy network［C］ // 2023 IEEE 12th Data Driven Control and Learning Systems Conference. Xiangtan， 2023：1377-1380.

[16] Hong X， Chen S. Modeling of complex-valued Wiener systems using B-spline neural network［J］. IEEE Transactions on Neural Networks， 2011， 22（5）： 818-825.

[17] 吴德会. 非线性动态系统的Wiener神经网络辨识法［J］. 控制理论与应用， 2009， 26（11）： 1192-1196.

[18] Wu De-hui. Identification method for nonlinear dynamic system using Wiener neural network［J］. Control Theory & Applications， 2009， 26（11）： 1192-1196.

[19] Fukushima K. Neocognitron： a self organizing neural network model for a mechanism of pattern recognition unaffected by shift in position［J］. Biological Cybernetics 1980， 36（4）： 193-202.

[20] Mohsen S， Mohammad F. Online identification of non-linear dynamic systems by Wiener system using subspace method and neural networks［J］. Transactions of the Institute of Measurement and Control， 2016， 40（2）： 666-674.

[21] Narendra K S， Annaswamy A M. Persistent excitation in adaptive systems［J］. International Journal of Control， 1987， 45（1）： 127-160.

[22] Bai E W. An optimal two-stage identification algorithm for Hammerstein-Wiener nonlinear systems［J］. Automatica， 1998， 34（3）： 333-338.

[23] Frandsen R， Poul E， Kristian J， et al. Unconstrained optimization［M］. Lyngby： Technical University of Denmark， 1999.

[24] Chiuso A. On the asymptotic properties of closed-loop CCA-type subspace algorithms： equivalence results and role of the future horizon［J］. IEEE Transactions on Automatic Control， 2010， 55（3）： 634-649.

[25] Loke M H， Dahlin T. A combined Gauss-Newton and quasi-Newton inversion method for the interpretation of apparent resistivity pseudo sections［C］// The 3rd Meeting of the European Association for Environmental and Engineering Geophysics. Aarhus， 1997： 139-142.

[26] Bansal J C， Singh P K， Saraswat M， et al. Inertia weight strategies in particle swarm optimization［C］// 2011 Third World Congress on Nature and Biologically Inspired Computing. Salamanca， 2011： 633-640.

[27] Larimore W. Canonical variate analysis in identification， filtering， and adaptive control［C］ // Proceedings of the 29th IEEE Conference on Decision and Control. Honolulu， 1990： 596-604.

[28] Gómez J C， Baeyens E. Subspace-based identification algorithms for Hammerstein and Wiener models［J］. European Journal of Control， 2005， 11（2）： 127-136.

[29] Chinthaned N， Sanposh P. Robust geometric control of a two-tank system［C］ //2016 13th International Conference on Electrical Engineering/Electronics， Computer， Telecommunications and Information Technology. Chiang Mai， 2016： 1-4.

基金资助

国家自然科学基金资助项目(62003082)

河北省自然科学基金资助项目(F2021501018)

河北省教育厅科学技术研究资助项目(ZD2022148)