定义非线性系统模型为：

式中，${y_k} \in R $是系统的实际输出；${\mathit{\boldsymbol{x}}_k} \in {R^{{n_x}}} $是回归向量；$ \mathit{\boldsymbol{p}} \in {R^n}$是待估计的参数向量；$f( \cdot, \cdot ) $是向量函数；$ {e_k} \in R$是不可测量误差。这里$ f({\mathit{\boldsymbol{x}}_k}, \mathit{\boldsymbol{p}})$对于参数向量p是非线性的。假设误差e_k有界，即：

式中，$ {\varepsilon _k}$是已知正数。

假设已获得数据集$\{ {\mathit{\boldsymbol{x}}_k}, {y_k}\} _{k = 1}^N $，由所有与模型式(1)、数据集$ \{ {\mathit{\boldsymbol{x}}_k}, {y_k}\} _{k = 1}^N$和有界误差假设式(2)相一致的参数向量所组成的集合可以表示为：

集合P_N称为参数可行集。式(3)可改写为：

式中，

定义加权无穷范数$\left\| \mathit{\boldsymbol{u}} \right\|_\infty ^{\bf{E}} $：$\left\| \mathit{\boldsymbol{u}} \right\|_\infty ^{\mathit{\boldsymbol{E}}} = {\max _i}\left| {{u_i}/{\varepsilon _i}} \right| $。可行集P_N的边界可以表达为：

e_k包含建模误差和测量噪声等，随机方法假设误差e_k的统计特性已知或部分已知。但是该假设在实际应用中可能会出现问题，原因有：当观测数据十分有限时，不易确定误差e_k的统计特性；在某些情况下误差e_k可能是确定性的。集员参数估计采用了更接近实际的有界误差假设式(2)，该方法的目的是有效描述参数可行集P_N。当系统是参数非线性时，可行集P_N比较复杂，一般很难精确描述。近似描述可行集的方法有4种，已经在前面介绍。需要说明的是，这里假设可行集P_N是有界且连通的。可行集是否连通涉及到可辨识性问题，该内容已超出本文的讨论范围，关于可辨识性请参考文献[16]及其引用文献。在假设可行集P_N是有界且连通的前提下，下面给出一种逼近可行集边界B_PN的方法。

根据几何与拓扑理论，可行集边界B_PN与n-1-sphere S^n-1同胚，它们之间存在同胚映射。因此，如果找到可行集边界B_PN与n-1-sphere S^n-1的同胚映射$ \phi ( \cdot )$，那么B_PN就可由S^n-1变换得到。但是由于可行集P_N形状复杂，$ \phi ( \cdot )$也就很难得到。为此，本文构造一个映射$\hat \phi ( \cdot ) $并使其能逼近映射$ \phi ( \cdot )$。使用$\hat \phi ( \cdot ) $可以把n-1-sphere S^n-1映射为可行集边界B_PN的近似$ {\hat B_{{P_N}}}$。构造映射$\hat \phi ( \cdot ) $包括以下3个步骤，这里用图 1辅助说明，但该方法不一定局限于二维：

图  1  映射$\hat \phi ( \cdot ) $的构造

1) 在可行集边界B_PN上均匀采样得到数据集$\mathit{\boldsymbol{\theta }}{}_i\left( {i = 1, 2, \cdot \cdot \cdot, l} \right) $，并使样本数据足够多。用Isomap将数据集$\{ {\mathit{\boldsymbol{\theta }}_i}\} _{i = 1}^l \subset {B_{{P_N}}} $映射为数据集$\{ {\mathit{\boldsymbol{\rho }}_i}\} _{i = 1}^l \subset {R^n} $；

2) 将数据集$\{ {\mathit{\boldsymbol{\rho }}_i}\} _{i = 1}^l \subset {R^n} $映射为数据集$\{ {\mathit{\boldsymbol{\eta }}_i}\} _{i = 1}^l \subset {S^{n-1}} $；

3) 基于数据集$\{ {\mathit{\boldsymbol{\theta }}_i}\} _{i = 1}^l $和$ \{ {\mathit{\boldsymbol{\eta }}_i}\} _{i = 1}^l$构造映射$ \hat \phi ( \cdot )$，以逼近可行集边界$ {B_{{P_N}}}$与n-1-sphere ${S^{n-1}} $之间的同胚映射$\phi ( \cdot ) $。

Isomap是一种概念上简单但是很实用的流形学习方法，它可以将输入数据映射到低维空间，并且在低维空间里数据的本征几何属性得到很好地保留。Isomap的优点在于较高的计算效率、较少的自由参数和非迭代全局优化。

实际上，可行集边界${B_{{P_N}}} $是一个$ n-1$维流形，但是很难通过Isomap将采样数据${\mathit{\boldsymbol{\theta }}_i}(i = 1, 2, \cdot \cdot \cdot, l) $嵌入到n-1维，因为可行集边界$ {B_{{P_N}}}$是“环形”的^[17]。所以，这里数据被嵌入到n维。

为了达到比较好的逼近效果，需要在可行集边界${B_{{P_N}}} $上采集足够多的样本点${\mathit{\boldsymbol{\theta }}_i} $。文献[11]提出了一种获取可行集边界的点的方法，但是这种方法受限于它的低效率。为了不丢失可行集边界的信息，它需要计算相当多的点，这就加大了Isomap算法的运行成本，并且该方法所获得的采样数据也不是均匀分布的，使Isomap的计算结果变差。为了克服这些缺点，文献[18]提出了一种简单却有效的采样方法，这里用图 2辅助说明：

图  2  可行集边界采样

1) 通过先验知识或者最优化方法得到一个包含可行集的盒子；

2) 定义一个均匀的网格以覆盖这个盒子；

3) 用无导数线搜索方法找到可行集边界与网格中盒子的边的交点。

Isomap将数据对$ {\mathit{\boldsymbol{\theta }}_i}$与$ {\mathit{\boldsymbol{\theta }}_j}$之间的欧式距离$ d(i, j)$作为输入。用Isomap将数据集$ \{ {\mathit{\boldsymbol{\theta }}_i}\} _{i = 1}^l$映射为数据集$\{ {\mathit{\boldsymbol{\rho }}_i}\} _{i = 1}^l $包含以下步骤：

1) 构建邻域图G。当$ {\mathit{\boldsymbol{\theta }}_i}$是$ {\mathit{\boldsymbol{\theta }}_j}$的K个最邻近点中的一个时，则将$ {\mathit{\boldsymbol{\theta }}_i}$和$ {\mathit{\boldsymbol{\theta }}_j}$相连，从而构成一个图G，设定路径的长度为$ d(i, j)$。

2) 计算最短路径。若$ {\mathit{\boldsymbol{\theta }}_i}$与$ {\mathit{\boldsymbol{\theta }}_j}$相连，则初始化最短路径$ {d_G}(i, j) = d(i, j)$；否则，令${d_G}(i, j) = + \infty $，对$ s = 1, 2, \cdot \cdot \cdot, l$，更新所有元素$ {d_G}(i, j)$为$ \min \{ {d_G}(i, j), {d_G}(i, s) + {d_G}(s, j)\} $。更新过程完成之后，可以得到最短距离矩阵${\mathit{\boldsymbol{D}}_G} = \{ {d_G}(i, j)\} $。

3) 计算嵌入。令${\lambda _s} $是矩阵$\mathit{\boldsymbol{\tau }}({\mathit{\boldsymbol{D}}_G}) $的第s个特征值(降序)，算子τ定义为$ \mathit{\boldsymbol{\tau }}(\mathit{\boldsymbol{D}}) =-\mathit{\boldsymbol{HLH}}/2$，其中矩阵L是距离的平方$\{ {L_{ij}} = D_{ij}^2\} $，H是中心矩阵$ \{ {H_{ij}} = {\delta _{ij}}-1/l\} $，$ {\delta _{ij}}$是Kronecker δ函数。$ \nu _s^i$是第s个特征向量的第i个分量，然后令向量${\mathit{\boldsymbol{\rho }}_i} $的第s个分量为$\sqrt {{\lambda _s}} \nu _s^i $。

运行Isomap所获得的嵌入保持了流形的本征几何属性，但是它的形状是不规则的。不难发现嵌入数据具有零均值和对角协方差。为了得到有规则形状的数据集，需将数据集$ \{ {\mathit{\boldsymbol{\rho }}_i}\} _{i = 1}^l$映射为数据集$\{ {\mathit{\boldsymbol{\eta }}_i}\} _{i = 1}^l \subset {S^{n-1}}:{\mathit{\boldsymbol{\eta }}_i} = {\mathit{\boldsymbol{\rho }}_i}/{\left\| {{\mathit{\boldsymbol{\rho }}_i}} \right\|_2}, i = 1, 2, \cdot \cdot \cdot, l $。

在完成前两个步骤后，每一个向量${\mathit{\boldsymbol{\theta }}_i} \in {B_{{P_N}}} $都有一个与之对应的向量${\mathit{\boldsymbol{\eta }}_i} \in {S^{n-1}} $。但是为了构造映射$ \hat \phi ( \cdot )$，需要计算与新的向量$\mathit{\boldsymbol{\eta }} \in {S^{n-1}} $相对应的向量$\mathit{\boldsymbol{\theta }} $。为了解决这个问题，这里用到了一种基于数据集$ \{ {\mathit{\boldsymbol{\theta }}_i}\} _{i = 1}^l$和$\{ {\mathit{\boldsymbol{\eta }}_i}\} _{i = 1}^l$的非参数方法^[19]。该方法包括以下步骤：

1) 在向量$ {\mathit{\boldsymbol{\eta }}_i}(i = 1, 2, \cdot \cdot \cdot, l)$中找到向量$ \mathit{\boldsymbol{\eta }}$的$K' $个最近邻；

2) 通过最小化代价函数$\sigma (\mathit{\boldsymbol{w}}) = {\left\| {\mathit{\boldsymbol{\eta }}-\sum\limits_{i = 1}^l {{w_i}{\mathit{\boldsymbol{\eta }}_i}} } \right\|_2} $(约束条件为：$ \sum\limits_{i = 1}^l {{w_i} = 1} $且若$ {\mathit{\boldsymbol{\eta }}_i}$不是向量$ \mathit{\boldsymbol{\eta }}$的最近邻，则$ {w_i} = 0$)计算最优线性重构权重$ {w_i}$；

3) 计算向量$\mathit{\boldsymbol{\theta }} = \sum\limits_{i = 1}^l {{w_i}{\mathit{\boldsymbol{\theta }}_i}} $。

用映射$ \hat \phi ( \cdot )$可求得与向量$ \mathit{\boldsymbol{\eta }} \in {S^{n-1}}$对应的向量$ \mathit{\boldsymbol{\theta }}$，而且不难看出$\hat \phi ( \cdot ) $是可行集边界$ {B_{{P_N}}}$与n-1-sphere $ {S^{n-1}}$之间同胚映射$\phi ( \cdot ) $的近似。一旦构造出$ \hat \phi ( \cdot )$，就可以将${S^{n-1}} $映射为$ {B_{{P_N}}}$的近似。

要构造映射$ \hat \phi ( \cdot )$需要确定两个自由参数：使用Isomap将数据集$\{ {\mathit{\boldsymbol{\theta }}_i}\} _{i = 1}^l $映射为数据集$\{ {\mathit{\boldsymbol{\rho }}_i}\} _{i = 1}^l $时所涉及的近邻个数K和计算与新的向量$ \mathit{\boldsymbol{\eta }} \in {S^{n-1}}$相对应的向量$ \mathit{\boldsymbol{\theta }}$时所涉及的近邻个数K'。根据文献[19]可确定K'。一旦K'选定，参数K就决定了近似边界的性能。本文提出一种能够确定参数K的方法。

令$ {\hat \phi _K}( \cdot )$表示完成上述过程后得到的映射(参数K'已经选定)。在n-1-sphere ${S^{n-1}} $上均匀采样得到向量$ {\mathit{\boldsymbol{\bar \eta }}_i}, i = 1, 2, \cdot \cdot \cdot, \bar l$，并假设数据足够多。引入变量$\sigma (K) $表示边界近似误差，有：

式中，

$\sigma (K) $越小，边界近似效果越好。因此，最优近邻个数$ {K^ * }$应该是$\arg {\min _K}\sigma (K) $。最后，所构造的映射$ \hat \phi ( \cdot )$定为${\hat \phi _{{K^ * }}}( \cdot ) $。

考虑如下非线性系统模型^[2]：

式中，$ {y_k}$和${u_k} $分别是系统输出和输入；${p_1} $和$ {p_2}$是参数。可以看出，回归向量${\mathit{\boldsymbol{x}}_k} = {[{y_{k-1}}, {u_{k-1}}]^{\rm{T}}} $。实验过程中，设定$ {p_1} = 0.5$、$ {p_2} =-0.2$、$ {y_1} = 0$，并且设计输入$ {u_k} = \cos (k/3) + \sin (k) + 1$。假设误差${e_k}$在区间$[-1.5, 1.5] $上服从均匀分布，因而，误差界为1.5。实验获得数据集$ \{ {\mathit{\boldsymbol{x}}_k}, {y_k}\} _{k = {\rm{2}}}^{{\rm{45}}}$，误差序列$\{ {e_k}\} _{k = 2}^{45} $={0.628 4, −0.213 3, −0.586 1, −0.931 0, −0.919 7, 0.546 7, −0.591 7, 0.125 0, −1.047 4, 0.593 7, −0.364 9, 1.080 0, 1.061 0, 0.280 7, −0.010 3, 1.199 3, 0.964 9, 0.434 7, 0.953 9, 0.480 7, −0.474 1, −0.630 8, −0.476 4, 0.102 2, 0.681 3, -0.572 1, 1.015 5, 0.204 2, −0.388 8, 0.608 2, 0.139 7, −0.165 4, 0.583 7, 0.363 9, 0.884 5, 1.370 5, 0.067 8, 1.140 4, −0.981 1, 1.439 2, −0.685 7, −0.743 0, 1.127 2, 0.711 9}。

借鉴文献[1]的思想，采用最优化方法得到一个包含可行集的盒子$ [{\rm{0}}{\rm{.2}}, 1] \times [-{\rm{0}}{\rm{.4}}, 0.{\rm{2}}]$，定义一个尺寸为$ {\rm{41}} \times {\rm{31}}$的均匀网格以覆盖这个盒子。计算可行集边界与网格中盒子的与$ {p_{\rm{2}}}$轴平行的边的交点，共获取76个。图 3显示了这些可行集边界采样点，参数$K' $设定为2。在1-sphere ${S^{\rm{1}}} $上均匀取400个点$ (\bar l = 400)$代入式(9)和式(10)，以确定最优近邻个数K。由于$ \arg {\min _K}\sigma (K) = 28$，则设定K为28。图 4所示为σ与K的对应关系。图 5a和图 5b分别表示用本文方法和SVM方法^[12]得到的可行集近似边界。对于后者，使用了在$[{\rm{0}}{\rm{.2}}, 1] \times [-{\rm{0}}{\rm{.4}}, 0.{\rm{2}}] $上均匀采样得到的130个点，及RBF核宽度为$\sqrt {4 \times {{10}^{-3}}} $的LS-SVM来获得最优结果。图中实线和虚线分别表示精确边界和近似边界，加号表示参数向量的真值。从图中可以看出，本文方法比SVM方法具有更好的逼近精度。同时也可以看出，用单一凸集合作为可行集外界的方法所得结果要比本文方法保守性高。

图  3  可行集边界采样

图  4  σ与K的对应关系

图  5  可行集近似边界

考虑经典的单室开放一级吸收模型^[11]：

式中， ${y_k}$是在 ${x_k}$时刻观察到的药物浓度； ${K_a}$是吸收速率常数；${K_e} $是消除速率常数；Dose是药物剂量；V是分布容积。假设常数${K_a}$和$ {K_e}$是待估计量，实验过程中，设定${K_a} = 1.{\rm{63}}\;{\rm{h}}{{\rm{r}}^{ - 1}}$、${K_e} = 0.{\rm{38}}\;{\rm{h}}{{\rm{r}}^{ - 1}}$、$V = {\rm{6}}0\;{\rm{L}}$。口服药物1 100 mg后，分别在0.25、0.5、1、1.5、2、4、6、8、10、12、18、24 h后测量血浆中的药物浓度，从而获得一个输入输出数据集。假设误差${e_k}$服从均值为0 mg/L和标准差为0.25/3 mg/L的截断正态分布，误差序列$\{ {e_k}\} _{k = {\rm{1}}}^{{\rm{12}}}$={−0.180 9, −0.004 9, −0.084 2, 0.051 2, 0.042 3, 0.141 0, 0.049 3, −0.053 6, 0.031 7, −0.084 1, −0.001 6, −0.004 0}。因而，误差界为0.25 mg/L。

借鉴文献[1]的思想，采用最优化方法得到一个包含可行集的盒子$[1.{\rm{54}}, 1.{\rm{66}}] \times [0.{\rm{355}}, 0.{\rm{395}}]$。定义一个尺寸为${\rm{31}} \times {\rm{21}} $的均匀网格以覆盖这个盒子。计算可行集边界与网格中盒子的边的交点，共获取84个。图 6显示了这些可行集边界采样点。参数K'设定为6。在1-sphere $ {S^{\rm{1}}}$上均匀取400个点$ (\bar l = 400)$代入式(9)和式(10)，以确定最优近邻个数K。由于$\arg {\min _K}\sigma (K) = 25 $，则设定K为25。图 7示出了σ与K的对应关系。图 8a和8b分别为用本文方法和SVM方法^[12]得到的可行集近似边界。对于后者，使用了在$ [1.{\rm{54}}, 1.{\rm{66}}] \times [0.{\rm{355}}, 0.{\rm{395}}]$上均匀采样得到的147个点，及RBF核宽度为$\sqrt {2 \times {{10}^{-4}}} $的LS-SVM来获得最优结果。图中实线和虚线分别表示精确边界和近似边界，加号表示参数向量的真值。从图中可以看出，本文方法比SVM方法具有更好的逼近精度。

图  6  可行集边界采样

图  7  σ与K的对应关系

图  8  可行集近似边界

本文给出一种新的非线性系统集员参数估计方法，可以处理针对复杂模型较难构建最小包含函数的问题。该方法假设可行集是有界且连通的，寻找可行集边界与n-1-sphere之间同胚映射的近似。映射构造成功后，就可以将n-1-sphere映射为可行集的近似边界。本文通过两个例子说明该方法的优势。接下来的工作将研究如何将方法延伸到用于有界非连通可行集的边界逼近。