1. 这次要做的事情用MATLAB把最小二乘递推辨识跑通做控制系统或者信号处理的朋友早晚都会碰到一个场景系统结构已知但里面的参数拿不准。比如一个直流电机传递函数结构是K / (Ts 1)K和T到底是多少你手头只有输入输出数据没有内部细节。这类问题就是系统辨识而最小二乘估计和它的递推算法RLS是最经典也最常用的解决工具之一。我在本科做课程设计、后来带学生做项目时经常看到有人一上来就抄一段递推最小二乘代码跑完发现参数不收敛或者收敛到错误的值然后就开始怀疑算法有问题。绝大多数情况下问题出在输入信号设计上——你拿什么信号去激励系统直接决定了辨识结果的品质。所以这次我用一个完整的案例把M序列生成→系统建模→递推最小二乘估计→误差分析这条链路全部走一遍代码直接给全每个关键步骤都讲清楚为什么这么做。这篇内容适合三类人正在做系统辨识课程设计的学生、需要在线估计参数来调试控制器的工程师、以及想搞清楚RLS和M序列到底怎么回事的入门者。不需要太多前置知识会基本的MATLAB语法就能跟上。2. 为什么偏偏是M序列和递推最小二乘2.1 最小二乘估计的核心逻辑要理解递推最小二乘先把非递推版本说清楚。假设系统是单输入单输出SISO的离散时间模型y(k) -a1*y(k-1) - a2*y(k-2) - ... - an*y(k-n) b1*u(k-1) b2*u(k-2) ... bm*u(k-m) e(k)这里的未知参数是a1, a2, ..., an, b1, b2, ..., bm。如果我们把模型写成矩阵形式Y Φ * θ其中 Φ 是由输入输出历史数据构成的矩阵θ 是未知参数向量。最小二乘的估计式是θ̂ (ΦᵀΦ)⁻¹ ΦᵀY这个式子看着简单但它需要把所有数据一次性装进 Φ 矩阵然后做矩阵求逆。数据量一大内存占用高、计算量大而且没法在线更新。实际应用中数据是不断采集进来的系统的参数还可能是时变的这时候就需要递推形式每来一组新数据用上一时刻的估计值加上一个修正量得到新估计值。递推算法的每一次更新只涉及向量和矩阵乘法复杂度固定这才是它可以长期在线运行的关键。2.2 遗忘因子的作用RLS基本递推公式网上到处都是但很多人不知道里面的遗忘因子 λ 是从哪里冒出来的更不知道它为什么要存在。标准最小二乘假设所有历史数据权重相同这适合参数不变的时不变系统。但真实系统往往有慢漂移比如电机温度升高导致电阻变大、机械磨损导致摩擦系数变化。如果旧数据和当前数据权重一样旧数据会拖新数据的后腿参数就追不上真实变化。遗忘因子的思想是离当前时刻越远的数据权重越小最新的数据起主导作用。数学实现就是在协方差矩阵的递推更新里乘一个 λ通常取0.95~0.99让历史信息随时间指数衰减。λ怎么选我自己的经验是λ越接近1算法越平稳但跟踪变化越慢λ越小跟踪越快但噪声敏感度也越高。对绝大多数课程设计和一般性工程测试λ取0.98左右比较稳健既有不错的跟踪能力又不会因为噪声剧烈抖动。2.3 为什么参数辨识需要M序列这是很多人忽略的关键点。RLS能不能有效辨识参数取决于矩阵ΦᵀΦ是否可逆、是否“足够良态”。如果输入信号设计得不好比如一直用常数信号、正弦信号或者简单的阶跃信号ΦᵀΦ可能亏秩或者条件数极大导致参数辨识结果严重失真甚至完全无法收敛。数学上有个说法叫持续激励条件persistently exciting输入信号必须包含足够丰富的频率成分才能充分激发系统的所有模态。M序列也就是最长线性反馈移位寄存器序列m序列本质上是一个二值伪随机信号它在一个周期内的自相关函数接近理想的脉冲函数——这正好满足持续激励条件。可以这么理解滤波器的频率覆盖足够宽才能把系统各个频段的响应都看清。M序列的频率分量均匀分布在较宽频带内用通俗的话说就是“白得很”所以它天然适合做辨识激励信号。3. M序列生成从原理到MATLAB实现3.1 m序列的生成原理m序列由一个n级线性反馈移位寄存器LFSR产生。结构上就是n个寄存器串联根据反馈多项式决定哪几级寄存器的值做异或运算后反馈到第一级。这个反馈多项式必须是本原多项式才能保证序列周期最大即2ⁿ - 1。比如4级寄存器本原多项式f(x) x⁴ x 1对应寄存器第3级和第4级做异或反馈。寄存器初始状态不能全0否则序列永远出不来。以[1 0 0 0]为初始状态每来一个时钟脉冲寄存器移位一次状态变化周期是15。注意n4时周期是2⁴-115n8时周期是255n12时是4095。实际做辨识时M序列的长度一般要覆盖系统过渡过程的好几倍否则激励不够充分约束条件无法满足。3.2 用移位寄存器实现M序列下面给出最贴近硬件原理的生成代码。我用的是MATLAB但逻辑完全可以用C/Verilog搬走。% 参数设置 n 6; % 寄存器级数周期为 63 N 2000; % 生成点数根据实验时间决定要大于多个周期 init_state [1 0 0 0 0 0];% 初始状态不能全零 % 本原多项式反馈抽头十进制表示n6 常用 [6 5] 表示 x^6 x^5 1 % 这里用抽头位置下标反馈位为第5级和第6级 tap1 5; % 第一个反馈抽头位置 tap2 6; % 第二个反馈抽头位置 reg init_state; m_seq zeros(1, N); for k 1:N out reg(end); % 输出最后一级 m_seq(k) out; feedback xor(reg(tap1), reg(tap2)); % 按本原多项式做异或 reg [feedback, reg(1:end-1)]; % 右移一位反馈值进入第一级 end运行之后得到的m_seq是0/1序列。注意它的均值不是0而是(2ⁿ⁻¹ - 1)/(2ⁿ - 1)直流分量有点大。这个直流分量不影响参数辨识但在后面把输入送入系统时为了贴近零均值激励条件通常把0/1映射成-1/1或-a/a。从0/1映射到双极性信号u 2 * m_seq - 1; % 0 - -1, 1 - 1如果想要幅值可调的激励信号再加一个系数u_amp a * u。3.3 直接调用MATLAB自带函数的方法MATLAB通信工具箱里有个现成函数comm.PNSequence可以一步生成m序列。之前上面的移位寄存器代码太原始了实际用自带函数更省事、也更不容易出错尤其是查本原多项式时交给工具箱处理就行。% 创建 PN 序列生成器对象 pnGen comm.PNSequence(... Polynomial, [6 5 0], ... % 本原多项式系数[6 5 0] 表示 x^6 x^5 1 InitialConditions, [1 0 0 0 0 0], ... SamplesPerFrame, N); m_seq pnGen(); % 生成 N 点序列用自带函数可以省去手动查本原多项式的麻烦。查多项式的时候我习惯直接记住几个常用的3级用[3 2 0]4级用[4 3 0]5级用[5 3 0]6级用[6 5 0]7级用[7 6 0]。这些配对经过验证周期分别是7、15、31、63、127。如果只有基本的MATLAB环境没有通信工具箱就用之前的移位寄存器手写代码。3.4 M序列参数怎么选这里给一个我自己反复验证过的选择逻辑参数选择原则我的建议寄存器级数 n周期需要大于系统过渡过程长度的3~5倍n6~10对应周期63~1023采样时间 Ts约为系统最快时间常数的1/10~1/5Ts T_min/5总点数 N至少为周期的5~10倍N 2000~10000连续采集信号幅值 u0线性工作区范围太大会激励非线性太小信噪比差先仿真试验取系统稳态输出的10%~30%一句话总结周期要长采样要快幅值要适中。这三个条件缺一不可不然后面的辨识效果肯定打折扣。注意M序列本身是周期的直接做估计时如果数据长度正好是周期的整数倍矩阵的性质会很好但也不要强制做成整数倍多取一段数据反而能降低边界效应。4. 递推最小二乘算法实现4.1 从批处理到递推的数学跳转批处理最小二乘的估计式是θ̂ (ΦᵀΦ)⁻¹ ΦᵀY递推的思路是在 k 时刻利用 k-1 时刻的估计值和新到的输入输出数据{u(k), y(k)}修正参数。标准RLS的三条关键递推方程是增益向量 K(k) P(k-1)·φ(k) / [λ φᵀ(k)·P(k-1)·φ(k)] 协方差矩阵更新 P(k) [I - K(k)·φᵀ(k)]·P(k-1) / λ 参数更新 θ̂(k) θ̂(k-1) K(k)·[y(k) - φᵀ(k)·θ̂(k-1)]其中φ(k)是回归向量由当前的输入输出历史数据构成。这三条公式不用死记理解逻辑更重要预测误差 实际输出 - 模型预测输出修正量 增益向量 × 预测误差。增益向量和协方差矩阵直接关联决定参数往哪个方向调整、调整幅度多大。关于初始化θ̂(0)一般取零向量或很小随机数P(0)取c·Ic 是较大常数如100~1000。协方差矩阵初始值反映了初始估计的不确定性c取得越大代表对初始参数越没信心算法前几步调整幅度就越大收敛更快c取得太小收敛会变慢甚至卡在初始值附近。4.2 完整MATLAB实现代码下面是我在一个二阶系统辨识任务里完整跑过的代码框架结构清晰可以直接改参数复用%% 参数配置 % 真实系统参数模拟对象 a1_true -1.5; a2_true 0.7; b1_true 0.1; b2_true 0.05; % 采样时间和信号 Ts 0.1; % 采样时间 N 3000; % 数据点数 % 遗忘因子 lambda 0.98; %% 生成M序列激励信号 pnGen comm.PNSequence(... Polynomial, [6 5 0], ... InitialConditions, [1 0 0 0 0 0], ... SamplesPerFrame, N); m_seq pnGen(); u 2 * double(m_seq) - 1; % 0/1 映射到 -1/1 %% 模拟真实系统输出加入噪声 % 系统模型: y(k) a1*y(k-1) a2*y(k-2) b1*u(k-1) b2*u(k-2) e(k) y zeros(1, N); e 0.05 * randn(1, N); % 高斯白噪声标准差0.05 for k 3:N y(k) -a1_true * y(k-1) - a2_true * y(k-2) ... b1_true * u(k-1) b2_true * u(k-2) e(k); end %% RLS递推辨识 theta_hat [0; 0; 0; 0]; % 参数初始估计 P 100 * eye(4); % 协方差矩阵初始值 theta_history zeros(4, N); % 保存每次迭代的参数估计值 for k 3:N % 构造回归向量[-y(k-1), -y(k-2), u(k-1), u(k-2)] phi [-y(k-1); -y(k-2); u(k-1); u(k-2)]; % RLS核心三步 K P * phi / (lambda phi * P * phi); error y(k) - phi * theta_hat; theta_hat theta_hat K * error; P (eye(4) - K * phi) * P / lambda; % 记录历史 theta_history(:, k) theta_hat; end %% 绘图参数估计轨迹 figure; plot(1:N, theta_history(1,:), b-, LineWidth, 1.2); hold on; plot(1:N, theta_history(2,:), r-, LineWidth, 1.2); plot(1:N, theta_history(3,:), g-, LineWidth, 1.2); plot(1:N, theta_history(4,:), m-, LineWidth, 1.2); yline(a1_true, b--); yline(a2_true, r--); yline(b1_true, g--); yline(b2_true, m--); xlabel(迭代次数 k); ylabel(参数估计值); legend(a1估计, a2估计, b1估计, b2估计, 真值); grid on; title(RLS参数估计收敛过程);这段代码跑完你能清楚地看到参数从初始值逐步逼近真值的过程。M序列在这里起到了关键的激励作用如果换成常数信号theta_history会明显不正常要么不动要么发散。4.3 参数辨识结果怎么评估光看参数轨迹还不够还需要估计误差和输出拟合两个指标来判断辨识质量。参数估计误差定义为误差 ||θ̂ - θ_true|| / ||θ_true||MATLAB里可以这样算% 参数估计误差相对2-范数 theta_true [a1_true; a2_true; b1_true; b2_true]; est_error sqrt(sum((theta_hat - theta_true).^2)) / sqrt(sum(theta_true.^2)); fprintf(参数估计相对误差: %.4f\n, est_error); % 用估计参数仿真模型输出与真实输出对比 y_hat zeros(1, N); for k 3:N y_hat(k) -theta_hat(1)*y(k-1) - theta_hat(2)*y(k-2) ... theta_hat(3)*u(k-1) theta_hat(4)*u(k-2); end figure; plot(y, k-); hold on; plot(y_hat, r--, LineWidth, 1.2); xlabel(k); ylabel(y); legend(实际输出, 辨识模型输出); grid on;预测误差曲线也要画出来就是每次迭代的error变量。如果算法收敛这个误差应该逐步减小并进入平稳状态。我习惯把误差画成半对数图纵轴用semilogy收敛趋势看得更直观。再强调一点估计误差和预测误差不是一回事。预测误差小不代表参数可靠因为系统可能过拟合或者信号激励不充分时预测误差也可以很小但参数远离真值。所以参数误差才是最终要看的指标。5. 踩过的坑和排查思路5.1 参数不收敛先怀疑激励信号如果你跑完RLS发现参数曲线一直振荡或者缓慢漂移、始终到不了真值附近第一件事不是调遗忘因子也不是改协方差初始化而是回到激励信号上检查。经验法则是计算输入信号的频谱确认其在系统带宽内有足够的功率覆盖。用pwelch或fft快速看频谱。M序列的频谱是离散线谱谱线间隔是1/(周期×Ts)覆盖带宽到约0.443/Ts等效带宽。如果系统带宽接近或超过这个覆盖范围就加大n或者减小Ts。我之前就栽过一回一个固有频率较高的二阶系统采样时间取太大M序列的频谱覆盖不到系统谐振点结果参数辨识出来的自然频率偏差超过30%。后来把Ts减小到原来的1/5参数估计精度马上提升一个量级。5.2 协方差矩阵变病态或爆炸RLS在高信噪比下收敛很快但协方差矩阵P(k)会越来越小最后可能出现数值病态导致增益向量更新异常。另一种情况是遗忘因子过小或噪声较大P矩阵可能突然增大参数估计剧烈跳变。排查手段有三个观察P矩阵对角元的变化如果数值持续下降后突然出现NaN或Inf说明数值稳定性出问题了对P矩阵施加上下限约束比如保证对角元不小于某个小阈值如1e-10采用平方根滤波QR分解形式数值稳定性更好但代码复杂度会上升。课程设计阶段我建议先用最简单的设置把流程跑通不要一上来就追求带遗忘因子平方根滤波的组合否则一个问题叠一个问题很难排查。5.3 噪声大小对辨识结果的影响我在前面的仿真里噪声标准差设了0.05这是个比较理想的值。实际系统里噪声更大比如信噪比低于10dB时RLS的估计方差会明显变大。有两个应对方向一是用增广最小二乘RELS把噪声模型也纳入参数向量对有色噪声特别有效二是增加数据长度因为最小二乘有一致性数据越长噪声的影响越能平均掉。还有一个细节仿真时加的噪声最好用固定种子的randn比如rng(2024)保证结果可复现。不然每次跑结果差异很大你都不知道是算法问题还是随机性问题。5.4 初始值设置不当导致前期振荡很多资料建议θ(0)0、P(0)1000*I这个组合在小噪声情况下收敛很快但有一个副作用最开始几步的修正量特别大参数轨迹会有明显跳变。如果系统本身非线性比较强这种大跳变可能让模型暂时偏离线性区。我的做法是θ(0)用一个粗略的估计值比如根据物理规律估算P(0)取10~100倍的适度值。这样收敛速度略微下降但全过程更平滑尤其适合输出信号幅值约束较强的场合。5.5 关于递推算法的效率问题RLS每次迭代的计算量约O(n²)n是参数个数。对4参数模型3000点数据秒算完MATLAB毫无压力。但如果参数个数上百、数据量百万级别就要考虑快速RLSFast RLS或LMS类算法了。实际项目中只要不是嵌入式实时环境RLS完全够用。真要上单片机注意矩阵运算用浮点P矩阵维度大的话内存和计算量可能吃紧建议做定点化改造。6. 怎样把这个实验扩展到自己的项目中6.1 扩展到带遗忘因子的时变系统跟踪如果系统参数随时间缓慢变化用固定λ1的RLS参数永远追不上真值。我建议先仿真一个参数突变场景验证遗忘因子的效果% 在第1500个点处让a1从-1.5突变到-1.2 for k 1500:N a1_true_now -1.2; % 实际代码里用一个随时间变化的变量 ... end对比 λ1 和 λ0.96 两种情况你会发现带遗忘因子的算法经过一段过渡期后能重新收敛到新参数λ1则牢牢卡在旧参数附近不动。这就是遗忘因子最直观的价值。λ的取值还可以自适应系统参数变化越快λ取越小噪声大则λ要适当提高。有些文献用变遗忘因子算法核心思想是每当检测到显著预测误差时自动降低λ。课程设计不用做这么复杂但理解这个逻辑对后面读文献很有帮助。6.2 扩展到有色噪声场景前面假设系统噪声是高斯白噪声这在实际中太理想了。如果噪声是相关的比如电机换向引起的周期性扰动标准RLS估计会存在偏差。这时候要把噪声模型也纳入辨识模型用增广递推最小二乘RELS。做法是把回归向量扩展为φ(k) [-y(k-1), ..., u(k-1), ..., e_hat(k-1), e_hat(k-2), ...]其中e_hat是预测误差的历史值。这样估计的参数向量里就多了噪声模型的参数虽然让问题维度上升但辨识精度提升非常明显。MATLAB实现改动很小主要就是扩展phi的构造和theta_hat的维度。6.3 从辨识到自适应控制的衔接参数辨识不是终点辨识出来的模型可以直接用于自校正控制或自适应控制。一个常见的做法是每个采样周期先用RLS更新参数再用当前参数计算控制器参数比如极点配置、最小方差控制实现“辨识-控制”一体化。我上学时做过一个电机转速控制实验用RLS在线辨识电机的传递函数每步更新PID参数成功实现了模型参数漂移情况下的稳定控制。相比固定PID自适应控制在参数变化20%以内时表现明显更好。这个方向很适合作为课程设计的延伸课题。6.4 可视化技巧实验报告里几个图特别加分M序列的时域波形频谱图让评审同学一眼看到激励信号特性参数估计轨迹和真值比较图能直观看到收敛速度预测误差曲线用半对数坐标展示收敛趋势最后跑一组“辨识模型预测 vs 实际输出”对比图验证模型有效性。四张图下来整个辨识实验的逻辑链就完整了。我习惯把这几张图画在一个figure的subplot里面报告排版时非常方便。7. 最后的一点实际操作心得这套流程我前前后后带过不少学生跑通最想强调的还是那句话辨识结果不好先别急着怪算法回头看看输入信号和采样时间选对没有。M序列本身已经是很优秀的激励信号了但用错了参数照样白搭。我自己的排查顺序永远是信号频谱够不够宽 → 数据长度够不够长 → 噪声影响大不大 → 遗忘因子合不合适 → 最后才动算法本身。另一个容易被忽略的点是数据类型。MATLAB里comm.PNSequence输出默认是logical类型直接参与乘法运算时要注意转换。我在初版代码里就因为这个吃过亏结果算出来的输出全是0排查了半天才发现是类型问题。代码里用double(m_seq)转一下就好。如果你只是想拿到一个能跑通的实验直接复制整合我上面的代码块改改系统参数就能出结果。如果你想把这里面的门道吃透建议把M序列换成伪随机二进制序列PRBS、把白噪声换成有色噪声、把批处理最小二乘和递推算法做一个对比实验——这些都是很好的深入方向。做工程和做研究很多时候差别不在于会不会套公式而是在于能不能折腾清楚边界条件在哪。把M序列和RLS这套组合折腾明白后面接触卡尔曼滤波、子空间辨识这些进阶方法时你会发现自己理解起来快很多。