电机控制算法综述03-滑模控制

文章目录
  1. 1. 滑模控制原理
  2. 2. 滑模面的构建
    1. 2.1. 确定系统的模型
    2. 2.2. 明确控制目标
    3. 2.3. 滑模面的确定原则
      1. 2.3.1. 跟踪误差阶次匹配
      2. 2.3.2. 收敛速度控制
      3. 2.3.3. 物理可实现性
      4. 2.3.4. 鲁棒性优化
  3. 3. 趋近率的确定
    1. 3.1. 常用的趋近率
    2. 3.2. 稳定性证明(李雅普诺夫方法)
    3. 3.3. 高频抖振抑制
  4. 4. 基于滑模控制的PID控制器Matlab代码

这一专栏将介绍电机控制算法,为接下来的课题研究提供参考。
本文首先介绍滑模控制原理,然后介绍采用滑模控制优化PID参数的电机控制系统Matlab仿真代码。

滑模控制原理

滑模控制(Sliding Mode Control, SMC)是一种非线性控制方法,通过设计特定滑动模态,使系统状态轨迹在有限时间内到达并保持在预设滑模面上,具有强鲁棒性和抗扰动特性。
滑模面(Sliding Surface),也称为切换面,是一个由设计者定义的超平面。滑模面是由状态变量及其导数组成的数学表达式,它决定了系统的行为模式。一旦系统的状态轨迹到达这个表面,理论上应该沿此表面滑向期望的目标点。
滑动模态(Sliding Mode)指的是当系统状态进入并保持在滑模面上时所表现出的一种特殊动态行为。在这种状态下,尽管外部扰动或参数变化可能会影响系统性能,但由于采用了特定的设计策略,使得系统能够维持一种相对稳定的运动趋势。

滑模面的构建

要想建立滑动模态,首先要构建合适的滑模面,而构建合适滑模面的前提是明确系统模型与控制目标。

确定系统的模型

电机的动力学系统(转子机械运动方程式)是一个二阶系统,其表达式如下。
\(T_e=J\frac{d^2 \theta}{dt^2}+B\frac{d\theta}{dt}+T_L=J\frac{d \omega}{dt}+B\omega+T_L\)
式中转速\(\omega=\frac{d\theta}{dt}\)
取转角\(\theta=x_1\)、转速\(\omega=x_2\),可以将其改写为状态方程的形式如下:
\(\left\{ {\begin{array}{*{20}{c}} { {\dot x_1} = x_2} \\ { {\dot x_2} = -\frac{B}{J}x_2 + \frac{1}{J}u - \frac{1}{J}d(t)} \end{array}} \right.\)
式中\(J\)为转动惯量,\(B\)为粘性摩擦系数。
\(u=T_e\)为电磁转矩,也是滑模控制器的输出量,被控对象(电机动力学系统)的输入量。
\(d(t)=T_L\)为负载转矩,也可以理解成系统的扰动。
这个二阶系统当然也可以写为:\(J{\ddot \theta}+B{\dot \theta}=u-d(t)\)

明确控制目标

这里要实现的目标是转速闭环控制,即转速误差为零。
定义转速误差\(e=\omega_{ref}-\omega\),式中\(\omega_{ref}\)是参考转速,也可以说是给定转速。
确定控制目标为\(e=0,{\dot e}=0\),在误差为零的同时确保误差加速度也为零,这样做可以提高系统抗扰动能力。

滑模面的确定原则

根据上述系统模型与控制目标,确定滑模面为:
\(s={\dot e}+\lambda e \quad (\lambda > 0)\)
下面介绍确定滑模面所需要遵循的原则。

跟踪误差阶次匹配

若系统阶次为\(n\),则滑模面的微分方程阶次为\((n-1)\)
这样可以保证滑动模态方程降阶为\((n-1)\)维,避免零动态不稳定问题,确保滑模运动与系统阶次匹配。
因此对于二阶系统,设计的是一阶滑模面。

收敛速度控制

式中\(\lambda\)决定误差\(e\)的收敛速度。
\(\lambda\)越大,系统收敛越快,但需满足\(\lambda_{max}<\frac{1}{\tau_s}\)
式中\(\tau_s\)为系统最小时间常数
工程上通常取\(\lambda=(3 \sim 5)/t_{set}\),式中\(t_{set}\)为期望调节时间。

注:关于时间常数
时间常数\(\tau\)是描述系统惯性特性的核心参数,表示系统响应达到稳态值63.2%所需的时间。
对于包含惯性、阻尼和刚度(弹性)的完整二阶机械系统,其动力学方程为
\(J{\ddot \theta}+B{\dot \theta}+K \theta=T\)
其时间常数根据阻尼比及其响应特性的不同可分为三种情况(此处略去不表)。
对于上文提及的二阶电机动力学系统,由于其不含刚度项(忽略弹性形变), 该系统的传递函数\(G(s)=\frac{1}{s(Js+B)}\)
其特征方程为\(s(Js+B)=0\),对应能求出两个极点。
一个是零极点(\(s=0\)),对应积分环节,反映位置累积特性,无时间常数。
另一个是惯性极点(\(s=-B/J\)),对应一阶惯性环节,时间常数\(\tau=\frac{J}{B}\)
综上所述,不含刚度项的二阶机械系统的时间常数与一阶机械系统相同,为\(\tau=\frac{J}{B}\)
对于典型一阶电气系统,其时间常数\(\tau=\frac{L}{R}\)
对于电机控制系统,其通常会包含多个子系统,其时间常数通常取各个子系统时间常数的最小值,即最小时间常数。

物理可实现性

滑模面表达式必须满足如下条件
1.仅包含可测量状态量;
2.不引入不可测高阶导数项;
3.满足传感器精度限制。

鲁棒性优化

滑模控制的鲁棒性优化旨在设计控制律,使得系统在存在外部扰动模型不确定性参数摄动时,仍能保持稳定并实现期望的动态性能。
其本质是通过滑模切换项的增益设计,抵消系统的不确定性影响,确保滑模运动始终收敛到预设的滑模面。

扰动上界\(D\):系统所受扰动(包括外部干扰、未建模动态、参数误差等)的最大幅值。
要求总扰动\(|d(t)|<D \quad t \ge 0\)
安全裕度\(\Delta\):为应对扰动上界估计误差或未预见干扰而额外增加的增益裕度。
\(\Delta = k_D \cdot D \quad k_D \in [0.1,0.3]\)
\(k_D\)过小,可能导致鲁棒性不足;若\(k_D\)过大,则会加剧高频抖振,增加能耗。
切换增益\(\eta\):控制律中切换项(符号函数项)的增益系数,用于补偿扰动。
要求\(\eta \ge D + \Delta\)

趋近率的确定

趋近律(Reaching Law)是滑模控制的核心组成部分,用于描述系统状态从初始位置向滑模面趋近的动态过程。其设计直接影响系统的收敛速度、抗干扰能力及控制输入的平滑性。趋近律通过定义状态变量\(s\)(即滑模面)的微分方程,确保系统轨迹在有限时间内到达并保持在滑模面上。

常用的趋近率

等速趋近律\({\dot s}=-\eta \cdot sign(s) \quad \eta>0\)
(线性收敛,速度较慢)
指数趋近律\({\dot s}=-ks -\eta \cdot sign(s) \quad \eta>0,k>0\)
(线性收敛,速度较慢。增大\(k\)可加速趋近,但需避免过大导致抖振加剧)
幂次趋近率\({\dot s}=-k|s|^\alpha sign(s) \quad k>0,1>\alpha>0\)
(有限时间收敛,速度最快,典型值取\(\alpha=0.5\)

稳定性证明(李雅普诺夫方法)

李雅普诺夫第二方法基本原理:若存在标量函数\(V(x)\)对状态空间中所有非零状态点\(x\)满足\(V(x)>0\)(正定)且\({\dot V}(x)<0\)(负定),则系统在原点渐近稳定。
基于上述原理分析滑模控制的稳定性,首先要选取合适的标量函数\(V(x)\)
所构造的标量函数要求其对于\(x\)具有连续一阶偏导数,且\(V(0)=0\)
这里取\(V=\frac{1}{2}s^2\)
显然所取的函数符合条件\(V(0)=0\)\(V(s)>0\),下面对其是否符合\({\dot V}(s)<0\)做进一步验证。
\({\dot V}=s{\dot s}=s[\eta sign(s)-ks] \le -\eta |s| \quad (\eta>0)\)
(这里采用指数趋近率)
由此可知所取的标量函数符合要求\({\dot V}(s)<0\),进而证明了所构造滑模控制系统的稳定性。

高频抖振抑制

为抑制滑模控制所存在的抖振问题,可以将趋近率中的符号函数替换为饱和函数。
\(sat(s/\Phi)=\left\{ {\begin{array}{*{20}{c}} {s/\Phi}&{|s|\le \Phi} \\ {sign(s)}&{|s|>\Phi} \end{array}} \right.\) 式中\(\Phi\)为边界层厚度,其选取需兼顾跟踪精度与平滑性。

另外,需要合理选择趋近率中的\(\eta\)以兼顾抗扰能力与平滑性。
\(\eta\)过大,则会使切换项剧烈变化,执行器高频颤振;若\(\eta\)过小,则会无法克服扰动,造成跟踪误差发散。

基于滑模控制的PID控制器Matlab代码

使用DeepSeek生成如下所示的Matlab代码:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
% 滑模控制优化PID参数的直流电机转速控制
clc; clear; close all;

% 电机参数
J = 0.0027; % 转动惯量
B = 0.0004924; % 阻尼系数

% 滑模控制参数
lambda = 80; % 滑模面系数(增大可加快收敛)【50->80】
eta = 20; % 切换项增益(增强抗扰动能力)【15->20】

% 自适应率系数
alpha = 100; % Kp调整系数【500->800】
beta = 10; % Ki调整系数【300】
gamma = 4; % Kd调整系数【100->150】

% 参考转速信号 (阶跃信号)
dtheta_ref = @(t) 10.0*(t >= 0.5); % 目标转速10rad/s
ddtheta_ref = @(t) 0; % 参考加速度为零

% 初始状态 [theta, theta_dot, integral_e, Kp, Ki, Kd]
x0 = [0; 0; 0; 30; 2; 2]; % 调整PID初始参数

% 仿真时间
tspan = [0 2];

% 解微分方程
[t, x] = ode45(@(t,x) motorSystem(t, x, J, B, lambda, eta, alpha, beta, gamma, dtheta_ref, ddtheta_ref), tspan, x0);

% 提取结果
theta_dot = x(:,2); % 电机实际转速
Kp = x(:,4);
Ki = x(:,5);
Kd = x(:,6);

% 绘制结果
figure;
subplot(2,1,1);
plot(t, theta_dot, 'b', t, dtheta_ref(t), 'r--');
title('转速跟踪');
legend('实际转速', '参考转速');
xlabel('时间(s)'); ylabel('转速(rad/s)');
grid on;

subplot(2,1,2);
plot(t, Kp, 'r', t, Ki, 'g', t, Kd, 'b');
title('PID参数自适应调整');
legend('Kp', 'Ki', 'Kd');
xlabel('时间(s)'); ylabel('参数值');
grid on;

% 系统微分方程定义(转速控制版本)
function dxdt = motorSystem(t, x, J, B, lambda, eta, alpha, beta, gamma, dtheta_ref, ddtheta_ref)
theta = x(1);
theta_dot = x(2);
integral_e = x(3);
Kp = max(x(4), 0); % 强制非负
Ki = max(x(5), 0);
Kd = max(x(6), 0);

% 计算转速跟踪误差
e = dtheta_ref(t) - theta_dot;
e_dot = ddtheta_ref(t) - (theta_dot*B/J); % 等价于 -theta_ddot

% 滑模面设计
s = e_dot + lambda*e;

% PID控制输出
u_pid = Kp*e + Ki*integral_e + Kd*e_dot;

% 滑模控制切换项(带饱和函数)
sat_thickness = 0.5;
if abs(s) <= sat_thickness
u_sw = eta*(s/sat_thickness);
else
u_sw = eta*sign(s);
end

% 总控制输入
u = u_pid + u_sw;

% 改进的自适应律(带死区和非负约束)
dead_zone = 0.05;
if abs(e) > dead_zone
dKpdt = alpha*e*s;
dKidt = beta*integral_e*s;
dKddt = gamma*e_dot*s;
else
dKpdt = 0;
dKidt = 0;
dKddt = 0;
end

% 系统动态方程
theta_ddot = (u - B*theta_dot)/J;

dxdt = [theta_dot;
theta_ddot;
e;
dKpdt;
dKidt;
dKddt];
end

该程序采用Matlab自带的求解器obe45对电机动力学系统微分方程进行求解,其使用变步长进行连续求解,其状态变量\(x\)包括转角\(\theta\),转速\(\dot \theta\),误差\(e\)的积分以及PID控制器参数\(K_p,K_i,K_d\)
该函数标准格式为[t, Xt] = ode45(odefun, tspan, X0)
其中odefun代指所导入的函数,可以是函数句柄、函数文件名、匿名函数句柄或内联函数名。
在程序中写作@(t,x) motorSystem(t, x, J, B, lambda, eta, k, alpha, beta, gamma, dtheta_ref, ddtheta_ref),其中@(t,x)为匿名函数,motorSystem是匿名函数所指代的函数。
tspan为仿真的时间区间,可以连续也可以离散,程序中设为[0 2],表示仿真时长为2s。
X0为状态变量初始值向量,t为返回列向量的时间点,Xt为返回对应t的求解列向量。

程序前半部分对各个参数进行定义,其中令转动惯量J=0.0027kg·m^2,粘性摩擦系数B=0.0004924N·m·s。该值与DeepSeek最初生成值不同,我按照Simulink中的永磁同步电机仿真模型默认值对其修改,该参数更符合电机控制的实际情况。
对于滑模控制参数,定义了滑模面系数\(\lambda\)与切换增益\(\eta\)
还有影响\(K_p,K_i,K_d\)更新速率的自适应系数\(\alpha,\beta,\gamma\),若要使控制系统\(K_p,K_i,K_d\)参数一直保持初始值不变,即进行普通的PID控制,可以将自适应系数全部设为0,以便于控制系统初始PID参数调试。

接下来对转速参考信号\(\dot \theta_{ref}\)与转速的导数参考信号\(\ddot \theta_{ref}\)进行定义,这里用到了匿名函数。
转速参考信号最初为0,在0.5s时刻为阶跃信号,给定转速10rad/s。
理论上在0.5s时刻其导数值为无穷大(Dirac函数),但在实际仿真中,ode求解器可能无法处理这种情况,因此这里转速的导数一直保持为0。

接下来研究motorSystem函数,其输入量包括时间\(t\),状态变量\(x\),以及各个参数\(J,B,\lambda,\eta,\alpha,\beta,\gamma,{\dot \theta_{ref}},{\ddot \theta_{ref}}\)
输出量dtxt是输入量状态变量\(x\)的微分量(对时间\(t\)求导),包括\({\dot \theta},{\ddot \theta},e,\frac{\text{d}K_p}{\text{d}t},\frac{\text{d}K_i}{\text{d}t},\frac{\text{d}K_d}{\text{d}t}\)

在该函数中做出了下述定义:
转速误差\(e={\dot \theta_{ref}} - {\dot \theta}\)
转速导数的误差\({\dot e}={\ddot \theta_{ref}} - \frac{B}{J}{\dot \theta}\)
滑模面\(s={\dot e}+\lambda e \quad (\lambda > 0)\)
总控制输入\(u=u_{pid}+u_{sw}\)
等效控制项\(u_{pid}=K_p \cdot e + K_i \cdot \int e \text{d}t + K_d \cdot {\dot e}\),基于PID控制维持标称系统性能。
切换项\(u_{sw}=\eta \cdot sat(s/\Phi)\),用于补偿系统不确定性和扰动。
切换项基于前文所述的等速趋近律构建,并将符号函数替换为饱和函数以抑制滑模控制所存在的抖振问题。
\(K_p,K_i,K_d\)的自适应率 \(\left\{ {\begin{array}{*{20}{c}} {\frac{\text{d}K_p}{\text{d}t}=\alpha * e * s} \\ {\frac{\text{d}K_i}{\text{d}t}=\beta * \int e * s} \\ {\frac{\text{d}K_d}{\text{d}t}=\gamma * {\dot e} * s} \end{array}} \right.\)
程序中对\(K_p,K_i,K_d\)的值进行了非负约束,通过代码K=max(K,0),使之保持为正数。
此外还添加了死区约束,当转速误差小于一个较小的值dead_zone时,将\(\frac{\text{d}K_i}{\text{d}t},\frac{\text{d}K_d}{\text{d}t}\)设为0,使得\(K_p,K_i,K_d\)的值不再发生变化。
电机动力学系统动态方程\({\ddot \theta} = (u - B*{\dot \theta})/J\)
该动态方程根据前文所述系统模型移项改写得到。