190 likes | 470 Views
SSI 方法. 该方法是基于线性离散时域状态空间模型的方法,以几何投影理论为基础,利用矩阵的 QR 分解(正交三角分解)和 SVD 分解(奇异值分解)技术,其本质是把将来输出的行空间投影到过去输出的行空间上,投影的结果保留了过去的全部信息,并用此来预测未来。. 线性的离散状态空间方程. 加入噪声: (过程噪声)和 (测量噪声)。 得下式:. 假定:白噪声, 和 满足. 对于白噪声激励的线性系统,其随机状态空间 模型可表示为:. 状态变量. 第 l 个测点在第 k 个采样间隔的输出向量. l 是输出的个数, n 是系统的阶数.
E N D
SSI方法 该方法是基于线性离散时域状态空间模型的方法,以几何投影理论为基础,利用矩阵的QR分解(正交三角分解)和SVD分解(奇异值分解)技术,其本质是把将来输出的行空间投影到过去输出的行空间上,投影的结果保留了过去的全部信息,并用此来预测未来。
线性的离散状态空间方程 加入噪声: (过程噪声)和 (测量噪声)。 得下式: 假定:白噪声, 和 满足
对于白噪声激励的线性系统,其随机状态空间 模型可表示为: 状态变量 第l个测点在第k个采样间隔的输出向量 l是输出的个数,n是系统的阶数 为输入噪声; 为观测噪声 A:系统状态矩阵 C:系统输出矩阵
将hankel矩阵的第i块行划分为 为非参考点输出,则
把将来输出的行空间投影到过去输出的行空间上:把将来输出的行空间投影到过去输出的行空间上:
投影矩阵 可分解为观测矩阵 和Kalman滤波状态序列的 的乘积
对投影矩阵进行奇异值分解得到: 随机子空间识别法通过取不同的加权矩阵将产生不同的识别法’ ① UPC算法②CVA 算法③PC 算法 由以上可得: ( ) 得到A、C最小二乘解: ( )
随机子空间方法步骤: 1、由实验输出数据构造系统的分块Hankel矩阵 2、计算特定分块Hankel矩阵的行空间投影即对 Hankel进行QR分解 3、对该投影进行奇异值分解从而得到广义的能观 矩阵 和状态序列Xi的非平稳Kalman滤波估计 4、由广义能观矩阵和非平稳Kalman滤波估计 确定系统矩阵A和输出矩阵C
Reference-based covariance-driven stochastic realization 构造Toeplitz矩阵 ( ) 由 所以
对Toeplitz矩阵进行奇异值分解 从而 又 故而
基于协方差的SSI • T0=Yf*Yp'; T1=Yf1*Yp'; • [U,S,V]=svd(T0); • U1=U(:,1:rank_svd);S1=S(1:rank_svd,1:rank_svd);V1=V(:,1:rank_svd); • Oi=U1*S1^0.5; C=S1^0.5*V1'; • %method 1 • Oi_1=Oi(1:(size(Oi,1)-size(data_sel,2)),:); • A=pinv(Oi_1)*Oi((size(data_sel,2)+1):end,:); • %method 2 • A1=pinv(Oi)*T1*pinv(C);
基于数据的SSI(基于参考点和一般的SSI共用一个程序)基于数据的SSI(基于参考点和一般的SSI共用一个程序) H=[Yp;Yf]; [Q,R]=qr(H'); Q=Q‘;R=R’;%QR分解 ---------------------------------------------------------------------------- • [U,S,V]=svd(Projection_now);%对Pi奇异值分解 • O_now=U1*S1^0.5; %Oi • X_now=pinv(O_now)*Projection_now; %Xi • O_truncated=O_now(1:(hankel_block_row-1)*sensor_num,:); %Oi-1 • X_next=pinv(O_truncated)*Projection_next;%Xi+1 ————————————————————————————————————— Projection_now=[R21;R31;R41]*Q1t; Projection_pre=[R41,R42]*[Q1t;Q2t];%两个投影矩阵 ————————————————————————————————————— A=X_next*pinv(X_now);
问题: • 1、投影分解为观测矩阵和kalman滤波状态序列的乘积的过程,投影的式子的来历 !
2、Hankel矩阵前的j(即Hankel矩阵的列数)的取值,j的底线?j与hankel矩阵的行数的关系,即j/i大小的影响。2、Hankel矩阵前的j(即Hankel矩阵的列数)的取值,j的底线?j与hankel矩阵的行数的关系,即j/i大小的影响。 3、运用协方差的方法中 ?
%method 1 Oi_1=Oi(1:(size(Oi,1)-size(data_sel,2)),:); A=pinv(Oi_1)*Oi((size(data_sel,2)+1):end,:); %method 2 A1=pinv(Oi)*T1*pinv(C); • 4、两种方法的比较 5、FFT幅值的问题 Y = fft(y,512); Pyy = Y.* conj(Y) / 512; f = 1000*(0:256)/512; plot(f,Pyy(1:257)) title('Frequency content of y') xlabel('frequency (Hz)')
程序界面的问题:文件保存的路径(mat和del的),绘图菜单下的命令程序界面的问题:文件保存的路径(mat和del的),绘图菜单下的命令