MATLAB - 模型预测控制（MPC）使用 ADMM 求解器四分之一汽车悬架悬挂系统动力学控制

系列文章目录

前言

一、四分车悬架模型

二、道路干扰剖面

三、设计模型预测控制器

四、设置优化求解器

五、辅助函数

前言

本例展示了如何为四分之一汽车悬架系统设计模型预测控制器 (MPC)，采用乘法交替方向法 (ADMM) 求解器来控制主动悬架的动态。MPC 控制器利用预测模型来预测和减轻路面不规则性和负载偏移的影响，从而确保最佳的车辆稳定性和乘客舒适性。

一、四分车悬架模型

四分之一汽车悬架系统模型是一个双自由度系统，通过关注一个车轮和四分之一车身质量的垂直动态来简化汽车悬架系统。该模型包括两个质量：弹簧质量 ms（对应车身）和非弹簧质量 mu（代表车轮组件）。悬架的弹簧和阻尼器（分别具有常数 ks 和 bs）将两个质量连接起来，而轮胎则被模拟为具有常数 kt 的弹簧。主动悬架力通过控制输入 u(t) 引入，路面效应通过激励 v(t) 模拟。系统动力学方程为

$\begin{cases}m_s\ddot{x}_s=F\\m_u\ddot{x}_u=-F-k_t(z_u-\nu(t))&\end{cases}$ ， $F=-k_{s}(z_{s}-z_{u})-b_{s}\left(\dot{z}_{s}-\dot{z}_{u}\right)+u.$

将状态向量定义为

$z(t)=\begin{bmatrix}\mathbf{x}\\\dot{\mathbf{x}}\end{bmatrix}$ ， $\mathbf{x}=\begin{bmatrix}x_s\\\\x_u\end{bmatrix}$ ， $\dot{\mathbf{x}}=\begin{bmatrix}\dot{x}_s\\\\\dot{x}_u\end{bmatrix}$

四分之一汽车悬架系统模型是通过线性状态空间系统实现的。

% Physical parameters
ms = 300;    % kg
mu = 60;     % kg
bs = 1000;   % N/m/s
ks = 16000;  % N/m
kt = 190000; % N/m% Define the mass matrix M.
M = [ms 0; 0 mu];
% Define the stiffness matrix F.
F = [-ks ks; ks -(ks + kt)];
% Define the damping matrix B.
B = [-bs bs; bs -bs];
% Define the input matrix q.
q = [1; -1];
% Define the disturbance matrix d.
d = [0; kt];% Define the state space matrices.
A = [zeros(2,2), eye(2,2); M\F, M\B];
B = [zeros(2,2); M\q M\d];
C = [1 -1 0 0];
D = 0;
plant = ss(A,B,C,D);

二、道路干扰剖面

在本例中，道路干扰被模拟为频率为 f Hz 的正弦波颠簸。扰动仅在 2/f 秒和 3/f 秒之间有效，以仿真车辆遇到的单次颠簸。扰动的振幅按系数（本例中为 0.025）缩放，以表示颠簸的严重程度。

% Define the frequency of the road disturbance.
f = 4; % Hz
Tstop = 2;% Define the sampling time.
Ts = 0.01;% Create a time vector from 0 to Tstop.
% with increments of Ts.
T = (0:Ts:Tstop)';% Define the road disturbance profile.
v = 0.025*(1 - cos(2*pi*f*T));% Activate the disturbance only between 2/f and 3/f seconds.
v(T < 2*1/f) = 0;
v(T > 3*1/f) = 0;

用列向量 [0 v] 表示道路干扰剖面，进行仿真。

Unn = zeros(length(T),2);
Unn(:,2) = v;

在没有主动控制输入的情况下仿真悬架系统，以提供系统在受到规定路面干扰时的基准响应。

% Use lsim to simulate the uncontrolled system response.
[Yuc,Tuc,Xuc] = lsim(plant,Unn,T);

绘制路面扰动对四分之一汽车悬架系统的影响以及悬架挠度随时间变化的曲线。

figure
title("Response of Uncontrolled Suspension to Road Disturbance")
plot(T, v, "DisplayName", "Road Disturbance"); 
hold on
plot(Tuc, Yuc, "DisplayName", "Suspension Deflection"); 
hold off
legend
ylabel("Meters")
xlabel("Time (seconds)")

三、设计模型预测控制器

默认情况下，被控对象状态空间模型的所有输入信号均为可控输入变量，所有输出均为可观测输出。将模型的第二个输入信号定义为与道路效应相对应的测量扰动。

plant = setmpcsignals(plant,"MD",2);

-->Assuming unspecified input signals are manipulated variables.

创建一个预测范围为 40 步、控制范围为 4 步、采样时间为 0.01 秒的 MPC 控制器。

Ts = 1e-2;
p = 40;
m = 4;
mpcobj = mpc(plant,Ts,p,m);

-->"Weights.ManipulatedVariables" is empty. Assuming default 0.00000.
-->"Weights.ManipulatedVariablesRate" is empty. Assuming default 0.10000.
-->"Weights.OutputVariables" is empty. Assuming default 1.00000.

将输出变量限制在 [-0.05, 0.05] 范围内。

mpcobj.OutputVariables(1).Min = -0.05;
mpcobj.OutputVariables(1).Max = 0.05;

设置 MPC 控制参数，以定义比例系数和权重。比例系数对输出变量进行归一化处理，权重则指定每个变量在成本函数中的相对重要性。在这种情况下，通过为可控变量的速率设置更高的权重，对控制输入变量的变化给予更多的重视（即更多的成本）。

mpcobj.OutputVariables(1).ScaleFactor = 1e-5;
mpcobj.Weights.ManipulatedVariables = 1;
mpcobj.Weights.ManipulatedVariablesRate = 10;
mpcobj.Weights.OutputVariables = 5;

审查 MPC 控制器设计，确保其按预期配置。审查功能可生成 MPC 控制器的设计审查。

review(mpcobj)

-->Converting model to discrete time.
-->Assuming output disturbance added to measured output #1 is integrated white noise.
-->"Model.Noise" is empty. Assuming white noise on each measured output.
Test report has been saved to:/tmp/Bdoc24b_2855429_3392907/index.html

将参考信号 r 定义为零向量，代表无悬架偏转的理想输出。为仿真定义路面扰动 v，其比例代表颠簸的严重程度，且仅在 2/f 至 3/f 秒的窗口内有效。然后设置 MPC 控制器的仿真选项，指定控制器不预测未来的可观测干扰。

% Define the reference signal as a vector of zeros. 
r = zeros(length(T),1);% Define the road disturbance "v" for the simulation.
v = 0.025 * (1 - cos(2*pi*f*T));
v(T < 2*1/f) = 0;
v(T > 3*1/f) = 0;% Set simulation options for the MPC controller. 
params = mpcsimopt;
params.MDLookAhead = "off";

四、设置优化求解器

对于只包含连续可控输入变量的 MPC 问题，“active-set ”算法是默认的优化方法。mpc 对象的优化器属性通常具有这种结构。

mpcobj.Optimizer

ans = struct with fields:OptimizationType: 'QP'Solver: 'active-set'SolverOptions: [1x1 mpc.solver.options.qp.ActiveSet]MinOutputECR: 0UseSuboptimalSolution: 0CustomSolver: 0CustomSolverCodeGen: 0

因为可控输入变量是连续的，所以优化问题是一个二次规划（QP）问题。如果有有限（离散）控制集，问题就变成了混合整数二次规划 (MIQP) 问题。更多信息，请参阅有限控制集 MPC。

在本例中，请指定 ADMM 求解器，它非常适合大规模和分布式优化问题。

mpcobj.Optimizer.Solver = "admm";

要修改求解器参数并使优化过程符合特定要求，可以通过 MPC 对象中优化器属性的 SolverOptions 属性访问和调整设置。通过 SolverOptions 属性，可以微调求解器行为的各个方面，如收敛公差和迭代限制选项。

mpcobj.Optimizer.SolverOptions

ans = ADMM with properties:AbsoluteTolerance: 1.0000e-04RelativeTolerance: 1.0000e-04MaxIterations: 1000RelaxationParameter: 1PenaltyParameter: 1.6000StepSize: 1.0000e-06

使用 “admm ”求解器仿真闭环系统。

[Yad, Tad, Uad, XPad, XCad] = sim(mpcobj,Tstop/Ts + 1,r,v,params);

-->Converting model to discrete time.
-->Assuming output disturbance added to measured output #1 is integrated white noise.
-->"Model.Noise" is empty. Assuming white noise on each measured output.

然后使用 “主动集 ”求解器仿真系统。该算法是在每个控制区间求解 QP 问题的默认选项。

mpcobj.Optimizer.Solver = "active-set";
[Yas, Tas, Uas, XPas, XCas] = sim(mpcobj,Tstop/Ts + 1,r,v,params);

-->Converting model to discrete time.
-->Assuming output disturbance added to measured output #1 is integrated white noise.
-->"Model.Noise" is empty. Assuming white noise on each measured output.

最后，使用 “内点 ”求解器仿真系统。这种算法是求解 QP 问题的另一种选择，与 “主动集 ”算法相比，可能具有不同的性能特点。

mpcobj.Optimizer.Solver = "interior-point";
[Yip, Tip, Uip, XPip, XCip] = sim(mpcobj,Tstop/Ts + 1,r,v,params);

-->Converting model to discrete time.
-->Assuming output disturbance added to measured output #1 is integrated white noise.
-->"Model.Noise" is empty. Assuming white noise on each measured output.

绘制三次仿真的结果。有关 plotResults 辅助函数的代码，请参阅辅助函数。

TT = {Tas, Tip, Tad, Tuc};
XX = {Yas, Yip, Yad, Yuc};
UU = {Uas, Uip, Uad, Tuc*0};plotResults(TT,XX,UU)
subplot(211); 
legend("active-set","interior-point","admm","uncontrolled")

五、辅助函数

这段代码定义了 plotResults 辅助函数。该函数绘制仿真结果，以比较它们的性能。图中显示了随时间变化的系统输出（悬挂偏转）和控制输入（主动悬挂力）。

function plotResults(TT,XX,UU)subplot(211); 
for i = 1:length(TT)plot(TT{i}, XX{i})hold on;
end
ylabel("Meters")
xlabel("Time, s")subplot(212); 
for i = 1:length(TT)plot(TT{i}, UU{i})hold on;
end
ylabel("N")
xlabel("Time, s")
end