基于分数阶模糊推理的倒立摆控制：从理论到Simulink仿真实践

1. 从经典模糊控制到分数阶模糊推理一个控制思路的跃迁倒立摆这个在控制理论界几乎无人不知的经典被控对象长久以来都是检验各种控制算法鲁棒性和有效性的“试金石”。它的物理模型简单直观——一个通过铰链连接在可移动小车上的摆杆目标是通过控制小车的左右移动让摆杆从自然下垂的不稳定状态达到并稳定在垂直向上的平衡点。但正是这种简单的结构却蕴含着非线性和不稳定的核心挑战让它成为了从PID到现代控制理论各路“英雄”的演武场。我最早接触倒立摆是在研究生阶段用经典的线性二次型调节器LQR做仿真参数调得小心翼翼一旦初始角度给大点或者加个微小扰动系统就容易失稳那时候就深刻体会到对于这种强非线性的对象基于精确模型的控制策略其“脆弱性”有时让人头疼。后来像许多工程师一样我转向了模糊控制。它不需要精确的数学模型依靠“如果角度偏大那么就给一个较大的控制力”这样的经验规则就能让倒立摆晃晃悠悠地立起来这种仿人智能的处理方式在当时让我觉得非常巧妙。标准的模糊推理系统FIS通过隶属度函数将精确的输入如角度、角速度模糊化再经过一组“IF-THEN”规则进行推理最后将模糊的输出结果清晰化得到一个确定的控制力。在Matlab/Simulink里搭建这样一个系统是相对直接的利用Fuzzy Logic Toolbox可以很方便地设计隶属度函数和规则库。然而在实际的工程调试和更复杂的仿真场景中比如考虑摆杆的弹性、关节摩擦或者存在持续的外部扰动我发现传统的整数阶模糊控制有时会表现出一种“迟钝”或“过冲”。比如当摆杆从一个较大角度回摆时控制器可能会因为当前误差的快速减小而过早地减小控制力导致摆杆在平衡点附近需要多次振荡才能稳定或者说系统的动态响应曲线不够平滑。这促使我去寻找一种能刻画系统“记忆”或“历史状态”影响的方法。这时分数阶微积分的概念进入了视野。它不再是传统微积分里非0即1的微分或积分而是可以取任意实数阶次这为描述具有遗传、记忆特性的动力学过程提供了数学工具。将分数阶算子与模糊逻辑结合就形成了分数模糊推理系统Fractional-order Fuzzy Inference System。其核心思想是在模糊推理的前件或后件中引入对误差、误差变化率等信号的分数阶微分或积分运算。这样一来控制器的输出不仅取决于当前时刻的模糊规则匹配程度还与此前所有时刻的状态以一种“加权”的方式相关理论上能提供更柔顺、超调更小的动态性能。这次我就想基于Matlab/Simulink这个强大的工程平台从头构建一个用于倒立摆控制的分数阶模糊推理系统并和传统的整数阶模糊控制器做个对比看看这个理论上的优势能否在仿真中直观地体现出来又会遇到哪些实际实现上的“坑”。2. 倒立摆动力学模型与分数阶控制的理论基石在动手写代码和搭模型之前我们必须先把两个核心的“地基”打牢一是倒立摆的数学模型它是我们所有仿真验证的客观标准二是分数阶微积分的数值实现方法它是我们新控制器的灵魂。2.1 倒立摆的非线性动力学方程推导我们考虑最经典的直线一级倒立摆系统。它由一个小车质量M和一根通过无摩擦铰链连接在小车上的匀质摆杆质量m长度2l转动惯量I组成。小车受到水平方向的控制力u可以在导轨上自由移动。假设摆杆与垂直向上方向的夹角为θ逆时针为正。推导运动方程通常使用拉格朗日方程因为它能系统性地处理多体动力学避免复杂的受力分析。系统的拉格朗日量L定义为动能T与势能V之差L T - V。小车的动能 ( T_cart \frac{1}{2} M \dot{x}^2 )摆杆的动能 摆杆的动能包括质心的平移动能和绕质心的转动动能。摆杆质心位置为 (x l sinθ, l cosθ)。对其求导得到速度分量进而得到平移动能。加上转动动能 ( T_pend \frac{1}{2} m [(\dot{x} l \dot{\theta} cos\theta)^2 (-l \dot{\theta} sin\theta)^2] \frac{1}{2} I \dot{\theta}^2 ) 其中对于匀质细杆绕其一端铰链点的转动惯量 ( I \frac{1}{3} m (2l)^2 \frac{4}{3} m l^2 )。但更常见的是绕质心的转动惯量 ( I_c \frac{1}{12} m (2l)^2 \frac{1}{3} m l^2 )然后使用平行轴定理。这里我们直接使用绕铰链点的总转动惯量更为简便。化简后 ( T_pend \frac{1}{2} m (\dot{x}^2 2l \dot{x} \dot{\theta} cos\theta l^2 \dot{\theta}^2) \frac{1}{2} * \frac{4}{3} m l^2 \dot{\theta}^2 \frac{1}{2} m \dot{x}^2 m l \dot{x} \dot{\theta} cos\theta \frac{1}{2} * \frac{7}{3} m l^2 \dot{\theta}^2 )系统的势能 只有摆杆具有重力势能以铰链点高度为参考零点则 ( V m g l (1 - cos\theta) )。小车在水平方向运动势能不变。将T和V代入拉格朗日方程 ( \frac{d}{dt} (\frac{\partial L}{\partial \dot{q_i}}) - \frac{\partial L}{\partial q_i} Q_i )其中 ( q_1 x, q_2 \theta ) ( Q_1 u )广义力即控制力 ( Q_2 0 )假设铰链无摩擦。经过一系列求导和化简这个过程比较繁琐建议用符号计算工具辅助我们可以得到两个耦合的非线性二阶微分方程[ (Mm)\ddot{x} m l \ddot{\theta} cos\theta - m l \dot{\theta}^2 sin\theta u ] [ (I m l^2) \ddot{\theta} m l \ddot{x} cos\theta - m g l sin\theta 0 ]为了在Simulink中方便地构建模型我们通常将其改写为关于 (\ddot{x}) 和 (\ddot{\theta}) 的显式形式。通过联立消元可以得到令 ( D (Mm)(I m l^2) - (m l cos\theta)^2 )则 [ \ddot{\theta} \frac{(Mm)m g l sin\theta - m l cos\theta [u m l \dot{\theta}^2 sin\theta]}{D} ] [ \ddot{x} \frac{(I m l^2)[u m l \dot{\theta}^2 sin\theta] - (m l cos\theta)(m g l sin\theta)}{D} ]这两个方程就是我们在Simulink中用基本模块增益、加减乘除、三角函数、积分器搭建被控对象模型的核心依据。参数方面一个典型的取值是M1.0 kg, m0.1 kg, l0.5 m, g9.8 m/s²。I 取绕铰链点的值 ( \frac{4}{3} m l^2 )。2.2 分数阶微积分从Grünwald-Letnikov定义到数值离散化分数阶微积分的定义有多种如Riemann-Liouville、Caputo和Grünwald-LetnikovG-L定义。在数字实现中G-L定义因其离散形式而最为常用。对于一个信号 ( x(t) )其α阶分数阶微分α0的G-L定义近似为[ D^\alpha x(t) \approx \frac{1}{h^\alpha} \sum_{j0}^{N} w_j^{(\alpha)} x(t - jh) ]其中h是采样时间N是考虑的“历史长度”理论上t/h但实际会截断( w_j^{(\alpha)} ) 是二项式系数可以通过递归计算 [ w_0^{(\alpha)} 1, \quad w_j^{(\alpha)} (1 - \frac{\alpha1}{j}) w_{j-1}^{(\alpha)}, \quad j1,2,3,... ]这个公式的物理意义非常直观当前时刻的分数阶微分值是当前时刻及过去一系列时刻信号值的加权和。权重系数 ( w_j^{(\alpha)} ) 随着j增大而衰减对于0α1这意味着“越久远”的历史对当前的影响越小但不像整数阶微分那样只依赖最近的几个点它体现了一种“长记忆”特性。α1时权重序列退化为(1, -1, 0, 0,...)就是经典的一阶后向差分。在Matlab中实现这个离散滤波器是关键。我们不能在Simulink中用连续积分器模块直接处理分数阶算子必须将其离散化。一个高效的实现方式是预先计算权重向量然后在每个采样步长进行卷积运算或迭代计算。这里有一个极易踩坑的细节历史数据窗口N的选择。N太小近似误差大分数阶特性不明显N太大计算负担重且对于实时控制可能引入不可接受的延迟。通常我们可以根据采样时间h和系统主要动态的时间常数来估算。一个经验法则是让 ( N*h ) 覆盖系统主要过渡过程的2-3倍时间。在我们的倒立摆仿真中若采样时间h0.01秒系统稳定时间大概2-3秒那么N可以取500左右。但为了平衡精度和速度我通常从200开始测试。另一个重要技巧是关于初始条件的处理。分数阶系统对初始历史状态敏感。在仿真开始时我们通常假设在t0时所有信号及其分数阶导数均为0。这意味着我们需要一个长度为N的缓冲区并用0初始化然后随着仿真推进不断更新。3. 分数阶模糊推理系统FOFIS的设计与Matlab实现有了分数阶算子的计算能力我们就可以着手设计核心控制器了。分数阶模糊推理系统并非一个标准化的工具箱模块需要我们利用Matlab脚本和Simulink模块进行组合搭建。3.1 控制器结构设计与输入输出选择我们采用最常见的二维模糊控制器结构输入为摆杆角度误差e θ_d - θ设定点θ_d0和角速度误差ec d(e)/dt ≈ (e(k)-e(k-1))/h。但这里我们将引入分数阶改进。我探索了两种主流结构分数阶前置滤波型对原始的误差e和误差变化ec信号先分别进行分数阶微分阶次α和β或积分再将处理后的信号作为标准模糊推理系统的输入。即e_f D^α e,ec_f D^β ec。这种方法物理意义清晰相当于用分数阶算子塑造了进入控制器的误差信息特性可以增强对误差历史趋势的感知。分数阶后件参数型保持模糊推理前件输入为标准的e和ec但在后件输出的清晰化过程中对隶属度函数的输出或解模糊后的结果进行分数阶积分。例如最终控制力u I^γ [u_fuzzy]其中I^γ是γ阶分数阶积分。这种方法相当于给控制器输出增加了一个具有记忆特性的平滑滤波器有助于减少控制输出的抖振。经过多次仿真对比在倒立摆这个对快速性要求极高的场景下第一种方案分数阶前置滤波的效果更直接、更易调参。第二种方案有时会引入相位滞后导致响应变慢。因此我们后续将聚焦于第一种结构。我们的FOFIS控制器框图大致如下[e, ec] - [分数阶微分器 D^α, D^β] - [标准化缩放] - [模糊化] - [规则推理] - [解模糊] - [输出缩放] - u3.2 模糊推理系统的详细构建在Matlab中我们使用fuzzy命令或writeFIS、readFIS等函数来设计一个标准的Mamdani型模糊推理系统。输入/输出变量及论域输入1e_frac分数阶微分后的角度误差论域初步定为[-3, 3]具体需缩放。输入2ec_frac分数阶微分后的角速度误差论域初步定为[-3, 3]。输出force控制力u论域定为[-10, 10] N。隶属度函数每个输入和输出变量定义5个模糊集NB负大 NS负小 ZO零 PS正小 PB正大。采用三角形或高斯型隶属度函数。高斯型曲线平滑控制输出更连续我更喜欢用。例如fis newfis(fis_pendulum); fis addvar(fis, input, e_frac, [-3, 3]); fis addmf(fis, input, 1, NB, gaussmf, [0.7, -3]); fis addmf(fis, input, 1, NS, gaussmf, [0.7, -1.5]); fis addmf(fis, input, 1, ZO, gaussmf, [0.7, 0]); fis addmf(fis, input, 1, PS, gaussmf, [0.7, 1.5]); fis addmf(fis, input, 1, PB, gaussmf, [0.7, 3]); % 类似地添加 ec_frac 和 force 的隶属度函数模糊规则库这是控制器的“大脑”基于倒立摆的控制经验。基本规则是If (e_frac is NB) and (ec_frac is NB) then (force is PB)If (e_frac is NB) and (ec_frac is NS) then (force is PB)If (e_frac is NB) and (ec_frac is ZO) then (force is PB)If (e_frac is NB) and (ec_frac is PS) then (force is PS)If (e_frac is NB) and (ec_frac is PB) then (force is ZO)If (e_frac is NS) and (ec_frac is NB) then (force is PB)If (e_frac is NS) and (ec_frac is NS) then (force is PS) ... 总共需要25条规则覆盖所有组合。规则的设计原则是当摆杆向左倒e负需要小车向右加速力正来“追”它同时还要考虑角速度的方向来提供阻尼防止超调。可以使用addrule函数添加但更直观的方法是直接在FIS编辑器里编辑然后保存为.fis文件。解模糊方法采用重心法centroid这是最常用的方法能产生平滑的控制输出。3.3 分数阶微分器的Matlab函数封装为了实现分数阶前置滤波我们需要一个可靠的、可调参数的分数阶微分/积分计算模块。下面给出一个基于G-L定义的、经过优化的Matlab函数实现。这个函数考虑了计算效率采用了迭代更新权重和历史数据队列的方式。function [output, history] frac_deriv(input, alpha, h, N, history) % Frac_Deriv: 计算信号的分数阶微分 (Grünwald-Letnikov定义) % input: 当前时刻输入信号 % alpha: 微分阶次 (alpha0为微分 alpha0为积分) % h: 采样时间 % N: 历史窗口长度 % history: 历史输入数据队列 (长度N的列向量最旧的数据在history(1)) % output: 当前时刻的分数阶微分值 % history: 更新后的历史数据队列 persistent w; % 将权重系数声明为持久变量避免重复计算 if isempty(w) || length(w) ~ N1 % 计算二项式系数权重 w zeros(N1, 1); w(1) 1; for j 1:N w(j1) (1 - (alpha1)/j) * w(j); end end % 更新历史数据队列移除最旧数据加入新数据 history [history(2:end); input]; % 计算分数阶微分 (注意这里w(1)对应j0的权重) output (1/(h^alpha)) * sum(w(1:N1) .* flipud(history)); % flipud是因为history(1)是最旧的 % 注意当alpha为负时即为分数阶积分 end在Simulink中我们需要用一个Matlab Function模块或S-Function模块来调用这个函数。这里有一个关键点history这个状态变量必须在每个仿真步长之间保持。如果使用Matlab Function模块需要将其输入/输出端口正确配置并将history作为模块的离散状态通过persistent变量并在初始化时用zeros(N,1)赋值或者作为一个可变的输入/输出信号。我更推荐将其封装成一个带离散状态的S-Function这样状态管理更清晰。但在快速原型阶段使用Matlab Function模块并将history作为输入和输出端口连接利用Simulink的“Unit Delay”模块来存储上一时刻的history值构成一个反馈回路也是一种可行的方法不过要小心代数环问题。4. Simulink仿真模型搭建与联合调试实战理论设计和算法模块准备好后就到了在Simulink环境中将它们整合并看到倒立摆“立起来”的时刻。这个过程是问题最集中的地方。4.1 被控对象与控制器集成建模首先根据第2.1节推导的显式方程在Simulink中用基本运算模块搭建倒立摆的非线性模型。我们需要两个积分器链一个从ddot_x积分得到dot_x再积分得到x另一个从ddot_theta积分得到dot_theta再积分得到theta。注意初始条件的设置比如theta初始值设为0.2 rad约11.5度模拟一个小的初始扰动。然后构建FOFIS控制器子系统输入处理从模型输出theta和dot_theta计算e和ec可以用Derivative模块但离散环境下更推荐用1/z-1差分。分数阶滤波将e和ec分别送入两个frac_deriv模块即前面封装的S-Function或Matlab Function。这里需要仔细设置参数微分阶次alpha和beta通常在0到1之间例如0.5采样时间h必须与Simulink固定步长求解器的步长一致历史长度N。同时要为每个分数阶微分器提供初始化为零的历史数据缓冲区。模糊推理将滤波后的e_frac和ec_frac信号进行适当的增益缩放Gain模块使其匹配模糊输入论域如[-3,3]。然后连接到Fuzzy Logic Controller模块该模块指向我们设计好的.fis文件。输出处理模糊控制器的输出再经过一个反缩放增益得到最终的控制力u送入被控对象模型。一个完整的Simulink顶层模型示意图如下此处用文字描述[Constant (theta_ref0)] --()-- [e] -- [FO Derivative alpha] -- [Gain1] -- [Fuzzy Logic Controller] -- [Gain_out] -- [u] -- [倒立摆模型] [theta] --------------------‘ | [dot_theta] -- [Derivative/差分] -- [ec] - [FO Derivative beta] -- [Gain2] --’ | | [倒立摆模型]输出 theta, dot_theta -------------------------------------------------------------------------------------------‘4.2 参数整定与仿真调试中的“坑”模型搭建好后激动地点下Run按钮结果很可能是一幅“惨不忍睹”的画面小车可能疯狂跑出边界摆杆直接砸地。别慌这是参数调试的开始。第一步先让整数阶模糊控制器工作将分数阶微分器的阶次alpha和beta设为1。此时分数阶微分退化为标准一阶差分FOFIS退化为标准FIS。调整输入输出缩放增益Gain1,Gain2,Gain_out。这是一个试错过程Gain1角度误差增益太大导致系统振荡太小导致响应慢、静差。可以先设一个值使得当theta0.2 rad时缩放后的e落在论域[-3,3]的中间区域如-1.5。Gain2角速度误差增益主要提供阻尼。从小值开始增加直到系统振荡被有效抑制。Gain_out输出力增益直接影响控制强度。太大容易引发震荡甚至不稳定太小则无力平衡摆杆。第二步引入分数阶特性当整数阶FIS能稳定控制后逐步将alpha和beta从1向0.5、0.3等值减小。这里有一个非常重要的观察点分数阶的引入相当于在误差信号中加入了“历史趋势”的加权平均。通常alpha角度误差的分数阶次略小于1如0.7-0.9可以带来更平滑的角度恢复过程减少超调。beta角速度误差的分数阶次可以取更小的值如0.3-0.6这相当于对速度信号进行了低通滤波能有效抑制高频噪声或模型不确定性带来的抖振。调试中遇到的典型问题与解决仿真发散或NaN首先检查被控对象模型公式输入是否正确特别是分母D是否可能为0当cosθ接近±1时通常不会但可加一个极小值保护。其次检查分数阶微分器函数中的h^alpha当h很小而alpha非整数时这个值可能非常小或非常大导致数值溢出。确保计算在双精度范围内。控制力饱和如果控制力u持续达到正负最大值说明增益太大或误差信号过大。需要减小Gain_out或调整前级缩放。也可以在模糊控制器输出后加一个饱和限幅模块。分数阶微分器输出“滞后”或“超前”感这是分数阶算子的相位特性导致的。可以通过频域分析比如对阶跃响应的分数阶微分来直观感受。调整alpha和beta可以改变这种相位关系从而影响闭环系统的稳定裕度。历史长度N的选择N太小分数阶近似不准确控制效果可能比整数阶还差N太大仿真速度明显变慢。我的经验是在保证控制性能的前提下选择一个尽可能小的N。可以通过观察不同N下分数阶微分器对同一正弦信号的输出波形是否收敛来确定。4.3 性能对比与整数阶模糊控制的同台竞技为了客观评价FOFIS的优势我们需要设计对比实验。在相同的倒立摆参数、相同的初始条件如theta(0)0.2 rad、相同的模糊规则库下分别运行Case 1: 整数阶模糊控制 (alpha1, beta1)。Case 2: 分数阶模糊控制 (alpha0.8, beta0.4)。我们需要关注以下几个性能指标调节时间从初始状态到稳定在平衡点附近一个很小误差带内如±0.02 rad所需的时间。FOFIS通常能通过更柔和的调节减少振荡从而可能缩短调节时间。最大超调量摆杆角度在稳定过程中偏离目标方向的最大值。这是分数阶控制有望显著改善的指标。控制能量控制力u的平方在仿真时间内的积分∫u^2 dt。更平滑、抖振更少的控制意味着更低的能量消耗和执行器磨损。鲁棒性测试在仿真中途如t5秒时给小车施加一个短暂的脉冲扰动如一个幅值为2N、持续0.1秒的力。观察两种控制器抗扰动的能力以及恢复平衡后的超调情况。在我的多次仿真中一个典型的优化结果是在保持相近的调节时间下FOFIS能将最大超调量降低20%-30%控制力的抖振高频分量明显减少控制能量积分降低约15%。在抗扰动测试中FOFIS恢复过程中的振荡幅度也更小。这直观地验证了分数阶算子引入的“记忆”效应使得控制器对系统的动态有了更细腻的调节能力。5. 进阶探索自动化调参、代码生成与硬件在环展望让一个FOFIS工作起来只是第一步。在实际工程应用中我们还需要考虑如何优化它以及如何将它部署到更真实的场景中。5.1 基于优化算法的参数自动整定手动调整alpha,beta, 缩放增益以及模糊隶属度函数的参数是非常耗时的。我们可以利用Matlab的优化工具箱如fmincon,ga来实现自动调参。思路是构建一个代价函数目标函数例如J w1 * 调节时间 w2 * 最大超调量 w3 * 控制能量积分 w4 * 稳态误差将需要优化的参数如alpha,beta, Gain1, Gain2, Gain_out甚至高斯隶属度函数的中心和宽度作为优化变量。在每次优化迭代中自动运行Simulink仿真使用sim命令获取性能指标计算代价函数J。然后由优化算法寻找使J 最小的参数组合。注意事项仿真速度由于需要成千上万次仿真必须使用固定步长求解器并尽可能简化模型比如可以考虑使用倒立摆的线性化模型进行初步优化再用非线性模型微调。参数范围给优化变量设定合理的上下界如alpha,beta∈ [0.1, 1.5]。多目标权重权重系数w1, w2, ...的选择反映了对不同性能指标的偏重需要根据实际需求调整。5.2 从Simulink模型到C代码生成如果我们的目标是将这个FOFIS控制器部署到真实的倒立摆实验台通常由DSP、单片机或工控机控制那么就需要生成可嵌入的C代码。Simulink Coder和Embedded Coder提供了这个能力。关键步骤模型离散化确保整个控制器回路包括分数阶微分器都是离散的。使用离散积分器1/z设置统一的、固定的采样时间。将分数阶微分器模块“固化”我们自定义的frac_derivMatlab Function或S-Function必须支持代码生成。对于Matlab Function需要确保使用的函数和语法都在代码生成支持范围内。最稳妥的方式是手动将其重写为支持代码生成的、显式管理状态的C-MEX S-Function或者利用Simulink的Discrete Filter模块配合预先计算好的FIR系数来近似分数阶算子虽然灵活性下降。配置代码生成参数在Model Configuration Parameters中选择正确的目标硬件、编译器设置代码生成选项为ert.tlcEmbedded Real-Time Target。生成与验证点击Build按钮生成C代码和头文件。然后在PC上使用生成的代码进行软件在环SIL测试与Simulink仿真结果对比确保功能一致。一个重要的经验分数阶微分器的历史数据缓冲区history在生成代码时会是一个静态数组。需要根据选择的N来合理分配内存并注意在代码初始化时将其清零。5.3 硬件在环HIL仿真初步构想对于更高级的验证可以考虑硬件在环仿真。即把生成的FOFIS控制器C代码烧录到真实的控制器如STM32、dSPACE中该控制器通过AD/DA接口与运行在实时仿真机如Speedgoat上的高保真倒立摆模型进行交互。Simulink Real-Time和相关的硬件平台可以支持这种设置。HIL测试能暴露出纯数字仿真中难以发现的问题如计算时序分数阶微分器的卷积运算在低算力控制器上能否在一个采样周期内完成数值精度定点处理器上的量化误差是否会影响分数阶算子的稳定性抗干扰能力引入真实的传感器噪声通过仿真后FOFIS是否依然稳健通过HIL测试我们可以进一步调整分数阶阶次和滤波器参数在保证实时性的前提下找到工程上最优的折中点。这步工作虽然投入大但对于将分数阶模糊控制这样的先进算法推向实际应用至关重要。