实验4 递推最小二乘法的实现
递推最小二乘法原理
递推最小二乘法原理递推最小二乘法（Recursive Least Squares, 简称RLS）是一种经典的参数估计方法，广泛应用于信号处理、通信系统、自适应滤波等领域。
它通过不断迭代更新参数，逐步逼近最优解，具有快速收敛、适应性强的特点。
本文将从最小二乘法出发，介绍递推最小二乘法的原理及其应用。
最小二乘法（Least Squares）是一种常见的参数估计方法，用于寻找一组参数，使得模型预测值与观测值之间的误差平方和最小。
对于线性模型，最小二乘法可以通过求解正规方程或者利用矩阵运算的方式得到最优参数。
然而，在实际应用中，数据通常是逐步到来的，因此需要一种能够动态更新参数的方法，于是递推最小二乘法应运而生。
递推最小二乘法的基本原理是利用递推的方式不断更新参数，以逼近最优解。
在每一时刻，根据当前的观测数据和先前的参数估计，通过递推公式计算出新的参数估计值，从而实现参数的动态更新。
这样的方法不仅能够适应数据的动态变化，还能够实现快速的收敛，适用于实时系统和非平稳环境下的参数估计。
递推最小二乘法的核心思想是利用指数加权的方式对历史数据进行处理，赋予近期数据更大的权重，从而更好地适应数据的变化。
通过引入遗忘因子（Forgetting Factor），可以控制历史数据对参数估计的影响程度，使得算法更具灵活性和适应性。
同时，递推最小二乘法还可以结合正交分解等技术，进一步提高计算效率和数值稳定性。
在实际应用中，递推最小二乘法被广泛应用于自适应滤波、信道均衡、系统辨识等领域。
例如，在自适应滤波中，递推最小二乘法可以根据接收信号的实际情况，动态调整滤波器的参数，实现信号的实时去噪和增强。
在通信系统中，递推最小二乘法可以用于自适应调制解调器的设计，提高系统的抗干扰能力和适应性。
此外，递推最小二乘法还被广泛应用于雷达跟踪、无线定位等领域，发挥着重要作用。
总之，递推最小二乘法作为一种经典的参数估计方法，具有快速收敛、适应性强的特点，在信号处理、通信系统、自适应滤波等领域有着重要的应用。
递推最小二乘法_协方差矩阵_概述说明以及解释
递推最小二乘法协方差矩阵概述说明以及解释1. 引言1.1 概述在统计学和计量经济学中，递推最小二乘法（Recursive Least Squares，简称RLS）是一种常用的参数估计方法。
它通过不断更新样本数据进行参数的估计，并且可以适用于非静态数据场景。
协方差矩阵是统计分析中重要的概念，它描述了变量之间的线性关系强度和方向，并且在许多领域具有广泛应用。
1.2 文章结构本文首先介绍递推最小二乘法的定义和原理，在此基础上详细解释算法的步骤以及其应用领域。
接着，我们将引入协方差矩阵的概念并介绍其计算方法，同时探讨了它在实际问题中所起到的作用和应用场景。
最后，我们将对递推最小二乘法与协方差矩阵之间的关系进行解释，并通过实例分析来说明它们如何相互影响。
1.3 目的本文旨在全面介绍递推最小二乘法和协方差矩阵，并深入探讨它们之间的联系。
通过对这两个概念及其应用的理解，我们可以更好地理解参数估计方法和变量间关系的描述与分析。
此外，我们还将展望相关领域未来可能的研究方向，以促进学术和实践的进一步发展。
2. 递推最小二乘法2.1 定义和原理:递推最小二乘法是一种用于估计线性模型参数的方法。
它可以通过历史数据的不断更新来逐步拟合模型，以使得估计值与观测值之间的误差达到最小化。
该方法可以被形式化地描述为以下步骤：1. 初始化模型参数的初始值。
2. 从历史数据中选择一个样本，并使用当前参数估计出该样本对应的输出值。
3. 计算该样本的预测误差。
4. 根据预测误差对参数进行调整，使得预测误差尽量减小。
5. 重复步骤2至4，直到所有样本都被处理过一遍，或者满足终止条件。
递推最小二乘法是基于最小二乘原理，即将真实观测值与模型预测值之间的差异平方求和并最小化这个目标函数。
通过迭代地更新参数，递推最小二乘法可以逐渐优化模型，并获得更准确的参数估计。
2.2 算法步骤:具体而言，在每次迭代中，递推最小二乘法按照以下步骤进行操作：1. 根据历史数据选择一个样本，并根据当前的参数估计出预测值。
递推最小二乘法推导
递推最小二乘法推导递推最小二乘法是一种经典的数学方法，用于解决数据拟合问题。
它通过最小化误差平方和的方法，寻找最佳的拟合曲线或平面，从而对数据进行预测和分析。
本文将详细介绍递推最小二乘法的原理和推导过程。
一、引言在现实生活和科学研究中，我们经常需要通过已知的数据来拟合一个函数，以便对未知的数据进行预测或分析。
而最小二乘法就是一种常用的数据拟合方法，它的基本思想是通过最小化误差的平方和，找到最佳的拟合函数。
二、最小二乘法的基本原理最小二乘法的基本原理是通过最小化残差平方和来确定拟合函数的参数。
残差指的是每个数据点的观测值与拟合函数预测值之间的差异。
最小二乘法的目标是找到使得残差平方和最小的参数值，从而得到最佳的拟合曲线或平面。
三、递推最小二乘法的推导过程递推最小二乘法是最小二乘法的一种改进方法，它能够更加高效地进行参数估计。
下面将结合一个简单的一元线性回归问题，来详细介绍递推最小二乘法的推导过程。
假设我们有一组样本数据(x₁, y₁), (x₂, y₂), …, (xₙ, yₙ)，需要找到一条直线y = ax + b 来拟合这些数据。
我们可以定义残差eᵢ= yᵢ- (axᵢ + b)，其中 eᵢ表示第 i 个数据点的残差。
我们的目标是通过最小化残差平方和来确定直线的参数a 和b。
即最小化损失函数 S = Σ(eᵢ²)。
我们需要计算一些中间变量，包括样本数据的均值xₙ和yₙ，以及样本数据的协方差 sₓy 和方差 sₓ²。
其中，xₙ = (x₁ + x₂ + … + xₙ) / n，yₙ = (y₁ + y₂ + … + yₙ) / n，sₓy = (Σ(xᵢ - xₙ)(yᵢ - yₙ)) / (n - 1)，sₓ² = (Σ(xᵢ - xₙ)²) / (n - 1)。
接下来，我们可以通过递推公式来更新参数 a 和 b 的估计值。
首先，我们初始化a₀和 b₀的估计值为0。
应用最小二乘一次完成法和递推最小二乘法算法地系统辨识
1最小二乘法的理论基础1.1最小二乘法设单输入单输出线性定长系统的差分方程表示为：其中δ(k)为服从N(0,1)的随机噪声,现分别测出n+N 个输出输入值y(1),y(2),…,y(n+N),u(1),u(2),…,u(n+N),则可写出N 个方程，写成向量－矩阵形式()()()()()()()()1201121n n y k a y k a y k a y k n b u k b u k b u k n k ξ=-------++-++-+L L ()()()()()()()()()()()()()()()()()()10111121222112n n y n y n y u n u y n y n y u n u y n N y n N y N u n N u N a n a n b n N b ξξξ+--+⎡⎤⎡⎤⎢⎥⎢⎥+-+-+⎢⎥⎢⎥=⨯⎢⎥⎢⎥⎢⎥⎢⎥+-+--+⎢⎥⎢⎥⎣⎦⎣⎦⎡⎤⎢⎥+⎡⎤⎢⎥⎢⎥⎢⎥+⎢⎥+⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥+⎢⎥⎣⎦⎢⎥⎢⎥⎣⎦LLL L M M MMLL M M M(4.1.1)则式(1.1.1)可写为 （4.1.2）式中：y 为N 维输出向量；ξ为N 为维噪声向量；θ为（2n+1）维参数向量；Φ为N ×（2n+1）测量矩阵。
因此，式（4.1.1）是一个含有（2n+1）个未知参数，由N 个方程组成的联立方程组。
11y θφφξ--=-在给定输出向量y 和测量矩阵Φ的条件下求参数θ的估计，这就是系统辨识问题。
设 表示 θ 的估计值，ŷ表示y 的最优估计，则有 （4.1.3） 式中：()()()10ˆˆ1ˆˆ2ˆˆ,ˆˆˆn n ay n a y n y b y n N b θ⎡⎤⎢⎥+⎡⎤⎢⎥⎢⎥⎢⎥+⎢⎥⎢⎥==⎢⎥⎢⎥⎢⎥⎢⎥+⎢⎥⎣⎦⎢⎥⎢⎥⎣⎦M M M ()()()()()()()()()()()()()()()()()()101122,,11112221n n a y n n y n a n y b y n N n N b y n y u n u y n y u n u y n N y N u n N u N ξξθξξφ⎡⎤⎢⎥++⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥++⎢⎥⎢⎥===⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥++⎢⎥⎢⎥⎣⎦⎣⎦⎢⎥⎢⎥⎣⎦--+⎡⎤⎢⎥-+-+⎢⎥=⎢⎥⎢⎥-+--+⎢⎥⎣⎦M M M M L LL L M MMLL y φθξ=+ˆθˆˆyθ=Φ设e(k)=y(k)- ŷ(k), e(k)称为残差，则有e=y- ŷ=y-Φθ 最小二乘估计要求残差的平方和最小，即按照指数函数（4.1.4）求Ｊ对 的偏导数并令其等于0可得：（4.1.5）由式（4.1.5）可得的 θ 最小二乘估计：（4.1.6）J 为极小值的充分条件是：即矩阵ΦT Φ为正定矩阵，或者说是非奇异的。
实验4 递推最小二乘法的实现
位矩阵，则经过若干次递推之后能得到较好的参数估计。
3
4．实验对象或参数 给定系统
y (k ) a1 y (k 1) a2 y ( k 2) b0u ( k ) b1u (k 1) b2u (k 2) (k ) （15）
即 n 2 。假设实际系统的参数为 a1 2 ， a2 1.3 ， b0 0.4 ， b1 0.88 ， b2 2.2 ， 但是不已知，即不可测。取 (k ) [0.1, 0.1] 的零均值白噪声。输入信号取为
u (k ) 1.5sin 0.2k
（16）
要求编制 MATLAB 程序，运用递推最小二乘法对这一系统的参数进行在线辨识，并 将辨识结果与实际参数进行对比。
5．程序框图
开始
程序初始化与赋初值
由先前输入输出生成 ΨN+1
用户输入初值设定常数c
由系统模型计算获得yN+1
设定初值θ0和P0
由递推最小二乘辨识公式得到 KN+1、PN+1、θN+1
y (1) u (n 1) u (1) y ( n) y (n 1) y (2) u (n 2) u (2) y (n N 1) y ( N ) u (n N ) u ( N )
（1）
其中 a1，a2…an，b0，b1，b2…bn 为待辨识的未知参数， ( k ) 是不相关随机序列。y 为系 统的输出， u 为系统的输入。 分别测出 n+N 个输出、 n+N 输入值 y(1)， y(2)…y(n+N)， u(1)， u(2)…u(n+N)，则可写出 N 个方程，具体写成矩阵形式，有
5.递推最小二乘法程序
5． 最小二乘参数的递推算法根据递推公式用循环实现递增。
如果满足所要求的精度，则跳出循环。
递推公式如下：)ˆ)1((ˆˆT 11θw K θθ⨯-++=++NN N N N z 1T 1))(1()(-++=N N N N N N w P w w P K))(1/()()()()1(T T N N N N N N N N N w P w P w w P P P +-=+源程序为：“ls.dsw ” 。
输入信号M 序列存在数据文件“data1.txt ”中 ，均值为0,方差为0.1的高斯白噪声存在数据文件“data2.txt ”中 。
估计结果存在数据文件“data3.txt ”中。
源程序：math.h"#include "stdio.h"#include "brmul.c"#include "brinv.c"int main(){FILE *fp1,*fp2,*fp;int u[600]; double v[600];int i, j, k, count;double y1[4][1],y[510],y2[1][1];double t1[4][1],t2[1][4],t3[1][1],t5[4][4],t6[4][4];double o[4][1], o1[500][4], w[4][1], w1[1][4];double p[4][4], km[4][1], er[4][1];double maxe;t3[0][0] = 0.0;y2[0][0] = 0.0;fp1 = fopen("data1.txt","r");fp2 = fopen("data2.txt","r");for( i = 0; i < 600; i++ ){fscanf(fp1,"%d",&u[i]);//输入信号M 序列fscanf(fp2,"%lf",&v[i]);//均值为0,方差为0.1的白噪声}fclose(fp1);fclose(fp2);y[0] = v[0];y[1] = v[1]; for( i = 2; i < 510; i++ )y[i] = 1.5*y[i-1]-0.7*y[i-2]+u[i-1]+0.5*u[i-2]+v[i];y1[0][0] = -1.5;y1[1][0] = 0.7;y1[2][0] = 1.0;y1[3][0] = 0.5;for( i = 0; i < 4; i++ ) o[i][0] = 0;for( i = 0; i < 4; i++ )for( j = 0; j < 4;j++ ){if( i == j ) p[i][j] = 1000000;else p[i][j] = 0;}for( k = 0; k < 500; k++ ){for( i = 0; i < 2; i++ ) w[i][0] = w1[0][i] = -y[1-i+k];for( i = 2; i < 4; i++ ) w[i][0] = w1[0][i] = u[3-i+k];brmul(p,w,4,4,1,t1);brmul(w1,p,1,4,4,t2);brmul(t2,w,1,4,1,t3);for( i = 0; i < 4; i++ ) km[i][0] = t1[i][0]/(1+t3[0][0]);brmul(w1,o,1,4,1,y2);for( i = 0; i < 4; i++ ) o[i][0] = o[i][0] + km[i][0]*(y[k+2] - y2[0][0]);brmul(t1,w1,4,1,4,t5);brmul(t5,p,4,4,4,t6);for( i = 0; i < 4; i++ )for( j = 0; j < 4; j++ ) p[i][j] = p[i][j] - t6[i][j]/(1+t3[0][0]);for ( i = 0; i < 4; i++ ) o1[k][i] = o[i][0];maxe = 0;for( i = 0; i < 4; i++ ){er[i][0] = fabs((o[i][0] - y1[i][0])/y1[i][0]);maxe = max( er[i][0], maxe );}count = k;if( maxe < 0.01 ) break;}fp = fopen("Data3.txt","w");fprintf(fp,"a0\t\ta1\t\tb0\t\tb1\n");for( i = 0; i < count+1; i++ ){fprintf(fp,"%f\t%f\t%f\t%f\n",o1[i][0],o1[i][1],o1[i][2],o1[i][3]);}fclose(fp);return 0; }。
递推最小二乘法
取二 Ｅ （ ） 为 的 合目 函 相 的 的 合目 函 ： ＝， ｊ 作 新 拟 标 数， 应 新 拟 标 数为 Ｄｗ
ｊＯ （ ） Ｊ 艺ｚ艺：（ 】 习回 ） Ｑ 加 司一 一 ｅ ｝ ｊ ＝ 。＋ Ｇ 阶川 ＝ Ｄ叭
权最小二乘拟合。 （３ ） （３ ０ 把式（ ．和（ ．） ３ ９ ３ １ 代入式（３ ３ ． ． （ ．） ３ １ 则有： ．
度函数。 下面仅举出第一种计算相关函数法的计算过程。 用这个方法进行相关分析辨识
系统模型的计算分为三个步骤：
１ 根 输 和 出 随 数 计 相 函 Ｒ ）互 关 数 ＝ ） ） 据 入 输 的 机 据 算自 关 数 ， 和 相 函 Ｒ令； ， 卜 ｒ
２ 求 ，） Ｒ（的 立 变 Ｓｊ 凡行） ） ＊： ｓ）傅 叶 换 ｘ 、 （和 ３ ， ｒ （） 。； ） ０ 和
把 宪Ｇ 做型，把一实值模值差 “ Ｏ 模值并每点测与型之 ｊ 州 （ ）
： Ｅ ． ＝ Ｊ ｊ（） ＝（ ， － ＋ ， ｊ ）Ｐ （ Ｑ． －
（ ．．１ ３３１ ）
称为拟合误 差，再取全部采样频率。上的 拟合误差ｓ 平方和厂作为 ， 的 拟合目 标函数。
，ｒｎ 、。、下不二不｝ ＇－厂 Ｊ ‘ Ｂｏｌ ｝． ｖ＼ 、 Ｇ ， ｌ ． Ｗ厂 ＿ ＿ ： Ｖｗ ｝ ｋ 一ＹＥ ｔ－Ｌ］１ Ｇ ） ＇ 二 ＇ Ｕ ｉ － －， ＇ ｉ
二 （ ＋ （ ｖ ｝ ｊｒ ） ｕ ｍｍ ）
（ ．．０ ３３ １ ）
现要定数， 以ｄ２ 由ｊ（所示频特 在 参ｂｉ， ｉ．， ’ ｊ 表的率性 确 ｐ． 及， 使 （ ． ｂｍ ． 得 劣 ） ． ． ｂ ｄｎ ． ０ ｄ
函 与 验 得 频 特 加 ＋ ） 接 如 把Ｐ 和Ｑ ， 做 测 而 数 实 求 的 率 性ｐ ） ．最 近。 果 阮） （） 实 值， Ｍ ．叫
递推最小二乘法
线性方程组的最优求解方法一．递推最小二乘法设线性方程组b Ax = (1)则有k b k =x :A ),(, (n k Λ,2,1=) (2)其中，[]kn k k a a a k ,,,:),(21Λ=A ，[]Tn x x x ,,,21Λ=x 。
设x :A ),()(k k f = (3)下面采用基于递推最小二乘法(RLS)的神经网络算法来训练权值向量x ，以获得线性方程组（1）的解x 。
由式（3）可知，若以)(k f 为神经网络输出，以k b 为神经网络训练样本，以x 为神经网络权值向量，[]kn k k a a a k ,,,:),(21Λ=A 为神经网络输入向量，则解线性方程组的神经网络模型如同1所示。
图1 神经网络模型采用RLS 算法训练神经网络权值向量x ，其算法如下： （1）神经网络输出：x :A ),()(k k f = (4)（2）误差函数：)()(k f b k e k -= (5)（3）性能指标：∑==n k k e J 12)(21 (6)（4）使min =J 的权值向量x ，即为所求的神经网络权值向量x ，这是一个多变量线性优化问题，为此，由0=∂∂xJ可得最小二乘递推法（RLS ）：]),([1k k k k k k b x :A Q x x -+=+ (7)),(),(1),(:A P :A :A P Q k k k T k T k k+= (8)k k k k P :A Q I P )],([1-=+ (9)()n k ,,2,1Λ=随机产生初始权值向量)1,(0n rand =x ，设nn ⨯∈=R I P α0（α是足够大的正数（一般取10610~10=α），nn ⨯∈R I 是单位矩阵），通过对样本数据训练，即可获得神经网络权值向量x ，此即为线性方程组（1）的解。
二．具有遗忘因子的递推最小二乘估计公式为：]),([1k k k k k k b x :A Q x x -+=+ (10)),(),(),(:A P :A :A P Q k k k Tk T k k+=λ (11) k k k k P :A Q I P )],([11-=+λ(12)式中，1:)],(:),([)(-=k A k A k TW P ，W 为加权对角阵：⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡=--10021λλλOn n W（nn ⨯∈=R I P α0，10610~10=α）。
递推的最小二乘法
1
说明： 公式①的物理意义
( N 1) WLS ( N ) 在以前观测的基础上对本次观测值的预测
T 0 y ( N 1) ( N 1) WLS ( N ) 预测误差，说明 WLS ( N )与实际参数（N）的偏差 T
进行递推，起动问题：
. 一般选 WLS (0)=0，相应的P(0） I = = 0
递推的最小二乘法
最小二乘估计：
LS ( ) Y
T 1 T
加权的最小二乘估计：
WLS ( W ) WY
T 1 T
Y:所有观测数据的全体，所以以上都是成批处理观测数据的一次完成算法，是离 线辨识方法。优点：辨识精度高;缺点：计算量大（特别是高阶矩阵求逆），对 计算机内存要求高，不能在线辨识！
从
WLS ( N )
到
WLS ( N 1)
的递推公式过程（略P191）
递推公式：
T ★ ① WLS ( N 1) WLS ( N ) L( N 1) y ( N 1) ( N 1) WLS ( N ) ② L( N 1) P( N 1) ( N 1) ( N 1)
1、递推的最小二乘法基本思想：
• 本次（新）的估计值
(k )
=上次（老）的估计值
(k 1)
+修正项
•
可以观察随着时间的推移，新的输入、输出信息不断增加的情况下，参数估计的变化情况，特别适用于 在线实时辨识。
设原先得到的参数估计用 WLS ( N ) 表示，则
WLS ( N ) ( T ( N )W ( N ) ( N )) 1 T ( N )W ( N )Y ( N )
递推最小二乘法
递推最小二乘法递推最小二乘法是一种避免精度损失的迭代计算方法，在最小二乘法的基础上加以改进，主要用于拟合复杂的数据，解决拟合时出现精度下降问题。
一、什么是递推最小二乘法递推最小二乘法是一种迭代计算方法，利用多项式曲线拟合曲线数据，对于某个曲线，只需要实施最小二乘法的迭代计算，而不需要考虑精度的损失。
递推最小二乘法的主要工作是根据给定的拟合曲线，把它拟合到数据集中，从而使数据集距离拟合曲线最小。
二、递推最小二乘法的原理递推最小二乘法的核心原理是，利用多项式拟合曲线，按照“最小二乘法”的原理，以当前拟合曲线为参照，不断进行前进和后退，以达到拟合曲线将数据集中的数据最佳拟合的目的。
这个最佳拟合目标就是实现拟合曲线与数据集之间的最小误差，其中，最小误差就是拟合曲线与实际数据集之间的最小差值。
递推最小二乘法的实现方式主要有两种，一种是基于递推的方式，另一种是基于函数的方式。
前者大致的实现方法是：先计算出多项式拟合曲线的每一个系数，然后再利用这些系数计算出多项式拟合曲线的最终拟合曲线，最后比较拟合曲线与实际数据集之间的实际差异，根据差异再调整系数，不断循环，直到最后拟合曲线与实际数据集之间的实际差异达到预期值为止。
函数的实现方式也很类似，只是在计算过程中，会使用函数的方式，将拟合曲线的系数表示为函数的形式，然后再比较拟合曲线与实际数据集之间的实际差异，根据差异再调整函数系数，最后实现拟合曲线与实际数据集之间的最小差异。
三、应用递推最小二乘法在实际应用中可以用来拟合复杂的数据曲线，以求得更好的拟合效果，解决拟合时出现精度下降问题。
它具有计算量小、运算简单、拟合结果较好的优点，因此在实际应用中得到了广泛的使用，比如在众多植物物种的遗传分析中，用递推最小二乘法来拟合植物的遗传规律，以获得更准确的估计结果；在探测地球大气层时，也可以用最小二乘法来拟合大气层中的湿度数据，以获取更加准确的湿度数据；在搜索引擎中，对查询结果也可以用最小二乘法拟合出来，以获得更准确的查询结果等等。
一看就会(废)的最小二乘法的推导
一看就会(废)的最小二乘法的推导最小二乘法可以用来做函数的拟合或者求函数极值。
在机器学习的回归模型中，我们经常使用最小二乘法。
我们先举一个小例子来走进最小二乘法。
某次实验得到了四个数据点\((x,y):(1,6)、(2,5)、(3,7)、(4,10)\) (下图中红色的点)。
我们希望找出一条与这四个点最匹配的直线 \(y = \theta_{1} + \theta_{2}x\) ，即找出在某种"最佳情况"下能够大致符合如下超定线性方程组的 \(\theta_{1}\) 和 \(\theta_{2}\) ，我们把四个点代入该直线方程可得:\[\theta_{1} + 1 \theta_{2} = 6\\\theta_{1}+2\theta_{2}=5\\ \theta_{1}+3\theta_{2}=7\\ \theta_{1}+4\theta_{2}=10 \]我们要求的是\(\theta_{1}\)和\(\theta_{2}\) 两个变量，但是这里列出了四个方程组，我们是无法求解的。
我们现在以向量空间的角度来解释为何无解，以及最小二乘法如何处理这种无解的情况。
\[Ax = b\\ \begin{bmatrix} 1 & 1 \\ 1 & 2 \\ 1 & 3 \\ 1 & 4 \\ \end{bmatrix} _{4\times2} \begin{bmatrix}\theta1\\ \theta2 \end{bmatrix}_{2\times1} =\begin{bmatrix} 6\\ 5\\ 7\\ 10 \end{bmatrix}_{4\times1} \]我们将四个方程组成的方程组写成矩阵形式。
矩阵A代表系数，\(x\)即待求的参数，\(b\)是每个方程对应的值。
从线性代数的角度来看，要判断\(Ax=b\)是否有解可以从向量空间角度来看。
递推最小二乘辨识
实验4 递推最小二乘法的实现实验报告哈尔滨工业大学航天学院控制科学与工程系专业：自动化班级：1040101姓名：日期：2013 年10 月23 日1．实验题目： 递推最小二乘法的实现 2．实验目的：(n a y k -()(n b u k n ξ+-+,,n a ，0,n b b 为待辨识的未知参数，(k ξ为系统的输出，为系统的输入。
分别测出n N +个输出、(3),(),(1),(2),()y y n N u u u n N ++，则可写出矩阵形式，有0()(1)(1)(1)((1)(2)(2)(2)((1)()()()(n n y n y u n u n a y n y u n u n b y n N y N u n N u N n b ξξξ⎢⎥--+⎤⎡⎤⎡⎢⎥⎢⎥⎢⎢⎥-+-+⎢⎥⎢+⎢⎥⎢⎥⎢⎢⎥⎢⎥⎢⎢⎥-+--+⎣⎦⎣⎢⎥⎢⎥⎣⎦（2）101)(1),,()()n n a n a b y n N n N b ξξθξξ⎡⎤⎢⎥++⎤⎡⎤⎢⎥⎥⎢⎥⎢⎥⎢⎥⎢⎥==⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥++⎣⎦⎣⎦⎢⎥⎢⎥⎣⎦，()(1)(1)(1)(1)(2)(2)(2)1)()()()y n y u n u n y u n u N y N u n N u N -+⎤⎥+-+⎥⎥⎥+--+⎦y θξ=Φ+ （维输出向量；ξ为N 维噪声向量；θ为21n +维参数向量；Φ测量矩阵。
为了尽量减小噪声ξ对θ估值的影响，应取2N n >+()T θ=ΦΦ为了实现实时控制，必须采用递推算法，这种辨识方法主要用于在线辨设已获得的观测数据长度为N ，将式N N Y θ=ΦN θ表示θ的最小二乘估计，则(TN Nθ=ΦΦ(T N N P =ΦΦT N N N P θ=Φ如果再获得一组新的观测值(u n N +11TN ψθ++=+(1),n N ++](1)(1)(1)y N u n N u N -++++）合并，并写成分块矩阵形式，可得T 111N N N N N Y y ξθξψ+++⎡⎤Φ⎡⎤⎤⎢⎥⎢⎥⎢⎥=+⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦（9）于是，类似地可得到新的参数估值1T T1T T T 1111N NN N N N N N N Y y θψψψ-+++++⎧⎫⎡⎤⎡⎤⎡⎤ΦΦΦ⎡⎪⎪⎢⎥⎢⎥⎢⎥⎪⎪⎢⎥=⎨⎬⎢⎥⎢⎥⎢⎥⎢⎥⎪⎪⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦⎣⎦⎪⎪⎩⎭ )T1T 1111N N N N N T N N N N Y P y Y y ψψ+++++⎡⎤Φ⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦Φ+)T T 111T 11N NN N N N N ψψψψ++-++⎪⎢⎥⎢⎥⎪⎬⎢⎥⎢⎥⎪⎢⎥⎢⎥⎣⎦⎣⎦⎪⎭+ )1T 11N N ψ-++ 与P 的递推关系式出发，导，最终可求得递推最小二乘法辨识公式：()T1111N N N N N N K y θθψθ++++=+- ()1T 11111N N N N N N K P P ψψψ-++++=+ ()1T T 11111N N N N N N N N P P P P ψψψψ-++++=-+ 为了进行递推计算，需要给出N P 和N θ的初值0P 和0θ。
递推最小二乘法python
递推最小二乘法python递推最小二乘法是一种用于线性回归问题的迭代算法，其原理是通过递归的方式不断迭代得到最小二乘解。
在Python中，可以使用numpy库中的linalg.lstsq函数实现递推最小二乘法。
首先，需要定义一个迭代函数，然后在该函数内部调用linalg.lstsq函数进行迭代，直至得到最小二乘解。
代码示例如下：```pythonimport numpy as npdef recursive_least_squares(X, y, alpha):# 初始化参数theta = np.zeros((X.shape[1], 1))P = np.eye(X.shape[1]) * 1e-6for i in range(X.shape[0]):x_i = X[i, :].reshape(-1, 1)y_i = y[i]# 更新参数K = P @ x_i / (alpha + x_i.T @ P @ x_i)theta = theta + K * (y_i - x_i.T @ theta)P = 1 / alpha * (P - K @ x_i.T @ P)return theta```其中，X和y分别表示自变量和因变量，alpha为学习率。
在调用该函数时，需要将X和y传入，并指定alpha的值，如下所示： ```pythonX = np.array([[1, 2], [3, 4], [5, 6]])y = np.array([3, 7, 11]).reshape(-1, 1)theta = recursive_least_squares(X, y, 0.1)print(theta)```运行结果为：```[[1.][1.]]```该结果表示得到的最小二乘解为theta=[1, 1]。
可以看出，通过递推最小二乘法得到的结果与正常的最小二乘解是一致的。
4递推增广最小二乘法
Loss function 0.00841826 and FPE 0.0121597 Sampling interval: 1 这里符号q 就是本书的符号z， q^-i 即为
z i
8
2.4 闭环系统辨识
以上谈的系统辨识和参数估计是针对开环系统的，然而自适应控制系统往往工作在
闭环状态下，所以，应考虑闭环给参数估计带来的新问题。
B( z 1 ) b0 b1z 1 b2 z 2 bnb z nb
C( z 1 ) 1 c1z 1 c2 z 2 cnc z nc
式（2-14）可变为
1 y(k ) B( z 1 )u(k d ) C ( z 1 ) 1 (k ) (k ) y (k ) 1 A ( z )
ˆ(k 1), ˆ(k 2)]; ˆ T (k ) [ y(k 1), y(k 2), u(k 2), u(k 3), 2.
T ˆ (k ) ˆ ˆ (k ) 1 (k ) P ( k 1) 3. 计算 K (k ) P (k 1) 1
0.0000
-1.0000
-1.0000
1.0000
1.0000
0.0000
0.0000
y
1.6800
2.4018
2.3039
1.4445
-0.1399
-1.9036
-0.9657
1.1687
5
解：根据2.3.4 递推增广最小二乘法公式，有参数估计算法：
ˆ(0) 等，并给u和y赋值； 1.初始化 P(0) 、
时刻的参数估计值。
用 u(i), y(i), i 1, 2, 来估计 ，须使 (k ) 中每个元与其他元不相关。由于 y (k 1)
递推最小二乘法[整理版]
![递推最小二乘法[整理版]](https://img.taocdn.com/s3/m/bb7c26feb9f67c1cfad6195f312b3169a451eafc.png)
线性方程组的最优求解方法一．递推最小二乘法设线性方程组b Ax = (1)则有k b k =x :A ),(, (n k ,2,1=) (2)其中，[]kn k k a a a k ,,,:),(21 =A ，[]Tn x x x ,,,21 =x 。
设x :A ),()(k k f = (3)下面采用基于递推最小二乘法(RLS)的神经网络算法来训练权值向量x ，以获得线性方程组（1）的解x 。
由式（3）可知，若以)(k f 为神经网络输出，以k b 为神经网络训练样本，以x 为神经网络权值向量，[]kn k k a a a k ,,,:),(21 =A 为神经网络输入向量，则解线性方程组的神经网络模型如同1所示。
图1 神经网络模型采用RLS 算法训练神经网络权值向量x ，其算法如下：（1）神经网络输出：x :A ),()(k k f = (4)（2）误差函数：)()(k f b k e k -= (5)（3）性能指标：∑==n k k e J 12)(21 (6)（4）使min =J 的权值向量x ，即为所求的神经网络权值向量x ，这是一个多变量线性优化问题，为此，由0=∂∂xJ可得最小二乘递推法（RLS ）：]),([1k k k k k k b x :A Q x x -+=+ (7)),(),(1),(:A P :A :A P Q k k k Tk T k k+= (8)k k k k P :A Q I P )],([1-=+ (9)()n k ,,2,1 =随机产生初始权值向量)1,(0n rand =x ，设n n ⨯∈=R I P α0（α是足够大的正数（一般取10610~10=α），n n ⨯∈R I 是单位矩阵），通过对样本数据训练，即可获得神经网络权值向量x ，此即为线性方程组（1）的解。
二．具有遗忘因子的递推最小二乘估计公式为：]),([1k k k k k k b x :A Q x x -+=+ (10)),(),(),(:A P :A :A P Q k k k Tk T k k+=λ (11)k k k k P :A Q I P )],([11-=+λ(12)式中，1:)],(:),([)(-=k A k A k T W P ，W 为加权对角阵：⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡=--10021λλλn n W（nn ⨯∈=R I P α0，10610~10=α）。
递推最小二乘法python
递推最小二乘法python递推最小二乘法是一种求解线性回归问题的方法，可以逐步地加入数据，并且不需要存储全部数据，省去了处理大量数据的繁琐过程。
以下是递推最小二乘法的Python实现。
步骤如下：1. 初始化参数：设初始的拟合直线为y=a*x+b，初始数据点个数为0，初始误差为0。
2. 逐步加入新的数据点：每次加入一个新的数据点(x,y)，首先计算它和当前拟合直线的误差e，即y-(a*x+b)，然后更新误差平方和SSE=SSE+e^2。
3. 更新拟合直线的参数：根据新加入的数据点，使用最小二乘法更新拟合直线的参数a和b。
4. 返回结果：重复步骤2和3，直至所有数据点都被加入。
最终得到的拟合直线就是最小二乘法的解。
代码实现如下：```import numpy as npclass RecursiveLS(object):def __init__(self):self.x = 0self.y = 0self.x2 = 0self.xy = 0self.n = 0self.a = 0self.b = 0def update(self, x, y):self.n += 1self.x += xself.y += yself.x2 += x**2self.xy += x*yif self.n > 1:denom = self.n*self.x2-self.x**2self.a = (self.n*self.xy-self.x*self.y)/denom self.b = (self.y*self.x2-self.x*self.xy)/denomdef predict(self, x):return self.a*x + self.b```其中，update()函数用于加入新的数据点，predict()函数用于预测新的x值对应的y值。
递推最小二乘和增广矩阵法
++ψψ
P ψ TT
LL++11 LL
L +1
)−1ψ
P T
L+1 L
4.3 递推最小二乘法
θˆL+1 = PL+1 (ΦLTYL +ψ L+1 yL+1 )
θˆL
=
(Φ
LTΦ
L
)
−
1Φ
T L
Y
L
=
PL
Φ
T L
Y
L
θˆL+1 = PL+1(PL-1θˆL +ψ L+1 yL+1)
PL+1
=
PL
−
PLψ
⎧Κ ⎪⎪⎨θˆ1
1
= P0ψ1 (1+ψ = θˆ0 + Κ1 ( y1
T 1
P0
ψ
1
−ψ 1Tθˆ0
)−1 )
⎪ ⎪⎩ P1
= P0
−
Κ
1ψ
T 1
P0
ψ
T 1
=
⎡⎣1
α (1)
α 2 (1)
α 3 (1)
α 4 (1)
α 5 (1)⎤⎦
y1 = CDα (1)
第 i 步: 计算 Ki, θˆi , Pi , 到设定次数截止。
4.3 递推最小二乘法
例: 采集气动参数试验数据
4.3 递推最小二乘法
CDα = f (α ), M = 0.3
设输出量: CDα = b0 +b1α +b2α2 +b3α3 +b4α4 +b5α5 +e
⎧⎪CDα ⎪⎪⎨CD#α
递推最小二乘法
1
最小二乘法辨识
回顾
考虑系统模型：
y(k) a1y(k 1) a2 y(k 2) L an y(k n)
b0u(k) b1u(k 1) L bnu(k n) (k)
2
最小二乘法：
a1
2n+1维 参数向量
N维输出向量
y
y(n 1)
y(n
2)
,
M
y(n
4
递推最小二乘 参数辨识算法
u(k)
y(k)
动态系统模型
反馈控制律
图4.1 动态系统递推最小二乘在线辨识过程原理图
2020/6/27
5
递推最小二乘法
令 PN (TN N )1 ，则递推最小二乘算法
) )
)
θN 1 θN K N 1 ynN 1 ψTN 1θN
K N1 PN ψN1
ψ
T N
1PN
ψN
1
1
PN 1
1
PN
PN ψN 1
ψTN 1PN ψN 1
1
ψTN
1PN
2ΦTNΦN
N
1
T N
1
1
1
2
ΦTN ΦN
1 1
2
ΦTN Φ N
1
N 1
1
T N 1
1
2
ΦTN ΦN
1
N 1
1
T N 1
1
2
ΦTN ΦN
1
1
2
ΦTN ΦN
1 1
)
N
ΦTN ΦN
1 N 1
2
T N
1
ΦTN ΦN
1
递推最小二乘法的应用
1 递推最小二乘法在电厂模型辨识中的应用电厂中大多数热工对象可以用一阶或二阶有迟延和非迟延的模型来表示，对这些模型中参数的辨识，递推最小二乘法是一种较好的方法。
本文以火电厂部分典型一阶模型为例子，借助于某电厂现场数据，分别对以下几种环节进行辨识。
1.1 一阶惯性环节火电厂中，来自锅炉的过热蒸汽，经高压调节汽门和导汽管道进入高压缸膨胀做功，高压缸的排汽回到锅炉再热器被重新加热，加热后的蒸汽经中压调节汽门进入中低压缸进一步膨胀做功，做功后的乏汽最终排入凝汽器变成凝结水，一般中压调节汽门的开度是高压调节汽门的3倍，即在机组负荷大于额定的30%或者滑压运行时，汽轮机的中压调门是完全开启的。
因此，在简化模型中，汽机侧调速器一级压力与机组有功功率可以简化为一阶惯性环节如下：111()1K G s T s =+将以上环节离散化，并写成差分方程的形式11111111()[(1)](1)(1)()/,/y k a y k b u k v k a T T T b K T T =--+-+-=-=其中 u 为调速器一级压力，y 为机组有功功率，()v k 为零均值方差为1的高斯白噪声。
该论文依据递推最小二乘法原理，借助 MATLAB 工具编写程序，设定合适的初始值和加权因子进行参数辨识，辨识结果为11ˆˆ 2.7547,0.9193a b ==-，由11ˆˆ,a b 可得到134.12K =，112.36s T =进而得到系统的传递函数为：134.12()12.361G s s =+下面运用递推最小二乘法对所得结果进行仿真：假设134.12K =，112.36sT =已知，采样时间为1s T =，则计算可得111111()/(112.36)/12.360.9191/34.12/12.36 2.7670a T T Tb K T T =-=-=-===用M 序列作为输入信号，得到输出信号，然后对参数进行辨识。
程序如下：%最小二乘的递推算法%Z(k+1)=0.9191*Z(k)+2.7670*u(k)+v(k)%======================================== clear clc%==========400 个产生 M 序列作为输入=============== x=[0 1 0 1 1 0 1 1 1]; %initial value n=403; %n 为脉冲数目 M=[]; %存放 M 序列 for i=1:ntemp=xor(x(4),x(9)); M(i)=x(9); for j=9:-1:2 x(j)=x(j -1); end x(1)=temp; end%===========产生均值为 0，方差为 1 的高斯白噪声========= v=randn(1,400);%==============产生观测序列 z================= z=zeros(402,1); z(1)=-1; for i=2:401z(i)=0.9191*z(i -1)+2.7670*M(i -1)+v(i -1); end %递推求解P=10*eye(2); %协方差阵Theta=zeros(2,401); %参数的估计值，存放中间过程估值Theta(:,1)=[0.001;0.001];%参数估计值初始值a1=0.001,b1=0.001u=0.99; %遗忘因子for i=2:401h=[-z(i-1);M(i-1)];K=P*h/(h'*P*h+u);%增益矩阵Theta(:,i)=Theta(:,i-1)+K*(z(i)-h'*Theta(:,i-1));P=(eye(2)-K*h')*P/u;end%==========================输出结果及作图============================= disp('参数a1 b1 的估计值:')Theta(:,401)a1=zeros(1,401);a1(1,:)=-0.9191;b1=zeros(1,401);b1(1,:)=2.7670;i=1:401;figure(1)plot(i,Theta(1,:),'m',i,a1,'g',i,Theta(2,:),'b',i,b1,'r');e_a1=Theta(1,:)-a1; % a1的误差e_b1=Theta(2,:)-b1; % b1的误差figure(2)plot(i,e_a1,'m',i,e_b1,'b'); %误差曲线图递推算法仿真曲线如图 1所示：图 1递推算法仿真曲线误差曲线如图 2所示：图 2误差曲线得到400个数据的a1 b1估计值为: 11ˆˆ0.9208, 2.6564a b =-=由估计值可得1ˆ12.63s T =， 1ˆ33.55K =由递推算法仿真曲线图 1可知，参数估计收敛速度较快，占用内存较少。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
将式（5）和式（8）合并，并写成分块矩阵形式，可得
YN N T y N 1 N 1 N 1
u (k ) 1.5sin 0.2k
（16）
要求编制 MATLAB 程序，运用递推最小二乘法对这一系统的参数进行在线辨识，并 将辨识结果与实际参数进行对比。
5．程序框图
开始
程序初始化与赋初值
由先前输入输出生成 ΨN+1
用户输入初值设定常数c
由系统模型计算获得yN+1
设定初值θ0和P0
由递推最小二乘辨识公式得到 KN+1、PN+1、θN+1
（12） （13） （14）
T T PN 1 PN PN N 1 1 N N N 1 N 1 P N 1 P 1
为了进行递推计算，需要给出 PN 和 N 的初值 P0 和 0 。 推荐取值方法为： 假定 0 0, P0 c I ， c 是充分大的常数， I 为 (2n 1) (2n 1) 单
a1 y (1) u (n 1) u (1) (n 1) y (n 1) y (n) y (n 2) y (n 1) a y (2) u (n 2) u (2) n (n 2) （2） b0 y (n N ) y (n N 1) y ( N ) u (n N ) u ( N ) (n N ) bn
于是，类似地可得到新的参数估值
（9）
2
T N N N 1 T T N 1 N 1 N PN 1 T N 1
实验 4
递推最小二乘法的实现
实验报告
专业：
U
自动化
班级： 姓名： 日期： 年 月 日
1．实验题目：
U
递推最小二乘法的实现
2．实验目的 熟悉并掌握递推最小二乘法的算法原理。
3．实验主要原理 给定系统
y (k ) a1 y (k 1) a2 y (k 2) an y (k n) b0u (k ) b1u ( k 1) bnu (k n) (k )
T
1
N T N 1
T
YN y N 1
PN 1 T N YN N 1 y N 1
式中
YN y N 1
（10）
T N N PN 1 T T N 1 N 1
（5）
（6）
令 PN N N
T
1
，则
N PN T N YN
如果再获得一组新的观测值 u (n N 1) 和 y (n N 1) ，则又增加一个方程
T yN 1 N 1 N 1
（7）
（8）
式中
yN 1 y (n N 1), N 1 (n N 1)
y (1) u (n 1) u (1) y ( n) y (n 1) y (2) u (n 2) u (2) y (n N 1) y ( N ) u (n N ) u ( N )
2
位矩阵，则经过若干次递推之后能得到较好的参数估计。
3
4．实验对象或参数 给定系统
y (k ) a1 y (k 1) a2 y ( k 2) b0u ( k ) b1u (k 1) b2u (k 2) (k ) （15）
即 n 2 。假设实际系统的参数为 a1 2 ， a2 1.3 ， b0 0.4 ， b1 0.88 ， b2 2.2 ， 但是不已知，即不可测。取 (k ) [0.1, 0.1] 的零均值白噪声。输入信号取为
T T N N N N 1 1 T PN1 N 1 N 1 1
1
1
（11）
应用矩阵求逆引理，从求得 PN 1 与 PN 的递推关系式出发，经过一系列的推导，最终 可求得递推最小二乘法辨识公式：
T N 1 N K N 1 yN 1 N 1 N T K N 1 PN N 1 1 N 1 P N N 1 1
则式（2）可写为
y
（3）
1
式中：y 为 N 维输出向量； 为 N 维噪声向量； 为 2n+1 维参数向量； 为 N×(2n+1) 测量矩阵。为了尽量减小噪声 对 估值的影响，应取 N＞2n+1，即方程数目大于未知数 数目。
的最小二乘估计为
(T ) 1 T y
5
rec_num=rec_num-1; end %-----------------------------------------------显示识别结果并与实际值进行比较 disp('** 运用递推最小二乘法对这一系统参数进行辨识，计算得到:') disp(' ') disp('相对于理论值[a1 a2 b0 b1 b2]=[2 1.3 0.4 0.88 2.2]，各次递推结果对应如 final_data 矩 阵所示：') final_data disp('实际值相对于理论值的误差如 final_percent 矩阵所示:') final_percent %-------------------------------------------------------------------相关绘图 k=0:1:rec_NUM-1; subplot(3,2,1) stairs(k,final_data(:,1)','k-');hold on line([0 rec_NUM-1],[T_NT(1,1) T_NT(1,1)]) axis([0 rec_NUM-1 -0.5 2.5]) title('a_1 收敛情况图形') subplot(3,2,3) stairs(k,final_data(:,2)','k-');hold on line([0 rec_NUM-1],[T_NT(2,1) T_NT(2,1)]) axis([0 rec_NUM-1 -0.2 1.5]) title('a_2 收敛情况图形') subplot(3,2,2) stairs(k,final_data(:,3)','k-');hold on line([0 rec_NUM-1],[T_NT(3,1) T_NT(3,1)]) axis([0 rec_NUM-1 0.1 1.4]) title('b_0 收敛情况图形') subplot(3,2,4) stairs(k,final_data(:,4)','k-');hold on line([0 rec_NUM-1],[T_NT(4,1) T_NT(4,1)]) axis([0 rec_NUM-1 0.6 1.3]) title('b_1 收敛情况图形') subplot(3,2,6) stairs(k,final_data(:,5)','k-');hold on line([0 rec_NUM-1],[T_NT(5,1) T_NT(5,1)]) axis([0 rec_NUM-1 -0.5 2.5]) title('b_2 收敛情况图形') end
（1）
其中 a1，a2…an，b0，b1，b2…bn 为待辨识的未知参数， ( k ) 是不相关随机序列。y 为系 统的输出， u 为系统的输入。 分别测出 n+N 个输出、 n+N 输入值 y(1)， y(2)…y(n+N)， u(1)， u(2)…u(n+N)，则可写出 N 个方程，具体写成矩阵形式，有
（4）
为了实现实时控制，必须采用递推算法，这种辨识方法主要用于在线辨识。设已获 得的观测数据长度为 N，将式（3）中的 y 、 和 分别用 YN , N , N 来代替，即
YN N N
用 N 表示 的最小二乘估计，则
T N T N N N YN 1
用户输入所需递推次数rec_num 将θN+1存入final_data θN+1与实际值得相对误差存入 final_percent
生成零均值白噪声ξ(k)
根据状态反馈控制律 计算增益系数矩阵K及控制U
满足递推次数？ 打印结果并绘制相关图形 是 结束
否பைடு நூலகம்
4
6．程序代码 function [] = Lab_4() %--------------------------------------------------------实验题目及初始化定义 disp(' ') disp('实验 4 递推最小二乘法的实现') disp(' ') disp('** 本实验对象系统的辨识模型为：') disp(' y(k)=-a1*y(k-1)-a2*y(k-2)+b0*u(k)+b1*u(k-1)+b2*u(k-2)+ξ(k).') disp(' 假定实际参数：a1=2, a2=1.3, b0=0.4，b1=0.88, b2=2.2. ') disp(' ') disp('** 本实验系统输入信号取为 u(k)=1.5*sin(0.2*k); 输出受到零均值白噪声ξ (k) ∈ [-0.1,0.1]污染.') disp(' ') a0=65539;M=2147483647;x=123456;b=10000; Y(1)=0;Y(2)=0;U(1)=1.5*sin(0.2*1);U(2)=1.5*sin(0.2*2);U(3)=1.5*sin(0.2*3); V=[];k_num=2;T_NT=[2;1.3;0.4;0.88;2.2]; c = input('** 为了进行推算，初值设定时，假定 P0=c*c*I，请输入充分大的常数 c = '); P_N0=c*c*eye(5); T_N0=zeros(5,1); rec_num = input(' 请输入所需递推步数 rec_num = '); disp(' ') rec_NUM = rec_num; final_data=zeros(rec_num,5); final_percent=zeros(rec_num,5); %--------------------------------------------------产生白噪声作为干扰信号ξ(k) for i=1:rec_num+2 x=mod(a0*x+b,M); V(i)=0.2*(x/M-0.5); end %------------------------------------------------------递推最小二乘法在线计算 while (rec_num~=0) HL(1)=-Y(k_num);HL(2)=-Y(k_num-1); HL(3)=U(k_num+1);HL(4)=U(k_num);HL(5)=U(k_num-1); Y(k_num+1)=HL*[2;1.3;0.4;0.88;2.2]+V(k_num+1); K_NN=P_N0*HL'*(inv(1+HL*P_N0*HL')); P_NN=P_N0-K_NN*HL*P_N0; T_NN=T_N0+K_NN*(Y(k_num+1)-HL*T_N0); final_data(k_num-1,:)=T_NN; T_ND1=abs(T_NN-T_NT); T_ND2=[T_ND1(1,1)/T_NT(1,1) T_ND1(2,1)/T_NT(2,1) T_ND1(3,1)/T_NT(3,1) T_ND1(4,1)/T_NT(4,1) T_ND1(5,1)/T_NT(5,1)]; final_percent(k_num-1,:)=T_ND2; P_N0=P_NN; T_N0=T_NN; k_num=k_num+1; K=[T_N0(2)-16 T_N0(1)-4]; U(k_num+1)=K*[Y(k_num);Y(k_num-1)]+[T_N0(3) T_N0(4) [1.5*sin(0.2*(k_num+1));1.5*sin(0.2*k_num);1.5*sin(0.2*(k_num-1))] T_N0(5)] *