受控自回归滑动平均模型（ARMAX）的系统辨识

受控自回归滑动平均模型 (Controled Auto Regression and Moving Average model, CARMA)，亦称带外部输入的自回归滑动平均模型 (Auto Regression and Moving Average model with eXogenous input, ARMAX)是应用非常广泛的线性系统模型，本文介绍该模型的一种系统辨识方法：最大似然法。

系统模型

$$\begin{gathered} y_k+a_1{y}_{k-1}+\cdots a_n{y}_{k-n}= \\ {b}_1{u}_{k-1}+b_n{u}_{k-n}+e_k+c_1{e}_{k-1}+\cdots +c_n{e}_{k-n} \end{gathered}$$
即
$$y_k=-\sum _{i=1}^n{a}_i{y}_{k-i}+\sum _{i=1}^n{b}_i{u}_{k-i}+\sum _{i=1}^n{c}_i{e}_{k-i}+e_k$$
该模型在 ARMA 的基础上考虑了外部输入 $u(k)$ 对输出的影响，为方便讨论，在此设定 $a,b,c$ 的下标都是从 $1$ 到 $n$，实际上不必如此。

似然函数

在数理统计学中，似然函数是一种关于统计模型中的参数的函数，表示模型参数中的似然性。
给定输出x时，关于参数θ的似然函数L(θ|x)（在数值上）等于给定参数θ后变量X的概率：

$$L(\theta |X)=P(X|\theta )$$

在统计学中，“似然性”和“概率”有明确的区分。概率用于在已知一些参数的情况下，预测接下来的观测所得到的结果，而似然性则是用于在已知某些观测所得到的结果时，对有关事物的性质的参数进行估计。

假定误差 $e(k)$ 服从 $0$ 均值、方差为 ${\sigma }^2$ 的高斯分布，则有似然函数为：
P ( Y N ∣ U N , Θ ) = P ( y N , … ∣ u N , … , Θ ) = P ( y 0 ) Π k = 1 N P ( y k ∣ Y k − 1 , U N , Θ ) = P ( y 0 ) Π k = 1 N P ( e k ) = P ( y 0 ) ( 2 π ) − N / 2 σ − N exp ⁡ [ − ( 1 / 2 σ 2 ) ∑ k = 1 N e k 2 ] (1)

$$\begin{array}{ll} P(Y_N | U_N,\Theta) &= P(y_N, \ldots | u_N, \ldots, \Theta) \\\\ &= P(y_0)\Pi _{k=1}^NP(y_k| Y_{k-1},U_N,\Theta) \\\\ &= P(y_0)\Pi _{k=1}^NP(e_k) \\\\ &= P(y_0)(2\pi)^{-N/2}\sigma^{-N} \exp \left[ -(1/2\sigma^2)\sum_{k=1}^N e_k^2\right] \tag{1} \end{array}$$ P(YN​∣UN​,Θ)​=P(yN​,…∣uN​,…,Θ)=P(y0​)Πk=1N​P(yk​∣Yk−1​,UN​,Θ)=P(y0​)Πk=1N​P(ek​)=P(y0​)(2π)−N/2σ−Nexp[−(1/2σ2)∑k=1N​ek2​]​(1)
其中
$$\begin{array}{ll}{e}_k& =y_k-{\varphi }_i\Theta \\ & \\ {\varphi }_k& =[-y_{k-1},\dots ,-y_{k-n}|u_{k-1},\dots ,u_{k-n}|e_{k-1},\dots ,e_{k-n}{]}^{\top }\\ & \\ \Theta & =[a_1,\dots ,a_n|b_i,\dots ,b_n|c_1,\dots ,c_n{]}^{\top }\end{array}$$
最大化似然(1)等价于最小化负对数似然(2)
$$J(\sigma ,\Theta )=N\ln \sigma +\frac{1}{2{\sigma }^2}\sum _k^N{e}_k^2\text{\hspace{1em}(2)}$$

辨识过程

在实际辨识过程中，由于参数 $\Theta$ 未知，误差 $e_k$ 不能精确获得，所以计算过程中用其估计值 ${\nu }_k\simeq e_k$ 来替代。

代价函数：
$$J({\sigma }_{\nu },\Theta )=N\ln {\sigma }_{\nu }+\frac{1}{2{\sigma }_{\nu }^2}\sum _k^N{\nu }_k^2\text{\hspace{1em}(2)}$$

• 采集 $N$ 组数据，估计一个初始 ${\Theta }_0$，比如先假定误差项系数 $c_i=0$，用最小二乘法求解$a_i,b_i$；
• $k=0$
• 固定 $\Theta$（即固定$\nu$），更新 ${\sigma }_{\nu }$:
  $${\sigma }_{\nu }=\underset{{\sigma }_{\nu }}{\mathrm{argmin}}J=\sqrt{\sum _k{\nu }_k^2/N}$$
• 固定${\sigma }_{\nu }$， 更新 $\Theta$，即最小化
  $$L=\sum _k{\nu }_k^2$$
  采用牛顿法更新参数：
  $${\Theta }_{t+1}={\Theta }_t-H^{-1}{\nabla }_{\Theta }L$$
• $k=k+1$，重复以上交替优化过程直到
  $$\frac{{\sigma }_t^2-{\sigma }_{t-1}^2}{{\sigma }_{t-1}^2}<10^{-4}$$

关于参数的梯度与海森矩阵

梯度

$$\frac{\partial L}{\partial \Theta }=2\sum {\nu }_k\frac{\partial {\nu }_k}{\partial \Theta }\text{\hspace{1em}(3)}$$
这里稍稍有一点复杂，因为：
$${\nu }_k=y_k+\sum _{i=1}^n{a}_i{y}_{k-i}-\sum _{i=1}^n{b}_i{u}_{k-i}-\sum _{i=1}^n{c}_i{\nu }_{k-i}=y_k-{\varphi }_k^{\top }\Theta$$
所以
$$\begin{gathered} \frac{\partial {\nu }_k}{\partial a_i}=y_{k-i}-\sum _{l=1}^n{c}_l\frac{\partial {\nu }_{k-l}}{\partial a_i} \\ \frac{\partial {\nu }_k}{\partial b_i}=-u_{k-i}-\sum _{l=1}^n{c}_l\frac{\partial {\nu }_{k-l}}{\partial b_i} \\ \frac{\partial {\nu }_k}{\partial c_i}=-{\nu }_{k-i}-\sum _{l=1}^n{c}_l\frac{\partial {\nu }_{k-l}}{\partial c_i} \end{gathered}$$
可以合并成
$$\frac{\partial {\nu }_k}{\partial \Theta }=-{\varphi }_k^{\top }-\sum _{l=1}^n{c}_l\frac{\partial {\nu }_{k-l}}{\partial \Theta }$$
所以$\partial L/\partial \Theta$ 已求得！

海森矩阵

假定 $\partial {\nu }_k/\partial \Theta$ 为行向量
$$\begin{gathered} H=\frac{{\partial }^{2}L}{\partial {\Theta }^2}=\frac{\partial }{\partial \Theta }\left(2\sum {\nu }_k\frac{\partial {\nu }_k}{\partial \Theta }\right) \\ =2\left(\sum _k{\left(\frac{\partial {\nu }_k}{\partial \Theta }\right)}^{\top }\left(\frac{\partial {\nu }_k}{\partial \Theta }\right)+{\nu }_k\frac{{\partial }^2{\nu }_k}{\partial {\Theta }^2}\right) \end{gathered}$$
主要需要考虑 ${\partial }^2\nu /\partial {\Theta }^2$:
$$\begin{gathered} \frac{{\partial }^2{\nu }_k}{\partial a_i\partial a_j}=-\sum _{l=1}^n{c}_l\frac{{\partial }^2{\nu }_{k-l}}{\partial a_i\partial a_j} \\ \frac{{\partial }^2{\nu }_k}{\partial a_i\partial b_j}=-\sum _{l=1}^n{c}_l\frac{{\partial }^2{\nu }_{k-l}}{\partial a_i\partial b_j} \\ \frac{{\partial }^2{\nu }_k}{\partial a_i\partial c_j}=-\frac{\partial {\nu }_{k-j}}{\partial a_i}-\sum _{l=1}^n{c}_l\frac{{\partial }^2{\nu }_{k-l}}{\partial a_i\partial c_j} \\ \frac{{\partial }^2{\nu }_k}{\partial b_i\partial c_j}=-\frac{\partial {\nu }_{k-j}}{\partial b_i}-\sum _{l=1}^n{c}_l\frac{{\partial }^2{\nu }_{k-l}}{\partial a_i\partial c_j} \\ \frac{{\partial }^2{\nu }_k}{\partial c_i\partial c_j}=-2\frac{\partial {\nu }_{k-j}}{\partial c_i}-\sum _{l=1}^n{c}_l\frac{{\partial }^2{\nu }_{k-l}}{\partial c_i\partial c_j} \end{gathered}$$
写成矩阵形式：
$$\begin{gathered} H_k=-\sum _{l=1}^n{c}_l{H}_{k-l}-G_k-G_k^{\top } \\ {G}_k={\left[0\dots 0|0\dots 0|{\left(\frac{\partial {\nu }_{k-1}}{\partial \Theta }\right)}^{\top }\dots {\left(\frac{\partial {\nu }_{k-n}}{\partial \Theta }\right)}^{\top }\right]}_{3n\times 3n} \end{gathered}$$
再次强调一下：$\partial {\nu }_k/\partial \Theta$ 是行向量，因为 $\Theta$ 是列向量！

到此，完整的计算过程应该心中有数了吧！