车辆稳定域划分，车辆稳定边界拟合 车辆稳定性相平面MATLAB程序绘制。 根据确定的简化魔术公式轮胎模型，建立车辆非线性二自由度运动微分方程，并进而对相平面图进行绘制。 包括横摆角速度与质心侧偏角的相平面，以及质心侧偏角速度与质心侧偏角的相平面。 具有Simulink和m脚本两种版本，而且可对稳定区域划分，鞍点位置与车速，路面附着系数的三维关系图等。

function Fy = magic_formula(alpha, mu)
    B = 0.3; C = 1.6; D = mu*1500;
    Fy = D*sin(C*atan(B*alpha));
end

这个非线性特性直接影响车辆横摆力矩，参数B控制曲线斜率，C决定峰值形状，D和路面摩擦系数mu直接相关。注意当alpha超过8度时，侧向力开始饱和。

二自由度模型的状态方程写成这样：

function dxdt = vehicle_model(t,x,u)
    beta = x(1);  % 质心侧偏角
    r = x(2);     % 横摆角速度
    v = u(1);     % 车速
    mu = u(2);    % 摩擦系数
    
    m = 1500;     % 质量
    Iz = 2500;    % 转动惯量
    a = 1.2; b = 1.5; % 轴距分配
    
    alpha_f = beta + a*r/v - delta;
    alpha_r = beta - b*r/v;
    
    Fyf = magic_formula(alpha_f, mu);
    Fyr = magic_formula(alpha_r, mu);
    
    dxdt = [ (Fyf + Fyr)/(m*v) - r;
             (a*Fyf - b*Fyr)/Iz ];
end

注意分母里的车速v可能导致奇点，实际仿真时要加个速度下限。这个微分方程的非线性全藏在alphaf和alphar里，相平面图的复杂行为就是从这里长出来的。

相平面绘制核心是暴力遍历初始条件：

v_range = 10:5:40;  % 车速扫描范围
beta_grid = linspace(-0.3,0.3,20);
r_grid = linspace(-1,1,20);
[beta0, r0] = meshgrid(beta_grid, r_grid);

figure;
for v = v_range
    trajectories = cell(size(beta0));
    for i = 1:numel(beta0)
        [~,x] = ode45(@(t,x) vehicle_model(t,x,[v,0.8]), [0 5], [beta0(i),r0(i)]);
        trajectories{i} = x;
    end
    quiver(beta0, r0, diff_x, diff_r); % 绘制向量场
    hold on
end

这个双重循环跑起来挺吃硬件，建议用parfor并行加速。稳定边界出现在向量场发散的位置，这时候轨迹会快速偏离平衡点。

鞍点搜索用牛顿迭代法：

function [beta_eq, r_eq] = find_saddle(v, mu)
    fun = @(x) vehicle_model(0,x,[v,mu]);
    options = optimoptions('fsolve','Display','off');
    equilibrium = fsolve(fun, [0;0], options);
    
    J = jacobian(fun, equilibrium); % 计算雅可比矩阵
    [V,D] = eig(J);
    if any(real(diag(D)) > 0) % 存在正实部特征值
        beta_eq = equilibrium(1);
        r_eq = equilibrium(2);
    else
        beta_eq = NaN;
        r_eq = NaN;
    end
end

特征值实部的符号决定平衡点性质，鞍点必然有一正一负两个实部。三维关系图用meshgrid生成参数空间：

[V,MU] = meshgrid(10:2:40, 0.1:0.1:1);
saddle_points = arrayfun(@(v,mu) find_saddle(v,mu), V, MU);
surf(V, MU, saddle_points);

这个曲面会清晰显示临界车速随摩擦系数变化的趋势，当曲面出现断崖式下跌时，说明系统在该区域失去稳定。

Simulink版本的优势在于实时调参——挂着滑块控件边拖参数边看相图变化，适合课堂演示。不过真要搞批量计算还是脚本利索。代码全放在GitHub了，需要自取，记得点个star。下次飙车前先跑个仿真看看相图稳不稳？