模糊控制——（四）模糊自适应整定PID控制

hvp

已于 2024-08-20 09:40:33 修改

阅读量6.1k

点赞数 48

文章标签：自动化

于 2024-06-25 14:38:22 首次发布

本文链接：https://blog.csdn.net/qq_40619699/article/details/139793640

版权

系列文章目录

模糊控制——（一）理论基础
 模糊控制——（二）设计流程
 模糊控制——（三）模糊洗衣机
 模糊控制——（四）模糊自适应整定PID控制
 模糊控制——（五）Sugeno模糊模型

关键字：模糊控制；PID；MATLAB；Simulink

文章目录

系列文章目录
前言
一、控制原理
二、设计模糊控制器
三、MATLAB实现
总结
参考文献

前言

提示：这里可以添加本文要记录的大概内容：

例如：随着人工智能的不断发展，机器学习这门技术也越来越重要，很多人都开启了学习机器学习，本文就介绍了机器学习的基础内容。

提示：以下是本篇文章正文内容，下面案例可供参考

一、控制原理

如今PID依然是在工业中使用最为广泛的基础控制方法，PID参数的整定方法很多，但大多数都以=对象特性为基础。随着计算机技术的发展，PID的参数也可以以操作人员的控制经验与知识为基础，通过计算机自动整定，这种控制方法就是专家PID控制，在之前的文章专家PID控制（一）也有简单的介绍，感兴趣的小伙伴可以开启传送门。不过由于操作者的经验不易精确描述，控制过程中各种信号量与评价指标不易定量表示，专家PID的控制方法也受到了限制。
模糊理论被认为是解决上述问题的一个有效途径，运用模糊数学的基本理论和方法，把规则的条件、操作用模糊集表示，并把这些模糊控制规则及有关信息（如评价指标、初始PID参数等）入知识库，计算机根据控制系统的实际响应（专家系统的输入条件），运用模糊推理整定PID控制器参数。
Alt

图1.1 模糊PID控制原理框图

自适应模糊控制器中多以误差 $e$ 、误差变化率 $ec$ 作为输入，注意这里的误差变化率 $ec=e\left( k \right)-e\left( k-1 \right)$ 并没有除采样时间，针对PI控制器上述图1.1中的结构可以简化为如图1.2的模型：
Alt

图1.2 模糊PI控制原理框图

其中 $r_{in}$ 、 $e$ 、 $k_p$ 、 $k_i$ 、 $k_d$ 、 $u$ 和 $y_{out}$ 分贝表示系统输入、误差、PID控制器的比例系数、PID控制器的积分系数、PID控制器的微分系数、控制器输入和系统输出。
想要理解模糊控制器调整参数的原理，首先需要对典型的二阶系统有初步的认识，这里以传递函数为：
$G\left( s \right)=\frac{10}{s^2+2s} \tag{1}$ 的二阶系统为例，观察其阶跃响应的动态特性。在Simulink中进行如图1.3的仿真建模
Alt

图1.3 典型二阶系统的节约响应Simulink建模

在MATLAB中利用m文件绘制系统输出 $c$ 、误差 $e$ 和误差的变化率 $ec$ 的图像，代码如下：

c=out.c;
e=out.e;
ec=ones(1001,1);
[i,~]=size(ec);
ec(1)=0;
for n=2:1:i
    ec(n)=e(n)-e(n-1);
end
figure(1);
plot(c);
xlim([0 1001]);
figure(2);
plot(e);
xlim([0 1001]);
figure(3);
plot(ec);
xlim([0 1001]);

注意，这里首先需要运行Simulink中的仿真将系统输出与误差导入工作区。最后分别得到：
Alt

图1.4 典型二阶系统阶跃响应的系统输出

Alt

图1.5 典型二阶系统阶跃响应的误差

Alt

图1.6 典型二阶系统阶跃响应的误差变化率

Alt

图1.7 典型二阶系统阶跃响应的系统输出、误差、误差变化率

将以上三个图像做如图1.7的处理，得到比例系数 $k_p$ 的整定规则：

在 $\sim t_1$ 时间段内 $e$ 从正值（+）向零值（0）变化，此时系统输出处于上升阶段。 $e$ 的值为正大（P），为了减小 $e$ ，增大比例系数 $k_p$ ，即 $\Delta k_p$ 取正大（P）
在 $t_1 \sim t_2$ 时间段内 $e$ 在0附近变化零值（Z），此时需要观察 $ec$ ， $ec$ 的值为负大（N），此时超调由不断增大的趋势，故减小比例系数 $k_p$ ，即 $\Delta k_p$ 取负大（N）
在 $t_2 \sim t_3$ 时间段内误差从0向负值（-）变化，此时系统输出处于超调阶段， $e$ 的值为负大（N），为了减小超调，减小比例系数 $k_p$ ，即 $\Delta k_p$ 取负大（N）
在 $t_3 \sim t_4$ 时间段内误差从负值（-）向0变化，此时系统输出处于超调阶段， $e$ 的值为负大（N），为了减小超调，减小比例系数 $k_p$ ，即 $\Delta k_p$ 取负大（N）
在 $t_4 \sim t_5$ 时间段内 $e$ 在0附近变化零值（Z），此时需要观察 $ec$ ， $ec$ 的值为正大（P），此时系统输出小于设定输出并有减小的趋势，故增大比例系数 $k_p$ ，即 $\Delta k_p$ 取正大（P）
在 $t_7$ 时刻 $e$ 在0附近变化零值（Z）， $ec$ 也在0附近变化零值（Z），此时系统输出处于极值，为了减小误差，增大比例系数 $k_p$ ，即 $\Delta k_p$ 取正大（P）

根据上述规则可以写出整定 $k_p$ 的模糊规则表：
Alt

图1.8 整定比例系数模糊规则表

对于比例系数 $k_i$ 的整定，选取积分分离策略——当被控量与设定值偏差较大时，取消积分作用；当被控量接近给定值时，引入积分控制，以消除静差，提高精度。故，整定 $k_i$ 的模糊规则表：
Alt

图1.9 整定积分系数模糊规则表

二、设计模糊控制器

1.工程环境

以下的案例基于：

PC操作系统：WIN11
MATLAB：2019a

2.控制器结构

模糊控制器的输入为误差 $e$ 、误差变化率 $ec$ 输出为比例系数变化量 $\Delta k_p$ 、积分系数变化量 $\Delta k_i$ ，故控制器为两输入两输出结构。

3.输入输出模糊集

其实在上文中，我们几乎已经确定了输入输出的模糊集：
$\begin{align} e &=\{N,Z,P\}=\{-1,0,1\} \\ ec &=\{N,Z,P\}=\{-1,0,1\} \\ \Delta k_p &=\{N,Z,P\}=\{-1,0,1\} \\ \Delta k_i &=\{N,Z,P\}=\{-1,0,1\} \end{align} \tag{2}$

4.输入输出隶属函数

首先给出被控对象的传递函数：
$G_p\left( s \right)=\frac{133}{s^2+25s} \tag{3}$ 我们要研究的是模糊整定PID在阶跃输入 $r\left( t \right)=1$ 情况下的系统输出特性。类似于上文提到的方法，先在Simulink中搭建模型：
Alt

图2.1 二阶传递函数建模

需要提醒的是以上变量导入到工作表中注意初始默认的数据类型是时间序列，如果需要数组数据类型则需要进行参数设置：
Alt

图2.2 以二维数组的形式导出数据到工作区

而且需要注意的是导出数据的长度是由Simulink中的步长与仿真时间计算得到的，例如博主这里的仿真时间1.5s，固定步长0.01s：
Alt

图2.3 设定仿真时间与步长

这样得到的数据长度就是 $1.5 \div 0.01 + 1 = 151$ ：
Alt

图2.4 导入工作区的数据长度

随后在MATLAB中绘图，代码如下：

c=out.c;
e=out.e;
ec=ones(151,1);
[i,~]=size(ec);
ec(1)=0;
ts=1;
for n=2:1:i
    ec(n)=(e(n)-e(n-1))/ts;
end
figure(1);
plot(c);
xlim([0 151]);
figure(2);
plot(e);
xlim([0 151]);
figure(3);
plot(ec);
xlim([0 151]);

得到误差 $e$ 和误差变化率 $ec$ 的图像：
Alt

图2.5 误差

Alt

图2.6 误差变化率

可以看到两者的变化范围，都在 $\left[-1,1\right]$ 以内，设定误差 $e$ 、误差变化率 $ec$ 、比例系数 $k_p$ 和积分系数的隶属函数：
Alt

图2.7 误差的隶属函数

Alt

图2.8 误差变化率的隶属函数

Alt

图2.9 比例系数的隶属函数

Alt

图2.10 积分系数的隶属函数

5.模糊规则

模糊规则可根据图1.7下的文字描述得到：
$\begin{align} &IF \quad e \quad is \quad N \quad and \quad ec \quad is \quad N \quad then \quad \Delta k_p \quad is \quad N \quad and \quad \Delta k_i \quad is \quad Z \\ &IF \quad e \quad is \quad N \quad and \quad ec \quad is \quad Z \quad then \quad \Delta k_p \quad is \quad N \quad and \quad \Delta k_i \quad is \quad Z \\ &IF \quad e \quad is \quad N \quad and \quad ec \quad is \quad P \quad then \quad \Delta k_p \quad is \quad N \quad and \quad \Delta k_i \quad is \quad Z \\ &IF \quad e \quad is \quad Z \quad and \quad ec \quad is \quad N \quad then \quad \Delta k_p \quad is \quad N \quad and \quad \Delta k_i \quad is \quad P \\ &IF \quad e \quad is \quad Z \quad and \quad ec \quad is \quad Z \quad then \quad \Delta k_p \quad is \quad P \quad and \quad \Delta k_i \quad is \quad P \\ &IF \quad e \quad is \quad Z \quad and \quad ec \quad is \quad P \quad then \quad \Delta k_p \quad is \quad P \quad and \quad \Delta k_i \quad is \quad P \\ &IF \quad e \quad is \quad P \quad and \quad ec \quad is \quad N \quad then \quad \Delta k_p \quad is \quad P \quad and \quad \Delta k_i \quad is \quad Z \\ &IF \quad e \quad is \quad P \quad and \quad ec \quad is \quad Z \quad then \quad \Delta k_p \quad is \quad P \quad and \quad \Delta k_i \quad is \quad Z \\ &IF \quad e \quad is \quad P \quad and \quad ec \quad is \quad P \quad then \quad \Delta k_p \quad is \quad P \quad and \quad \Delta k_i \quad is \quad Z \\ \end{align} \tag{4}$

6.模糊控制表

模糊控制表如图1.8与图1.9。

7.模糊推理

模糊推理包含了规则匹配、规则触发、规则前提推理、规则总可信度和规则总输出，具体步骤可召唤传送门模糊控制——（三）模糊洗衣机。

8.反模糊化

将得到的规则总输出利用反模糊函数进行计算，的到确切值的过程称为反模糊化，详情可见上一步的传送门。

三、MATLAB实现

首先定义一个m文件用于设计模糊控制器，并保存模糊接口系统文件

% Fuzzy Control Design for PID controller
% 为PID参数整定设计的模糊控制器
clc;
clear;
close all;

a=mamfis();
a.Name="FuzzyPI";
%% 
% 创建输入与输出模糊变量并设置隶属度函数
% 创建第一个输入模糊变量与隶属度函数
a.Input(1).Name='e';
a.Input(1).Range=[-1,1];
a=addMF(a,'e','zmf',  [-1,-1/3],   'Name','N');
a=addMF(a,'e','trimf',[-2/3,0,2/3],'Name','Z');
a=addMF(a,'e','smf',  [1/3,1],     'Name','P');
% 创建第二个输入模糊变量与隶属度函数
a.Input(2).Name='ec';
a.Input(2).Range=[-1,1];
a=addMF(a,'ec','zmf',  [-1,-1/3],   'Name','N');
a=addMF(a,'ec','trimf',[-2/3,0,2/3],'Name','Z');
a=addMF(a,'ec','smf',  [1/3,1],     'Name','P');
% 创建第一个输出模糊变量与隶属度函数
a.Output(1).Name='kp';
a.Output(1).Range=(1/3)*[-10,10];
a=addMF(a,'kp','zmf',  (1/3)*[-10,-3],'Name','N');
a=addMF(a,'kp','trimf',(1/3)*[-5,0,5],'Name','Z');
a=addMF(a,'kp','smf',  (1/3)*[3,10],  'Name','P');
% 创建第二个输出模糊变量与隶属度函数
a.Output(2).Name='ki';
a.Output(2).Range=(1/30)*[-3,3];
a=addMF(a,'ki','zmf',  (1/30)*[-3,-1], 'Name','N');
a=addMF(a,'ki','trimf',(1/30)*[-2,0,2],'Name','Z');
a=addMF(a,'ki','smf',  (1/30)*[1,3],   'Name','P');
%% 
% 创建并添加模糊规则
rulelist=[1 1 1 2 1 1;
          1 2 1 2 1 1;
          1 3 1 2 1 1;
          
          2 1 1 3 1 1;
          2 2 3 3 1 1;
          2 3 3 3 1 1;
          
          3 1 3 2 1 1;
          3 2 3 2 1 1;
          3 3 3 2 1 1;];
a=addRule(a,rulelist);                                                      % 添加模糊规则
showrule(a)                                                                 % 显示模糊规则
%% 
% 设定反模糊函数
a.defuzzification='centroid';
%% 
% 保存模糊文件
writeFIS(a,'fuzzyPI');
%%
% 用模糊工具箱进行仿真
a_new=readfis('fuzzyPI');                                              %读取模糊文件
disp('---------------------------------------------------------------');
disp('fuzzy controller  table:e=[-1,1], ec=[-1,1]');
disp('---------------------------------------------------------------');
fuzzy FuzzyPI;                                                             %模糊规则设计器
ruleview FuzzyPI;                                                          %规则查看器
u=evalfis([-1,0],a_new)
%%
% 绘图
figure(1);
plotfis(a_new);
figure(2);
plotmf(a_new,'input',1);
figure(3);
plotmf(a_new,'input',2);
figure(4);
plotmf(a_new,'output',1);
figure(5);
plotmf(a_new,'output',2);

随后在另外一个m文件中定义传递函数并利用保存的fis文件加载设计的模糊控制器，实现模糊整定PI控制

% Fuzzy PI controller
% 为洗衣机设计的模糊控制器
clc;
clear;
close all;
%%
% 加载模糊接口系统文件
a_new=readfis('FuzzyPI');
%%
% 离散函数
ts=0.001;
sys=tf(133,[1,25,0]);
% dsys=c2d(sys,ts,'zoh');
dsys=c2d(sys,ts,'z');                                                      %将连续时间传递函数转为离散时间传递函数，其中‘z’表示用“零阶保持法“离散
[num,den]=tfdata(dsys,'v');                                                %返回离散传递函数的分母与分子系数，其中‘v’表示结果用vector向量表示，如果不用‘v’则结果是ceil元胞数组
%%
% 初始化参数
kp0=0;
ki0=0;
u_1=0;u_2=0;
c_1=0;c_2=0;
e_1=0;
ec_1=0;
ei_1=0;


kp=ones(1,1000);
ki=ones(1,1000);
time=ones(1,1000);
r=ones(1,1000);
u=ones(1,1000);
c=ones(1,1000);
e=ones(1,1000);
%%
% 在迭代中更参数
for k=1:1:1000
    time(k)=k*ts;
    r(k)=1;
    k_pid=evalfis([e_1,ec_1],a_new);
    kp(k)=kp0+k_pid(1);
    ki(k)=ki0+k_pid(2);
    u(k)=kp(k)*e_1+ki(k)*ei_1;
    c(k)=-den(2)*c_1-den(3)*c_2+num(2)*u_1+num(3)*u_2;
    e(k)=r(k)-c(k);
    u_2=u_1;u_1=u(k);
    c_2=c_1;c_1=c(k);
    ec_1=e(k)-e_1;
    e_1=e(k);
    ei_1=ei_1+e(k)*ts;
end

figure(1);
plot(time,r,'r',time,c,'b:','linewidth',2);
xlabel('time/s');ylabel('r,c');
legend('Ideal postion','Practical postion');
figure(2);
subplot(211);
plot(time,kp,'r','linewidth',2);
xlabel('time/s');ylabel('kp');
subplot(212);
plot(time,ki,'r','linewidth',2);
xlabel('time/s');ylabel('ki');
figure(3);
plot(time,u,'r','linewidth',2);
xlabel('time/s');ylabel('u');