悬架模型自抗扰控制：1/4车主动悬架控制Simulink模型资料，包含Plant_Passive、Plant_Act、NLSEF、ESO、ADRC_QuaterCar2018a、ADRC_QuaterCar等文件，模型自带PID、自抗扰控制，支持ADRC算法，仿真结果有车身加速度、悬架动挠度对比，适合做悬架控制、车辆动力学、控制算法研究。

针对你提出的 1/4 车主动悬架自抗扰控制（ADRC）完整的 Simulink + MATLAB 仿真方案。这套方案涵盖了你提到的所有关键模块：被动悬架（Plant_Passive）、主动悬架（Plant_Act）、扩张状态观测器（ESO）、非线性状态误差反馈（NLSEF）以及完整的 ADRC 控制器。

以下是详细的模型搭建逻辑和核心代码。

系统数学模型

1/4 车模型包含两个自由度：车身质量 m_s（簧上质量）和车轮质量 m_u（簧下质量）。

运动微分方程

设车身位移为 z_s，车轮位移为 z_u，路面输入为 z_r，控制力为 u。

m_s ddot{z}_s = -k_s(z_s - z_u) - c_s(dot{z}_s - dot{z}_u) + u

m_u ddot{z}_u = k_s(z_s - z_u) + c_s(dot{z}_s - dot{z}_u) - k_t(z_u - z_r) - u

其中：
k_s, c_s：悬架弹簧刚度和阻尼系数。
k_t：轮胎刚度。

Simulink 模型结构

在 Simulink 中，你需要构建以下主要子系统：

Plant_Act：主动悬架模型（包含控制输入 u）。
Plant_Passive：被动悬架模型（u=0，用于对比）。
ESO：扩张状态观测器（观测总扰动 f）。
NLSEF：非线性状态误差反馈（计算控制量 u_0）。
ADRC Controller：顶层封装，包含 ESO 和 NLSEF。

核心 MATLAB Function 代码

你需要将以下代码放入 Simulink 的 MATLAB Function Block 中。

扩张状态观测器（ESO）

这是 ADRC 的核心，用于实时估计“总扰动”（包括路面干扰和参数摄动）。

function [z1, z2, z3, f] = ESO(z_s, u, b0, beta1, beta2, beta3, x1, x2, x3, Ts)
% ESO: 扩张状态观测器
% 输入: z_s (车身加速度/位移), u (控制量), b0 (控制增益估计值)
% beta1/2/3 (观测器带宽参数), x1/2/3 (上一时刻状态), Ts (采样时间)
% 输出: z1, z2, z3 (状态估计), f (总扰动估计)

% 状态更新 (使用欧拉法)
e = z_s - x1; % 误差

% 观测器律
x1 = x1 + Ts * (x2 + beta1 * e);
x2 = x2 + Ts * (x3 + beta2 * e + b0 * u);
x3 = x3 + Ts * (beta3 * e);

z1 = x1;
z2 = x2;
z3 = x3;
f = x3; % 总扰动估计
end

非线性状态误差反馈（NLSEF）

这是 ADRC 的控制律，取代传统的 PID。

function u = NLSEF(e, edot, w, r, h, z1, z2, z3, f, w_c, b0)
% NLSEF: 非线性状态误差反馈
% 输入: e (误差), edot (误差微分), w (带宽), r (误差权重), h (采样时间)
% z1/z2/z3 (ESO状态), f (总扰动), b0 (控制增益)

% 非线性误差组合 (使用符号函数或饱和函数)
beta1 = 5 * w; % 通常设为 5~10倍带宽
beta2 = 10 * w^2;

% 简化的线性状态误差反馈 (易于工程实现)
u0 = -beta1 * e - beta2 * edot;

% 补偿总扰动
u = (u0 - f) / b0;

% 限幅处理 (防止控制量过大)
u_max = 2000; % 根据执行器能力设定
u_min = -2000;
if u > u_max
u = u_max;
elseif u < u_min
u = u_min;
end
end

完整的 ADRC 控制器代码（顶层）

function [u, z1, z2, z3, f] = ADRC_Controller(z_s_ref, z_s_act, u_prev, b0, beta, w_c, Ts, …)
persistent x1 x2 x3;

if isempty(x1)
    x1 = 0; x2 = 0; x3 = 0;
end

% 1. 计算误差
e = z_s_ref - z_s_act;
edot = (e - e_prev) / Ts;
e_prev = e;

% 2. 调用 ESO
[z1, z2, z3, f] = ESO(z_s_act, u_prev, b0, beta(1), beta(2), beta(3), x1, x2, x3, Ts);
x1 = z1; x2 = z2; x3 = z3; % 更新状态

% 3. 调用 NLSEF
u = NLSEF(e, edot, w_c, 1, Ts, z1, z2, z3, f, b0);

end

参数设置建议

在 MATLAB 工作区中定义以下参数：

% 车辆参数
ms = 290; % 簧上质量 (kg)
mu = 50; % 簧下质量 (kg)
ks = 16812; % 悬架刚度 (N/m)
cs = 1000; % 悬架阻尼 (N*s/m)
kt = 190000;% 轮胎刚度 (N/m)

% ADRC 参数
b0 = 1/ms; % 控制增益估计值
beta = [100, 2500, 1e5]; % ESO 带宽参数 (通常设为 w0^1, w0^2, w0^3)
w_c = 50; % 控制器带宽 (rad/s)

仿真结果分析

运行模型后，你应该对比以下曲线：

车身加速度 (ddot{z}_s)：ADRC 应显著降低峰值，提高乘坐舒适性。
悬架动挠度 (z_s - z_u)：ADRC 应限制在物理行程范围内。
轮胎动载荷：确保接地性。

Simulink 中的 S-Function 或 MATLAB Function 模块。

被控对象代码（Plant_Act）

这是主动悬架的物理模型，包含簧上质量 m_s 和簧下质量 m_u。

S-Function 模板代码

define S_FUNCTION_NAME plant_active_suspension
define S_FUNCTION_LEVEL 2

include “simstruc.h”

// 参数索引
define MS ssGetSFcnParam(S,0) // 簧上质量
define MU ssGetSFcnParam(S,1) // 簧下质量
define KS ssGetSFcnParam(S,2) // 悬架刚度
define CS ssGetSFcnParam(S,3) // 悬架阻尼
define KT ssGetSFcnParam(S,4) // 轮胎刚度

define U_PLANT 0 // 输入：控制力 u
define U_ROAD 1 // 输入：路面位移 zr

static void mdlInitializeSizes(SimStruct *S)
{
ssSetNumSFcnParams(S, 5); // 5 个参数
if (ssGetNumSFcnParams(S) != ssGetSFcnParamsCount(S)) return;

// 设置输入端口：控制力 u + 路面输入 zr
ssSetNumInputPorts(S, 2);
ssSetInputPortWidth(S, 0, 1);
ssSetInputPortWidth(S, 1, 1);
ssSetInputPortDirectFeedThrough(S, 0, 0);

// 设置输出端口：车身加速度、悬架动挠度等
ssSetNumOutputPorts(S, 4);
ssSetOutputPortWidth(S, 0, 1); // 车身加速度
ssSetOutputPortWidth(S, 1, 1); // 悬架动挠度
ssSetOutputPortWidth(S, 2, 1); // 车轮加速度
ssSetOutputPortWidth(S, 3, 1); // 轮胎动载荷

// 设置连续状态：4个状态 [zs, zs_dot, zu, zu_dot]
ssSetNumContStates(S, 4);

ssSetNumSampleTimes(S, 1);
ssSetOptions(S, SS_OPTION_EXCEPTION_FREE_CODE);

}

static void mdlOutputs(SimStruct *S, int_T tid)
{
real_T *y = ssGetOutputPortSignal(S,0);
real_T *x = ssGetContStates(S);
real_T *u = ssGetInputPortSignal(S, U_PLANT);
real_T *zr = ssGetInputPortSignal(S, U_ROAD);

// 状态变量
real_T zs = x[0];     // 车身位移
real_T dzs = x[1];    // 车身速度
real_T zu = x[2];     // 车轮位移
real_T dzu = x[3];    // 车轮速度

// 参数
real_T ms = *mxGetPr(MS);
real_T mu = *mxGetPr(MU);
real_T ks = *mxGetPr(KS);
real_T cs = *mxGetPr(CS);
real_T kt = *mxGetPr(KT);

// 计算输出
// 1. 车身加速度 (用于评价舒适性)
y[0] = ( -k(zs-zu) - cs(dzs-dzu) + u[0] ) / ms;

// 2. 悬架动挠度 (zs - zu)
y[1] = zs - zu;

// 3. 车轮加速度
y[2] = ( k(zs-zu) + cs(dzs-dzu) - kt*(zu-zr[0]) - u[0] ) / mu;

// 4. 轮胎动载荷 (归一化或实际力)
y[3] = kt * (zu - zr[0]);

}

static void mdlUpdate(SimStruct *S, int_T tid) {}

static void mdlDerivatives(SimStruct *S)
{
real_T *x = ssGetContStates(S);
real_T *dx = ssGetdX(S);
real_T *u = ssGetInputPortSignal(S, U_PLANT);
real_T *zr = ssGetInputPortSignal(S, U_ROAD);

// 状态变量
real_T zs = x[0];
real_T dzs = x[1];
real_T zu = x[2];
real_T dzu = x[3];

// 参数
real_T ms = *mxGetPr(MS);
real_T mu = *mxGetPr(MU);
real_T ks = *mxGetPr(KS);
real_T cs = *mxGetPr(CS);
real_T kt = *mxGetPr(KT);

// 状态方程
dx[0] = dzs;
dx[1] = ( -k(zs-zu) - cs(dzs-dzu) + u[0] ) / ms;
dx[2] = dzu;
dx[3] = ( k(zs-zu) + cs(dzs-dzu) - kt*(zu-zr[0]) - u[0] ) / mu;

}

ifdef MATLAB_MEX_FILE
include “simulink.c”
else
include “cg_sfun.h”
endif

ADRC 控制器核心算法（ESO + NLSEF）

这部分通常放在 MATLAB Function 模块中，接收车身加速度（或位移）作为反馈。

function u = fcn(zs_measured, zr, Ts, …
b0, w_o, w_c, …
ms, mu, ks, cs, kt)

% ADRC Controller for Active Suspension
% 输入: zs_measured - 测量的车身加速度 (或位移，取决于观测器设计)
% 输出: u - 控制力

persistent z1 z2 z3 u0

if isempty(z1)
z1 = 0; % 观测状态1 (估计的车身速度)
z2 = 0; % 观测状态2 (估计的总扰动 f)
z3 = 0; % 观测状态3 (估计的路面干扰相关项，可选)
end

% — 1. 扩张状态观测器 (ESO) —
% 这里假设观测器观测 [z_s_dot, f, …]
% 为了简化，这里展示二阶ESO形式
% f 代表总扰动: f = k(zs-zu) + cs(dzs-dzu) - u_disturbance

% 通常需要将加速度信号积分一次得到速度，或者直接观测加速度
% 这里假设 zs_measured 是车身加速度 y = zs_dot_dot

% 简化的ESO更新逻辑 (伪代码，需根据实际微分方程调整)
% beta1, beta2 为观测器带宽参数
beta1 = 2 * w_o;
beta2 = w_o^2;

% 观测器误差
e = z1 - zs_measured;

% 状态更新 (欧拉法)
z1 = z1 + Ts * (z2 + beta1 * e);
z2 = z2 + Ts * (z3 + beta2 * e + b0 * u0);
% 注意：这里 b0*u0 是控制输入对状态的影响

% — 2. 非线性状态误差反馈 (NLSEF) —
% 设定目标：车身加速度为 0 (跟随平坦路面)

e = 0 - z1; % 误差

% 简单的线性反馈或非线性反馈 (fal函数)
% u0 = -k1 * e - k2 * z2;
% 其中 k1, k2 与 w_c 相关

k1 = 2 * w_c;
k2 = w_c^2;

u0 = -k1 * e - k2 * z2;

% — 3. 计算最终控制量 u —
% 根据观测到的总扰动 z2 进行补偿
% u = (u0 - z2) / b0;
% 注意：b0 是控制输入通道的增益

% 这里的 b0 通常与 ms 有关，假设 b0 = 1/ms
% 实际中 b0 需要通过系统辨识或理论推导得出

u = (u0 - z2) / b0;

% 限幅处理 (防止控制力过大)
u_max = 2000; % 例如 ±2000N
u_min = -2000;
u = max(min(u, u_max), u_min);

PID 对比代码

作为基准，PID 控制器非常简单：

function u = fcn(zs_measured, zr, Ts, Kp, Ki, Kd)

persistent e_prev integral

if isempty(e_prev)
e_prev = 0;
integral = 0;
end

% 误差 (假设目标是车身静止，路面干扰为zr)
e = 0 - zs_measured;

% 积分项
integral = integral + e * Ts;

% 微分项
de = (e - e_prev) / Ts;

% PID 输出
u = Kp * e + Ki * integral + Kd * de;

e_prev = e;

参数设置建议

在 Simulink 的工作区（Workspace）中定义以下参数：

% 物理参数
ms = 290; % 簧上质量 (kg)
mu = 40; % 簧下质量 (kg)
ks = 16800; % 悬架刚度 (N/m)
cs = 1000; % 悬架阻尼 (N.s/m)
kt = 190000;% 轮胎刚度 (N/m)

% ADRC 参数 (需整定)
b0 = 1/ms; % 控制增益近似
w_o = 1000; % 观测器带宽 (通常取系统带宽的5-10倍)
w_c = 100; % 控制器带宽

以上代码和参数配置为你构建了完整的主动悬架控制仿真基础。

黑色线 (Passive)：被动悬架（无控制），作为基准。
棕色线 (Active)：理想主动悬架（通常指天钩或最优控制参考）。
浅黄线 (Estimate)：ADRC 中 ESO（扩张状态观测器）对总扰动的估计值。
紫色线 (PID)：传统 PID 控制器的结果。

Simulink 环境，包含 1/4 车模型、路面输入模块、PID 控制器 和 ADRC 控制器。

以下是复现该曲线的核心代码，主要包含 ADRC 控制器的 MATLAB Function 代码 和 PID 控制器的代码。

ADRC 控制器代码 (MATLAB Function)

这是实现“Estimate”和主动控制的核心。

function [u, z1, z2, z3, f] = fcn_adrc(zs, d_zs, u, b0, beta1, beta2, beta3, weso, h)
% ADRC 控制器：包含 ESO 和 NLSEF
% zs: 车身位移
% d_zs: 车身速度
% u: 控制力
% b0: 控制量增益
% beta1, beta2, beta3: ESO 观测器增益
% weso: 带宽
% h: 采样时间

% 参数初始化
persistent x1 x2 x3
if isempty(x1)
x1 = 0;
x2 = 0;
x3 = 0;
end

% — 扩张状态观测器 (ESO) —
% 状态更新
x1 = x1 + h * (x2 - beta1 * (x1 - zs));
x2 = x2 + h * (x3 - beta2 * (x1 - zs) + b0 * u);
x3 = x3 + h * (-beta3 * (x1 - zs));

% 输出估计值
z1 = x1; % 估计的位移
z2 = x2; % 估计的速度
z3 = x3; % 估计的总扰动 f

% — 非线性状态误差反馈 (NLSEF) —
% 定义误差
e1 = 0 - z1; % 位移误差 (设定值为0)
e2 = 0 - z2; % 速度误差 (设定值为0)

% 非线性组合 (最简单的线性组合，也可以用fal函数)
beta1_nl = weso^2;
beta2_nl = 2 * weso;

% 计算控制量
u0 = (beta1_nl * e1 + beta2_nl * e2 - z3) / b0;

% 饱和限幅 (防止控制力过大)
u_max = 2000; % 根据实际情况调整
u_min = -2000;
u = max(min(u0, u_max), u_min);

% 输出扰动估计用于绘图
f = z3;

PID 控制器代码 (MATLAB Function)

这是实现“PID”曲线的核心。

function u = fcn_pid(zs, d_zs, Kp, Ki, Kd, h)
% PID 控制器
% zs: 车身位移
% d_zs: 车身速度
% Kp, Ki, Kd: PID 参数

persistent integral
if isempty(integral)
integral = 0;
end

% 计算误差
e = 0 - zs; % 设定位移为0

% 积分项
integral = integral + e * h;

% PID 输出
u = Kp * e + Ki * integral + Kd * (-d_zs); % 速度反馈为负

% 限幅
u_max = 2000;
u_min = -2000;
u = max(min(u, u_max), u_min);

1/4 车模型参数 (Model Parameters)

在 Simulink 的 Model Workspace 中设置以下参数：

% 悬架参数
ms = 240; % 簧上质量 (kg)
mu = 32; % 簧下质量 (kg)
ks = 16000; % 悬架刚度 (N/m)
cs = 980; % 悬架阻尼 (N.s/m)
kt = 190000;% 轮胎刚度 (N/m)

% 控制参数
b0 = 1/ms; % ADRC 控制增益

% ESO 观测器增益 (通常根据带宽 weso 设计)
weso = 50; % 观测器带宽
beta1 = 3weso;
beta2 = 3weso^2;
beta3 = weso^3;

% PID 参数 (需要根据 Ziegler-Nichols 或其他方法整定)
Kp = 2000;
Ki = 100;
Kd = 500;

% 仿真参数
Ts = 0.001; % 采样时间 (s)

路面输入代码 (Road Input)

为了复现图中的随机路面，使用白噪声通过滤波器生成：

% 路面不平度滤波器传递函数
% Gq = (pin0)^2 * q0 * V / (s^2 + (pin0*V)s)
n0 = 0.127; % 参考空间频率 (m^-1)
V = 20; % 车速 (m/s)
q0 = 256e-6;% 路面功率谱密度 (m^3)

num = [(pin0)^2 * q0 * V];
den = [1, pin0*V, 0];
road_filter = tf(num, den);

% 在 Simulink 中使用 Band-Limited White Noise 输入到这个 Transfer Fcn

如何连接

Plant_Act：使用你之前的 S-Function 或 State-Space 模块，输入为 u 和 zr，输出为 zs, zu, ds。
ADRC Block：输入 zs, ds, u，输出控制力 u 和估计扰动 f。
PID Block：输入 zs, ds，输出控制力 u。
Scope：连接 zs-zu（悬架动挠度）到示波器，并设置标签为 Passive, Active, Estimate, PID。