简介:分数阶滑模控制(FOSMC)算法是分数阶微积分与滑模控制的结合,能有效应对非线性、时变及不确定性系统。压缩包提供了该算法在MATLAB/Simulink环境下的完整实现,面向控制理论与工程应用的研究者、研究生及高年级本科生。共24个文件,包含12个MATLAB脚本(.m)、10个Simulink仿真图像(.png)、1个.mdl和1个.slx格式的模型文件,整体大小仅315KB,结构精简便于直接加载和运行。资源已有811人学习,控制器设计涉及运动学模型、动力学模型、轨迹规划与控制误差计算等环节,从分数阶微分方程建模、滑动表面设计到切换函数优化均有对应的脚本和模型演示。读者可通过仿真图像观察位置误差、速度误差及控制输入的变化,并借助模型文件调整参数,深入理解FOSMC的鲁棒性与抖振抑制机制,适合作为科研预研或课程实战的参考资料。
1. 分数阶滑模控制到底解决什么问题:从整数阶滑模的抖振说起
滑模控制以对参数摄动和外部扰动的强鲁棒性著称,但把第一个整数阶滑模控制器放进 Simulink 跑闭环时,看到的往往是控制量在平衡点附近高频来回切——这就是抖振。导师扫了一眼波形就撂下一句话:“这送给做电机的人,人家直接退货。”分数阶滑模控制算法正是在这个痛点上流行起来的:把滑模面和趋近律里的整数阶微积分推广到 0~1 之间的分数阶,让切换动作从瞬时翻转变成带记忆的渐进过渡,在保住鲁棒性的前提下把抖振幅值降一个数量级以上。这篇笔记围绕 Matlab/Simulink 实现,从分数阶算子近似、控制器搭建、参数调试到典型翻车场景排查完整走一遍,适合做四旋翼、电机驱动、双向储能变换器等非线性对象控制的工程师和学生。
2. 分数阶微积分与滑模控制的结合点:为什么分数阶能把抖振压下来
2.1 三种分数阶定义与 Simulink 里的落地选择
分数阶微积分不是“一个”定义,工程上常见的就有 Riemann-Liouville(RL)、Caputo 和 Grunwald-Letnikov(GL)三种,它们的区别主要体现在初值处理和离散化方式上。
- RL 定义适合理论推导,但其分数阶导数初值没有直观物理含义,在控制问题里对“初始状态”的解释很别扭。
- Caputo 定义对整数阶初值友好,拉普拉斯变换后初始条件可以直接用 x(0)、ẋ(0),所以多数控制论文写稳定性证明时习惯用 Caputo。
- GL 定义本质是差分加权和的历史累积,天然面向数值计算,做离散仿真和嵌入式实现最方便。
在 Simulink 闭环仿真里,绝大多数人不会去解析算分数阶算子,而是用频域近似。最经典的是 Oustaloup 滤波器:把 s^λ 在一个频段 [w_b, w_h] 内拟合成整数阶有理传递函数,然后直接拖一个传递函数模块或 LTI System 块进模型。这样做的代价是近似只在一个频段里成立,滤波器阶次 N 越大近似越准,但相位滞后和计算量也跟着涨。我的落地原则是:理论推导用 Caputo,连续仿真用 Oustaloup 近似,离散化验证用 GL 递推。三者混着用没有关系,只要边界清楚即可。
2.2 分数阶滑模面的构造:从 c·e + ė 到带记忆的滑模面
整数阶滑模面通常写成 s = c·e + ė,控制目标是把误差轨迹逼到 s = 0 这条流形上。分数阶化之后,最常见的形式是分数阶积分滑模面:
s = ė + c·e + μ·D^{−λ}e
这里 λ ∈ (0,1),μ 是分数积分项的权重系数,D^{−λ}e 表示对误差做分数阶积分。对比整数阶写法,多出来的 D^{−λ}e 把误差的历史信息引入了滑模面,相当于给系统加了一个长时记忆。当状态接近滑模面时,切换项 sign(s) 的翻转频率会因为这个分数阶项的惯性作用而明显降低,表现出来就是控制量不再像锯齿一样高频跳变。
这里稍微泼一点冷水:不是把任何整数阶次换成小数就能消抖。如果 λ 取值太靠近 0.9 以上,滑模面响应变慢,跟踪误差上升;太靠近 0.1 以下,分数阶项接近普通积分,和整数阶 PI 型滑模区别不大。后面第四章会给一组我在仿真里反复验证过的范围。
2.3 趋近律的分数阶改造与稳定性判断
有了滑模面还要设计趋近律。整数阶常用等速趋近律 ṡ = −ε·sign(s) 或指数趋近律 ṡ = −ε·sign(s) − k·s。分数阶版本可以直接把趋近律左侧换成分数阶导数:
D^{1−λ}s = −ε·sign(s) − k·s
这个形式在理论上跟整数阶一致:用李雅普诺夫候选函数 V = ½s²,对其求分数阶导数后可以推出到达条件,收敛条件是切换增益满足 ε > |d(t)|,d(t) 为集总扰动上界。这个条件和整数阶滑模形式统一,但正因为左侧是分数阶算子,s 到达零点的过程是渐进收敛而不是瞬间穿越,控制量的高频抖动因此被削掉一大截。
我在仿真里验证一个分数阶控制器是否真的“更优”,从来不看论文结论,而是做一组同对象对比:整数阶和分数阶分别跑同样的正弦扰动,保存控制量波形和跟踪误差。分数阶的切换频率应明显稀疏,而稳态误差量级不变。如果误差变大,先怀疑参数而不是算法本身。
3. 用 Matlab/Simulink 搭建分数阶滑模控制:从被控对象到闭环仿真
3.1 被控对象建模:以带非线性和扰动的一阶倒立摆简化模型为例
分数阶滑模控制适用于一类二阶非线性系统,写成一般形式:
ẋ₁ = x₂
ẋ₂ = f(x₁,x₂) + b·u + d(t)
这里 f(x₁,x₂) 代表模型已知的动力学部分,b 是控制增益,d(t) 是外部扰动。为了后面仿真好复现,用一个简化的一阶倒立摆摆角模型:
f(x₁,x₂) = −a·sin(x₁) − c_p·x₂
b = 1
d(t) = 0.5·sin(2π·t)
对应参数:a = 9.8,c_p = 0.1。目标轨迹取 xd(t) = 0.5·sin(t),让控制器在正弦跟踪工况下暴露动态性能差异。这个模型可以替换成四旋翼俯仰通道、电机位置环或者双向储能变换器的电压环,分数阶滑模控制器的结构不变,变的只是 f(x)、b 和扰动上界。
在 Simulink 里给这个对象建模,我习惯用积分器模块而不是 MATLAB Function 现算,因为后面要接 Oustaloup 滤波器和 S 函数,纯积分链路最容易排查代数环。模型里放两个 Integrator,第一个积分结果就是 x₁,第二个积分结果作为 x₂ 反馈到 f 的计算端;外部扰动用一个 Sine Wave 模块叠加在 x₂ 的导数上。所有状态量通过 Goto/From 或者 Bus 引出,避免线缆交叉。
3.2 生成分数阶算子近似:Oustaloup 滤波器的 MATLAB 实现
Oustaloup 滤波器的核心函数不长,我自己维护了一个版本,输入频段边界、滤波器递归次数和分数阶次,输出一个 zpk 对象,可以直接喂给 Simulink 的 LTI System 模块。
function G = oustaloup_sf(wb, wh, N, lambda) % oustaloup_sf 生成分数阶算子 s^lambda 的有理近似 % 输入: % wb - 拟合频段下限 (rad/s),通常取 0.001~0.01 % wh - 拟合频段上限 (rad/s),通常取 1000~10000 % N - 递归次数,滤波器实际阶数为 2N+1 % lambda - 分数阶次,0<lambda<1 表示微分,-1<lambda<0 表示积分 % 输出: % G - zpk 对象,可直接传给 tf() 或 Simulink LTI System 块 if lambda >= 1 || lambda <= -1 error('lambda 必须落在 (-1,0) 或 (0,1) 区间'); end k = 1:N; % 递归求解零极点频率 w_k = wb * (wh/wb).^((k + N + 0.5*(1-lambda))/(2*N+1)); w_kp = wb * (wh/wb).^((k + N + 0.5*(1+lambda))/(2*N+1)); % 增益对齐,保证滤波器中频增益接近 s^lambda 的理论值 K = (wh/wb)^(-lambda/2) * prod(w_kp) / prod(w_k); % 零极点形式,所有零极点都在左半平面 G = zpk(-w_kp, -w_k, K); end这段代码的关键参数有三个:lambda 决定近似的是微分还是积分,N 决定滤波器阶数,wb/wh 决定拟合频段。拟合频段必须覆盖系统的闭环带宽。比如系统机械带宽 10 rad/s,那就至少取 wb=0.01、wh=1000。频段太小,分数阶特性会在实际工作频率上失真;频段太大,滤波器阶数不变的情况下两端的匹配误差会增大。生成完用 bode(G) 扫一眼,看中频段的幅频斜率是否接近 20·λ dB/dec,这是一个非常快的体检方法。
3.3 控制器实现:用 Simulink 模块搭分数阶滑模控制器
控制器部分我不写复杂的 S 函数,直接在 Simulink 里用 Sum、Gain 和 MATLAB Function 模块搭出来,理由是好调参、好加 Scope 观察中间变量。控制律基于第二章的分数阶积分滑模面:
s = ė + c·e + μ·D^{−λ}e
u = (1/b)·(−f + ẍd − c·ė − μ·D^{1−λ}e − ε·sign(s) − k·s)
D^{−λ}e 和 D^{1−λ}e 各用一个 Oustaloup 滤波器生成。前者 lambda 取 −0.6,后者 lambda 取 0.4,注意 1−λ=0.4 正好是互补阶次。控制器的核心计算放到一个 MATLAB Function 块里:
function u = fsmc_controller(e, edot, frac_der, frac_int, f_est, b_est) % 分数阶滑模控制器核心计算 % 输入: % e - 跟踪误差 xd - x1 % edot - 误差导数 xd_dot - x2 % frac_der - 外部滤波器输出的 D^(1-lambda) e % frac_int - 外部滤波器输出的 D^(-lambda) e % f_est - 被控对象模型已知部分 f(x1,x2) % b_est - 控制增益估计值 % 输出: % u - 控制量 % 控制器参数,实际调试时改成全局变量或从Mask读入 c = 8; % 滑模面比例增益 mu = 1.5; % 分数积分项增益 eps = 0.5; % 切换增益,必须大于扰动上界 k = 20; % 趋近律线性增益 % 分数阶滑模面 s = edot + c*e + mu*frac_int; % 控制律:等效控制 + 切换控制 u = (1/b_est) * (-f_est + 0.5*cos(0.5) - c*edot - mu*frac_der ... - eps*sign(s) - k*s); end这里有个容易踩坑的地方:代码里的目标加速度项 xd_ddot 如果是随时间变化的解析函数,不要写在控制器函数里固定值,应该从外部 Signal 引入。上面示例里目标轨迹 xd=0.5·sin(t),加速度是 −0.5·sin(t),我直接把这个信号作为控制器的一个输入端接进 MATLAB Function,而不是把解析表达式写死在函数体里。好处是后续换目标轨迹不用改控制器代码,只改信号源。
控制器参数 c、mu、eps、k 我建议全放在 MATLAB 工作区里,模型里用变量名引用。这样跑参数扫描时只需要写一个 for 循环改工作区变量再 sim(),不需要打开模型手工改 Gain 模块。
3.4 完整闭环模型搭建与求解器设置
把对象模型、两个 Oustaloup 滤波器、MATLAB Function 控制器按下列顺序连接:
设定点信号 xd 经过一个减法器得到 e,xd 的导数用 Derivative 模块得到 xd_dot,再与 x2 相减得 edot。e 同时送入两个 Oustaloup 滤波器,分别输出 frac_der 和 frac_int。控制器输出的 u 直接接到对象模型的 u 输入端,对象模型反馈 x1、x2 回减法器。为了方便观察,用 Scope 接四处信号:x1 和 xd 的对比、s 滑模面、u 控制量、e 跟踪误差。
求解器配置是这套仿真最容易翻车的环节。我建议把 Solver 设为 ode45,最大步长不要超过 0.001s。Oustaloup 滤波器是一个 2N+1 阶的高阶传递函数,如果你用默认的变步长再加一个宽松的最大步长,Simulink 会为了维持精度把步长压到 1e-6 甚至更小,仿真速度慢到让人怀疑电脑坏了。反之,如果你为了追求速度把最大步长设成 0.01s,高频段的滤波器动态会被步长截断,仿真结果看起来像另一个系统。先固定最大步长 0.001s 跑通,再根据控制量波形决定是否放宽。
4. 分数阶滑模控制器必调的三个参数:λ、趋近律系数与滤波器阶次
4.1 分数阶次 λ 的取值范围与效果对照
λ 是整个控制器里最“玄学”的参数。它的取值直接影响滑模面的记忆特性和收敛速度,以下是几组仿真中反复验证的典型区间:
| λ 取值范围 | 系统表现 | 适合场景 |
|---|---|---|
| 0.2~0.4 | 滑模面接近整数阶,记忆效应弱,响应快但抖振压制有限 | 对响应速度要求高的场合 |
| 0.4~0.7 | 记忆效应适中,抖振和跟踪精度平衡最好 | 大多数二阶非线性系统首选 |
| 0.7~0.9 | 强记忆效应,抖振很小但响应变慢,误差收敛拖尾 | 扰动幅度大、控制量受限的场合 |
我一般把 λ 的初值设在 0.5,跑通后以 0.05 为步长向两侧扫描。判断标准不是看误差,而是看控制量 u 的开关次数。取一段稳定时间窗,统计 sign(s) 翻转次数,翻转次数下降 50% 以上且跟踪误差 RMS 上升不超过 20%,这个 λ 就是可用的。
4.2 趋近律参数 c、eps、k 的调参顺序
这三个参数和 λ 不是平级关系。我的调参顺序是先 c 再 k 最后 eps,因为 c 决定滑模面的“斜率”,k 决定趋近过程的线性收敛速度,eps 决定稳态附近抵抗扰动的最小切换强度。
c 从 1 开始,逐渐加大直到系统阶跃响应不出现超调震荡。c 太大,滑模面过陡,等效控制量放大,容易激发未建模动态;c 太小,误差收敛慢,滑模面趋于纯分数项。
k 的调节看趋近过程:给系统一个阶跃目标,观察 s 从初始值降到零区间的耗时。k 加大这个耗时缩短,但 k 过大一样会引起控制量尖峰。一般取 λ 确定后系统最快收敛值的 60%~80%。
eps 是最后一个碰的参数。它的下限由扰动上界决定,理论条件是 eps > |d(t)|_max。多数人犯错是把 eps 调得过大,觉得“切换大力出奇迹”,结果抖振又回来了。我的做法是先把 eps 设成扰动幅值的 1.5 倍,然后逐步减小,直到控制量波形出现可察的高频毛刺就回退 20%。
4.3 Oustaloup 滤波器阶次 N 与拟合频段的折衷
滤波器参数对控制器性能的影响常被忽略,但它直接决定了分数阶算子的“纯度”。N 取 3~5 在多数连续仿真里就足够,N 取到 8 以上除了拖慢仿真速度,相位滞后还会在高频段扭曲滑模面的真实形态。
拟合频段的上下界比 N 重要得多。一个实用的经验公式:wb 取系统闭环带宽的 1/100,wh 取闭环带宽的 100 倍。如果系统闭环带宽约 10 rad/s,那 wb=0.1、wh=1000。这样保证分数阶特性覆盖主要的能量频段,同时避免极端频段带来的数值病态。
可以用一个简单方法检查滤波器是否合格:把 Oustaloup 滤波器和一个理想分数阶积分器分别在相同正弦输入下跑,比较输出的幅值和相位。幅值误差超过 5% 或相位偏差超过 3°,就调整频段或增加 N。
5. 分数阶滑模仿真避坑指南:5 个典型翻车现场与排查方法
5.1 现象:仿真步长越跑越小,十分钟才走 0.1 秒
原因几乎都出在 Oustaloup 滤波器和变步长求解器的组合上。滤波器阶次 N 高、转折频率密集,ode45 的误差控制器会认为系统“刚性”程度高而不断压缩步长。这不是模型发散,是数值刚度导致的效率灾难。
解决:先检查 wb 是否太小,如果 wb 比系统工作频率低了 4 个数量级以上,滤波器在低频段引入了极慢极点,直接影响步长。把 wb 上调到闭环带宽的 1/50~1/100。另外把最大步长锁定到 0.001s,禁止求解器无限压小步长。如果还慢,把滤波器 N 降到 3,绝大多数仿真模型的等效精度差异在可接受范围。
5.2 现象:控制量波形在平衡点附近高频振荡,切换频率比整数阶滑模还高
这种情况一般是 eps 设得过大,切换项 sign(s) 主导了控制量。分数阶滑模的优势在于降低开关频率,但前提是切换增益刚好压过扰动。eps 过大时,分数阶项的缓冲作用被过强切换淹没,抖振不减反增。
解决:把 eps 从当前值逐步减半观察控制量波形。如果减到扰动幅值附近波形还没有明显高频毛刺,说明还有余量。另外检查 mu 的取值:mu 是滑模面分数积分项的权重,如果 mu 太小,分数阶记忆效应在滑模面里占比不足,D^{−λ}e 项几乎是摆设。
5.3 现象:闭环仿真直接发散,误差跑到 1e10 量级
先按顺序查三件事:滑模面 s 是否出现计算溢出的中间量,Oustaloup 滤波器是否因传递函数实现形式产生了非最小相位零点,被控对象参数 b_est 是否错误地取了负号。最常见的翻车点是把 b_est 设成了真实 b 的相反数,等效控制项的正反馈直接吹飞系统。
另一个隐蔽原因是目标轨迹的加速度信号没有接入控制器。某些实现里我把 xd_ddot 写死在函数体里,一旦目标轨迹换成阶跃或者分段信号,固定值的加速度项会失控。排查方法是在控制器输出端加一个 Data Type Conversion 和饱和限幅,先把异常控制量截住,再逐项检查输入信号。
5.4 现象:仿真结果看起来不错,但换成另一台电脑或 MATLAB 版本后波形变了
这个属于数值实现细节。不同 MATLAB 版本对 zpk 对象的内部实现方式有差异,Oustaloup 滤波器从零极点形式转换成传递函数时,极点和零点的排序变化会影响数值精度。2023b 和 2021a 跑同一模型,高频段波形可能有细微差别。
解决:把滤波器对象在仿真前用 tf() 显式展开,并用 balred 做一次降阶,显式指定状态空间实现。这样滤波器在模型里是一个确定性的状态空间矩阵,而不是依赖版本内部算法的 zpk 对象。再用 ss() 转成状态空间后保存为 .mat,模型里直接从工作区加载。
5.5 现象:分数阶滑模面的初值对结果影响巨大,换一个初值系统就换一个性能
分数阶算子自带历史记忆,仿真起点附近的一段“预热期”会持续影响后续滑模面的形态。如果系统初始误差大,D^{−λ}e 在起点附近累积的负值会造成滑模面初始位置偏差,控制器要花更长时间归位。
解决:给积分器设置合理的初始状态。Simulink 里的 Integrator 模块可以直接填 Initial Condition,把 x1、x2 的初值设成目标轨迹的初值,而不是默认的 0。这样误差起始为零,分数阶积分项的初始记忆也是零,滑模面从原点出发。实机上则要先跑一段“零点校准”,确认传感器反馈等于零再切闭环。
6. 从 Simulink 仿真走向实机:分数阶算子的离散化与鲁棒性验证
实机部署前,Oustaloup 连续滤波器不能直接生成 C 代码,因为它的零极点数量在嵌入式平台上会造成过重负担。我的做法是把控制律改写成离散 GL 实现:在 MATLAB Function 里维护一个误差历史缓冲区,每次采样周期按 GL 系数加权求和近似分数阶算子。系数的计算放到一个初始化脚本里预先算好,运行时只做乘加运算,MCU 上的开销基本可忽略。
鲁棒性验证我习惯用蒙特卡洛批量测试,而不是只看一两条波形。脚本循环 50 次,每次给扰动幅值和相位随机扰动,记录跟踪误差 RMS 和控制量最大幅值,最后用箱形图检查离群点。如果某些参数组合下误差明显变大,先用 5.3 节的方法排查,再考虑调整 λ。正常结果应该是误差 RMS 分散度在 20% 以内且没有离群值。
另外不要相信连续仿真的时间响应,实机第一个版本务必用固定步长离散仿真验证一次:求解器设为 ode4、步长 1ms,与后续 C 代码生成的采样周期保持一致。这一步能提前暴露离散化带来的相位滞后,避免在实机上被“万物皆低通”的教育。我踩过最大的坑就是连续仿真调好的 eps 在固定步长下完全不够,扰动照样穿过滑模面。记住一句教训:分数阶滑模控制器的参数不是调一次就完事,更换仿真步长、滤波器频段或目标轨迹后,都要回归测试一遍 λ 和趋近律参数的组合。
希望这套从原理到排错再到实机验证的路径能帮你在分数阶滑模控制上少走几个月弯路。
本文还有配套的精品资源,点击获取