主要包含根轨迹+频率分析+状态空间+数字控制器+PID。

拿到物理模型→系统分析→列出物理运动方程式→拉普拉斯变换→传递函数+状态空间方程→引入控制器→系统输出响应波形即可达到期望要求

一.物理建模得被控对象

1.传递函数

% 传递函数
s = tf('s');
P_pitch = (1.151*s + 0.1774)/(s^3 + 0.739*s^2 + 0.921*s);

2.状态空间方程

做法一：直接构建 A、B、C、D 矩阵

% 状态空间方程
A = [-0.313, 56.7, 0;
     -0.0139, -0.426, 0;
      0, 56.7, 0];
B = [0.232;
     0.0203;
     0];
C = [0, 0, 1];
D = [0];
pitch_ss = ss(A, B, C, D);

做法二：由传递函数得到状态方程

(1) 自己写出来传递函数分子、分母

% 自己写传递函数分子分母
num = [1.151, 0.1774];
den = [1, 0.739, 0.921, 0];
[A, B, C, D] = tf2ss(num, den); % num 是开环传递函数分子，den 是分母
pitch_ss_1 = ss(A, B, C, D);

(2) 直接由传递函数转换

pitch_ss_2 = ss(P_pitch); % 其中 P_pitch 是开环传递函数

注意：此命令得到状态方程的 A、B、C、D 并不会直接保存在工作区

二.系统分析

1.设计要求

• 超调量

• 上升时间

• 稳定时间

• 稳态误差

2.在阶跃响应下开环传递函数输出波形

做法一：带 figure

% 看看开环传递函数输出响应波形
t = 0:0.01:10;     % 时间从 0→10s，中间每间隔 0.01 就记录一下数值大小，即采样
figure(1);         % 告诉软件，我要画一个图，先把图方框准备好
step(P_pitch, t);  % 阶跃响应 step 下，开环传递函数的输出响应波形
grid on;           % 显示网格
xlabel('时间');    % 横坐标的名字
ylabel('飞机的俯仰角'); % 纵坐标的名字
title('开环阶跃响应曲线');  % 标题

做法二：直接输出波形图

% 直接输出波形图
step(P_pitch), grid;  % 不加分号直接显示输出
axis([0 10 0 5]);     % 限制坐标范围。横坐标范围 0-10，纵坐标范围 0-5
xlabel('时间');
ylabel('飞机的俯仰角');
title('开环阶跃响应曲线');

做法三：线性系统分析器

% 启动线性系统分析器，输入波形是 step，研究对象是开环传递函数，时间为 0→10
linearSystemAnalyzer('step', P_pitch, 0:0.01:10);
% 可在分析器里看到输出阶跃响应波形

注意：线性系统分析器功能强大，例如

• 可以更改输入响应。右键→线性仿真→设计信号（或者在命令语句将 step 改为正弦波等等）

• 右键可以更改想看的波形，可以改看零极点分布图、Bode 图、根轨迹图等等

% 可以同时输入多个传递函数，观察波形分析
rP_pitch = 0.1/(0.5*s + 1); % 随便设的
linearSystemAnalyzer('step', rP_pitch, P_pitch, 0:0.1:10);

3.看看开环传递函数零极点

做法一：pole、zero 指令

fprintf('开环传递函数极点为');
pole(P_pitch)    % 开环极点，注意没有分号
zero(P_pitch)    % 开环零点

做法二：pzmap 指令

% 同时看零极点
pzmap(P_pitch);
axis([-1 1 -1 1]);

4.得到闭环传递函数

fprintf('闭环传递函数为');
sys_cl = feedback(P_pitch, 1);

5.对输出Y(s) 进行因式分解以求得拉普拉斯反变换

（若无特殊要求，此步省略）

输出Y(s) = \text{传递函数} G(s) \times \text{输入} R(s)。此处输入R(s) = \dfrac{1}{s}。

R = 1/s;
Y = zpk(sys_cl * R);

①即Y(s) = \dfrac{1.151 (s + 0.1541)}{s (s + 0.08805) (s^2 + 0.6509s + 2.015)}

②当前目标是因式分解成如下形式

③现在求A、 B、 C、 D（也可以手动计算），若使用 MATLAB 计算，相关命令为 residue

[r, p, k] = residue([1.151, 0.1774], [1, 0.739, 2.072, 0.1774, 0]);

得到

④因此A = 1， B = -0.0881，现在求C、 D

[num, den] = residue(r(1:2), p(1:2), k);
tf(num, den);

⑤可得

⑥因此C = -0.5605， D = -0.4036， A、 B、 C、 D 全部得到，可以进行拉普拉斯反变换求解

syms s
A = 1;                  
B = -0.0881;
C = -0.5605;
D = -0.4036;
F = A/s + B/(s + 0.08805) + (C*s + D)/(s^2 + 0.6509*s + 2.015);
ilaplace(F);

⑦这样y(t) 就得到了

三.控制方法与分析

1.根轨迹法

先找补偿器增益K_p，看看单纯用K_p 能不能满足要求

（1）做法一：

% 打开控制系统设计器
controlSystemDesigner('rlocus', P_pitch);
% 命令中，VIEWS 必须为以下一个或多个项: 'rlocus'、'bode'、'nichols' 和 'filter'

注意：控制系统设计器功能强大，例如

• 可同时查看根轨迹与输出响应，并且根轨迹的零极点可拖动，移动根轨迹上的零极点从而改变K，使得输出响应实时变化

• 针对不同的回路结构，可以在编辑架构里选择不同的回路

• 点击控制器G，可实时查看当前回路增益K_p 与开环传递函数

• 根轨迹中，右键→设计需求→新建，可添加要求。其中上升时间与固有频率有响应转换关系，要求添加完毕后，非阴影区（白色区）即是期望区域

• 输出响应也可以右键添加指标，显示当前的稳定值、上升时间等

（2）做法二：

% 绘制根轨迹图
rlocus(P_pitch);
axis([-0.6 0 -0.6 0.6]); % 坐标范围
sgrid(0.6, 0.9);         % 第一个是阻尼比，第二个是固有频率（有公式算出来的）

根轨迹图有了，接下来找K_p。找K_p 就必须先在根轨迹上选一个点。

[Kp, pole] = rlocfind(P_pitch);  % 选择极点后对应的 K 值

命令执行后会让你在根轨迹上选一个点，选择此点后会显示①此点对应的K_p 值 + ②离这个点最近的、实轴上的极点坐标

以上两个方法都得到了K_p 值，把这个K_p 值乘进开环传递函数看看输出响应波形，如果单纯的K_p 无能为力，则还有超前/滞后补偿器可以抢救。

超前补偿器

适用条件：当前的根轨迹与所要求的非阴影区（期望区域、要求区域）不沾边，具体来说当前根轨迹范围完全在期望区域的右侧，就不得不向左移动根轨迹。（换句话说，当我们希望闭环的根往左移动（更稳定时），引入超前补偿即可）

传递函数：

（1）做法一：

继续在控制系统设计器中调节

• 点击右上角预设项→选项→零点/极点/增益→确定

• 右击空白处→编辑器补偿器→右击→添加零极点→超前→修改实极点、实零点的值（选择有技巧）（此例的零点选择为-0.9，极点选择为-3）

• 随后新根轨迹图就有了，拖动极点以达到期望的输出响应波形

（2）做法二：

% 带超前补偿器的根轨迹图
z0 = -0.9;
p0 = -3;
s = tf('s');
C_lead = (s - z0)/(s - p0);  % 构建超前补偿器
rlocus(C_lead * P_pitch);     % 看根轨迹图
axis([-0.6 0 -0.6 0.6]);
sgrid(0.6, 0.36);
% 重新选 Kp
[Kp, pole] = rlocfind(C_lead * P_pitch);

再用此K_p 去求输出波形

Kp = 200;
sys_cl = feedback(Kp * C_lead * P_pitch, 1);
t = 0:0.1:20;
step(sys_cl, t);

（3）做法三：

可以尝试不同类型的动态补偿器（添加零点和极点）

滞后补偿器

适用条件：响应速度都挺好，就是稳态误差太大，稳定值跟期望值差老远。因此滞后补偿器可以在减小稳态误差的同时，仍能满足暂态要求（不会过大的改变系统原有的特性）。

传递函数：

完全同超前补偿器，唯独不同之处在于z 值与p 值的选取，超前补偿器|z| < |p|，滞后补偿器|z| > |p|

例如 超前补偿器：

滞后补偿器：

其他步骤一模一样。

特殊情况：

在某些特殊情况下，例如什么补偿器都加了，共轭极点已经拉到期望范围的边缘了，但是输出波形还是不满足要求，有可能是因为实轴上的、靠近虚轴的第三极点作为主导极点，不仅减缓了系统响应，而且压制了共轭极点的影响，这个情况下，可以将共轭极点拉到期望范围外，依然会满足要求。

2.频率分析法

先看开环 Bode 图

（1）做法一：

% 在控制系统设计器里看开环传函 Bode 图
controlSystemDesigner('bode', P_pitch);

（2）做法二：

% 绘制 Bode 图
bode(P_pitch), grid;

（3）做法三：

% 查看增益裕度与相位裕度
margin(P_pitch), grid;
% 若增益裕度与相位裕度都为正，则系统稳定

设计比例控制器K_p

如果当前 Bode 图与输出波形不符合要求，则有以下两种方法找到增益K_p

（1）做法一：比较适用于一阶系统

①稳态误差公式： \text{error} = \dfrac{1}{1 + M_{w \rightarrow 0}} \times 100

②在 Bode 图上找到w \rightarrow 0 时对应的幅值 dB（例如-34 dB），带入公式计算得到稳态误差，同时也可以再看看闭环输出响应来验证一下（$K_p$ 取1 即可）

sys_cl = feedback(Kp * P_pitch, 1);
step(sys_cl);

③减小稳态误差，即抬高稳态值，例如要求稳态误差小于2\%，带入公式可以计算得出期望幅值=33.8 dB，因此从原本的-34 dB → 期望的33.8 dB，一次性抬高67.8 dB。再根据幅值与比例增益K_p 之间的关系，可得K_p = 2455。

④随后便可以将K_p = 2455 乘入传递函数，观察开环传函 Bode 图以及闭环输出波形

⑤特殊情况，由于这种做法K_p 往往取的很大，在某些特殊控制情况下（例如汽车巡航控制，输出速度响应波形变化太快，现实里根本无法实现）（例如K_p 过大导致系统有振荡，相位裕度太小，系统不稳定等），就不得不降低K_p，可以引入滞后补偿器。

（2）做法二：比较适用于多阶系统

①直接考虑相位裕度，一般为了得到满意的性能，相位裕度一般为30^\circ~ 60^\circ，这里直接拉满令P_m = 60^\circ。 180^\circ - \text{相位}120^\circ = \text{相位裕度}60^\circ，因此直接找到相位120^\circ 所对应的点

②在原始开环 Bode 图中找相位120^\circ 所对应的频率（例如=9.17 rad/s），将此频率作为新的截止频率

③查看频率为9.17 所对应的mag（幅值图的横坐标）与相角phase

[mag, phase, w] = bode(P_pitch, 9.17);

④得到数据（举例）$mag = 0.0157$，再将mag 换算成幅值

20 * log10(0.0157);

⑤得到幅值= -36.0820 dB，为了使当前选中频率9.17 rad/s 为截止频率，就必须将此幅值-36.082 dB →0 dB，因此需要抬高36.082 dB，换算成K_p 值 ≈63。

⑥利用此K_p 值再去看新 Bode 图（ P_m \approx 60^\circ）与输出阶跃响应。如果输出波形不满意，引入补偿器

超前补偿器

特别注意：使用补偿器的时候要从系统的开环传递函数从头设计，以上得到的K_p 值全部扔掉。

使用条件： 加快截止频率且增加相位裕度。它为系统增加了正相位，这样可以增加相位裕度，从而增加阻尼。与此同时，它还可以在高频率下增加开环频率响应的幅度，进而增强截止频率和整体速度（稳定时间降低）。

传递函数：

步骤：

①先只引入一个K，先初步决定K = 10，看看根轨迹 + 输出波形

K = 10;
margin(K * P_pitch), grid;
figure;
sys_cl = feedback(K * P_pitch, 1);
step(sys_cl), grid;

②在本例中，增加K = 10 会使得输出有振荡，这是因为相位裕度变小了（本例P_m = 10.4^\circ），不过没关系，超前补偿器可以增大阻尼、消除振荡。

③现在选择\alpha， \alpha 需要由最大相位确定。

\alpha 与最大相位的关系：

④现在找最大相位，超调量 → 阻尼比 → 相位裕度三者也有关系。举例超调量小于10\% → 阻尼比\zeta 就要大于0.6 → 相位裕度要大于60^\circ。当前的相位裕度=10.4^\circ，要提到60^\circ，再来个50^\circ 就够了。再由于超前补偿器的特性，要额外多加几度，最终选定添加55^\circ 的相位超前。因此最大相位\phi_m = 55^\circ。

⑤根据公式可计算得到\alpha \approx 0.1

⑥根据下面关系式计算出超前补偿器所造成的幅度增加量

⑦从（ K = 10 的）Bode 图可得知，截止频率为3.49 rad/s，幅值= -10 dB 的频率为6.1 rad/s，因此引入\alpha 后新截止频率=6.1 rad/s，利用这个新截止频率来计算T

⑧至此超前补偿器的三个参数全部得到K = 10， \alpha = 0.1， T = 0.52。看看波形图。

% 重新将开环传递函数写一遍
s = tf('s');
P_pitch = (1.151*s + 0.1774)/(s^3 + 0.739*s^2 + 0.921*s);
 
% 引入超前补偿器
K = 10;
alpha = 0.1;
T = 0.52;
C_lead = K * (T*s + 1)/(alpha*T*s + 1);
margin(C_lead * P_pitch), grid;
 
% 闭环输出响应
sys_cl = feedback(C_lead * P_pitch, 1);
step(sys_cl), grid;
title('加入了 K、α、T 后的拥有超前补偿器的输出阶跃响应图');

⑨输出波形图很难直观看到超调量等数据，采用 stepinfo 命令，可看到全部的性能指标值。

stepinfo(sys_cl);

⑩如果某些指标还是差一些，可以微调K、 \alpha 的值，例如本例超调量依然略大，取\alpha = 0.0718（添加60^\circ 的相位超前）（ T = 0.5654），就可以解决问题。

滞后补偿器

与根轨迹里的滞后补偿器完全相同

适用条件：滞后补偿器可以在保持带宽频率不变的情况下增加低频段的幅值增益，这样就可以抬高稳定值，并且带宽频率不变

传递函数：

零极点的选择：通常来说，与根轨迹法一样，滞后补偿器的零点要与系统的极点放在一起；至于极点选个靠近虚轴的就可以，例如p = 0.01 即可

% 找到回路增益 Kp，滞后补偿器 C_lag 后，命令如下（步骤同根轨迹）
bode(Kp * C_lag * P_pitch), grid;
% 在控制系统设计器中也可以实现（步骤同根轨迹）
controlSystemDesigner('bode', P_pitch);

3.状态空间法

观察能控性（全状态反馈控制器）

先给出能控性矩阵M = [B\ AB\ A^2B\ \ldots\ A^{n-1}B]

A = [-0.313, 56.7, 0;
     -0.0139, -0.426, 0;
      0, 56.7, 0];
B = [0.232;
     0.0203;
     0];
C = [0, 0, 1];
D = [0];
 
sys_order = order(pitch_ss);   % 找到最小的 MIMO 系统
 
mo = ctrb(A, B);      % 构建能控性矩阵 M = [B AB ...]
sys_rank = rank(mo);  % 找到矩阵的秩

若 sys_order == sys_rank，则矩阵满秩，则能控

观察能观性（状态观测器）

no = obsv(A, C);   % 构建能观性矩阵 N = [C; CA; ...]'
sys_rank = rank(no);

设计全状态反馈控制器

经过对原状态空间进行全状态反馈的变换，可以得新状态空间表达式

（1）做法一：极点放置法

适用于标准的一阶二阶系统（不适用于本例）

①从\det(sI - (A - BK)) 的行列式（矩阵的根）中找到闭环极点，根据阶数可获得几个可任意放置的极点（一阶可放一个，二阶可放两个）

②先找到极点要放的位置，例如p_1 = -5 + i， p_2 = -5 - i（举例），调用命令得到K 值

% 开环传递函数（举例）
J = 0.01; % 转子的转动惯量
b = 0.1;  % 电机粘性摩擦常数
K = 0.01; % 电动势常数、电机转矩常数
R = 1;    % 电阻
L = 0.5;  % 电感
 
% 状态方程
A = [-b/J,    K/J;
      K/L, -R/L];
B = [0;
     1/L];
C = [1, 0];
D = 0;
motor_ss = ss(A, B, C, D);
 
% 放极点，求 Kc
p1 = -5 + 1i;
p2 = -5 - 1i;
Kc = place(A, B, [p1, p2]);

③K_c = [13.0100, -1.0000]，将K 带回新状态方程，看输出响应波形

motor_ss = ss(A, B, C, D);
 
motor_ss_new = ss(A - B * Kc, B, C, D);
t = 0:0.01:3;
step(motor_ss_new, t);
title('状态控制器的输出响应');

④通常来说，放置极点不可能一次成功，一般都需要引入预补偿器N_{\text{bar}}，目的是缩放输入，抬高稳定值。

Nbar = rscale(motor_ss, Kc); % 利用原状态方程算 Nbar
 
t = 0:0.01:10;
step(motor_ss_new * Nbar, t); % 将新状态方程乘以 Nbar
grid;
title('状态控制器 + 预补偿器的输出响应');

注意：rscale 函数（查找比例因子以消除稳态误差函数）内容可以在 MATLAB 官方帮助中找到，复制粘贴到 MATLAB 当前目录中

（2）做法二：使用 LQR（线性二次调节器）来生成最佳增益矩阵K（适用本例）

①若要使用 LQR，需要定义两个参数：状态成本加权矩阵Q 和控制权重矩阵R，为图简单，选择R = 1，$Q = p \cdot C^T \cdot C$

②LQR 参数全了，接下来得到反馈控制矩阵K，先令加权因子p = 2

% 状态空间方程
A = [-0.313, 56.7, 0;
     -0.0139, -0.426, 0;
      0, 56.7, 0];
B = [0.232;
     0.0203;
     0];
C = [0, 0, 1];
D = [0];
pitch_ss = ss(A, B, C, D);
 
% 设计 LQR
p = 2;
Q = p * C' * C;
R = 1;
[K] = lqr(A, B, Q, R);

③可得K = [-0.5034, 52.8645, 1.4142]，根据此K 得到新状态方程，且查看输出响应

% 新状态方程
sys_cl = ss(A - B * K, B, C, D);
 
% 查看输出响应
step(sys_cl), grid;
title('LQR 下的闭环阶跃响应');

④如果输出波形响应太慢，则增大 LQR 里的加权因子p 值（本例将加权因子增大到p = 50），这样可以得到更大的K，加快系统响应，但是与此同时也会降低稳态值。如果输出波形稳态值过小，稳态误差过大，则再调节 LQR 也没救了，需要再同样引入预补偿器N_{\text{bar}}

⑤同样使用 rscale（查找比例因子以消除稳态误差函数）来得到N_{\text{bar}}

Nbar = rscale(A, B, C, D, K); % 在 LQR 的基础上继续做，因此 K 为之前由 LQR 得到的 K

⑥得到N_{\text{bar}} 看看输出波形

% 带 LQR + Nbar 的新状态空间方程
sys_cl = ss(A - B * K, B * Nbar, C, D); % 也可以同做法一，先求得 sys_cl_new，再乘以 Nbar，是一样的
 
step(sys_cl), grid;
title('LQR + Nbar 下的闭环阶跃响应');

预补偿器N_{\text{bar}}

由于上述两种做法都可能用到N_{\text{bar}}，因此做一些补充

• 添加的预补偿器N_{\text{bar}} 并不位于反馈回路，这就意味着如果有未知干扰，预补偿器对此部分的偏差无法消除

• 换句话说，预补偿器的作用就是将输入放大或者缩小，导致了输出稳态值的提升，但是稳态误差说到底还是没真正被消除。根本原因就是N_{\text{bar}} 不在反馈回路，它不能真正改变输出的波形，也就是说不会真正减小稳态误差

• 有一种方法是将预补偿器与积分控制相结合，既可以缩放输出波形，又可以一定程度减小稳态误差

4.PID

PID 传递函数：

% 打开控制系统设计器
s = tf('s');
P_pitch = (1.151*s + 0.1774)/(s^3 + 0.739*s^2 + 0.921*s);
 
controlSystemDesigner(P_pitch);

左上角→调节方法→PID 调节→选择鲁棒响应时间/经典设计公式

• 鲁棒响应时间就是自动谐调 PID 参数，可以为所有类型的 PID 调节参数，可用于设计稳定、不稳定、积分的系统

• 经典设计公式算法需要稳定的或带有积分的设备，无法调整导数滤波器，但是会有很多的公式方法

此例选择鲁棒响应时间算法→设计模式时域表示→参数调调调→设计成功

5.数字控制器

（1）做法一：状态方程 + LQR

①先得到离散状态空间方程，使用 c2d 命令将连续时间模型→离散时间模型，c2d 需要三个参数：系统模型 + 采样时间 + 保持电路的类型（通常用零阶保持器 zoh）

②采样时间：采样频率至少比闭环 Bode 图的带宽频率（闭环截止频率）（不是开环截止频率）大 30 倍。根据闭环波特图可得，带宽频率为 2 rad/sec，因此得到 0.01s 的采样时间

③使用 c2d 命令，可以得到原离散状态空间模型

A = [-0.313, 56.7, 0;
     -0.0139, -0.426, 0;
      0, 56.7, 0];
B = [0.232;
     0.0203;
     0];
C = [0, 0, 1];
D = [0];
sys_ss = ss(A, B, C, D);
Ts = 0.01;
sys_d = c2d(sys_ss, Ts, 'zoh');

④验证能控性

co = ctrb(sys_d);
Controllability = rank(co);

⑤离散全状态反馈控制系统，得到新离散状态空间方程

⑥依旧利用 LQR，不过是离散版本

A = sys_d.A;
B = sys_d.B;
C = sys_d.C;
D = sys_d.D;
p = 50;
Q = p * C' * C;
R = 1;
[K] = dlqr(A, B, Q, R);

⑦得到K，就可以得到离散状态空间方程了，且看看输出波形

time = 0:Ts:10;
theta_des = ones(size(time));
sys_cl = ss(A - B * K, B, C, D, Ts);
[y, t] = lsim(sys_cl, theta_des, time);
% 画图
stairs(t, y);
title('DLQR 下的闭环阶跃响应');

⑧可能仍需要添加预补偿器N_{\text{bar}}，但是离散版本不可以用 rscale 自动求了，需要手动调调改改，最后选定N_{\text{bar}} = 6.95。重新得离散状态空间方程 + 输出波形

Nbar = 6.95;
sys_cl = ss(A - B * K, B * Nbar, C, D, Ts); % 多乘了 Nbar
[y, t] = lsim(sys_cl, theta_des, time);
% 画图
stairs(t, y);
title('DLQR 下的闭环阶跃响应');

（2）做法二：传递函数 + 根轨迹

①将连续传函转换成离散传函

Ts = 1/50;
dp_pitch = c2d(P_pitch, Ts, 'zoh');
zpk(dp_pitch); % 变成因式分解形式，更容易看零极点

②利用 zgrid 命令找到离散的可接受符合要求的根轨迹区域，然后在该区域里找增益K。而使用 zgrid 命令需要两个参数：自然频率W_n + 阻尼比\zeta。

假设上升时间为 5 秒，超调量为10\%，因此固有频率大于 0.36，阻尼比大于 0.6。但由于 zgrid 命令的固有频率单位是 rad/sample，因此换算成 0.0072 rad/sample

Wn = 0.0072;
zeta = 0.6;
 
% 画根轨迹图
rlocus(dp_pitch);
zgrid(zeta, Wn);
axis([0.97 1 -0.05 0.05]);

后面步骤与代码同根轨迹法

③选极点 → 得K_p 值 → 画闭环输出波形 → 不满意 → 乘上超前/滞后补偿器 → 重新选极点 → 得新K_p 值 →K_p × 超前/滞后补偿器 × 原离散开环传递函数 → 得闭环传递函数 → 画闭环输出波形

（3）做法三：传递函数 + PID + 根轨迹

①将连续传函转换成离散传函

Ts = 1/50;
dp_pitch = c2d(P_pitch, Ts, 'zoh');
zpk(dp_pitch);

②把连续 PID 转换成离散 PID

% 离散 PID
Kp = 100;
Ki = 200;
Kd = 10;
 
C = Kp + Ki/s + Kd*s;
dC = c2d(C, Ts, 'tustin');

③用离散 PID 去控制离散传递函数，看看输出响应波形

sys_cl = feedback(dC * dp_pitch, 1);
[x2, t] = step(sys_cl, 12);  % 12 的意思是时间范围
% 画图
stairs(t, x2);

如果发现输出响应不稳定，但该 PID 参数在连续域可以得到很好的调节结果（在离散域就调节不好），需要进根轨迹找原因

④看看离散系统 + 离散 PID 的根轨迹

% 不稳定，看看根轨迹
rlocus(dC * dp_pitch);
axis([-1.5 1.5 -1 1]);
title('离散系统带 PID 的根轨迹');

⑤经过很复杂的、我不会的分析后发现，PID 控制器有个我们不想要的极点，因此就需要改变 PID 控制器。因此给离散 PID 增加了一个-0.82 的极点（随便举例）

% 给离散 PID 再增加一个极点
z = tf('z', Ts);
dC_new = dC / (z + 0.82);
% 看看新根轨迹
rlocus(dC_new * dp_pitch);

⑥在新根轨迹上找个极点，看看对应的K 值。此时得到K 值

⑦在此K 值下 + 新离散 PID + 原离散开环传递函数的输出波形如何？

sys_cl = feedback(K * dC_new * dp_pitch, 1); % 其中 K 值为上一步选择极点得到的值
[x3, t] = step(sys_cl, 8);   % 8 的意思是时间范围
 
% 画输出波形图
stairs(t, x3);
title('离散域下开传函 + PID + 根轨迹');

⑧如果不行就重新选K，一直达到设计要求