这一周，我在学习递推最小二乘法RLS、递推增广最小二乘法RELS以及具有遗忘因子的递推最小二乘法FFRLS，在matlab中实现了相应的代码，并写了适用于三者的S函数，在Simulink中得到验证。

以下内容仅用于自我总结与分享交流，内容较多。另外，如有不妥，可在下方指出交流，多谢。

文章目录

摘要正文内容1. RELS算法简介1.1 基本公式1.2 计算步骤 2. RELS的代码实现2.1 产生数据2.2 RELS函数2.3 参数估计过程可视化 3. RELS的S函数及Simulink仿真3.1 Simulink仿真模型3.2 Random Number模块的问题3.3 S函数形式的RELS 4. 其他说明

摘要

仅详细介绍RELS，递推最小二乘法RLS和具有遗忘因子的FFRLS可以由RELS得到，前者可以看作是RELS的一类特例，后者则是在RLS的基础上加入了简单的遗忘因子。

RELS算法简介Matlab代码实现RELS的S函数及Simulink仿真其他说明

正文内容

1. RELS算法简介

1.1 基本公式

了解的同学这块内容可以直接跳过哈。

自校正控制的被控对象可用离散时间随机线性模型来描述，离散时间随机线性模型具有如下形式：

A ( z − 1 ) y ( k ) = z − d B ( z − 1 ) u ( k ) + C ( z − 1 ) ξ (1) A(z^{-1})y(k)=z^{-d}B(z^{-1})u(k)+C(z^{-1})\xi\tag{1} A(z−1)y(k)=z−dB(z−1)u(k)+C(z−1)ξ(1)考虑一般情况ARMA模型， C ( z − 1 ) ≠ 1 C(z^{-1})\neq1 C(z−1)​=1。式中， A ( z − 1 ) = 1 + a 1 z − 1 + a 2 z − 2 + ⋯ + a n a z − n a A(z^{-1})=1+a_1z^{-1}+a_2z^{-2}+\cdots+a_{n_a}z^{-n_a} A(z−1)=1+a1​z−1+a2​z−2+⋯+ana​​z−na​， B ( z − 1 ) = b 0 + b 1 z − 1 + b 2 z − 2 + ⋯ + b n b z − n b B(z^{-1})=b_0+b_1z^{-1}+b_2z^{-2}+\cdots+b_{n_b}z^{-n_b} B(z−1)=b0​+b1​z−1+b2​z−2+⋯+bnb​​z−nb​， C ( z − 1 ) = 1 + c 1 z − 1 + c 2 z − 2 + ⋯ + c n c z − n c C(z^{-1})=1+c_1z^{-1}+c_2z^{-2}+\cdots+c_{n_c}z^{-n_c} C(z−1)=1+c1​z−1+c2​z−2+⋯+cnc​​z−nc​。式(1)可变为

y ( k ) = [ 1 − A ( z − 1 ) ] y ( k ) + B ( z − 1 ) u ( k − d ) + [ C ( z − 1 ) − 1 ] ξ ( k ) + ξ ( k ) y(k)=[1-A(z^{-1})]y(k) + B(z^{-1})u(k-d)+[C(z^{-1})-1]\xi(k)+\xi(k) y(k)=[1−A(z−1)]y(k)+B(z−1)u(k−d)+[C(z−1)−1]ξ(k)+ξ(k)将其表示为

y ( k ) = φ T θ + ξ ( k ) (2) y(k)=\varphi^T\theta+\xi(k)\tag{2} y(k)=φTθ+ξ(k)(2)式中，

φ ( k ) = [ − y ( k − 1 ) , − y ( k − 2 ) , ⋯ , − y ( k − n a ) , u ( k − d ) , u ( k − d − 1 ) , ⋯ , u ( k − d − n b ) , ξ ( k − 1 ) , ξ ( k − 1 ) , ξ ( k − 2 ) , ⋯ , ξ ( k − n c ) ] (3) \varphi(k) = [-y(k-1), -y(k-2),\cdots,-y(k-n_a),u(k-d),u(k-d-1),\\ \cdots,u(k-d-n_b),\xi(k-1),\xi(k-1),\xi(k-2),\cdots,\xi(k-n_c)]\tag{3} φ(k)=[−y(k−1),−y(k−2),⋯,−y(k−na​),u(k−d),u(k−d−1),⋯,u(k−d−nb​),ξ(k−1),ξ(k−1),ξ(k−2),⋯,ξ(k−nc​)](3) θ T = [ a 1 , a 2 , ⋯ , a n a , b 0 , b 1 , ⋯ , b n b , c 1 , c 2 , ⋯ , c n c ] (4) \theta^T=[a_1,a_2,\cdots,a_{n_a},b_0,b_1,\cdots,b_{n_b},c_1,c_2,\cdots,c_{n_c}]\tag{4} θT=[a1​,a2​,⋯,ana​​,b0​,b1​,⋯,bnb​​,c1​,c2​,⋯,cnc​​](4)在 φ T ( k ) \varphi^T(k) φT(k)中，由于 ξ \xi ξ不可测，所以只能用其估计值 ξ ^ \hat{\xi} ξ^​来替换， ξ ^ ( k ) = y ( k ) − y ^ ( k ) = y ( k ) − φ T ^ ( k ) θ ^ (5) \hat\xi(k)=y(k)-\hat{y}(k)=y(k)-\hat{\varphi^T}(k)\hat{\theta}\tag{5} ξ^​(k)=y(k)−y^​(k)=y(k)−φT^​(k)θ^(5)基于递推最小二乘法RLS，RELS的递推公式如下

θ ^ ( k ) = θ ^ ( k − 1 ) + K ( k ) [ y ( k ) − φ T ^ ( k ) θ ^ ( k − 1 ) ] K ( k ) = P ( k − 1 ) φ ^ ( k ) [ 1 + φ T ^ ( k ) P ( k − 1 ) φ ^ ( k ) ] − 1 P ( k ) = [ I − K ( k ) φ T ^ ( k ) ] P ( k − 1 ) \hat{\theta}(k)=\hat{\theta}(k-1)+K(k)[y(k)-\hat{\varphi^T}(k)\hat{\theta}(k-1)]\\ K(k)=P(k-1)\hat{\varphi}(k)[1+\hat{\varphi^T}(k)P(k-1)\hat{\varphi}(k)]^{-1}\\ P(k)=[I-K(k)\hat{\varphi^T}(k)]P(k-1) θ^(k)=θ^(k−1)+K(k)[y(k)−φT^​(k)θ^(k−1)]K(k)=P(k−1)φ^​(k)[1+φT^​(k)P(k−1)φ^​(k)]−1P(k)=[I−K(k)φT^​(k)]P(k−1) ξ ( k ) \xi(k) ξ(k)的估计可用下列方法之一来进行：

预测形式为 ξ ^ ( k ) = y ( k ) − φ T ^ ( k ) θ ^ ( k − 1 ) \hat{\xi}(k)=y(k)-\hat{\varphi^T}(k)\hat{\theta}(k-1) ξ^​(k)=y(k)−φT^​(k)θ^(k−1)滤波形式为 ξ ^ ( k ) = y ( k ) − φ T ^ ( k ) θ ^ ( k ) \hat{\xi}(k)=y(k)-\hat{\varphi^T}(k)\hat{\theta}(k) ξ^​(k)=y(k)−φT^​(k)θ^(k)

并且当 k ≤ 1 k\leq1 k≤1时，有 ξ ^ ( k ) = 0 \hat{\xi}(k)=0 ξ^​(k)=0。

以上内容来自徐湘元的《自适应控制与预测控制》第二章系统辨识。 ξ ( k ) \xi(k) ξ(k)的计算是一个很纠结的地方。最开始我是看的柴天佑的《自适应控制》，这本书上面用的公式只有预测形式的，而仿真结果是噪声参数c收敛到0，我以为是代码问题，调试了两天，还是解决不了问题。最终，还是看到《自适应控制与预测控制》书上推荐使用滤波形式的公式，参数估计就对了。

1.2 计算步骤

递推增广最小二乘法RELS的计算步骤如下所示：

置初值 θ ^ ( 0 ) = 0 \hat{\theta}(0)=0 θ^(0)=0， P ( 0 ) = 0 P(0)=0 P(0)=0， ξ ^ ( 0 ) , ξ ^ ( 1 ) , ⋯ , ξ ^ ( 1 − n c ) \hat{\xi}(0), \hat{\xi}(1),\cdots, \hat{\xi}(1-n_c) ξ^​(0),ξ^​(1),⋯,ξ^​(1−nc​);按照式(3)构造增广数据向量 φ ^ ( k − 1 ) \hat{\varphi}(k-1) φ^​(k−1);根据RELS的递推公式，计算 K ( k ) K(k) K(k);计算 θ ( k ) \theta(k) θ(k)。计算 P ( k ) P(k) P(k)。计算滤波形式的 ξ ^ ( k ) \hat{\xi}(k) ξ^​(k)；递推一步 k → k + 1 k\to{k+1} k→k+1，并返回步骤2，构造新的数据向量。

2. RELS的代码实现

Matlab代码主要分为以下几个模块：

产生数据（包括输入输出及干扰）RELS函数参数估计过程可视化

2.1 产生数据

采用如下所示的带有有色噪声的被控对象模型

y ( k ) + 1.5 y ( k − 1 ) + 0.6 y ( k − 2 ) = 2 u ( k − 3 ) − 1.4 u ( k − 4 ) + ξ ( k ) + 1.2 ξ ( k − 1 ) + 0.85 ξ ( k − 2 ) y(k)+1.5y(k-1)+0.6y(k-2)=2u(k-3)-1.4u(k-4)+\xi(k)+1.2\xi(k-1)+0.85\xi(k-2) y(k)+1.5y(k−1)+0.6y(k−2)=2u(k−3)−1.4u(k−4)+ξ(k)+1.2ξ(k−1)+0.85ξ(k−2)

式中， ξ ( k ) \xi(k) ξ(k)为均值为0、方差为0.8的正态白噪声， u ( k ) u(k) u(k)为方差为1的正态白噪声。

产生输出数据y的代码主要参考了博文: 自适应控制——RELS、MVC(参数已知)、MVSTC(参数未知)

其实我最初是用simulink搭建模型去产生输出数据的，但是后来我发现这样不好，至于原因后文会讲到。为了在matlab中根据差分方程产生其中的y，借用一下这位朋友的代码片段了，另外我添加了数据可视化以及输入数据和干扰数据的保存。

%递推增广最小二乘法% 产生输入、干扰及输出数据L = 150;a=[1.5 0.6];b=[2,-1.4];c=[1,1.2,0.85];%系统输入输出表达式% 随机数生成器设置rng('default')s = rng(0);e = 0.8*randn(L,1);%精确干扰u = randn(L,1); % 输入 y = zeros(L,1); % 输出seeta_stand = [a,b,c(2:3)];e0 = zeros(1,2);u0 = zeros(1,4);y0 = zeros(1,2);fa = zeros(1,6);for i = 1:Lfa = [-y0,u0(3:4),e0]';y(i) = fa'*seeta_stand' + e(i);for ie = length(e0):-1:2e0(ie) = e0(ie-1);ende0(1) = e(i);for iu = length(u0):-1:2u0(iu) = u0(iu-1);endu0(1) = u(i);for iy = length(y0):-1:2y0(iy) = y0(iy-1);endy0(1) = y(i);endt = 0:0.01:(L-1)*0.01;%输入、干扰、输出可视化figuresubplot(211)plot(t,u,'b',t,e,'g')![在这里插入图片描述](https://img-/1115115234652.PNG?x-oss-process=image/watermark,type_ZmFuZ3poZW5naGVpdGk,shadow_10,text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L0x1Y2t5X0h1c3Q=,size_16,color_FFFFFF,t_70)legend('u','e')subplot(212)plot(t,y)legend('y')% 保存数据到mat文件t = 0:0.01:1.49;input_source = [t;u'];ksai_source = [t;e'];save('input.mat','input_source');save('ksai.mat', 'ksai_source');

输入输出数据以及干扰数据如下图所示

2.2 RELS函数

RELS函数主要接收的数据有阶次 n a n_a na​， n b n_b nb​， n c n_c nc​，输入 u u u，输出 y y y以及迭代次数 N N N。

function [theta_a, theta_b, theta_c] = myRELS( na, nb, nc, d, u, y, N)%递推增广最小二乘法% na：输出部分的阶次% nb：输入部分的阶次% nc：噪声部分的阶次% d：被控对象的时延% u：输入信号% y：输出信号% K：迭代次数% 设置初值len_b = nb + 1;len = na + len_b + nc;theta = zeros(len, 1);%B中第一个元素也要加入P = eye(len) * 10^6;% 以na=2 nb=1 nc=2 d=3为例% 构造数据向量的子数组y_part = zeros(na, 1); % [y(k-1); y(k-2)]u_part = zeros(len_b+d-1, 1); % [u(k-1); u(k-2); u(k-3); u(k-4)]%多出d-1用来实现右移功能，且第一个元素能保存当前时刻的ue_part = zeros(nc, 1); % 噪声估计初值 [e(k-1), e(k-2)]res = [];k = 1; %k是迭代/采样次数，从1开始，1对应t=0时刻的初始状态; y、u在零时刻有值while k <= N% 认为t=0(k=1)时刻就是第一个采样点时刻fai = [-y_part; u_part(d:end); e_part]; % 数据向量y_k = y(k); % 采样% 计算K(k)K = P * fai / (1 + fai' * P * fai); % 计算theta的最小二乘估计theta = theta + K * (y_k - fai' * theta); % theta定义时就是列向量res = [res; theta'];% 更新PP = ( eye(len) - K * fai') * P;% 计算误差，书上的theta应该为k，fai实际上的表述应该为fai(k)，而不是k-1ksai = y_k - fai' * theta;% 更新a部分，当前的值到下一时刻就是k-1时刻for i = length(y_part):-1:2y_part(i) = y_part(i-1);endy_part(1) = y_k;% 更新b部分for i = length(u_part):-1:2u_part(i) = u_part(i-1);endu_part(1) = u(k);% 更新ksai(c)部分for i = length(e_part):-1:2e_part(i) = e_part(i-1);ende_part(1) = ksai;% 更新迭代次数kk = k + 1;endtheta_a = res(:, 1:na);theta_b = res(:, na+1:end-(nc-1)-1);theta_c = res(:, end-(nc-1):end);end

2.3 参数估计过程可视化

通过RELS函数得到的参数估计过程可以通过plot绘制得到，代码我就不贴了，如图所示。注意，输入和干扰的随机序列不同时，参数估计也会不同，相应的曲线也会不同。

3. RELS的S函数及Simulink仿真

网上关于递推最小二乘法RLS、增广最小二乘法RELS的博客及代码有很多，都是在脚本中写的，几乎没有用S函数写的。使用脚本写的RELS仿真有这么几种方式：

生成输入 u u u和干扰 ξ \xi ξ的N个数据，其中N为迭代次数；然后在每次循环迭代中，计算输出 y ( k ) y(k) y(k)以及RELS的 K ( k ) 、 P ( k ) 、 θ ( k ) 、 ξ ( k ) K(k)、P(k)、\theta(k)、\xi(k) K(k)、P(k)、θ(k)、ξ(k)。生成输入 u u u、干扰 ξ \xi ξ、输出 y y y的N个数据；然后在每次循环迭代中，只计算RELS的 K ( k ) 、 P ( k ) 、 θ ( k ) 、 ξ ( k ) K(k)、P(k)、\theta(k)、\xi(k) K(k)、P(k)、θ(k)、ξ(k)。

其实，这两者在本质上没什么差别。实际情况下，在信号采集的中间时间点 t k ≤ t N t_k\leq{t_N} tk​≤tN​，我们是不可能获得输入 u u u、干扰 ξ \xi ξ、输出 y y y在整个过程的N个数据。所以在.m脚本中的仿真往往差那么一点在线递推的感觉。

然而，Simulink刚好能给你这种在线的感觉，在simulink中仿真更符合现实情况的工作方式。为此，我特地学了一下S函数，并且用S函数实现了“三合一”版本的RELS，所谓的“三合一”是指，只要给定合适的参数，该模块就可以执行RLS、RELS、FFRLS的功能。

3.1 Simulink仿真模型

这里，仍然使用前面的带有有色噪声的被控对象模型

y ( k ) + 1.5 y ( k − 1 ) + 0.6 y ( k − 2 ) = 2 u ( k − 3 ) − 1.4 u ( k − 4 ) + ξ ( k ) + 1.2 ξ ( k − 1 ) + 0.85 ξ ( k − 2 ) y(k)+1.5y(k-1)+0.6y(k-2)=2u(k-3)-1.4u(k-4)+\xi(k)+1.2\xi(k-1)+0.85\xi(k-2) y(k)+1.5y(k−1)+0.6y(k−2)=2u(k−3)−1.4u(k−4)+ξ(k)+1.2ξ(k−1)+0.85ξ(k−2)

先看一下搭建的Simulink仿真模型，如下图所示。输入及干扰从2.1中生成的 i n p u t . m a t input.mat input.mat和 k s a i . m a t ksai.mat ksai.mat文件中导入（干扰 ξ \xi ξ的读音不是ksai，而是xi，包涵一下这个小问题）。另外，设置求解器为固定步长的离散类型。由于输入和输出是150个数据，固定步长为0.01，所以仿真时间调整为1.49s。

看到这里，可能有朋友会问，你怎么不用Simulink自带的Random Number模块呢？那样不是更方便嘛。

其实，我最开始就是用的这个模块，用它仿真后得到的噪声参数 c c c的估计有问题，如下图所示。左图是simulink仿真出来的；右图是把simulink产生的输入干扰和输出数据导入到Matlab的工作空间，再利用2.2的myRELS函数得到的曲线。

从图中可以看到曲线相同，说明S函数my_online_RLS实现了函数my_RELS的功能。那么，问题究竟出在了哪里呢？（这个问题让我挣扎了好久，好久。。。）

3.2 Random Number模块的问题

突然，灵光一现，会不会是Simulink有问题呢？通过查看Simulink Random Number的文档，发现Random Number使用的随机数生成器是Matlab 4.0版本的。而且，要在Simulink中使用其他随机数生成器，需要在matlab中生成随机数，然后导入到Simulink中。那么，难道4.0版本的生成器就不可以吗？

The generator algorithm is identical to the one used in MATLAB Version 4.0 by the rand and randn functions. For details on the mcg16807 algorithm, see Choosing a Random Number Generator (MATLAB) in the MATLAB® documentation.To use other algorithms supported by MATLAB in a Simulink model, generate a stream of random numbers in MATLAB, and store the output as a .mat file. Use this .mat file as the random number input for your simulation.

然后，查看随机数的帮助文档。最终，发现在Matlab帮助文档Replace Discouraged Syntaxes of rand and randn一页指出：v4和v5的生成器不再推荐使用；所有的生成器除了’twister’以外都是有问题的。

The v4 and v5 generators are no longer recommended unless you are trying to exactly reproduce the random numbers generated in earlier versions of MATLAB.All of the generators except ‘twister’ are flawed.

因此，放弃使用Random Number！从而Simulink仿真模型就是3.1图示的情况了。

3.3 S函数形式的RELS

S函数形式的RELS函数需要输入的参数有阶次 n a n_a na​， n b n_b nb​， n c n_c nc​， d d d， α \alpha α。其中，遗忘因子 α = 1 \alpha=1 α=1就成了一般的RELS， n c = 0 n_c=0 nc​=0就成了一般的RLS。至于S函数的代码细节，我就不讲了，相信大家都能看懂。

function [sys,x0,str,ts,simStateCompliance] = my_online_RLS(t,x,u,flag,na,nb,nc,d,alpha)% 在线递推最小二乘法，可以考虑白噪声和有色噪声% u ：输入，包括系统的输入和输出% na：差分方程的A的阶次% nb：差分方程的B的阶次% nc：差分方程的C的阶次% d ：差分方程输入的时延% alpha：遗忘因子n = [na, nb, nc]; % 将三个阶次保存在数组中switch flagcase 0[sys,x0,str,ts,simStateCompliance]=mdlInitializeSizes(n,d);case 1sys=mdlDerivatives(t,x,u);case 2sys=mdlUpdate(t,x,u,n,d,alpha);case 3sys=mdlOutputs(t,x,u);case 4sys=mdlGetTimeOfNextVarHit(t,x,u);case 9sys=mdlTerminate(t,x,u);otherwiseDAStudio.error('Simulink:blocks:unhandledFlag', num2str(flag));endfunction [sys,x0,str,ts,simStateCompliance]=mdlInitializeSizes(n,d)global P K len %后面Update算法中要用到，定义为全局变量len = sum(n) + 1; % 参数theta的长度sizes = simsizes;sizes.NumContStates = 0;% na多一个保存当前时刻的位置；nb多一个常数项以及保存d个时刻的位置；nc多一个保存当前时刻的位置sizes.NumDiscStates = len * 2 + d - 1; %前面len + d - 1为fai；后面len为thetasizes.NumOutputs= len; % thetasizes.NumInputs= 2; % 最外层输入信号和输出结果sizes.DirFeedthrough = 0;sizes.NumSampleTimes = 1; % at least one sample time is neededsys = simsizes(sizes);x0 = zeros(sizes.NumDiscStates,1);str = [];ts = [-1 0];P = eye(len) * 1e6;% 初始化P，矩阵阶次同faiK = 0;simStateCompliance = 'UnknownSimState';function sys=mdlDerivatives(t,x,u)sys = [];function sys=mdlUpdate(t,x,u,n,d,alpha)na = n(1); nb = n(2); nc = n(3);input = u(1); y = u(2);global P K len% 这里的x就是数据向量fai，也是初始化时的x0y_part_rb = na; % y_part的右边界y_part = x(1:y_part_rb); % na个u_part_rb = y_part_rb + (nb+1) + d - 1; % u_part的右边界u_part = x(y_part_rb+1 : u_part_rb); % nb + 1 + d - 1个ksai_part_rb = u_part_rb + nc;ksai_part = x(u_part_rb+1 : ksai_part_rb); % 共nc个theta = x(ksai_part_rb + 1:end); %共len=na+nb+1+nc个 end-(len-1) = ksai_part_rb% 构造数据向量if nc > 0 fai = [y_part; u_part(d:end); ksai_part];elsefai = [y_part;u_part(d:end)];end% 计算当前时刻的K(k)，计算中用到的变量都是上一时刻的K = (P * fai) / (alpha + fai' * P * fai);% 计算当前时刻的参数估计值thetatheta = theta + K * (y - fai' * theta); % 这里计算的theta和后面的fai要一起更新到状态量当中% 计算当前时刻的PP = (eye(len) - K * fai') * P / alpha;% 计算噪声估计值ksai = y - fai' * theta;% 更新状态量中输出y的部分% t = 0.01时(k=1)，y(k-1)=y(0)进入fai(k-1),y(1)保存在y_part(0)for i = length(y_part):-1:2y_part(i) = y_part(i-1);endy_part(1) = -y; % 记录当前时刻的-y% 更新状态量中输入u的部分for i = length(u_part):-1:2u_part(i) = u_part(i-1);endu_part(1) = input; % 记录当前时刻的ufai_u = u_part(d:end);% 更新状态量中噪声ksai的部分， 计算当前时刻的噪声ksai，if nc > 0for i = length(ksai_part):-1:2ksai_part(i) = ksai_part(i-1);endksai_part(1) = ksai; % 记录当前时刻的ksai,下一时刻用end% fai，theta更新到状态量中x = [y_part;u_part;ksai_part;theta]; %[fai;theta]是错的，x是状态量，包括所有延时的量sys = x;function sys=mdlOutputs(t,x,u)global lentheta = x(end-(len-1):end);sys = theta;function sys=mdlGetTimeOfNextVarHit(t,x,u)sys = [];function sys=mdlTerminate(t,x,u)sys = [];

最终，Simulink仿真得出的参数估计与脚本得出的结果相同。

4. 其他说明

终于码完了，第一次写这么长的文章，把我这一周的工作做个总结，跟大家一起分享一下，希望能对你们有点帮助。

最后，附上所有文件的下载链接，其中包括所有上述文件以及myplot.m绘图函数。另外，所有文件都是在MatlabRa版本编写的。