MATLAB中搭建一个四分之一二自由度的半主动悬架模型并为它一个滑模控制器，并用前向欧拉法离散化，采样时间取0.01

以下为MATLAB代码：

% 定义系统参数
m = 250; % 质量
k1 = 20000; % 弹簧刚度
k2 = 20000; % 弹簧刚度
c1 = 500; % 阻尼系数
c2 = 500; % 阻尼系数
h = 0.01; % 采样时间

% 定义状态向量
% x1: 车身垂直位移
% x2: 车身垂直速度
% x3: 车轮垂直位移
% x4: 车轮垂直速度
x = zeros(4, 1);

% 定义输入向量
% u1: 车身加速度
% u2: 半主动阻尼器力
u = zeros(2, 1);

% 定义输出向量
y = x(1);

% 定义滑模控制器参数
K = 200; % 滑模控制器增益
epsilon = 0.1; % 滑模面宽度

% 定义离散化矩阵
A = [1 h 0 0; (-k1-k2)/m -c1/m c2/m 0; 0 0 1 h; k2/m c2/m (-k2*u(2))/m -c2/m];
B = [0; 1/m; 0; 0];

% 定义初始状态
x0 = [0; 0; 0; 0];

% 定义仿真时间
T = 10;
t = 0:h:T;

% 定义存储变量
X = zeros(4, length(t));
U = zeros(2, length(t));
Y = zeros(1, length(t));

% 滑模控制器
s = x(2) - x(4);
if abs(s) < epsilon
    sdot = 0;
else
    sdot = sign(s);
end
u(2) = (-K*s - c2*x(4) + c2*x(2)) / k2;

% 循环仿真
for i = 1:length(t)
    % 更新状态和输出
    xdot = A*x + B*u;
    x = x + h*xdot;
    y = x(1);
    
    % 存储状态和输入
    X(:, i) = x;
    U(:, i) = u;
    Y(i) = y;
    
    % 更新滑模控制器
    s = x(2) - x(4);
    if abs(s) < epsilon
        sdot = 0;
    else
        sdot = sign(s);
    end
    u(2) = (-K*s - c2*x(4) + c2*x(2)) / k2;
    
    % 更新输入
    u(1) = 0.5*sin(2*pi*t(i));
end

% 绘图
subplot(2, 1, 1);
plot(t, Y);
xlabel('Time (s)');
ylabel('Displacement (m)');
title('Quarter-Car Semi-Active Suspension');

subplot(2, 1, 2);
plot(t, U(2,:));
xlabel('Time (s)');
ylabel('Semi-Active Damper Force (N)');