题目简述
控制一个三阶传递函数的阶跃响应，采样周期为1ms，传递函数如下：
$$
G(s)=\frac{523500}{s^3+87.35s^2+1047s}
$$
- 编程实现上述系统的单位阶跃响应，并画出波形图
- 编程实现用传统PID对上述系统进行校正，要求设置10组参数，给出结果图
- 编程实现用专家PID对上述系统进行校正，给出结果图
理论知识
传统PID校正
加入了反馈后的PID控制：
$$
e_{ss}=\lim_{t\to \infty}e(t)=\lim_{s\to 0}sE(s)=\lim_{s\to 0}\frac{sR(s)}{1+PID(s)\times G(s)}
$$
其中$PID(s)$是$PID$的传递函数，$R(s)$是输入，当$sE(s)$在收敛域内时，一般来说可以看到当信号稳定后，误差$e$是趋近于0的。
专家PID校正
首先先将连续信号离散化，专家控制系统有5大规则：
$|e(k)|>M_1$，此时控制输出信号为常数，相当于进行开环控制。
$e(k)\Delta e(k)>0$或者$\Delta e(k)=0$
$|e(k)|>M_1$，此时控制PID输出为：
$$
u(k)=u(k-1)+k_1{k_p[e(k)-e(k-1)]+k_ie(k)+k_d[e(k)-2e(k-1)+e(k-2)]}
$$$|e(k)|<M_2$，此时控制PID输出为：
$$
u(k)=u(k-1)+k_p[e(k)-e(k-1)]+k_ie(k)+k_d[e(k)-2e(k-1)+e(k-2)]
$$
$e(k)\Delta e(k)<0$，$\Delta e(k)\Delta e(k-1)>0$或者$\Delta e(k)=0$,此时控制PID输出如下：
$$
u(k)=u(k-1)
$$$e(k)\Delta e(k)<0$，$\Delta e(k)\Delta e(k-1)<0$
$|e(k)|>M_2$，此时控制PID输出为：
$$
u(k)=u(k-1)+k_{1}k_{p}e(k)
$$$|e(k)|<M_2$，此时控制PID输出为：
$$
u(k)=u(k-1)+k_{2}k_{p}e(k)
$$
$|e(k)<\epsilon|$，此时采用PI控制，输出为：
$$
u(k)=u(k-1)+k_p[e(k)-e(k-1)]+k_ie(k)
$$
问题求解
求解单位阶跃响应
由于所求阶跃响应是时域响应，所以需要将s域的响应变换到时域，基本思路如下：
先利用式子
$$
Y(s)=G(s)\times{\frac{1}{s}}
$$
求出该系统在输入为单位阶跃信号时s域的响应再利用
matlab中的ilaplace函数将响应$Y(s)$变换到时域，即可求得该系统的单位阶跃响应
代码如下：
clear;clc;
syms s t;
% 传递函数
G=523500/(s^3+87.35*s^2+10470*s);
% 传递函数乘以1/s得到s域的阶跃响应
Y=G*(1/s);
% 阶跃响应变换到时域
y_t=ilaplace(Y,t);
disp(y_t);
将得到的时域的$y(t)$写成函数形式，以便于画图调用：
%%%将y_t函数改写成如下的函数
function fval=y(t)
%%%%求时域下的阶跃响应，输入时间求出响应值
for i=1:length(t)
fval(i)=50*t(i)+(1747*exp(-(1747*t(i))/40)* (cos((13699991^(1/2)*t(i))/40) -(5323991*13699991^(1/2)*sin((13699991^(1/2)*t(i))/40))/23933884277))/4188 - 1747/4188;
end
end
设置采样间隔为1ms，持续时间为0.5s，画出波形图：
%%%画图
t=0:0.001:0.5;
figure (1)
plot(t,y(t),'linewidth',2)
xlabel("时间/s");
ylabel("阶跃响应值");
title("系统G(s)的阶跃响应图")
得到的波形图如下：
传统PID控制校正
结果展示
选取不同的PID参数，画出图像，效果如下：
结果分析
从上面的10张图像可以得出一些结论：
- 如果kp很大，输出曲线会出现振荡的情况，见图（1）（3）（4）（5）
- 如果kp很小，输出曲线回收敛得很慢，见图（9）
- 如果kp适中，那么收敛速度会较快，也不会出现很强烈的振荡现象，见图（6）（7）（10）
- 如果ki很大，输出曲线的最大值会超过1，ki太小输出曲线最大值回小于1，见图（8）（9）
- 如果kd很大，输出曲线会出现值变得很大的情况，见图（2）
综上：当三个参数都选取恰当时，会使输出曲线像一个一阶系统，见图（6）（10）
专家PID控制校正
在传统PID的基础上加上五大规则，即构成了专家PID控制
结果展示
当选取kp、ki、kd三个参数选取合适时，改变五大规则里的参数，比较同一组PID参数下不同规则的结果（由于规则参数较多，在这就不详细说明了），展示如下：
结果分析
由此可见，在同一组kp、ki、kd参数的条件下，专家PID校正需要一定的专业经验知识才能调参更好
代码展示
求阶跃响应
clear;clc;
syms s t;
% 传递函数
G=523500/(s^3+87.35*s^2+10470*s);
% 传递函数乘以1/s得到s域的阶跃响应
Y=G*(1/s);
% 阶跃响应变换到时域
y_t=ilaplace(Y,t);
disp(y_t);
%%%画图
t=0:0.001:0.5;
figure (1)
plot(t,y(t),'linewidth',2)
xlabel("时间/s");
ylabel("阶跃响应值");
title("系统G(s)的阶跃响应图")
%%%将a中的函数改写成如下的时域函数
function fval=y(t)
%%%%求时域下的阶跃响应，输入时间求出响应值
for i=1:length(t)
fval(i)=50*t(i)+(1747*exp(-(1747*t(i))/40)*(cos((13699991^(1/2)*t(i))/40) -(5323991*13699991^(1/2)*sin((13699991^(1/2)*t(i))/40))/23933884277))/4188 - 1747/4188;
end
end
传统PID
%%传统PID控制
clear;clc;
ts = 0.001; %采样时间
sys = tf(5.235e005,[1,87.35,1.047e004,0]); %传递函数
dsys = c2d(sys,ts,'z'); %连续模型离散化
[num,den] = tfdata(dsys,'v'); %获得分子分母
u_1=0;u_2=0;u_3=0;
y_1=0;y_2=0;y_3=0;
x = [0,0,0]';
x2_1 = 0;
%比例积分微分参数
kp = 0.5;
ki = 0.02;
kd =0.2;
error_1 = 0;
for k = 1:1:500 %运行时间0.5s
time(k)=k*ts;
r(k) = 1.0; %第k次控制器输出
u(k) = kp*x(1) +kd*x(2) + ki*x(3);
%线性变换 Z变换
y(k) = -den(2)*y_1-den(3)*y_2-den(4)*y_3+...
num(1)*u(k)+num(2)*u_1+num(3)*u_2+num(4)*u_3;
error(k) = r(k) - y(k);
%Return of parameter
u_3=u_2;u_2=u_1;u_1=u(k);
y_3=y_2;y_2=y_1;y_1=y(k);
x(1)=error(k); %计算P
x2_1=x(2);
x(2)=(error(k)-error_1)/ts; %计算D
x(3)= x(3)+error(k)*ts; %计算 I
error_1=error(k);
end
figure(1);
plot(time,r,'b',time,y,'r','linewidth',2);
xlabel('时间(s)');ylabel('输出');
grid on
title({'传统PID控制阶跃响应曲线',['kp=',num2str(kp),',ki=',num2str(ki),...
',kd=',num2str(kd)]})
% figure(2);
% plot(time,r-y,'r','linewidth',2);
% xlabel('times(s)');ylabel('error');
% grid on
% title('误差响应曲线')
专家PID
展示效果最好的那一组专家规则，即第一组规则参数
%%专家PID控制
clear;clc;
ts = 0.001; %采样时间
sys = tf(5.235e005,[1,87.35,1.047e004,0]); %传递函数
dsys = c2d(sys,ts,'z'); %连续模型离散化
[num,den] = tfdata(dsys,'v'); %获得分子分母
u_1=0;u_2=0;u_3=0;
y_1=0;y_2=0;y_3=0;
x = [0,0,0]';
x2_1 = 0;
%比例积分微分参数
kp = 0.6;
ki = 0.03;
kd =0.002;
error_1 = 0;
for k = 1:1:500 %运行时间0.5s
time(k)=k*ts;
r(k) = 1.0; %第k次控制器输出
u(k) = kp*x(1) +kd*x(2) + ki*x(3);
%规则1：控制误差输出的绝对值，避免超调
if abs(x(1))>0.8
u(k) = 0.45;
elseif abs(x(1))>0.4
u(k) = 0.40;
elseif abs(x(1))>0.2
u(k) = 0.12;
elseif abs(x(1))>0.01
u(k) = 0.10;
end
%规则2：误差为某一参数，未发生变化
if x(1)*x(2)>0 || (x(2)==0)
if abs(x(1))>=0.05
u(k)=u_1+ 2*kp*x(1);
else
u(k)=u_1+ 0.4*kp*x(1);
end
end
%规则3：误差的绝对值朝小的方向变化
if (x(1)*x(2)<0 && x(2)*x2_1>0)||(x(1)==0)
u(k)=u(k);
end
%规则4：误差处于极值状态
if x(1)*x(2)<0 && x(2)*x2_1<0
if abs(x(1))>=0.05
u(k)=u_1+2*kp*error_1;
else
u(k)=u_1+0.6*kp*error_1;
end
end
%规则5：误差的绝对值很小
if abs(x(1))<=0.001 %Pi控制
u(k)=0.5*x(1)+0.01*x(3);
end
if u(k)<=-10
u(k)=-10;
end
if u(k)>=10
u(k)=10;
end
%线性变换 Z变换
y(k) = -den(2)*y_1-den(3)*y_2-den(4)*y_3+num(1)*u(k)+num(2)*u_1+...
num(3)*u_2+num(4)*u_3;
error(k) = r(k) - y(k);
%Return of parameter
u_3=u_2;u_2=u_1;u_1=u(k);
y_3=y_2;y_2=y_1;y_1=y(k);
x(1)=error(k); %计算P值
x2_1=x(2);
x(2)=(error(k)-error_1)/ts; %计算D值
x(3)= x(3)+error(k)*ts; %计算 I值
error_1=error(k);
end
figure(1);
p1=plot(time,r,'b','linewidth',2);
hold on;
p2=plot(time,y,'r','linewidth',2);
xlabel('时间(s)');ylabel('输出');
grid on
title({'专家PID控制阶跃响应曲线',['kp=',num2str(kp),',ki=',num2str(ki),...
',kd=',num2str(kd)]})
hold on;
p3=plot(time,r-y,'k','linewidth',2);
legend([p2,p3],'输出值','error','location','best')
