H无穷示例——MATLAB

四分之一汽车悬挂模型（Quarter-Car Suspension Model）

以四分之一汽车悬架为模型，将轮胎简化为弹簧(弹性系数为${\mathbf{k}}_{\mathbf{t}}$），轮胎组件简化为质量${\mathbf{m}}_{\mathbf{w}}$，汽车被动悬架简化为弹簧（弹性系数为${\mathbf{k}}_{\mathbf{s}}$）和阻尼器（阻尼系数为${\mathbf{b}}_{\mathbf{s}}$），液压执行器简化为输入力${\mathbf{f}}_{\mathbf{s}}$，车身则简化为质量${\mathbf{m}}_{\mathbf{b}}$。
​

图源：https://www.youtube.com/watch?v=kRt7H0k8A4k

简化后的模型如下(MATLAB示例截图):
​
其中，$x_b$表示车身竖直方向位移，$x_w$表示轮胎竖直方向位移，$r$则为路面干扰（崎岖不平等等)。

  根据牛顿第二定律可列下式
$$\{\begin{array}{ccc}{m}_b{a}_b& =& -k_s(x_b-x_w)-b_s({\dot{x}}_b-{\dot{x}}_w)+10^3{f}_s\\ m_w{a}_w& =& -k_t(x_w-r)+k_s(x_b-x_w)+b_s({\dot{x}}_b-{\dot{x}}_w)-10^3{f}_s\end{array}$$

1.模型的状态空间方程及对应代码

  定义状态变量为$x=(x_1,x_2,x_3,x_4{)}^T=(x_b,{\dot{x}}_b,x_w,{\dot{x}}_w{)}^T$,可得状态空间方程如下:
$$\begin{array}{ccc}{\dot{x}}_1& =& x_2\\ {\dot{x}}_2& =& \frac{-k_s(x_b-x_w)-b_s({\dot{x}}_b-{\dot{x}}_w)+10^3{f}_s}{m_b}\\ \dot{x_3}& =& x_4\\ {\dot{x}}_4& =& \frac{-k_t(x_w-r)+k_s(x_b-x_w)+b_s({\dot{x}}_b-{\dot{x}}_w)-10^3{f}_s}{m_w}\end{array}$$
  根据状态空间方程的一般形式
$$\{\begin{array}{ccc}\frac{dz(t)}{dt}& =& Az(t)+Bu(t)\\ y(t)& =& Cz(t)+Du(t)\end{array}$$
输出量$y(t)=(x_1,s_d,{\dot{x}}_2)$,其中$s_d$为悬架位移，由此可得

$\left[\begin{array}{c}{\dot{x}}_1\\ {\dot{x}}_2\\ {\dot{x}}_3\\ {\dot{x}}_4\end{array}\right]=\left[\begin{array}{cccc}0& 1& 0& 0\\ \frac{1}{m_b}[-k_s& -b_s& k_s& b_s]\\ 0& 0& 0& 1\\ \frac{1}{m_w}[k_s& b_s& -k_s-k_t& -b_s]\end{array}\right]\left[\begin{array}{c}{x}_1\\ x_2\\ x_3\\ x_4\end{array}\right]+\left[\begin{array}{cc}0& 0\\ \frac{10^3}{m_b}& 0\\ 0& 0\\ -\frac{10^3}{m_w}& \frac{k_t}{m_w}\end{array}\right]\left[\begin{array}{c}{f}_s\\ r\end{array}\right]$

$\left[\begin{array}{c}{x}_1\\ s_d\\ {\dot{x}}_2\end{array}\right]=\left[\begin{array}{cccc}1& 0& 0& 0\\ 1& 0& -1& 0\\ \frac{1}{m_b}[-k_s& -b_s& k_s& b_s]\end{array}\right]\left[\begin{array}{c}{x}_1\\ x_2\\ x_3\\ x_4\end{array}\right]+\left[\begin{array}{cc}0& 0\\ 0& 0\\ -\frac{10^3}{m_b}& 0\end{array}\right]\left[\begin{array}{c}{f}_s\\ r\end{array}\right]$

则

$A=\left[\begin{array}{cccc}0& 1& 0& 0\\ \frac{1}{m_b}[-k_s& -b_s& k_s& b_s]\\ 0& 0& 0& 1\\ \frac{1}{m_w}[k_s& b_s& -k_s-k_t& -b_s]\end{array}\right]$, $B=\left[\begin{array}{cc}0& 0\\ \frac{-10^3}{m_b}& 0\\ 0& 0\\ \frac{-10^3}{m_w}& \frac{k_t}{m_w}\end{array}\right]$,
$C=\left[\begin{array}{cccc}1& 0& 0& 0\\ 1& 0& -1& 0\\ \frac{1}{m_b}[-k_s& -b_s& k_s& b_s]\end{array}\right]$, $D=\left[\begin{array}{cc}0& 0\\ 0& 0\\ -\frac{10^3}{m_b}& 0\end{array}\right]$
  此部分转换MATLAB代码如下¹：

% 模型参数
mb = 300;	%kg
mw = 60;	%kg
bs = 1000;	%N/m/s
ks = 16000; %N/m
kt = 190000; %N/m

%状态矩阵
A = [0 1 0 0;[-ks -bs ks bs]/mb;
	 0 0 0 1;[ks   bs -ks-kt -bs]/mw];
B = [0 0; 1e3/mb 0; 0 0; -1e3/mw kt/mw]; 
C = [1 0 0 0; 1 0 -1 0; A（2,:）];
D = [0 0; 0 0; B(2,:)];

qcar = ss(A,B,C,D);
qcar.StateName = {'body travel(m)';'body vel(m/s)';'wheel travel(m)';'wheel vel(m/s)'};
qcar.InputName = {'fs';'r'};
qcar.OutputName = {'xb';'sd';'ab'};

  示例中提到’The transfer function from actuator to body travel and acceleration has an imaginary-axis zero with natural frequency 56.27 rad/s. This is called the tire-hop frequency.’

  轮胎跳动频率”（tire-hop frequency）是指汽车悬架系统中的一种固有频率，在这个频率下，轮胎相对于地面会发生垂直方向的振动或跳动。这一现象通常与悬架系统的动力学特性有关，特别是涉及到轮胎和悬架之间的弹簧和阻尼元件的耦合。

  使用函数’tzero’，它可以接受一个传递函数或状态空间模型对象，并返回该系统的传递零点。

tzero(qcar{'xb','ab'},'fs');

  以上这段代码计算了从输入量‘执行器输出力$f_s$ ’到输出量‘车身位移$x_b$’和‘车身加速度$a_b$’的传递函数矩阵的传递零点。

‘Similarly, the transfer function from actuator to suspension deflection has an imaginary-axis zero with natural frequency 22.97 rad/s. This is called the rattlespace frequency.’

zero(qcar('sd','fs'));

  使用函数’zero’计算输入量‘执行器输出力$f_s$ ’到输出量‘悬架位移$s_d$’传递函数矩阵的传递零点。
  函数’tzero’与’zero’的区别²

  ‘Road disturbances influence the motion of the car and suspension. Passenger comfort is associated with small body acceleration. The allowable suspension travel is constrained by limits on the actuator displacement. Plot the open-loop gain from road disturbance and actuator force to body acceleration and suspension displacement.’
  示例中指出乘客舒适度主要与车身加速度有关，而悬架位移主要受执行器位移的限制。分别绘制在’道路干扰r’作为输入量、'执行器输出力$f_s$’作为输入量，'车身加速度$a_b$’和’悬架位移$s_d$’作为输出量下的开环伯德图。

bodemag(qcar({'ab','sd'};'r'),'b',qcar({'ab','sd'};'fs','r',{1,100}));
legend('Road disturbance(r)','Actuator force(fs)','location','SouthWest');
title({'Gain from road dist(r) and actuator force(fs)';
	'to body accel(ab) and suspension travel(sd)'});

运行代码，可得

  ‘Because of the imaginary-axis zeros, feedback control cannot improve the response from road disturbance to body acceleration at the tire-hop frequency, and from to suspension deflection at the rattlespace frequency.
  Moreover, because of the relationship and the fact that the wheel position roughly follows at low frequency (less than 5 rad/s), there is an inherent trade-off between passenger comfort and suspension deflection: any reduction of body travel at low frequency will result in an increase of suspension deflection.’
  从bode图中可以看出：输入为$f_s$，输出分别为$a_b$、$s_d$时，都在某个频率下出现了增益衰减的尖峰。则意味着系统在该频率下的增益显著下降，在这个频率下，系统对输入信号的响应非常弱，几乎没有放大作用。
  示例还讨论了在低频输入下，乘客舒适度和悬架挠度的固有权衡：任何车身位移的减少会导致悬架挠度的增加³。

2.液压驱动器模型

  使用一阶传递函数$\frac{1}{1+\frac{s}{60}}$表示名义液压驱动器，其最大位移为0.05m。
  示例中使用了一系列执行器模型来解释执行器和四分之一车模型的建模误差和可变性。这一系列模型由一个名义（标称）模型和频率依赖的不确定性组成。
  低频时（< 3 rad/s），模型可以在标称值上下40%变化；模型不确定度在15 rad/s时超过100%，在大约1000 rad/s时达到2000%。
  加权函数$W_{unc}$用于随频率调制不确定性的数量。
代码如下:(Live script)

ActNom = tf(1,[1/60 1]);
Wunc = makeweight(0.40,15,3);  %低频处幅值20lg0.4；高频处幅值为20lg3；在频率为15 rad/s时经过幅值20lg1
%(ultidyn)Uncertain linear time-invariant dynamics
unc = ultidyn('unc',[1,1],'SampleStateDim',5);  
Act = ActNom*(1 + unc * Wunc);
Act.InputName = 'u';
Act.OutputName = 'fs';
bodemag(unc);

bodemag(Act);

rng('default');		%设置随机数生成器为默认，确保 ultidyn 生成的不确定性模型（unc）在每次运行时是可重复的。
bode(Act,'b',Act.NominalValue,'r+',logspace(-1,3,120))

系统模型

1.系统框图（matlab示例）

• 将路面干扰r建模为信号d，为了模拟7cm的道路干扰带宽 ，使用恒定权重系数$W_{road}$=0.07。 (To model broadband road deflections of magnitude seven centimeters, we use the constant weight.)
• 前文提到过乘客舒适度与悬架之间的权衡，可以对悬架位移$s_d$、车身加速度$a_b$引入3组权重系数，分别对应’Comfort’、'Balance’和’Road Handling’模式。从’Comfort’到’Road Handling’模式，权重系数$W_{s_d}$逐渐增加，$W_{a_b}$逐渐减小。
• 将传感器的测量噪声建模为信号$d_2$和$d_3$，分别对应车身加速度传感器和悬架位移传感器，并分别使用权重系数$W_{d_2}$=0.01和$W_{d_3}$=0.5模拟噪声带宽。 （We use $W_{d_2}$=0.01 and $W_{d_3}$=0.5 to model broadband sensor noise of intensity 0.01 and 0.5, respectively.）
  

  控制目标可以被描述为干扰抑制目标：最小化干扰$d_1、d_2、d_3$对控制力u、悬架位移$s_d$和车身加速度$a_b$的加权组合影响，当使用范数$H_{\infty }$（峰值增益）来衡量 “影响” 时，这相当于设计一个控制器，该控制器最小化从干扰输入$d_1、d_2、d_3$到误差信号$e_1、e_2、e_3$的范数。
(The control objectives can be reinterpreted as a disturbance rejection goal: Minimize the impact of the disturbances on a weighted combination of control effort u, suspension travel $s_d$, and body acceleration $a_b$. When using the $H_{\infty }$norm (peak gain) to measure “impact”, this amounts to designing a controller that minimizes the norm from disturbance inputs $d_1,d_2,d_3$ to error signals $e_1,e_2,e_3$.)
  此外，示例中提到 ‘Use a high-pass filter for $W_{act}$ to penalize high-frequency content of the control signal and thus limit the control bandwidth.’ （对于高通滤波器限制控制信号的高频部分从而限制控制带宽很是不解。）
代码如下：

Wroad = ss(0.07); Wroad.u = 'd1'; Wroad.y = 'r';
Wact = 0.8*tf([1 50],[1 500]); Wact.u = 'u'; Wact.y = 'e1';
Wd2 = ss(0.01); Wd2.u = 'd2'; Wd2.y = 'Wd2';
Wd3 = ss(0.5); Wd3.u = 'd3'; Wd3.y = 'Wd3';

  建立’HandlingTarget’传递函数作为道路干扰到悬架位移输$s_d$出的闭环控制目标；建立’ComfortTarget’传递函数作为道路干扰到车身加速度$a_b$输出的闭环控制目标。‘Because of the actuator uncertainty and imaginary-axis zeros, only seek to attenuate disturbances below 10 rad/s.’

HandlingTarget = 0.04 * tf([1/8 1],[1/80 1]);
ComfortTarget = 0.4 * tf([1/0.45 1],[1/150 1]);

Targets = [HandlingTarget; ComfortTarget];
bodemag(qcar({'sd','ab'}, 'r')*Wroad, 'b', Targets, 'r--',{1,1000}),grid;
title('Response to road disturbance')
legend('Open-loop','Closed-loop target')

  定义权重$W_{s_d}$和$W_{a_b}$为’ComfortTarget’和’Handling Target’的倒数。为了对应三种模式，建立三套权重$(\beta W_{s_d},(1-\beta )W_{a_b})$ 对应不同模式：comfort($\beta =0.01$)，balanced($\beta =0.5$)、handling($\beta =0.99$)

beta = reshape([0.01 0.2 0.66],[1 1 3]);
Wsd = beta / HandlingTarget;
Wsd.u = 'sd'; Wsd.y = 'e3';
Wab = (1-beta) / ComfortTarget;
Wab.u = 'ab'; Wab.y = 'e2';

使用’connect()'连接系统

sdmeas = sumblk('y1 = sd + Wd2');
abmeas = sumblk('y2 = ab + Wd3');
ICinputs = {'d1';'d2';'d3';'u'};
ICoutputs = {'e1';'e2';'e3';'y1';'y2'};
qcaric = connect(qcar(2:3,:),Act,Wroad,Wact,Wab,Wsd,Wd2,Wd3,sdmeas,abmeas,ICinputs,ICoutputs)

其中，'qcar(2:3, : )'为以下状态空间方程：

H$\infty$设计

1.性能指标及对应控制器

使用’hinfsyn’为混合因子$\beta$的每个值计算H$\infty$控制器。

ncont = 1; % 一个控制信号，u
nmeas = 2; %两个测量信号，sd和ab
K = ss(zeros(ncont,nmeas,1));
gamma = zeros(3,1);
for i = 1:3
	[K(:, :, i),~,gamma(i)] = hinfsyn(qcaric(:, :, i),nmeas, ncont);
end

gamma

2.闭环模型

K.u = {'sd','ab'}; K.y = 'u';
CL = connect(qcar, Act.Nominal,K,'r',{'xb';'sd';'ab'});
qcar(:,'r');
bodemag(qcar(:, 'r'),'b',CL(:, :, 1), 'r-.', CL(:, :, 2), 'm-.',CL(:, :, 3), 'k-.',{1,140}),grid
legend('Open loop', 'Comfort', 'Balanced', 'Road Handling','location','Southeast');
title('Body travel, Body acceleration and suspension deflection due to road')

---

1. 与MATLAB原示例代码中不同的是，矩阵中输入量’fs’ 和’r’的前后顺序相反，qcar.InputName = {‘fs’;‘r’}，这会导致输入矩阵B和直接传递矩阵D的略微差异。 ↩︎

2. 主要区别总结（AI）
   计算内容：
   zero 函数计算的是传递函数的零点，适用于分析单个传递函数的行为。这些零点是传递函数分子多项式的根。
   tzero 函数计算的是系统的传递零点，它考虑了系统的所有输入和输出。传递零点是系统在特定输入输出组合下可能发生的系统耦合失效点，反映了系统结构的特定性质。
   适用范围：
   zero 函数可以用于分析单个传递函数，适合 SISO 和 MIMO 系统中的单一通道分析。
   tzero 函数更适合用于分析整个系统的传递零点，尤其是在 MIMO 系统中，它能够捕捉到系统输入输出间的复杂关系。
   系统分析中的应用：
   如果你仅需要分析某个输入输出通道的传递函数零点，可以使用 zero 函数。
   如果你需要分析系统的结构特性、奇异值和传递零点等全局性质，尤其是在 MIMO 系统中，那么 tzero 函数是更合适的选择。 ↩︎

3. (AI):
   **低频（<5 rad/s）：**在低频率下，车轮的位置通常与车身的位置变化趋势较为一致，这意味着路面的不平整会以较低的频率传递到车身上。
   **车身运动的减少：**为了提高乘客的舒适度，设计悬架系统时通常希望减少车身在低频率下的运动（即减少车身位移和加速度）。
   悬架挠度对乘客舒适度的影响主要体现在低频率下的车身运动和悬架挠度之间的权衡。如果悬架系统过于刚硬，虽然可以减少悬架挠度，但会将更多的振动传递到车身上，降低舒适度。反之，如果悬架过于柔软，尽管可以提高舒适度，但会导致更大的悬架挠度，影响车辆的操控性和稳定性。 ↩︎