简介一份面向Matlab/Simulink用户的分数阶滑模控制FOSMC算法资源包聚焦非线性、时变及不确定性系统的鲁棒控制设计适合控制理论与工程应用方向的研究生、工程师及进阶本科生。压缩包共24个文件包括12个m脚本、10个png曲线图、1个mdl模型与1个slx模型其中m脚本覆盖轨迹生成、运动学/动力学建模、控制器与误差计算等模块png图为各状态响应曲线可直接用于结果对比与论文写作mdl/slx模型支持在Simulink中运行和修改。已有811人学习该资源关注度较高。通过研读代码与模型可掌握分数阶微积分建模、滑动表面设计、切换函数平滑处理及参数整定方法还能根据具体对象快速替换模型扩展为自定义控制方案。资源包仅315KB内容精炼、结构清晰适合作为课程设计或课题初期的参考模板。1. 分数阶滑模控制算法在 Simulink 里落地别被 D^μ 吓住核心就三个问题分数阶滑模控制算法FOSMC听起来唬人但真正在 Simulink 里复现过一轮的人都知道控制律推导占三成分数阶算子怎么离散化占七成。最近我把一套基于 MATLAB Simulink 的 FOSMC 模型完整拆开重搭发现最值得注意的就三件事——用 Grünwald-LetnikovGL公式把分数阶积分写成可编程的累加器、把滑模面从 sėλe 改成 sėλD^{-μ}e、再用边界层饱和函数把抖振压到可接受范围。这套方案适合正在做电机、机械臂、四旋翼等二阶系统鲁棒控制的人不需要额外工具箱固定步长下跑得很稳。2. 控制律结构拆解GL 离散化、Oustaloup 逼近与参数初值怎么定2.1 三种分数阶算子实现路径FOTF、GL、Oustaloup 怎么选在 MATLAB 里算分数阶导数/积分常见的有三条路。第一条是直接用现成的 FOTF 工具箱Fractional Order Transfer Function国内最常用的版本是薛定宇老师维护的。它有 fotf 对象可以定义 s^α 形式的分数阶传递函数也可以直接对信号做分数阶求导和积分。好处是验证理论特别方便几行代码就能画出 D^0.5 sin(t) 的曲线。但缺点也很明显它是第三方工具箱不同 MATLAB 版本下兼容性不稳定想把它嵌进 Simulink 连续模型里还要自己转成离散滤波器或状态空间绕一圈回来不如直接写 GL。第二条是自己写 GLGrünwald-Letnikov公式的离散化。GL 定义是D^α f(t) ≈ h^{-α} Σ_{j0}^{N} (-1)^j C(α,j) f(t-jh)其中 h 是离散步长N 是记忆长度C(α,j) 是推广到实数阶的二项式系数。注意这个公式对 α0微分和 α0积分都成立一个函数就能把分数阶积分器和分数阶微分器都做了。它的好处是零依赖、逻辑透明、代码量小而且可以把记忆长度 N 当作一个显式参数来控制计算精度。代价是必须用固定步长仿真步长一变结果就乱。第三条是 Oustaloup 滤波器近似。它把 s^α 在一个频带 [ω_l, ω_h] 内近似成一个整数阶传递函数。这个近似的阶数是 2N1 阶比如 N4 就是 9 阶传递函数。好处是它可以像普通传递函数一样放进 Transfer Fcn 模块甚至在变步长连续仿真里也能跑。坏处是频带外会失真而且阶数高的时候数值敏感容易跑出一堆轻微的数值振荡。我一般只在需要连续域分析时才用它做步进仿真首选 GL。三者的取舍用一张表说清楚实现路径依赖固定步长要求适用场景主要风险FOTF 工具箱第三方工具箱不强制理论验证、快速画波形版本兼容差Simulink 集成绕GL 离散化无必须仿真、代码生成、嵌入式移植步长一变就发散Oustaloup 滤波器无不强制连续域设计、变步长仿真频带外失真、阶数高易震荡我给的建议是做学位论文的仿真验证用 GL要在连续系统框架里做频域分析用 Oustaloup只想在命令行里验证一下分数阶算子性质用 FOTF。下面整套模型都按 GL 来写。2.2 滑模面、趋近律和控制律推导为什么分数阶能压抖振先交代被控对象。假设一个带摩擦和外部扰动的二阶系统J θ̈ b θ̇ d(t) u其中 u 是控制输入d(t) 是外部扰动J 是转动惯量b 是粘滞摩擦系数。这个模型在电机、单关节机械臂、四旋翼姿态里到处可见所以拿它讲 FOSMC 最直观。设跟踪误差 e θ_d - θ。经典整数阶滑模面是 s ė λe。分数阶的做法是把误差的整数阶积分换成 μ 阶分数阶积分s ė λ D^{-μ} e, μ ∈ (0,1)这里 D^{-μ} 是 μ 阶分数阶积分算子。当 μ0 时 D^0 e e滑模面退化成 s ė λe所以这套设计天然兼容整数阶滑模做对比实验时只要把 μ 设成 0 就行不需要另搭一套模型。这一点后面会反复用到是个很省事的性质。趋近律取常见的指数趋近律加边界层ṡ -ε·sat(s/φ) - k·ssat 是饱和函数φ 是边界层厚度。对 s 求导代入系统方程解出等效控制整理后得到控制律u J θ̈_d b θ̇ J λ D^{1-μ} e J (ε sat(s/φ) k s)前三项是等效控制Jθ̈_d 补偿期望加速度bθ̇ 补偿摩擦JλD^{1-μ}e 把误差的分数阶动力学引入控制最后一项是切换控制负责抵抗扰动和模型误差。为什么分数阶版本通常比整数阶版抖振小一个常见的解释是μ 阶积分把误差的历史记忆揉进滑模面等效控制里多了一项 D^{1-μ}e。当 μ 在 0 到 1 之间时1-μ 也在 0 到 1 之间这是一个比一阶微分更“温和”的算子它不会像整数阶微分那样放大高频噪声所以切换项需要付出的补偿量变小抖振的底噪也就下来了。抖振这东西有点玄学理论上说得再好最终还是要看仿真和实验里的数据后面第 4 章会专门给量化指标。2.3 控制参数初值一张表说清 μ、λ、ε、k、φ 的试凑起点新手最容易踩的坑是拿到一篇文章的参数直接抄抄完发现模型发散就怀疑算法有问题。实际上 FOSMC 的一堆参数互相耦合得知道每个参数管的到底是哪一段。参数含义建议初值范围我常用的起点调参方向μ滑模面分数阶积分阶次0.2 ~ 0.80.5稳态误差压不住就增大 μλ滑模面比例系数5 ~ 5020收敛慢就增大过冲大就减小ε切换增益抵抗扰动1 ~ 10扰动上界的 1.2 倍抖振大就减小跟踪偏差大就增大k指数趋近律线性增益10 ~ 5015收敛慢就增大但别超过 ε 太多φ边界层厚度0.01 ~ 0.10.05抖振明显就增大稳态精度要求高就减小调参有个经验顺序先把 μ 定在 0.5λ 和 k 从中间值开始调出基本稳定的收敛曲线然后逐步加大扰动观察跟踪误差按需增大 ε最后再调 φ 来平衡抖振和稳态精度。不要一上来就同时动五个参数那样你根本不知道是哪个参数把系统搞坏的。注意这里的 λ 和 k 在控制律里的位置不同。λ 乘在分数阶积分项上影响等效控制的补偿k 乘在趋近律的线性项上直接决定 s 回到滑模面的速度。两者协同但各有侧重第 4 章对比实验里你会看到它们对结果的贡献不一样。3. Simulink 模型搭建MATLAB Function 块、积分器被控对象与数据采集3.1 模型级联结构与模块清单整个模型自上而下分四段期望轨迹 → 误差计算 → FOSMC 控制器 → 被控对象。扰动在对象内部注入不经过控制器这样才符合真实情况——控制器并不知道扰动的具体波形。用到的模块清单如下模块类型参数或说明theta_dSine Wave幅值 1频率 0.5 rad/s期望轨迹 θ_d sin(0.5t)eSumtheta_d - thetaControllerMATLAB FunctionFOSMC 控制律GL 算子输入 [t;e;thetadot]dSine Wave幅值 0.5频率 3 rad/s外部扰动 d(t)Plant两个积分器 Gain实现 Jθ̈ bθ̇ u - deout/uout/toutTo Workspace采集误差、控制量、时间关键设置两处一是求解器必须选固定步长我用 ode4步长 h0.001二是在 MATLAB Function 块的 Edit Data端口和数据管理器里把 mu、lambda、epsilon、k、phi 声明为 Parameter并绑定工作区变量名而不是写死在代码里。先把这一步做对后面能省大量返工。3.2 MATLAB Function 块里的 FOSMC 实现控制器是这个模型的心脏。下面的代码把 GL 分数阶算子和控制律实现在一个 MATLAB Function 块里function u fosmc(in, mu, lambda, epsilon, k, phi) % FOSMC 控制器GL 离散化分数阶算子 % in [t; e; thetadot]三个输入合成一个向量 % mu, lambda 等通过 MATLAB Function 块 Parameter 绑定工作区变量 t in(1); e in(2); thetadot in(3); persistent buf_mu buf_1mu coef_mu coef_1mu e_prev init_done N 100; % GL 记忆长度 h 0.001; % 与模型定步长保持一致 if isempty(init_done) % 首次调用时预计算二项式系数 init_done 1; e_prev e; buf_mu zeros(N1, 1); buf_1mu zeros(N1, 1); alpha_mu -mu; % 分数阶积分 alpha_1mu 1 - mu; % 分数阶微分 coef_mu zeros(N1, 1); coef_1mu zeros(N1, 1); coef_mu(1) 1; coef_1mu(1) 1; for j 1:N coef_mu(j1) coef_mu(j) * (alpha_mu - j 1) / j; coef_1mu(j1) coef_1mu(j) * (alpha_1mu - j 1) / j; end coef_mu coef_mu .* (-1).^(0:N); % 并入 (-1)^j coef_1mu coef_1mu .* (-1).^(0:N); end edot (e - e_prev) / h; % 误差导数一阶差分 e_prev e; buf_mu(2:end) buf_mu(1:end-1); % 历史窗口右移 buf_mu(1) e; buf_1mu(2:end) buf_1mu(1:end-1); buf_1mu(1) e; Dmume h^(-alpha_mu) * (coef_mu * buf_mu); % D^{-mu} e D1mu h^(-alpha_1mu) * (coef_1mu * buf_1mu); % D^{1-mu} e s edot lambda * Dmume; % 分数阶滑模面 sat max(min(s / phi, 1), -1); % 饱和函数替代 sign thetadd_d -0.25 * sin(0.5 * t); % 期望加速度解析式 u 0.5 * (thetadd_d lambda * D1mu) 0.1 * thetadot ... 0.5 * (epsilon * sat k * s); end逻辑说明整个函数的核心是把 GL 公式的有限记忆窗口实现成了循环移位数组。每一步采样把所有历史值右移一位当前误差放到窗口头部然后和二项式系数做点积乘上 h^{-α}就得到当前时刻的分数阶算子输出。D^{-μ}e 用于构造滑模面 sD^{1-μ}e 用于等效控制两者共用同一个误差信号只是 α 不同所以维护两个缓冲区。参数说明N 是记忆长度N 越大历史信息保留越多但计算量线性增加。步长 0.001、N100 时单步只做 101 次乘加开销可忽略对精度要求高可以取 200。h 必须和模型定步长一致不然 GL 公式的离散基准就错了这是最容易翻车的地方。控制律里的 0.5 和 0.1 分别是 J 和 b 的数值如果你换被控对象这两个数要跟着改。建议把 J 和 b 也绑定成工作区参数别学我图省事直接写死在代码里。3.3 被控对象与扰动注入用积分器搭二阶系统被控对象 Jθ̈ bθ̇ u - d 在 Simulink 里用两个积分器搭最直观输入 u-d → Gain 1/J → 第一个积分器 → 得到 θ̇ → Gain b 反馈到输入形成摩擦项 → 第二个积分器 → 得到 θ。这样搭的好处是每一路信号都能直接用 Scope 看排查问题时比黑匣子传递函数直观得多。扰动 d 用一个 Sine Wave 模块幅值 0.5、频率 3 rad/s加在对象输入端。注意扰动不要经过控制器否则控制器相当于预知了扰动路径鲁棒性验证就没意义了。采集部分用三个 To Workspace 模块分别输出 e、u、θ。变量名设成 eout、uout、thetaout格式选 Array。仿真时长 20 s。按这个结构搭完模型是期望轨迹模块 → 误差和模块 → Controller → 被控对象 → 反馈回误差同时 theta 和 thetadot 作为额外输入进 Controller。thetadot 直接从第一个积分器输出引线不需要再对 θ 求导这样能少一个数值噪声源。4. 脚本联调与对比实验参数扫描、μ0 退化校验和抖振量化4.1 主脚本工作区变量驱动模型重跑模型搭好之后真正干活靠的是主脚本。Simulink 模型里的 mu、lambda 这些参数必须绑定工作区变量在 MATLAB Function 块的 Edit Data 里把参数勾为 Parameter并绑定变量名脚本里给变量赋值再调用 sim()模型才会按新参数跑。很多人栽在这一步直接在脚本里写mu 0.5;但模型块里用的是常数模块脚本改多少都白搭。脚本框架如下%% 初始化 clear; clc; h 0.001; T 20; J 0.5; b 0.1; %% FOSMC 参数 mu 0.5; lambda 20; epsilon 2; k 15; phi 0.05; %% 跑分数阶滑模 sim(fosmc_demo.slx); e_f eout; u_f uout; t tout; %% 看基本曲线 figure(1); subplot(2,1,1); plot(t, e_f); title(FOSMC error); grid on; subplot(2,1,2); plot(t, u_f); title(FOSMC u); grid on;逻辑说明sim() 执行前工作区必须存在模型引用的所有变量。给 mu 赋值后模型里的 GL 系数初始化逻辑会按新 μ 重算。这是 Persistent 初始化机制带来的好处——每次仿真开始时 init_done 为空函数首次调用会重新预计算系数所以更换参数后不需要重启 MATLAB。参数说明h 是仿真步长也是 GL 公式里的离散步长两处必须一致。T 是仿真时长这里设 20 s 足够让 0.5 rad/s 的期望轨迹跑完大约三个周期能同时观察暂态和稳态。如果你想把扫描结果存档建议在脚本里加循环把每组参数的 ITAE 和抖振幅值写进结构体后面画对比图就方便了。4.2 整数阶对照把 μ 置零退化为标准滑模这套模型最有价值的一点是 μ0 时自动退化成整数阶滑模ISMC。因为 D^0 e e滑模面变成 s ė λe控制律跟标准 ISMC 完全一致。做对比实验不需要另建模型只需要把 μ 改成 0 再跑一遍%% 整数阶对照mu0 退化 mu 0; sim(fosmc_demo.slx); e_i eout; u_i uout;对照时注意两点一是保持 λ、ε、k、φ 完全不变只改 μ这样对比的才是分数阶和整数阶的差异而不是参数调优的差异二是 μ 从 0.5 改到 0 之后D^{1-μ}e 变成 D^1 e即一阶微分GL 公式里 alpha_1mu 1二项式系数退化为普通差分数值上同样可靠不会出现奇异情况。如果你发现 μ0 时结果异常优先检查是不是两套参数没同步。4.3 怎么读结果ITAE、抖振幅值与收敛时间对比只靠肉眼盯曲线不够我常用三个量化指标%% 指标计算 ITAE_f sum(abs(e_f) .* t) * h; % 积分时间绝对误差 ITAE_i sum(abs(e_i) .* t) * h; jitter_f std(diff(u_f)); % 控制量高频分量标准差 jitter_i std(diff(u_i)); t_settle_f t(find(abs(e_f) 0.01, 1)); % 首次进入稳态带 t_settle_i t(find(abs(e_i) 0.01, 1)); fprintf(FOSMC: ITAE%.3f, jitter%.3f, settle%.2fs\n, ITAE_f, jitter_f, t_settle_f); fprintf(ISMC : ITAE%.3f, jitter%.3f, settle%.2fs\n, ITAE_i, jitter_i, t_settle_i);逻辑说明ITAE 给误差乘上时间权重后期误差权重更大能反映长时间稳态跟踪品质jitter 用相邻两步控制量差值的标准差表示抖振越厉害这个值越大收敛时间取误差首次进入 ±0.01 rad 带的时刻。这三个指标分别对应跟踪精度、抖振强度、快速性刚好覆盖滑模控制最关心的三个维度。参数说明ITAE 的求和结果依赖步长 h脚本里必须乘 h不然换步长后指标不可比。jitter 对噪声敏感如果被控对象输出有测量噪声先加一阶低通再算 diff。还有一个细节t(find(...)) 在误差从来没进过 ±0.01 带时会返回空脚本会报错稳妥做法是先判断 isempty 再赋值。从我这组初值跑出来的典型结果是FOSMC 的 ITAE 比 ISMC 低 10%~20%jitter 低 30% 以上收敛时间两者接近。但注意这是固定参数下的对比不同的 λ、ε 组合下差距方向可能变化。做论文的话建议做参数扫描把 μ 从 0.1 扫到 0.9画出指标随 μ 变化的曲线比单独跑两三条曲线有说服力得多。5. 避坑与常见问题步长不一致、代数环、参数埋死这五关分数阶滑模在 Simulink 里翻车的地方高度集中下面五条是我实际踩过的按现象、原因、解决写清楚每条都是血泪经验。5.1 变步长跑 GL 算子结果复现不了还满天 NaN现象模型用 ode45 变步长跑同一组参数第一次跑和第二次跑结果不完全一致改一个小参数就 NaN仿真直接崩。原因GL 公式里 h 是固定离散步长而 ode45 会自适应改步长。MATLAB Function 块里的 h 写死 0.001实际求解步长却可能跳到 0.01 甚至更大GL 的高阶项误差被指数放大很快溢出成 NaN。解决把模型求解器改成固定步长。Solver settings 里 Type 选 Fixed-stepSolver 选 ode4Fixed-step size 填 0.001保证和 GL 的 h 严格一致。从那以后我只要看到谁用 GL 算子配变步长就知道离 NaN 不远了。5.2 代数环控制器和被控对象互相等待现象仿真一开始就弹 Algebraic loop detected有时候甚至直接卡死不更新。原因控制器输出 u 进被控对象被控对象输出 θ 反馈回控制器算误差误差又进控制器算 u形成闭合的代数依赖。GL 版本其实用的是历史数据理论上不存在本步依赖但 MATLAB Function 块输入输出直接相连时求解器仍会按连续路径检测代数环。解决在误差 e 进入 MATLAB Function 块之前加一个 Memory 模块把路径断成半拍延迟0.001 s 的延迟对跟踪影响可以忽略。如果你用的是 Oustaloup 滤波器方案代数环是真实存在的因为滤波器的状态更新需要当前步输入这时候要么把滤波器拆成显式状态空间要么干脆换 GL。5.3 参数埋死在块里脚本里改一万遍都没用现象脚本里写了mu 0.8; sim(fosmc_demo.slx);结果输出曲线和 mu0.5 时一模一样一点变化没有。原因模型里把 mu 写成了常数模块或直接写进 MATLAB Function 块代码里跟工作区变量没有绑定关系。脚本改的变量根本传不到模型里。解决在 MATLAB Function 块的 Edit Data 里把 mu、lambda、epsilon、k、phi 声明为 Parameter绑定工作区变量名。这样每次 sim() 前只要给工作区变量赋值模型就自动更新。这是个习惯问题我在第 3 章就强调过但几乎每个初学者都会在这里卡一次。5.4 用整数阶积分器串联假装分数阶频率特性压根不对现象有人用一串 Integrator 模块串联比如四个积分器想近似 D^{-1/2}结果系统发散或者稳态行为完全不符合理论。原因n 个整数阶积分器串联是 (1/s)^n幅频特性是 -n×20dB/dec 的直线斜率而分数阶积分 D^{-μ} 的斜率是 -μ×20dB/dec。μ 在 0 到 1 之间时串联整数阶积分器要么斜率太陡要么响应完全不是一回事没有任何近似意义。解决要么用 GL 离散化要么用 Oustaloup 滤波器搭一个真正的分数阶传递函数。整数阶模块拼不出分数阶算子这不是精度问题是频率特性本质不同。5.5 ε 和 φ 的配合抖振不减反增的疑难杂症现象扰动变大后把 ε 从 2 加到 5抖振反而更凶甚至出现极限环振荡把 φ 调大稳态误差又涨回来。原因ε 是切换增益它必须刚好盖过扰动上界多出来的部分全是抖振能量源。φ 是边界层厚度φ 太小s 在零附近高速穿越时饱和函数近似 sign等效于没消抖φ 太大等效控制被边界层吸收稳态精度下降。解决先给扰动一个上界估计ε 取值不要超过上界 1.2 倍然后从 φ0.1 往小调每调一次看一次 jitter 指标找到抖振和稳态误差的平衡点。我一般把 φ 和 ε 绑在一起调ε 固定后φ 从大到小递进直到 jitter 不再明显下降为止再回头微调 ε。这套流程很笨但有效比靠感觉试快得多。6. 进阶把 FOSMC 封装成自定义库再走一遍代码生成验证6.1 子系统封装与掩码参数模型验证没问题后我会把控制器整体选中右键 Create Subsystem 封装再用 Mask Editor 把 mu、lambda、epsilon、k、phi 暴露成掩码参数。这样双击模块就是一个参数面板不用打开 MATLAB Function 块去改代码。封装好的模块拖进自定义库保存后续新模型直接引用比复制粘贴稳得多。掩码参数校验在 Mask 的 Initialization 回调里做mu 限制在 0~1lambda 必须大于 0phi 不能为 0参数非法直接 error() 中断比仿真炸了再去猜原因强。6.2 代码生成的验证清单要往硬件在环走仿真通过还差一步。代码生成前的检查项检查项要求原因求解器固定步长连续求解器不能直接生成嵌入式代码代数环必须清零代码生成阶段会直接报错MATLAB Function勾选代码生成只用基本运算避免字符串、eval记忆长度 N100~200太大占内存MCU 吃不消之后做模型在环验证生成代码接回 Simulink 和原模型对比误差应该在 1e-6 量级。上 HIL 前再打点日志把控制参数和关键状态量记录好。这套 FOSMC 模型、脚本和参数模板已经整理成可以直接跑的资源拿到手之后建议先按第 4 章流程跑一遍 μ0 退化校验确认模型本身没问题再开始调参数。从那以后我拿到任何一套 FOSMC 模型第一件事就是跑这个校验再看 jitter 有没有异常这套流程走通了才敢往控制器里加自己的东西。希望帮到你。本文还有配套的精品资源点击获取
