电机控制算法综述01-1-PID控制器的Matlab仿真实现

文章目录
  1. 1. matlab代码实现PID控制
    1. 1.1. 自建离散PID函数,配合差分方程形式的被控对象
    2. 1.2. 自建离散PID函数,配合差分方程形式的电机被控对象
    3. 1.3. matlab内置PID函数,配合传递函数形式的被控对象
    4. 1.4. matlab内置函数ode45,配合微分方程形式的电机被控对象
  2. 2. Simulink实现PID控制
    1. 2.1. 采用最基本的模块构建PID模块
    2. 2.2. Simulink自带PID模块
    3. 2.3. S-function建立PID模块

这一专栏将介绍电机控制算法,为接下来的课题研究提供参考。
本文将总结各个PID控制器实现方法(纯matlab代码、Simulink自带模块、S函数),被控对象的不同建模方法(传递函数、微分方程、Simulink自带离散或连续模块),分析各个方法在不同方面的优劣。
后面的文章将结合《先进PID控制MATLAB仿真》书籍相关内容,对PID控制器参数的传统整定方法(Ziegler-Nichols方法等)进行介绍。

matlab代码实现PID控制

自建离散PID函数,配合差分方程形式的被控对象

这一章节采用单个matlab脚本文件,编写针对差分方程形式的被控对象进行离散PID控制的matlab程序,不涉及Simulink仿真和S函数。
文中程序改编自RBF神经网络PID控制器的设计与仿真
对于任何一个控制系统仿真实验,无论用的是什么软件,进行的是什么仿真,第一个要考虑的问题就是仿真步长和仿真时间。
这里取步长(采样时间)Ts=0.01s,仿真时间T=1s。与此同时,就可以计算出程序运行的总步数N=T/Ts=100。
接下来构建被控对象,该系统的核心公式为:
\(y(k+1) = a_1 \cdot y(k) + a_2 \cdot y(k-1) + b_1 \cdot u(k) + b_2 \cdot u(k-1)\)
这是一个二阶离散系统,式中y(k+1)就是下一时刻被控对象的输出量,如果被控对象是电机,那么y可以是转速;u(k)是PID控制器输出量(被控对象的输入量),如果被控对象是电机,u的物理意义根据控制系统实际情况会有所不同,可能是电压、电流、转矩等。
该系统的传递函数为\(G(z)=\frac{Y(z)}{U(z)}=\frac{b_1 z^{-1}+b_2 z^{-2}}{1-a_1 z^{-1}-a_2 z^{-2}}\)

接下来是PID控制器设计,采用的是位置式PID,其核心公式为:
\(u(k)=K_p \cdot e(k) + K_i \cdot \sum e(k) + K_d \cdot ec(k)\)
与PID参数相乘的三项分别是误差、误差累积和与误差变化量,它们的计算公式为:
\(\left\{ \begin{gathered} e(k)=r(k)-y(k) \\ \sum e(k) = \sum e(k-1) + T_s \\ ec(k)=e(k)-e(k-1) \end{gathered} \right.\)
注:一些文献中的误差变化量公式为\(ec(k)=\frac{e(k)-e(k-1)}{T_s}\)
然后采用max与min函数将控制量u限幅至-10至10之间:
\(u(k)=\max(\min(u(k), 10), -10)\)

离散被控对象+离散PID的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
% 自建离散PID函数,配合差分方程形式的被控对象
clc;clear;close all;
% 仿真参数
Ts=0.01;T=1;N=T/Ts;
% 被控对象(二阶离散系统)
a1 = 1.5; a2 = -0.7; b1 = 0.2; b2 = 0.1;
% PID参数
Kp = 0.2;Ki = 30;Kd = 1;
% 初始化系统变量
y = zeros(1, N+1); % 输出 1对应t=0s N对应t=0.99s N+1对应t=1s
u = zeros(1, N+1); % 控制量
e = zeros(1, N+1); % 误差
ec = 0;
sum_e = 0; % 误差积分
% 参考信号(阶跃信号)
r = zeros(1, N+1);
step_time = 0.5; % 阶跃时间0.5秒
step_index = round(step_time / Ts) + 1; % 转换为索引(round为四舍五入函数)
r(step_index:end) = 1; % 阶跃信号
% 主循环
for k = 2:N+1 %2:101共100个数据 T=0.01s到T=1s
% 计算误差和误差变化率
e(k) = r(k) - y(k);
ec = e(k) - e(k-1);
sum_e = sum_e + e(k) * Ts;
% 计算控制量(PID输出)
u(k) = Kp * e(k) + Ki * sum_e + Kd * ec;
u(k) = max(min(u(k), 10), -10); % 限幅
% 应用控制量到系统
y(k+1) = a1 * y(k) + a2 * y(k-1) + b1 * u(k) + b2 * u(k-1);
end
% 绘制结果
t = (1:N) * Ts; %T=0.01到T=1s共100个数据
figure;
subplot(3,1,1);
plot(t, r(2:N+1), 'r--', t, y(2:N+1), 'b-');%T=0.01s到T=1s共100个数据
xlabel('时间 (s)');ylabel('输出量');
xlim([0 1]);ylim([0 1.2]);
legend('给定值', '系统输出');
title('系统响应');

subplot(3,1,2);
plot(t, e(2:N+1));
xlabel('时间 (s)');ylabel('误差');
xlim([0 1]);ylim([-0.05 0.05]);
title('误差');

subplot(3,1,3);
plot(t, u(2:N+1));
xlabel('时间 (s)');ylabel('控制量');
xlim([0 1]);ylim([0 1.5]);
title('系统输入');
该仿真程序所描述的就是一个最简单的控制系统,如下图所示:
图1 最简单的控制系统
程序中取\(a_1=1.5,a_2=-0.7,b_1=0.2,b_2=0.1\),系统的输出公式为: \(y(k+1)=1.5y(k)-0.7y(k-1)+0.2u(k)+0.1u(k-1)\)
程序中的给定信号为阶跃信号,在0.5s时给定值为1。
程序后半部分绘制仿真所得的波形,包括给定值r与系统输出值y、误差量e、控制量u随时间变化波形,如下图所示:
图2 离散PID控制仿真波形

可以观察到给定值出现阶跃后,系统输出值在0.04s后便达到给定值附近。
而控制量u最终稳定在0.67左右。这是因为当系统稳定时,有y(k+1)=y(k)=y(k-1)=r,u(k)=u(k-1),带入到前面的系统输出公式可以得到y=1.5y-0.7y+0.2u+0.1u=0.8y+0.3u,解得u=2/3,与实验结果相符。
这一现象也说明,系统控制量的稳定值u(∞)只与被控对象有关,与所选取的PID参数无关,当然若PID参数选取的不好会导致系统不稳定,但如果系统稳定对应的控制量稳定值就是唯一的。

自建离散PID函数,配合差分方程形式的电机被控对象

上文仿真程序的被控对象是一个二阶离散系统,下面根据电机运动方程建立离散系统模型。
文中程序改编自电机控制算法综述01-2-模糊PID控制器的Matlab仿真
电机的运动方程如下:
\(T_e-T_L=J\frac{d^2 \theta}{dt^2}+B\frac{d\theta}{dt}=J\frac{d \omega}{dt}+B\omega\)
式中转速\(\omega=\frac{d\theta}{dt}\)\(T_e\)是电磁转矩,\(T_L\)是负载转矩,\(J\)为转动惯量,\(B\)为粘性摩擦系数(有的文献也用\(D\)表示)。
运动方程转换为差分方程形式,程序中的代码如下:
speed = (1 - B*Ts/J)*speed + (Kt*Ts/J)*u;
\(y(t+1)=(1-\frac{B\cdot Ts}{J})y(t)+\frac{K_t \cdot T_s}{J}u(t)\)
这是一个一阶离散系统。其中速度speed就是被控对象的输出量y;u就是控制量,它是PID控制器的输出量,也是被控对象出入量,这里的物理意义就是电压。Ts为采样时间。
对应的传递函数可以写作:\(G(z)=\frac{Y(z)}{U(z)}=\frac{K_t \cdot T_s/J z^{-1}}{1-(1-B\cdot T_s/J) z^{-1}}\)
该差分方程再转换回运动方程的形式,如下:
\(K_t \cdot u=J\frac{y(t+1)-y(t)}{T_s}+B \cdot y(t)\)
对于不同的电机,电压到转矩的转换过程各有不同,这取决于电机的动态特性(与电枢电阻、电感、磁链等相关)以及功率变换器件工作特性,但通常电压与转矩(这里指电磁转矩与负载转矩之差)呈正比关系。
程序中对转换过程进行简化,取\(T_e-T_L=K_t \cdot u\),式中\(K_t\)是转矩常数,也就是电压转换到转矩的比例系数。

离散电机被控对象+离散PID的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
% 离散PID函数,差分方程形式的电机被控对象
clear; clc;
J = 0.01; % 转动惯量
B = 0.1; % 阻尼系数
Kt = 0.15; % 转矩常数
Ts = 0.001; % 采样时间
T = 0:Ts:0.1; % 仿真时间
speed = 0;target_speed = 100; % 目标转速
Kp = 100;Ki = 40;Kd = 0.001;% PID参数
e_prev = 0;integral = 0;
speed_history = zeros(size(T));% 存储结果
u_history = zeros(size(T));

for k = 2:length(T) %2:101共100个数据 T=0.001s到T=0.1s
e = target_speed - speed;
ec = (e - e_prev) / Ts;
integral = integral + e * Ts;
u = Kp * e + Ki * integral + Kd * ec;%计算控制量
u = max(min(u, 220), -220);%电压限幅
% 更新系统状态(简化电机模型)
speed = (1 - B*Ts/J)*speed + (Kt*Ts/J)*u;
speed_history(k) = speed;
u_history(k) = u;
e_prev = e;
end

figure;
subplot(2,1,1);
plot(T, speed_history, 'b', T, target_speed*ones(size(T)), 'r--');
title('转速响应'); xlabel('时间(s)'); ylabel('转速(rpm)'); legend('实际', '目标');
xlim([0 0.1]);ylim([0 120]);
subplot(2,1,2);
plot(T, u_history, 'g');
title('控制量'); xlabel('时间(s)'); ylabel('控制电压(V)');
xlim([0 0.1]);ylim([0 250]);
离散PID和离散电机被控对象的仿真程序如上所示,设定仿真时长T=0.1s,仿真步长Ts=0.001s,目标转速为100r/min。
转动惯量J=0.01,摩擦系数B=0.1,转矩常数Kt=0.15。
PID控制器与前一个程序类似,唯一区别为误差变化量公式为\(ec(k)=\frac{e(k)-e(k-1)}{T_s}\)
控制量限幅设为220,也就是将电压限制在220V以内。
运行该程序,得到电机转速(被控对象输出)与电压(PID控制器输出,被控对象输入)波形如下图所示。
图3 离散PID电机控制仿真波形

可以看到电机在36ms时转速达到给定值100r/min附近。
电机运行可以分为两个阶段,第一个阶段时恒压升速阶段,电压维持在限幅值220V,转速可近似看作是匀速上升。
第二个阶段转速达到稳定值,电压也随之下降,最终稳定在66.67V左右。
这个电压(控制量)稳定值也可以按照之前的方法进行计算。
稳态下有y(∞)=y(t+1)=y(t),程序中的差分方程可转化为:
\(y(\infty)=(1-\frac{B\cdot Ts}{J})y(\infty)+\frac{K_t \cdot T_s}{J}u(\infty)\)
\(\Rightarrow \frac{B\cdot Ts}{J}y(\infty)=\frac{K_t \cdot T_s}{J}u(\infty) \Rightarrow u(\infty)=\frac{B}{K_t}y(\infty)=\frac{0.1}{0.15}\times 100\approx 66.67\)
计算值与实验值相符。

下面进一步研究恒压升速阶段,根据限幅电压与转矩常数,可以计算出这一阶段的电磁转矩与负载转矩之差\(T_e-T_L=0.15 \times 220=33N \cdot m\)
带回到运动方程公式,可得\(K_t \cdot u=33N \cdot m=J\frac{y(t+1)-y(t)}{T_s}+B \cdot y(t)\)
当电机刚刚起动时,转速为0,也就是y(t)=0
此时可计算出\(y(t+1)-y(t)=\frac{33\cdot T_s}{J}=3.3r/min\)
实验波形中t=0.001s时刻的转速为3.3r/min,与计算值相符。
当电机转速接近目标转速时,假设转速达到100r/min
此时可计算出\(y(t+1)-y(t)=\frac{(33-B \cdot 100)\cdot T_s}{J}=2.3r/min\)
实验波形中t=0.034s时刻转速为95.52r/min,t=0.035s时刻转速为97.86r/min,两者转速差为2.34r/min,与计算结果相近。
可见随着转速升高,转速差由3.3r/min降低至2.3r/min,也就是说转速加速度随着转速升高而降低。
这主要是受到了摩擦系数B的影响,相当于是为电机施加了一个额外的负载转矩,且该额外负载转矩大小与转速成正比。
当转速达到100r/min时,摩擦系数产生的等效额外负载转矩值为10N·m。
若取摩擦系数B=0,则电机恒压升速阶段的加速度保持恒定。

matlab内置PID函数,配合传递函数形式的被控对象

文中程序改编自电机控制器算法综述02-1-蚁群算法
前文采用离散传递函数的被控对象建立了离散PID控制系统程序,接下来利用matlab内置pid函数,编写连续传递函数被控对象的PID控制系统程序。
为了便于对比,首先将上一个程序中的差分方程形式的电机被控对象传递函数转换为连续传递函数。
离散传递函数为\(G(z)=\frac{Y(z)}{U(z)}=\frac{K_t \cdot T_s/J z^{-1}}{1-(1-B\cdot T_s/J) z^{-1}}=\frac{0.015z^{-1}}{1-0.99z^{-1}}\)
输入以下代码在matlab中建立传递函数模型,tf函数的三个输入分别为分子、分母和采样时间。
G = tf([0.015],[1 -0.99],0.001)
然后用d2c函数将离散传递函数转换为连续传递函数。
Gc=d2c(G)
得到的连续传递函数为\(G(s)=\frac{15.08}{s + 10.05}\)
连续电机被控对象+matlab内置PID函数的仿真代码(点击展开或收起)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
% matlab内置PID函数,传递函数形式的被控对象
clear; clc;
Kp = 100;Ki = 40;Kd = 0.001;% PID参数
Tf = 0.001; % 滤波器时间常数
controller = pid(Kp, Ki, Kd, Tf); % 使用带滤波器的PID控制器(PIDF)
s = tf('s');
motor = 15.08 / (s+10.05);% 被控对象模型
sys_u = controller / (1 + controller * motor);% 从r到u的传递函数
sys = feedback(controller * motor, 1); % 闭环系统(从r到y的传递函数)
t = 0:0.001:0.1; % 仿真时间
y = step(sys, t); % 系统输出响应
u = step(sys_u, t); % 控制器输出响应

figure;
subplot(2,1,1);
plot(t,y);
plot(t, y, 'b', t, ones(size(t)), 'r--');
title('转速响应'); xlabel('时间(s)'); ylabel('转速(rpm)');legend('实际', '目标');
xlim([0 0.1]);ylim([0 2]);
subplot(2,1,2);
plot(t,u);
title('控制量'); xlabel('时间(s)'); ylabel('控制电压(V)');
xlim([0 0.1]);ylim([0 110]);
该程序采用matlab内置的pid函数,其传递函数为\(C(s)=K_p+\frac{K_i}{s}+\frac{K_d s}{T_f s+1}\)
这里的微分项传统PID控制器多了一个滤波器,滤波系数为\(T_f\)
这是因为原微分项\(K_d s\)物理上不可实现,实际设备需要低通滤波器限制高频增益。
引入滤波器可以使系统成为真有理传递函数,避免分子阶次高于分母阶次的非因果问题,造成仿真报错。
feedback函数的作用是计算带有反馈的闭环系统的传递函数,对于文中的闭环系统,从输入r到输出y的传递函数计算公式为:
\(Y(s)=\frac{C(s)G(s)}{1+C(s)G(s)}\)
step函数的作用是计算该传递函数的阶跃响应。
另外,要想查看PID控制器输出量u随时间变化情况,需要得到从输入r到控制器输出u的传递函数,其公式为:
\(U(s)=\frac{C(s)}{1+C(s)G(s)}\)
同样用step函数得到该传递函数的阶跃响应。
图4 连续与离散系统波形对比

运行仿真程序,就能得到上图左侧的y与u随时间变化曲线。
该程序进行的是阶跃响应,相当于目标转速为1r/min,同时没有对控制器输出u进行限幅,对离散PID程序进行同样的调整,得到上图右侧的转速与电压随时间变化曲线。
可以观察到连续系统的转速与控制量波动更小,曲线更平滑
电机起动时的控制器输出u最大,随后迅速降低,其峰值约为100。
经过实验发现该值受到比例系数\(K_p\)影响,与\(K_p\)大小呈正比。
这是因为在电机起动阶段,转速误差值较大,PID控制器比例部分起主要作用,可以使得转速快速跟踪给定。
而随着转速误差减小,比例部分的作用也会不断变弱,到转速稳定阶段,积分部分开始起主要作用,用于消除稳态静差。

连续与离散系统转速稳定阶段的控制器输出u稳定值皆为0.66,这也与理论计算结果相符。

matlab内置函数ode45,配合微分方程形式的电机被控对象

文中程序改编自电机控制算法综述03-滑模控制
上文采用matlab内置的pid函数建立的连续传递函数被控对象的控制系统存在不能对pid函数进行限幅,只能进行单位阶跃响应的问题,下面采用matlab的微分方程(ordinary differential equation,ODE)求解器,建立微分方程形式的电机被控对象的PID控制程序。
为简化程序,考虑到当前电机模型下PID控制中的微分环节的作用不大,因此这里采用的是PI控制。(实际上是DeepSeek生成的程序中没微分项)

微分方程形式的电机电机被控对象PI控制的仿真代码(点击展开或收起)
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
J = 0.01;      % 转动惯量 (kg.m^2)
b = 0.1; % 阻尼系数 (N.m.s)
K = 0.15; % 电机转矩常数 (N.m/A)
Kp = 100;Ki = 40;
target_speed = 100; % 目标转速 (rad/s)
tspan = [0 0.1]; % 仿真时间范围
x0 = [0; 0]; % 初始状态 [转速; 误差积分]
[t, x] = ode45(@(t,x) motor_dynamics(t, x, target_speed, J, b, K, Kp, Ki), tspan, x0);
% 提取状态变量
speed = x(:, 1); % 电机实际转速
integral_error = x(:, 2); % 误差积分项
% 计算控制信号和误差
error = target_speed - speed;
control_signal = Kp*error + Ki*integral_error;
control_signal = max(min(control_signal, 220), -220);

figure;
subplot(2,1,1);
plot(t, speed, 'b', t, target_speed*ones(size(t)), 'r--');
title('转速响应'); xlabel('时间(s)'); ylabel('转速(rpm)');legend('实际', '目标');
xlim([0 0.1]);ylim([0 120]);
subplot(2,1,2);
plot(t, control_signal, 'g');
title('控制量'); xlabel('时间(s)'); ylabel('控制电压(V)');
xlim([0 0.1]);ylim([0 250]);

% 定义电机动态方程
function dxdt = motor_dynamics(t, x, target_speed, J, b, K, Kp, Ki)
speed = x(1); % 当前转速
integral_error = x(2); % 误差积分项
error = target_speed - speed;% 计算当前误差
u = Kp * error + Ki * integral_error; % 控制电压
u = max(min(u, 220), -220);%电压限幅
dspeed_dt = (K*u - b*speed) / J;% J*dw/dt = K*u - b*w
dxdt = [dspeed_dt; error];% 返回状态导数(误差积分的变化率=当前误差)
end

obe45使用变步长进行连续求解,该函数标准格式为[t, Xt] = ode45(odefun, tspan, X0)
函数的输入有三个,其中odefun代指所导入的函数,可以是函数句柄、函数文件名、匿名函数句柄或内联函数名。
在程序中写作@(t,x) motor_dynamics(t, x, target_speed, J, b, K, Kp, Ki),其中@(t,x)为匿名函数,motor_dynamics是匿名函数所指代的函数。
tspan为仿真的时间区间,可以连续也可以离散,程序中设为[0 0.1],表示仿真时长为0.1s。X0为状态变量初始值向量。
函数的输出有两个,其中t为返回列向量的时间点(由于是变步长,因此各个时间点的间隔不固定),Xt为返回对应t的求解列向量。

本程序中状态变量x包括转速和误差积分这两项。
motor_dynamics函数除了必须有的xt这两个参数外,还带入了另外一些参数,分别是目标转速target_speed、转动惯量J、摩擦系数b、转矩常数K、比例系数Kp与积分系数Ki
返回的求解列向量Xt就是状态变量的导数,或者说是微分值。程序中为转速的变化量dspeed_dt和转速误差值error
其中转速变化量由运动方程移项获得,即
\(T_e-T_L=K_t u=J\frac{d \omega}{dt}+B\omega \Rightarrow \frac{d \omega}{dt}=\frac{K_t u-B\omega}{J}\)
图5 微分方程形式电机PI控制仿真波形

转速与电压波形如上图所示,可见其与图3类似,转速也是在36ms左右达到给定转速附近。
但转速稳定阶段的电压波动相较于图3更大,这可能是变步长的缘故。

Simulink实现PID控制

这里介绍采用Simulink搭建电机PID控制系统的方法,主要分为三种,第一种是采用最基本的模块(Commonly Used Blocks中的Gain、Sum、Delay、Saturation等基础模块),搭建出一个PID控制模块;第二种是采用自带的PID模块;第三种是用S-function建立PID模块。

采用最基本的模块构建PID模块

图6 Simulink搭建PID控制系统

按照前文所述搭建出的电机PID控制系统如上图所示,其电压与转速波形也与图3相吻合。

Simulink自带PID模块

Simulink自带的PID模块可以实现连续时间和离散时间的PID控制算法,支持抗饱和、外部重置和信号跟踪等高级功能。
该模块具体使用方法参见:https://ww2.mathworks.cn/help/simulink/slref/pidcontroller.html
图7 Simulink内置PID模块

该模块时域可设置为连续时间或离散时间,其中连续时间公式为:
\(U(s)=P+I\frac{1}{s}+D\frac{N}{1+N\frac{1}{s}}=P+I\frac{1}{s}+D\frac{s}{\frac{1}{N}s+1}\)
在配置界面中,与I相乘的\(\frac{1}{s}\)和与D相乘的\(\frac{N}{1+N\frac{1}{s}}\)分别被称作积分器滤波器
其中N是滤波器系数,默认值为100,该值与前文matlab内置的PID函数的滤波器系数Tf的关系式为:
\(N=\frac{1}{T_f}\)
时域设置为离散时间的公式如下,式中Ts是采样周期:
\(U(z)=P+I\cdot T_s \frac{1}{z-1}+D\frac{N}{1+N\cdot T_s \frac{1}{z-1}}\)
该式中同样包含滤波器系数N,若取消勾选“使用滤波导数”,则公式为:
\(U(z)=P+I\cdot T_s \frac{1}{z-1}+D\frac{1}{T_s}\frac{z-1}{z}\)
该公式不含滤波系数,与D相乘的部分被称作微分器
PID控制器设置为离散形式时,还可以调整积分器和滤波器方法(连续函数中积分项与微分项中的s用什么样的含z表达式代替)。
默认设置为前向欧拉(Forward Euler):\(\frac{1}{s}=\frac{T_s}{z-1}\)
这种方法适用于采样时间较小的情况,此时奈奎斯特极限远大于控制器的带宽。对于较大的采样时间,使用前向欧拉方法可能会导致不稳定,即使是对一个在连续时间下稳定的系统进行离散化。
设置为后向欧拉(Backward Euler)时,\(\frac{1}{s}=\frac{T_s z}{z-1}\)
这种方法对稳定的连续时间系统进行离散化总能得到一个稳定的离散时间结果。
设置为梯形(Trapezoidal)时,\(\frac{1}{s}=\frac{T_s}{2}\frac{z+1}{z-1}\)
这种方法对稳定的连续时间系统进行离散化总能得到稳定的离散时间结果。在所有可用的积分方法中,梯形方法在离散化系统的频域特性和相应的连续时间系统之间提供了最接近的匹配。

打开PID控制器配置界面,可以看到五个选项卡。
“主要”界面可以设置比例、积分、微分系数和滤波器系数N。
可以设置源的形式,如果设置为外部源,那么上述参数就不是在配置界面中设置,这时PID控制器会多出四个输入端口,分别输入P、I、D、N。
另外,还可以选择自动调节PID参数的方法,实现PID参数自动整定,调节方法包括基于传递函数与基于频率响应两种。
参考链接(步骤清晰)利用MATLAB工具箱自动整定SIMULINK PID参数
最后还有一个选项为“启用过零检测”,过零检测(Zero Crossing Detection)只有采用变步长求解器时才有用。
过零检测的原理是根据采样过程中变量的变化率(前后两个采样点的变化量)来动态调整仿真步长,当前后两个采样点的值变化大时,则缩小采样步长;当前后两个采样点的值变化小时,则增大步长。
这样做可以增大变化率大的时候的采样频率,提高仿真精度;降低变化率小的时候的采样频率,提高仿真速度(效率)。

第二个选项卡“初始化”,可以设置积分器和滤波器的初始值。
外部重置(External reset)的作用是触发外部信号时将积分器和滤波器重置为前面设置的初始值,可选的外部信号有上升沿、下降沿、任一沿、电平信号。
跟踪模式(Tracking mode),这个没怎么看懂。 作用是使得模块输出跟随在TR端口提供的跟踪信号。启用跟踪模式后,跟踪信号与模块输出之间的差值通过增益Kt反馈到积分器输入。
跟踪模式的应用包括无冲击控制转换(Bumpless control transfer)和防止多回路控制结构中的模块饱和

第三个选项卡“饱和”,可以设置PID控制器输出u和积分器的上下限。
另外还可以设置抗饱和方法,包括反算(back-calculation)和钳位(clamping)
“反算”是指当模块输出饱和时,将饱和控制输出与未经饱和处理的原控制输出之间的差值加到积分器。反算系数Kb是加到积分器前要乘以的倍数,默认值为1。
“钳位”又称作有条件积分,当控制器输出超出限制且积分器输出与PID模块输入(误差e)具有相同符号时,积分停止。当控制器输出超出限制且积分器输出与模块输入具有相反符号时,积分恢复。
钳位适合具有相对较小死区时间的系统,对于具有较大死区时间的系统,瞬态性能会较差。

后面两个选项卡“数据类型”与“状态属性”与硬件代码生成相关,这里略。

S-function建立PID模块

S函数(S-function)是System-function系统函数的缩写,是指采用非图形化的方式(计算机语言)描述的一个功能块。
S函数的模块搭建与程序编写方法可以参考单神经元PID控制器的设计与仿真
图8 采用S函数搭建电机PID控制器
图9 S函数模块参数与封装配置

搭建的控制系统模型如上图所示,S函数一共有四个参数,分别是采样时间T和控制参数Kp,Ki,Kd。
状态变量x有三个,分别是误差e、误差累计和sum_e和误差变化量ec,这三个状态变量分别与Kp,Ki,Kd相乘得到PID控制器输出量。
输入量u有两个,分别是设定转速rin与实际转速y。
输出量除了PID控制器输出量外,还有e、sum_e、ec,连接到示波器上可以观察它们随时间的变化。

PID模块S函数代码(点击展开或收起)
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
function [sys,x0,str,ts] = PID_S(t,x,u,flag,T,Kp,Ki,Kd)
% PID控制器S函数实现
% 输入:u(1)=参考转速, u(2)=实际转速
% 输出:y(1)=控制量

switch flag
case 0 % 初始化
[sys,x0,str,ts] = mdlInitializeSizes(T);
case 2 % 离散状态更新
sys = mdlUpdate(t,x,u,T);
case 3 % 输出计算
sys = mdlOutputs(t,x,u,T,Kp,Ki,Kd);
case {1,4,9} % 未使用flag
sys = [];
otherwise
error(['未处理的flag = ',num2str(flag)]);
end

%% 初始化函数
function [sys,x0,str,ts] = mdlInitializeSizes(T)
sizes = simsizes;
sizes.NumContStates = 0; % 无连续状态
sizes.NumDiscStates = 3; % 状态变量x=[e,sum_e,ec]
sizes.NumOutputs = 4; % 输出量[u,e,sum_e,ec]
sizes.NumInputs = 2; % 输入量[r, y]
sizes.DirFeedthrough = 1; % 输入直接影响输出
sizes.NumSampleTimes = 1;
sys = simsizes(sizes);

% 状态变量初始化
x0 = [0,0,0];
str = [];
ts = [T 0]; % 采样时间

%% 状态更新函数
function sys = mdlUpdate(t,x,u,T)
e_prev = x(1);sum_e = x(2);
r = u(1);y = u(2);
e = r - y;
sum_e = sum_e + e*T;
if t == 0
ec = 0;
else
ec = (e - e_prev)/T;
end
sys = [e,sum_e,ec];

%% 输出计算函数
function sys = mdlOutputs(~,x,~,~,Kp,Ki,Kd)
e = x(1);sum_e = x(2);ec = x(3);
u_out = Kp*e + Ki*sum_e + Kd*ec;
sys = [u_out,e,sum_e,ec];
取采样时间为0.001s时,得到的电压与转速波形波动较大。
因此这里取采样时间T=1e-06s,得到电压与转速波形以及e、sum_e、ec随时间变化波形如下图所示。
图10 电压与转速波形
图11 e、sum_e、ec随时间变化波形