基于广义多项式混沌算法的Van der Pol方程参数辨识

1 广义随机多项式（Generalized Polynomial Chaos）

系统的动态特征可以表述为
$$\{\begin{array}{l}\dot{x}=F(t,x,u;\theta ),0\le t\le t_f,x(0)=x_0\\ y=C(\theta )x\end{array}\text{\hspace{1em}(1)}$$

其中，$x\in {\mathbb{R}}^{n_s}$，存在$n_s$个状态变量；$\theta \in {\mathbb{R}}^{n_p}$，存在$n_p$个未知参数；$C(\theta )$为线性矩阵。

根据gPC理论，将系统参数$\theta (\xi )$与状态变量$x(t,\xi )$利用正交多项式函数进行展开，第$i$个参数${\theta }_i$可以分解为随机变量${\xi }_i$表示的：
$${\theta }_i({\xi }_i)=\sum _{m=0}^{N_p}{\theta }_{i,m}{\varphi }_m({\xi }_i),i=1,...,n_p\text{\hspace{1em}(2)}$$

其中，任意$i$与$m$的${\theta }_{i,m}$这里认为已知（例如之前对数分布的$\theta (\xi )=e^{\mu +\sigma \xi }$关系），且${\varphi }_m({\xi }_i)$代表正交基函数。如果第$i$个参数${\theta }_i$对应的概率密度函数为${\rho }_i({\xi }_i)$，那么相应的正交关系为
$$\mathbb{E}[{\varphi }_m({\xi }_i){\varphi }_n({\xi }_i)]=\int {\varphi }_m({\xi }_i){\varphi }_n({\xi }_i){\rho }_i({\xi }_i)d{\xi }_i\text{\hspace{1em}(3)}$$

同时，第$j$个状态变量${\hat{x}}_j(t,\xi )$可以分解为随机变量$\xi$表示的：
$${\hat{x}}_j(t,\xi )=\sum _{|\alpha |=0}^{N_s}{x}_{j,\alpha }(t){\Phi }_{\alpha }(\xi ),j=1,...,n_s\text{\hspace{1em}(4)}$$

需要说明的是，标准的gPC理论将系统参数视为随机变量的函数，但是考虑到系统参数的不确定性会直接导致状态变量的不确定性，因而可以同样将状态变量视为随机变量的函数。由于${\xi }_i(i=1,...,n_p)$之间是相互独立的，从而有联合概率密度函数$\rho (\xi )$以及$n_p$变量多项式基函数${\Phi }_{\alpha }$：
$$\rho (\xi )=\prod _{i=1}^{n_p}{\rho }_i({\xi }_i),{\Phi }_{\alpha }(\xi )=\prod _{i=1}^{n_p}{\varphi }_i({\xi }_i),\xi =[{\xi }_1,...,{\xi }_i,...,{\xi }_{n_p}]\text{\hspace{1em}(5)}$$

相应的正交关系有
$$\mathbb{E}[{\Phi }_{\alpha }(\xi ){\Phi }_{\beta }(\xi )]=\int {\Phi }_{\alpha }(\xi ){\Phi }_{\beta }(\xi )\rho (\xi )d\xi \text{\hspace{1em}(6)}$$

值得注意的是，这里的$\alpha /\beta \in {\mathbb{N}}_0^{n_p}$属于向量形式的$n_p$维标号，即$\alpha =[{\alpha }_1,...,{\alpha }_{n_p}]$，其绝对值代表所有元素标号之和，即$|\alpha |={\alpha }_1+...+{\alpha }_{n_p}$。在此基础之上，第$j$个状态变量$x_j(t,\xi )$的分解系数$x_{j,\alpha }(t)$的总个数（式（4）的多项式总个数）为
$$r=\frac{(N_s+n_p)!}{N_s!n_p!}\text{\hspace{1em}(7)}$$

虽然系统参数的多项式分解系数是已知的，但是状态变量的多项式分解则需要单独计算。在$x_j(t,\xi )$未知的情况下，需要使用Galerkin方法或随机配置法求解；在$x_j(t,\xi )$已知的情况下，系数可以通过如下公式计算
$$x_{j,\mathbf{\alpha }}(t)=\frac{\mathbb{E}[x_j(t,\xi ){\Phi }_{\mathbf{\alpha }}(\xi )]}{\mathbb{E}[{\Phi }_{\mathbf{\alpha }}(\xi ){\Phi }_{\mathbf{\alpha }}(\xi )]}=\frac{\mathbb{E}[{\hat{x}}_j(t,\xi ){\Phi }_{\mathbf{\alpha }}(\xi )]}{\mathbb{E}[{\Phi }_{\mathbf{\alpha }}(\xi ){\Phi }_{\mathbf{\alpha }}(\xi )]}\text{\hspace{1em}(8)}$$

2 随机配执法（Stochastic Collection Method）

首先，依据全局概率密度函数$\rho (\xi )$选取一个确定的配置点集 { μ ( 1 ) , . . . , μ ( i ) , . . . , μ ( Q ) } \{\mu^{(1)},...,\mu^{(i)},...,\mu^{(Q)}\} {μ(1),...,μ(i),...,μ(Q)}，其中$Q\ge r$表示配置点的数量。需要特别注意的是，$Q$必须要大于等于$r$，否则会导致之后的$A^{\#}$矩阵不存在。全部点均存在于参数变量空间，即${\mu }^{(i)}\in {\mathbb{R}}^{n_p}$。利用这些点替代方程(1)中的随机变量$\xi$，有
$$\{\begin{array}{l}{\dot{z}}^{(i)}=F(t,z^{(i)}(t);\theta ({\mu }^{(i)})),i=1,...,Q\\ z^{(i)}(0)=z_0^{(i)}\end{array}\text{\hspace{1em}(9)}$$

其中，在第$i$个定点的第$j$个状态变量状态变量，有
$$z_j^{(i)}=\sum _{|\alpha |=0}^{N_s}{x}_{j,\alpha }(t){\Phi }_{\alpha }({\mu }^{(i)}),j=1,...,n_s\text{\hspace{1em}(10)}$$

这里，可以利用数值积分的方法求解(9)中每一个配置点下的状态变量，从而得到$Q$组解耦的集合，且每个集合包含$n_s$个状态变量，有
$$Z={\left[\begin{array}{c}{z}^{(1)}\\ ⋮\\ z^{(i)}\\ ⋮\\ z^{(Q)}\end{array}\right]}_{Q\cdot n_s\times 1}\in {\mathbb{R}}^{Q\times n_s},z^{(i)}={\left[\begin{array}{c}{z}_1^{(i)}\\ ⋮\\ z_j^{(i)}\\ ⋮\\ z_{n_s}^{(i)}\end{array}\right]}_{n_s\times 1}\in {\mathbb{R}}^{n_s}\text{\hspace{1em}(11)}$$

其次，同时考虑到式(4)可以表示为内积的形式，则$n_s$中第$j$个状态变量可以写为矩阵形式，有
$$\begin{gathered} {\hat{x}}_j(t,\xi )=\sum _{|\alpha |=0}^{N_s}{x}_{j,\alpha }(t){\Phi }_{\alpha }(\xi )=P(\xi {)}^T{X}_j(t),j=1,...,n_s \\ \text{\hspace{1em}(12)} \end{gathered}$$

其中，
$$P(\xi )={\left[\begin{array}{c}⋮\\ {\Phi }_{\alpha }(\xi )\\ ⋮\end{array}\right]}_{r\times 1}\in {\mathbb{R}}^r,X_j(t)={\left[\begin{array}{c}⋮\\ x_{j,\alpha }(t)\\ ⋮\end{array}\right]}_{r\times 1}\in {\mathbb{R}}^r\text{\hspace{1em}(13)}$$

进一步考虑全部的状态变量$\hat{x}(t,\xi )$，有
$$\hat{x}(t,\xi )={\left[\begin{array}{c}{\hat{x}}_1(t,\xi )\\ ⋮\\ {\hat{x}}_j(t,\xi )\\ ⋮\\ {\hat{x}}_{n_s}(t,\xi )\end{array}\right]}_{n_s\times 1}={\left[\begin{array}{ccc}{P}^T(\xi )& O& O\\ O& \ddots & O\\ O& O& P^T(\xi )\end{array}\right]}_{n_s\times r\cdot n_s}{\left[\begin{array}{c}{X}_1(t)\\ ⋮\\ X_j(t)\\ ⋮\\ X_{n_s}(t)\end{array}\right]}_{r\cdot n_s\times 1}=\mathbb{P}(\xi )X(t)\text{\hspace{1em}(14)}$$

之后，考虑$Q$个配置点，将式(11)与(14)联立，即$\hat{x}=z$可以得到
$$Z={\left[\begin{array}{c}{z}^{(1)}\\ ⋮\\ z^{(i)}\\ ⋮\\ z^{(Q)}\end{array}\right]}_{Q\cdot n_s\times 1}={\left[\begin{array}{c}\mathbb{P}({\mu }^{(1)})\\ ⋮\\ \mathbb{P}({\mu }^{(i)})\\ ⋮\\ \mathbb{P}({\mu }^{(Q)})\end{array}\right]}_{Q\cdot n_s\times r\cdot n_s}{\left[\begin{array}{c}{X}_1(t)\\ ⋮\\ X_j(t)\\ ⋮\\ X_{n_s}(t)\end{array}\right]}_{r\cdot n_s\times 1}=\mathbb{A}X(t)\text{\hspace{1em}(15)}$$

不难发现，向量$X(t)$由$X_j(t)$组成，而$X_j(t)$由$x_{j,\alpha }(t)$组成，也就是说$X(t)$中的$r\cdot n_s$个元素实际上对应$n_s$个状态变量的$r$个概率多项式分解系数。

最后，可以通过如下矩阵运算求解相应的分解系数
$$\begin{gathered} X(t)={\mathbb{A}}^{\#}Z,{\mathbb{A}}^{\#}=({\mathbb{A}}^T\mathbb{A}{)}^{-1}{\mathbb{A}}^T\in {\mathbb{R}}^{r\cdot n_s\times Q\cdot n_s} \\ \text{\hspace{1em}(16)} \end{gathered}$$

其中，矩阵运算${\mathbb{A}}^{\#}$表示Moore-Penrose伪逆，配置点的选择过程中需要确保这一伪逆存在（即$Q\ge r$）。

3 极大似然法（Maximum Likelihood Approach）

问题提出

下面研究在实际系统观测结果下，随机变量$\xi$最佳数值的迭代参数辨识方案。在此基础之上，基于式(2)即可实现未知系统参数$\theta (\xi )$的计算，这里重写如下
$${\theta }_i({\xi }_i)=\sum _{m=0}^{N_p}{\theta }_{i,m}{\varphi }_m({\xi }_i),i=1,...,n_p\text{\hspace{1em}(17)}$$

假设系统输出的噪声符合期望为零高斯分布，协方差矩阵为$R\in {\mathbb{R}}^{n_y\times n_y}$。同时，认为系统的每次采样间的观测结果是相互独立的，方程(1)中唯一的不确定来自于未知的参数，并且状态变量与系统参数的概率多项式分解存在。基于以上假设，可以将似然函数写为
$$\mathcal{L}(\xi |y_0,..,y_t)=\prod _{\tau =0}^t\rho (y(\tau )|\xi )\propto e^{-\frac{1}{2}{\sum }_{\tau =0}^t[y(\tau )-\hat{y}(\tau ,\xi ){]}^T{R}^{-1}[y(\tau )-\hat{y}(\tau ,\xi )]}\text{\hspace{1em}(18)}$$

其中，$\mathcal{L}(\xi |y_0,..,y_t)$是随机变量$\xi$在全时刻系统观测结果$y_0,..,y_t$条件下的似然函数，$\rho (y(t)|\xi )$是观测结果$y(t)$在参数$\xi$条件下的概率密度函数，$\hat{y}(t,\xi )$是随机模型的输出，即
$$\hat{y}=C(\theta (\xi ))\hat{x}(t,\xi )\text{\hspace{1em}(19)}$$

根据定义，极大似然辨识结果$\hat{\theta }$使得(18)具有最大值，等价于不含复数的指数部分具有最小值，即
$$\begin{gathered} \mathcal{J}(t,\xi )=\frac{1}{2}\sum _{\tau =0}^t[y(\tau )-\hat{y}(\tau ,\xi ){]}^T{R}^{-1}[y(\tau )-\hat{y}(\tau ,\xi )] \\ \hat{\xi }=argmin\mathcal{J}(t,\xi )\text{\hspace{1em}(20)} \end{gathered}$$

问题分析

$\mathcal{J}(t,\xi )$的迭代更新是该方法递归的重要条件，这里利用概率多项式的性质来分离方程中时间与未知参数两部分，从而实现这种递归方法。将(14)与(19)代入(20)，并进行代数化简，有
$$\mathcal{J}(t,\xi )=\frac{1}{2}\sum _{i=1}^{n_y}\sum _{j=1}^{n_y}[R^{-1}{]}_{i,j}[(\sum _{\tau =0}^t{y}_i{y}_j)-2C_i\mathbb{P}(\sum _{\tau =0}^{t}X{y}_j)+C_i\mathbb{P}(\sum _{\tau =0}^{t}X{X}^T)(C_i\mathbb{P}{)}^T]\text{\hspace{1em}(21)}$$

其中，$[R^{-1}{]}_{i,j}$表示协方差矩阵逆中第$i$行向量中的第$j$个元素，标量$y_k$表示观测结果向量$y$中$n_y$个元素中的第$k$个，$C_k$表示输出矩阵$C$中的第$k$行。

可以看到，求和项中仅包含与时间有关的$X$、$y$项，并包含系数$X$、$\mathbb{P}$，只有对状态变量使用概率多项式分解后才能够得到这样的结构。

为了方便实现这里考虑求和项的递归实现方案。对于任意一个矩阵$G(t)$定义如下形式

$$\begin{gathered} D_G^{(t)}=\sum _{\tau =0}^{t}G(\tau ) \\ {D}_G^{(t+\Delta t)}=\sum _{\tau =0}^{t+\Delta t}G(\tau ) \end{gathered}$$
即
$$\Rightarrow D_G^{(t+\Delta t)}=D_G^{(t)}+G(t+\Delta t)\text{\hspace{1em}(22)}$$

其中，$\Delta t$代表采样间隔。

定义$D_{y_i{y}_j}^{(t)}={\sum }_{\tau =0}^t{y}_i{y}_j\in \mathbb{R}$、$D_{X{y}_j}^{(t)}={\sum }_{\tau =0}^{t}X{y}_j\in {\mathbb{R}}^{r\cdot n_s}$、$D_{X{X}^T}^{(t)}={\sum }_{\tau =0}^{t}X{X}^T\in {\mathbb{R}}^{r\cdot n_s\times r\cdot n_s}$为固定的变量，并将(21)重写为
$$\mathcal{J}(t,\xi )=\frac{1}{2}\sum _{i=1}^{n_y}\sum _{j=1}^{n_y}[R^{-1}{]}_{i,j}[D_{y_i{y}_j}^{(t)}-2C_i\mathbb{P}{D}_{X{y}_j}^{(t)}+C_i\mathbb{P}{D}_{X{X}^T}^{(t)}(C_i\mathbb{P}{)}^T]\text{\hspace{1em}(23)}$$
其中定义的变量可以在固定内存中存储相关数据，并根据时间进行迭代求解。

梯度下降优化

基于梯度的参数更新规则，有
$${\hat{\xi }}_{t+1}={\hat{\xi }}_t-\Gamma \frac{\partial \mathcal{J}(t,\xi )}{\partial \xi }{|}_{\xi =\hat{\xi }}\text{\hspace{1em}(24)}$$

其中，${\hat{\xi }}_t$代表$t$时刻$\xi$的辨识结果，不同的$\Gamma$矩阵对应了不同类型的下降方法。

将(23)代入(24)，可以得到
$${\hat{\xi }}_{t+1}={\hat{\xi }}_t-\Gamma \sum _{i=1}^{n_y}\sum _{j=1}^{n_y}[R^{-1}{]}_{i,j}\times \left[\frac{\partial C_i\mathbb{P}}{\partial \xi }{|}_{\xi =\hat{\xi }}\left(D_{X{X}^T}^{(t)}(C_i\mathbb{P}{)}^T{|}_{\xi =\hat{\xi }}-D_{X{y}_j}^{(t)}\right)\right]\text{\hspace{1em}(25)}$$

注意到，上式中$D_{X{X}^T}^{(t)}$与$D_{X{y}_j}^{(t)}$可能会随时间无约束的增加，因此这里需要进行归一化，有
$${\hat{\xi }}_{t+1}={\hat{\xi }}_t-\Gamma \sum _{i=1}^{n_y}\sum _{j=1}^{n_y}[R^{-1}{]}_{i,j}\times \left[\frac{\partial C_i\mathbb{P}}{\partial \xi }{|}_{\xi =\hat{\xi }}\times \frac{D_{X{X}^T}^{(t)}(C_i\mathbb{P}{)}^T{|}_{\xi =\hat{\xi }}-D_{X{y}_j}^{(t)}}{1+||D_{X{y}_j}^{(t)}|{|}_2}\right]\text{\hspace{1em}(26)}$$

其中，标量$||D_{X{y}_j}^{(t)}|{|}_2=||{\sum }_{\tau =0}^{t}X{y}_j|{|}_2$。特别地，偏微分导数项有
$$\frac{\partial C_i\mathbb{P}}{\partial \xi }{|}_{\xi =\hat{\xi }}=\frac{\partial C_i}{\partial \xi }\mathbb{P}{|}_{\xi =\hat{\xi }}+C_i\frac{\partial \mathbb{P}}{\partial \xi }{|}_{\xi =\hat{\xi }}\text{\hspace{1em}(27)}$$

4 范德波尔方程参数辨识（Van der Pol）

Van der Pol方程可以描述为
$$\{\begin{array}{l}{\dot{x}}_1=x_2\\ {\dot{x}}_2=-x_1-\frac{1}{\theta }(x_1^2-1)x_2+\sin (50t)\\ y=x_1\\ x_1(0)=0.1,x_2(0)=-0.1\end{array}\text{\hspace{1em}(28)}$$

function xdot = mVanderPol(t,x)
global theta;
xdot = [x(2);   (1-x(1)^2)*x(2)/theta-x(1)+sin(50*t) ];
end

其中，未知系统参数$\theta \sim Uniform([5,11])$为均匀分布，因此可以令相关的随机变量为$\xi \sim Uniform([-1,1])$，对应的多项式分解基函数为$[-1,1]$区间上的Legendre多项式
$$\begin{gathered} {\varphi }_0(\xi )=1,{\varphi }_1(\xi )=\xi ,{\varphi }_2(\xi )=\frac{3}{2}{\xi }^2-\frac{1}{2},{\varphi }_3(\xi )=\frac{5}{2}{\xi }^3-\frac{3}{2}\xi , \\ {\varphi }_4(\xi )=\frac{35}{8}{\xi }^4-\frac{30}{8}{\xi }^2+\frac{3}{8},{\varphi }_5(\xi )=\frac{63}{8}{\xi }^5-\frac{70}{8}{\xi }^3+\frac{15}{8}\xi ,...\text{\hspace{1em}(29)} \end{gathered}$$

这里的近似属于强近似，因此1阶的参数展开即可实现平均分布的近似，有($n_p=1,N_p=1$)
$$\theta (\xi )=\sum _{m=0}^{N_p}{\theta }_m={\theta }_0+{\theta }_1\xi =8+3\xi \text{\hspace{1em}(30)}$$

状态变量同样展开为多项式分解的形式($\alpha$在$n_p=1$的情况下为标量，$n_s=1$，$N_s=5$，$r=6$)

$$公式修改为文末图\text{\hspace{1em}(31)}$$

上式中，因为只有1个参数变量，因此多项式基$\Phi =\varphi$，$x_{j,0}\sim x_{j,5}$为多项式未知系数，需要利用随机配置法求解。

根据随机配置法，随机变量$\xi$在$[-1,1]$内按照均匀分布的原则选取$Q$个点${\mu }^{(1)},...,{\mu }^{(Q)}$，因此有$P({\mu }^{(i)})$、$\mathbb P(\mu^{(i)}) $、$\mathbb A$进而可以计算得到$\mathbb A^#$。
$$\begin{gathered} P^T({\mu }^{(i)})=[{\varphi }_0({\mu }^{(i)}){\varphi }_1({\mu }^{(i)}){\varphi }_2({\mu }^{(i)}){\varphi }_3({\mu }^{(i)}){\varphi }_4({\mu }^{(i)}){\varphi }_5({\mu }^{(i)})]\in {\mathbb{R}}^{1\times 6} \\ \mathbb{P}({\mu }^{(i)})=\left[\begin{array}{cc}{P}^T({\mu }^{(i)})& O\\ O& P^T({\mu }^{(i)})\end{array}\right]\in {\mathbb{R}}^{2\times 12},\mathbb{A}=\left[\begin{array}{c}\mathbb{P}({\mu }^{(1)})\\ ⋮\\ \mathbb{P}({\mu }^{(Q)})\end{array}\right]\in {\mathbb{R}}^{2Q\times 12} \\ \text{\hspace{1em}(32)} \end{gathered}$$

在具体实现的过程中，接下来需要运行$Q$次数值计算，从而得到某一时刻$t$的$z^{(i)}(t)$（全部时长的结果是一条空间轨迹），进而得到$Z$。在此基础之上，基于${\mathbb{A}}^{\#}$与$Z$，可以得到每一时刻的$X(t)$（元素为分解系数${\hat{x}}_{j,1\sim 5}(t)$）。进一步地，代入(26)即可实现参数$\xi$的辨识，经由(30)可以得到系统参数$\theta$的辨识结果。仿真结果如下，相关代码见附录。

5 小结

这里仅仅对包含1个参数的情况进行了实施仿真。通过不同辨识参数的调整，得到如下结论：（1）改变测量噪声没有对收敛效果产生明显影响，但是增大采样率则有助于收敛速度的加快以及准确度的提升；（2）在进行随机配置过程中，将随机改为均匀分布能够得到同样的辨识效果；（3）该方法在系统参数发生动态跳变的过程中并没有很好的跟踪真实值，且修改采样步长无明显变化，应当对优化过程中的梯度算法进行调整。下一步的工作将着眼于多参数的辨识情况，以及优化算法中梯度下降能力的提升。

附录

%% 基本参数
global theta;   % 全局变量，vdp参数
real_theta = 1; esti_theta = 0;  % 真实值，辨识值   
e1 = 2; e2 = 1.5;
% 起始时间、计算步长、终止时间
t0 = 0; tf= 6; sample_stp = 0.001; 
% 总时长采样点数
t = t0:sample_stp:tf;  sample_num = length(t);   esti_xi = zeros(1,sample_num);
% 概率多项式展开参数
ns = 2; np = 1;Ns = 5; R = factorial(Ns+np)/(factorial(Ns)*factorial(np));  
% 随机变量选择Q（Q>R）个点 (-1,1)之间均匀分布
Gamma = 0.5;   Q = 15;   mu = [-1:2/Q:1];%rand(Q, 1)*2 - 1;
% Q个点对应的状态变量初始条件
x0 = zeros(ns,Q); x0(1,:) = 0.1; x0(2,:) = 10;  
%% 导入实际的采样输出
[ real_y,set_theta ] = getSampleData(x0,t0,sample_stp,sample_num,real_theta,15);
subplot(211);plot(t,real_y);hold on;xlabel('时间(s)');ylabel('电流(A)');grid on;legend('测量曲线');
%% 构造正交基矩阵，这部分与时间无关
% 勒让德多项式
xi = sym('xi'); phi = sym('phi', [R 1]); phi(1) = 1; phi(2)=xi;
for k=1:R-2  
    phi(k+1+1) = expand((2*k+1)/(k+1)*xi*phi(k+1) - k/(k+1)*phi(k));        
end
% 正交基矩阵 P
Pt = sym('Pt', [1 R]); 
for k=1:R
    Pt(k) = phi(k);
end
% 正交基矩阵 PP
PP = sym('PP1', [ns R*ns]);
for k=1:ns
    PP(k,:)=0;
    PP(k,(k-1)*R+1:k*R) = Pt;
end
% CP导数与CP矩阵计算
C = [1 0]; CPP = C*PP; dCPP = diff(CPP,xi,1);
% 
DXX = zeros(R*ns,R*ns); DXY = zeros(R*ns,1);
% 系数矩阵 AA 直接代入Q个随机变量
AA = sym('AA', [Q*ns R*ns]); % 正交基矩阵 AA
for k=1:Q
    xi = mu(k);
    AA((k-1)*ns+1:k*ns,:) = subs(PP);
end
AAp = double(inv(transpose(AA)*AA)*transpose(AA));
xi = sym('xi');
%% 迭代计算
tic;str = ['正在运行中'];
hwait = waitbar(0,'请等待>>>>>>>>');
for k=1:1:sample_num     
   t1 = t0 + sample_stp;      
   waitbar(k/sample_num,hwait,str);
   for q=1:Q      % 随机配置Q个随机变量点
       theta = e1 + e2 * mu(q); % 参数变化
       % 四阶Runge-Kutta计算
       K1(:,q) = mVanderPol(t0,x0(:,q));
       K2(:,q) = mVanderPol(t0+sample_stp/2,x0(:,q)+sample_stp*K1(:,q)/2);
       K3(:,q) = mVanderPol(t0+sample_stp/2,x0(:,q)+sample_stp*K2(:,q)/2);
       K4(:,q) = mVanderPol(t0+sample_stp,x0(:,q)+sample_stp*K2(:,q));
       x1(:,q) = x0(:,q) + sample_stp/6*(K1(:,q)+2*K2(:,q)+2*K3(:,q)+K4(:,q));           
   end  
   % 构造当前时刻X矩阵  列向量，每Q个为1组，ns组
   Z(:,k) = reshape(x0,Q*ns,1);
   Xfactor(:,k) = AAp*Z(:,k);  
   % 极大似然辨识
   xi = esti_xi(:,k); dCPP_num = double(subs(dCPP,xi)); CPP_num = double(subs(CPP,xi));
   DXX = DXX + Xfactor(:,k)*transpose(Xfactor(:,k));
   DXY = DXY + Xfactor(:,k)*real_y(:,k);
   esti_xi(:,k+1) = esti_xi(:,k) - 0.5*Gamma*(dCPP_num*(DXX*transpose(CPP_num) - DXY)/(norm(DXY)+1));
   % 为下一次迭代
   t0 = t1; x0 = x1; 
   Xcoll(:,:,k) = x0'; % 不同配置点的模拟轨迹
end
close(hwait); toc% 注意必须添加close函数
 xi = sym('xi');esti_xi(1)=[];
 subplot(212);plot(t,(e1 + e2 *esti_xi));hold on;plot(t,set_theta);
 xlabel('时间(s)');ylabel('参数\theta');grid on;legend('辨识结果');