控制具有长死区时间的过程：史密斯预估器

R2026b
打开脚本

此示例展示了 PI 控制对于具有长死区时间的过程的局限性，并说明了一种称为“史密斯预估器”的控制策略的优势。

此示例受以下文献启发：

A. Ingimundarson and T. Hagglund, "Robust Tuning Procedures of
Dead-Time Compensating Controllers," Control Engineering Practice,
9, 2001, pp. 1195-1208.

过程模型

该类过程的开环响应被建模为一阶加死区时间，其时间常数为 40.2 秒，时滞为 93.9 秒：

s = tf('s');
P = exp(-93.9*s) * 5.6/(40.2*s+1);
P.InputName = 'u';
P.OutputName = 'y';
P

P = 

  From input "u" to output "y":
                    5.6
  exp(-93.9*s) * ----------
                 40.2 s + 1
 
Continuous-time transfer function.

请注意，该时滞超出了时间常数的两倍。此模型在许多化工过程中都具有代表性。其阶跃响应如下所示。

step(P), grid on

PI 控制器

比例积分 (PI) 控制是过程控制中常用的一种方法。相应的控制架构如下所示。

补偿器 C 是一个标准形式的 PI 控制器，包含两个调节参数：比例增益 Kp 和积分时间 Ti。我们使用 PIDTUNE 命令来设计一个开环带宽为 0.006 弧度/秒的 PI 控制器：

Cpi = pidtune(P,pidstd(1,1),0.006);
Cpi

Cpi = 

             1      1 
  Kp * (1 + ---- * ---)
             Ti     s 

  with Kp = 0.0501, Ti = 47.3
 
Continuous-time PI controller in standard form

要评估 PI 控制器的性能，请关闭反馈环，并仿真对参考信号 ysp 和输出扰动信号 d 的阶跃变化的响应。由于反馈路径中存在延迟，因此需要使用 SS 命令将 P 或 Cpi 转换为状态空间表示：

Tpi = feedback([P*Cpi,1],1,1,1);  % closed-loop model [ysp;d]->y
Tpi.InputName = {'ysp' 'd'};

step(Tpi), grid on

闭环响应的超调在可接受范围内，但略显迟缓（约 600 秒后稳定）。增大比例增益 Kp 虽然可以加快响应速度，但也会显著增加超调，并很快导致不稳定：

Kp3 = [0.06;0.08;0.1];      % try three increasing values of Kp
Ti3 = repmat(Cpi.Ti,3,1);   % Ti remains the same
C3 = pidstd(Kp3,Ti3);       % corresponding three PI controllers
T3 = feedback(P*C3,1);
T3.InputName = 'ysp';

step(T3)
title('Loss of stability when increasing Kp')

长死区时间会使 PI 控制器的性能受到严重限制。这是因为 PI 控制器不知道死区时间的存在，当实际输出 y 与期望设定点 ysp 不匹配时，它的反应过于“急躁”。每个人在淋浴时都经历过类似的现象，水温调节存在较长的滞后时间。若因过于急躁而频繁调整，往往就会出现水温忽烫忽冷交替变化的情况。更明智的做法是先等待温度设定变化生效，再进行进一步调整。而且，一旦我们知道了旋钮调到什么位置能提供心仪的水温，便只需等待淋浴系统做出反应所需的时间，然后即可获得合适的水温。这种“最优”控制策略正是史密斯预估器方案背后的基本思想。

史密斯预估器

史密斯预估器的控制结构如下所示。

史密斯预估器使用内部模型 Gp 来预测过程的无延迟响应 yp（例如，给定旋钮设置将提供什么样的水温）。然后，它将该预测值 yp 与期望设定点 ysp 进行比较，以确定所需的调整（控制量 u）。为了防止系统漂移并抑制外部扰动，史密斯预估器还会将实际过程输出与考虑了死区时间的预测值 y1 进行比较。差值 dy=y-y1 通过滤波器 F 反馈，并构成总误差信号 e 的一部分。请注意，dy 相当于在等待足够长的时间让淋浴系统做出反应后所感知到的温度偏差。

部署史密斯预估器方案需要

  • 过程动态特性的模型 Gp 以及过程死区时间的估计值 tau

  • 补偿器和滤波器动态特性（C 和 F）的适当设置

基于过程模型，我们使用：

$$G_p(s) = {5.6 \over 1 + 40.2 s } , \tau = 93.9 $$

对于 F，使用时间常数为 20 秒的一阶滤波器来捕获低频扰动。

F = 1/(20*s+1);
F.InputName = 'dy';
F.OutputName = 'dp';

对于 C，我们重新设计 PI 控制器和 PI 控制器看到的整体被控对象，包括从 P、Gp、F 和死区时间获得的动态特性。借助史密斯预估器控制结构，我们能够增加开环带宽以实现更快的响应，同时增加相位裕度以减小超调。

% Process
P = exp(-93.9*s) * 5.6/(40.2*s+1);
P.InputName = 'u';
P.OutputName = 'y0';

% Prediction model
Gp = 5.6/(40.2*s+1);
Gp.InputName = 'u';
Gp.OutputName = 'yp';

Dp = exp(-93.9*s);
Dp.InputName = 'yp'; Dp.OutputName = 'y1';

% Overall plant
S1 = sumblk('ym = yp + dp');
S2 = sumblk('dy = y0 - y1');
Plant = connect(P,Gp,Dp,F,S1,S2,'u','ym');

% Design PI controller with 0.08 rad/s bandwidth and 90 degrees phase margin
Options = pidtuneOptions('PhaseMargin',90);
C = pidtune(Plant,pidstd(1,1),0.08,Options);
C.InputName = 'e';
C.OutputName = 'u';
C

C = 

             1      1 
  Kp * (1 + ---- * ---)
             Ti     s 

  with Kp = 0.574, Ti = 40.2
 
Continuous-time PI controller in standard form

PI 控制器与史密斯预估器的比较

为了比较两种设计的性能，首先推导史密斯预估器结构中从 ysp,d 到 y 的闭环传递函数。为了方便连接所有相关模块，请为这些模块的所有输入和输出通道命名，并通过 CONNECT 完成连线：

% Assemble closed-loop model from [y_sp,d] to y
Sum1 = sumblk('e = ysp - yp - dp');
Sum2 = sumblk('y = y0 + d');
Sum3 = sumblk('dy = y - y1');
T = connect(P,Gp,Dp,C,F,Sum1,Sum2,Sum3,{'ysp','d'},'y');

使用 STEP 比较史密斯预估器（蓝色）与 PI 控制器（红色）：

step(T,'b',Tpi,'r--')
grid on
legend('Smith Predictor','PI Controller')

史密斯预估器的响应速度要快得多，且无超调。通过绘制从 ysp 到 y 的闭环波特响应，也可以在频域中看出这种差异。请注意，史密斯预估器的带宽更高。

bode(T(1,1),'b',Tpi(1,1),'r--',{1e-3,1})
grid on
legend('Smith Predictor','PI Controller')

模型不匹配情况下的稳健性

在前面的分析中，内部模型

$$ G_p(s) e^{-\tau s} $$

与过程模型 P 完全匹配。在实际应用中，内部模型只是真实过程动态特性的近似模型，因此，了解史密斯预估器对于过程动态特性及死区时间不确定性的稳健性至关重要。

假设有两个能够代表过程参数不确定性范围的扰动被控对象模型：

P1 = exp(-90*s) * 5/(38*s+1);
P2 = exp(-100*s) * 6/(42*s+1);

bode(P,P1,P2), grid on
title('Nominal and Perturbed Process Models')

要分析稳健性，请将标称模型和扰动模型收集到一个过程模型数组中，重新构建 PI 控制器和史密斯预估器设计的闭环传递函数，并对闭环响应进行仿真：

Plants = stack(1,P,P1,P2);  % array of process models
T1 = connect(Plants,Gp,Dp,C,F,Sum1,Sum2,Sum3,{'ysp','d'},'y'); % Smith
Tpi = feedback([Plants*Cpi,1],1,1,1);   % PI

step(T1,'b',Tpi,'r--')
grid on
legend('Smith Predictor 1','PI Controller')

两种设计都对模型不匹配很敏感，闭环波特图证实了这一点：

bode(T1(1,1),'b',Tpi(1,1),'r--')
grid on
legend('Smith Predictor 1','PI Controller')

提高稳健性

要降低史密斯预估器对建模误差的敏感性，请检查内环和外环的稳定裕度。内环 C 的开环传递函数为 C*Gp，因此可通过以下方式获得其稳定裕度

margin(C * Gp)
title('Stability Margins for the Inner Loop (C)')

内环具有充足的增益裕度和相位裕度，因此接下来重点关注外环。使用 CONNECT 推导在内环闭合的情况下，从 ysp 到 dp 的开环传递函数 L：

Sum1o = sumblk('e = ysp - yp');  % open the loop at dp
L = connect(P,Gp,Dp,C,F,Sum1o,Sum2,Sum3,{'ysp','d'},'dp');

bodemag(L(1,1))

当过程模型和预测模型完全匹配时，该传递函数基本上应当为零。为了深入了解外环的稳定裕度，我们需要使用其中一个扰动过程模型，例如 P1：

H = connect(Plants(:,:,2),Gp,Dp,C,Sum1o,Sum2,Sum3,{'ysp','d'},'dy');
H = H(1,1);  % open-loop transfer ysp -> dy
L = F * H;

margin(L)
title('Stability Margins for the Outer Loop (F)')
grid on;
xlim([1e-2 1]);

此增益曲线在 0.04 弧度/秒附近有一个凸起，这会降低增益裕度，并增加闭环阶跃响应中的凸起。要解决此问题，请选择一个更早开始滚降且滚降更快的滤波器 F：

F = (1+10*s)/(1+100*s);
F.InputName = 'dy';
F.OutputName = 'dp';

验证增益裕度在 0.04 弧度/秒相位穿越频率附近是否有所改善：

L = F * H;
margin(L)
title('Stability Margins for the Outer Loop with Modified F')
grid on;
xlim([1e-2 1]);

最后，使用修改后的滤波器对闭环响应进行仿真：

T2 = connect(Plants,Gp,Dp,C,F,Sum1,Sum2,Sum3,{'ysp','d'},'y');

step(T2,'b',Tpi,'r--')
grid on
legend('Smith Predictor 2','PI Controller')

修改后的设计虽然会导致标称响应速度略微变慢，但可提供更稳定的性能。

改善抗扰

从 d 到 y 的闭环传递函数公式表明，F 的最优选择为

$$ F(s) = e^{\tau s} $$

其中 tau 是内部模型的死区时间。无论 P 与 Gp 之间的不匹配情况如何，此选择都能实现完美的抗扰。遗憾的是，这种“负延迟”不满足因果性，因此无法实现。在以下论文中：

  Huang, H.-P., et al., "A Modified Smith Predictor with an Approximate
  Inverse of Dead Time," AiChE Journal, 36 (1990), pp. 1025-1031

作者建议使用相位超前逼近：

$$ e^{\tau s} \approx { 1 + B(s) \over 1 + B(s) e^{-\tau s} }$$

其中 B 是一个低通滤波器，其时间常数与内部模型 Gp 相同。您可以按如下方式测试此方案：

定义 B(s) 和 F(s)。

B = 0.05/(40*s+1);
tau = totaldelay(Dp);
F = (1+B)/(1+B*exp(-tau*s));
F.InputName = 'dy';
F.OutputName = 'dp';

重新设计带宽降低的 PI 控制器。

Plant = connect(P,Gp,Dp,F,S1,S2,'u','ym');
C = pidtune(Plant,pidstd(1,1),0.02,pidtuneOptions('PhaseMargin',90));
C.InputName = 'e';
C.OutputName = 'u';
C

C = 

             1      1 
  Kp * (1 + ---- * ---)
             Ti     s 

  with Kp = 0.144, Ti = 40.2
 
Continuous-time PI controller in standard form

计算闭环模型 T3。

T3 = connect(Plants,Gp,Dp,C,F,Sum1,Sum2,Sum3,{'ysp','d'},'y');

比较 T3 与 T2 和 Tpi。

step(T2,'b',T3,'g',Tpi,'r--')
grid on
legend('Smith Predictor 2','Smith Predictor 3','PI Controller')

比较结果表明，我们的最新设计加快了抗扰响应，但代价是设定点跟踪变慢。