写在前面，最近在学自动控制原理，正巧我一个朋（er）友（zi），遇见了如下问题，起初以为只是把网上搜到的C代码重构一下能跑就行，没想到我的C语言基础不足以支持我修复这个漏洞百出的脚本，索性就在matlab重新写了。
以下是我的完成过程，自动控制菜鸟，欢迎交流
原问题：
原答案出处：https://www.docin.com/p-1079784869.html
完成步骤：
1.将G(s)离散化为G(z)(Z变换)。
2.将G(z)改写出离散差分方程。
3.step阶跃，绘制激励-响应曲线。
4.改变T，R0，和差分方式观察区别。
5.构造系统PID控制器输出序列。
6.迭代计算得到响应曲线。
7.根据绘图现象修改PID参数。
matlab代码如下：
%% tustin变换
% Author:HXF Time:2022.5.5 15:36
% syms K1 T1 Ts s z % 定义符号变量
% Gs = K1/(T1*s+1)^2
% GZ = subs(Gs,s,2*(z-1)/Ts/(z+1)) %离散化 tustin变换
% Gz = collect(GZ,z) % 将表达式G2中的以z为变量合并相同次幂；
%% 一阶向前差分
%
% syms K1 T1 Ts s z % 定义符号变量
% Gs = K1/(T1*s+1)^2
% GZ = subs(Gs,s,(z-1)/Ts) %离散化 一阶向前差分
% Gz = collect(GZ,z) % 将表达式G2中的以z为变量合并相同次幂；
%% 一阶向后差分
% % 本实验采用一阶向后差分
% syms K1 T1 Ts s z Yk Yk_1 Yk_2 Rk% 定义符号变量
% Gs = K1/(T1*s+1)^2
% GZ = subs(Gs,s,(1-1/z)/Ts) %离散化 一阶向后差分
% Gz = collect(GZ,z) % 将表达式G2中的以z为变量合并相同次幂；
% [num,den]=numden(Gz)
% a0 = 2*(T1/Ts)*(1+(T1/Ts))/(1+(T1/Ts)) ^2
% b0 = -(T1/Ts)^2/(1+(T1/Ts)) ^2
% c0 = K1/(1+(T1/Ts))^2
% Yk = a0*Yk_1+b0*Yk_2+c0*Rk % 差分方程
%% 对象动态特性实验
clear,clc;
% 通过赋值将符号计算转化为代数计算
K1=1.0;T1=10.0;R0=1;
% 绘制matlab阶跃响应
G1=tf([K1],[T1^2 2*T1 1])
hold on
step(G1) % 连续传递函数阶跃响应曲线 T趋于0（R0=1）
step(2*G1) % 连续传递函数阶跃响应曲线 T趋于0（R0=2）
% G21=c2d(G1,0.5,'zoh') % 零阶保持器离散化 一个脉冲序列变成一个连续的阶梯信号
% G21=c2d(G1,0.5,'foh') % 一阶保持器离散化 讯号重建为分段线性近似原始讯号的函数
G21=c2d(G1,0.5,'foh') % 一阶保持器离散化
step(G21) % 离散传递函数阶跃响应曲线 % T=0.5 R0=1
step(2*G21) % 离散传递函数阶跃响应曲线 % T=0.5 R0=2
G22=c2d(G1,2,'foh') % 一阶保持器离散化
step(G22) % 离散传递函数阶跃响应曲线 % T=2 R0=1
step(2*G22) % 离散传递函数阶跃响应曲线 % T=2 R0=2
G23=c2d(G1,5,'foh') % 一阶保持器离散化
step(G23) % 离散传递函数阶跃响应曲线 % T=5 R0=1
step(2*G23) % 离散传递函数阶跃响应曲线 % T=5 R0=2
hold off
legend('连续传递函数阶跃响应曲线 T趋于0 (R0=1)','连续传递函数阶跃响应曲线 T趋于0 (R0=2)', ...
'离散传递函数阶跃响应曲线(T=0.5s R0=1)','离散传递函数阶跃响应曲线(T=0.5s R0=2)', ...
'离散传递函数阶跃响应曲线(T=2s R0=1)','离散传递函数阶跃响应曲线(T=2s R0=2)', ...
'离散传递函数阶跃响应曲线(T=5s R0=1)','离散传递函数阶跃响应曲线(T=5s R0=2)');
title('对象动态特性曲线');
xlabel('t/s');
ylabel('响应Y');
%% 单回路PID控制实验
clear,clc;
% 可调整参数
L = 1; % 积分分离逻辑系数
T = 2; % 周期
Kp = 5;Ti = 14;Td = 0; % PID参数
K1=1.0;T1=10.0;R0=1;
% PID控制器参数
a = Kp*(1+(T/Ti)*L+(Td/T))
b = -Kp*(1+2*(Td/T))
c = Kp*(Td/T)
% 差分方程参数
a0 = 2*(T1/T)*(1+(T1/T))/(1+(T1/T)) ^2
b0 = -(T1/T)^2/(1+(T1/T)) ^2
c0 = K1/(1+(T1/T))^2
G1=tf([K1],[T1^2 2*T1 1])
G2=c2d(G1,T,'foh')
[num,den]=tfdata(G2,'v')
e=zeros(1,1000); %前一时刻的偏差
Ee=zeros(1,1000); %累积偏差
u=zeros(1,1000); %前一时刻的控制量
y=zeros(1,1000); %前一时刻的输出
r=ones(1,1000); %期望值
u=zeros(1,1000);%预先分配内存
time=zeros(1,1000);%时刻点（设定1000个）
for k=1:1:1000
% time(k)=(k)*T; %时间参数
if k==1
e(k)=r(1)-y(1);
else
e(k)=r(k)-y(k-1); %误差信号
end
if k==1
delta_u(k)=a*e(k);
elseif k==2
delta_u(k)=a*e(k)+b*e(k-1);
else
delta_u(k)=a*e(k)+b*e(k-1)+c*e(k-2); %系统PID控制器输出序列
end
if k==1
Ee(k)=e(k);
else
Ee(k)=Ee(k-1)+e(k); %误差的累加和
end
if k==1
u(k)=delta_u(k);
else
u(k)=u(k-1)+delta_u(k); %前一个的控制器输出值
end
if k==1
y(k)=c0*u(k);
elseif k==2
y(k)=a0*y(k-1)+c0*u(k);
else
y(k)=a0*y(k-1)+b0*y(k-2)+c0*u(k);%系统响应输出序列
end
end
hold on
plot(y);plot(r);
legend('阶跃响应曲线','期望值');
title('PID控制阶跃响应曲线');
xlabel('t/s');
ylabel('响应Y');
