相平面理解：

直观的比喻：碗里的玻璃球

想象一个底部平滑的碗（这就好比你固定的车速、路面和转角，它们决定了碗的形状）。

• 稳态：如果玻璃球正好在碗底，它就不动了。这就像你认为的“状态应该是固定的”。

• 相轨迹（动态演变）：如果我把玻璃球拿到碗壁的某个高处（这就是给了它一个初始的侧偏角和横摆角速度）然后松手。玻璃球会在碗里来回滚动，它的位置和速度时刻在变化，最终慢慢停在碗底。

相平面图，就是你把玻璃球放在碗壁上成百上千个不同的位置松手，然后用相机把球滚动的所有的轨迹全部重叠拍在同一张照片上。

怎么去实现呢？

如何联合simulink模型绘制相平面：

第一步：改造你的 Simulink 模型（关键！）

为了让 MATLAB 脚本能控制 Simulink，你需要让 Simulink 模型里的一些固定数值变成工作区变量。

1. 设置工况输入（设为变量）：

   在 Simulink 中，你的前轮转角、车速、纵向力、附着系数不要直接写死成数字。用 Constant（常数）模块，但里面的值分别填上变量名：delta_phase、V、Fxf、mu。

2. 设置初始状态（核心步骤）：

   找到你计算质心侧偏角 $\beta$ 和横摆角速度 $r$ 的那两个积分模块（Integrator, 图标是 1/s）。

   双击它们，把里面的 Initial condition（初始条件） 从 0 改成变量名：

   • $\beta$ 的积分模块初始条件设为：beta_init

   • $r$ 的积分模块初始条件设为：r_init

3. 设置输出提取（To Workspace）：

   在 $\beta$ 和 $r$ 的信号线上，分别接一个 To Workspace 模块。

   • 变量名分别设为 beta_out 和 r_out。

   • Save format（保存格式）强烈建议选择 Timeseries 或 Array（下面的脚本以 Timeseries 为例）。

不同工况设置：

将Fxf分别设置为-2kN、0kN和2kN,其中Fxf=-2kN表示前轮制动工况的漂移，Fxf=0kN表示后轮驱动车辆的漂移，Fxf=2kN表示前轮驱动工况的漂移。

用质心侧偏角β和横摆角速度r表征整个车辆系统,并将这两个状态量作为横纵坐标,建立相平面图,相平面中的每一个点代表一个质心侧偏角状态值和一个横摆角速度状态值的组合,当给定不同初始状态值时,系统就会随着时间的推移不断发生变化,对应的系统状态就会在相平面图上绘制出一条曲线,不同的初始条件会产生不同的曲线,多条曲线构成了相轨迹簇。

第一组：路面附着系数μ = 0.85和车辆速度V = 10m/s,δ = 0deg、δ = 5deg和δ = 15deg三种不同前轮转角下的相图

问题一：要想使用工作变量区中输出的变量（通常存在out里），需要在写脚本的时候先把他提出来。

%把速度从out里取出来

vel=out.y；%速度数值

time=out.t%时间

注意：注意在计算中设计的变量类型：
在新版 MATLAB/Simulink 中，out.xxx 默认导出的是 timeseries（时间序列）对象，而不是单纯的数值矩阵或标量。当你把一个 timeseries 对象直接塞进公式进行加减乘除时，结果依然是带有时间属性的复杂对象，而不是单纯的数字。最后这个奇怪的对象被送进了 atan 函数，atan 无法识别它，于是果断报错。

问题二：使用了现成的相平面绘制模型

他的模型相平面程序内嵌套了多个函数，较为复杂，迁移成本大，所以自己进行写相平面程序。

问题三;已经能跑出来了，但还存在很多问题，比如不流畅，线条不均匀不细致，找到平衡点或者鞍点的位置，有待解决。

2026.3.13 版本一

2026.3.16 版本二：

625次方针，存在些许毛刺，不够平顺，附代码如下：
对照正确代码：

% Simulink 单工况相平面高密度绘制脚本
clear; clc; close all;

%% 1. 车辆模型参数 (Vehicle Parameters)
m  = 1600;      % 整车质量 (kg)
Iz = 2800;      % 横摆转动惯量 (kg·m^2)
a  = 1.17;      % 质心到前轴距离 (m)
b  = 1.69;      % 质心到后轴距离 (m)
L  = a + b;     % 轴距
Cf = 148600;    % 前轮侧偏刚度 (N/rad)
Cr = 97600;     % 后轮侧偏刚度 (N/rad) 
g = 9.8;
h = 0.685;
mu = 0.85;

%% 2. 定义全局工况参数 (单工况)
mdl = 'three_model';      % Simulink 模型文件名
load_system(mdl);         % 预加载模型，加快循环速度

V = 10;                   % 车速 (m/s)
sim_time = '4';           % 仿真时长 (秒)

% [修改点1]：直接指定你要的具体工况：前轮纵向力 2000N，前轮转角 0度
Fxf_val = 2000;           
delta_deg = 0;            
delta_rad_val = delta_deg * pi / 180; % 转换为弧度

%% 3. 设定初始状态网格 (加密)
% [修改点2]：将网格点数从 7 提高到 15 (或者你可以改成 20，但注意仿真时间会变长)
% 15 * 15 = 225 个初始点，画出来的线簇会非常饱满
beta_grid = linspace(-0.4, 0.4, 25);     % 侧偏角初始值范围 (稍微拓宽了一点方便看全貌)
gamma_grid = linspace(-1.5, 1.5, 25);    % 横摆角速度初始值范围 

%% 4. 开始跑 Simulink 仿真并绘图
% 创建一个较大的单图窗口
figure('Name', 'Simulink Single Phase Plane', 'Position', [100, 100, 800, 600]);
hold on; box on; grid on;

total_sims = length(beta_grid) * length(gamma_grid);
current_sim = 0;
fprintf('开始仿真，总共需要运行 %d 次...\n', total_sims);

%遍历所有初始条件网格
for b0 = beta_grid
    for r0 = gamma_grid
        % 将当前循环的参数推送到 MATLAB 基础工作区
        assignin('base', 'beta_init', b0);
        assignin('base', 'gamma_init', r0);
        assignin('base', 'F_xf_phase', Fxf_val);
        assignin('base', 'delta_phase', delta_rad_val);
        assignin('base', 'V_phase', V);
        assignin('base', 'mu', mu);
       
        % 运行 Simulink 仿真
simout = sim(mdl, 'StopTime', sim_time, 'ReturnWorkspaceOutputs', 'on');
% 提取仿真结果
beta_traj = simout.beta_out.Data;
r_traj = simout.gamma_out.Data;
% 绘制当前起点的相轨迹
plot(beta_traj, r_traj, 'b-', 'LineWidth', 0.8);
% 画个黑点标记起始位置
plot(beta_traj(1), r_traj(1), 'k.', 'MarkerSize', 5);
% 打印进度 (可选，不想看可以注释掉)
current_sim = current_sim + 1;
if mod(current_sim, 20) == 0
fprintf('已完成: %d / %d 次仿真\n', current_sim, total_sims);
end
    end
end
 % for b0 = beta_grid
 %    for r0 = gamma_grid
 %        % [注意这里：把变量名统一为你在 Simulink 里实际填写的名字]
 %        assignin('base', 'beta_init', b0);
 %        assignin('base', 'gamma_init', r0);       % 修改为 r_init
 %        assignin('base', 'Fxf_phase', Fxf_val);     % 修改为 Fxf
 %        assignin('base', 'delta_phase', delta_rad_val);
 %        assignin('base', 'V', V);             % 修改为 V
 %        assignin('base', 'mu', mu);
 % 
 %        try
 %            % 运行 Simulink 仿真 (加入超时和错误捕捉)
 %            simout = sim(mdl, 'StopTime', sim_time, 'ReturnWorkspaceOutputs', 'on');
 % 
 %            % 提取仿真结果
 %            beta_traj = simout.beta_out.Data; 
 %            r_traj = simout.gamma_out.Data;       % 确保对应 Simulink To Workspace 里的 r_out
 % 
 %            % 绘制当前起点的相轨迹
 %            plot(beta_traj, r_traj, 'b-', 'LineWidth', 0.8);
 %            plot(beta_traj(1), r_traj(1), 'k.', 'MarkerSize', 5);
 % 
 %        catch ME
 %            % 如果仿真发散崩溃，打印警告并继续下一个点，不会中断整个程序
 %            warning('初始点 (beta=%.2f, r=%.2f) 仿真发散或失败。跳过此点。', b0, r0);
 %            % 可选：画个红叉标记发散点
 %            plot(b0, r0, 'rx', 'MarkerSize', 8, 'LineWidth', 1.5);
 %        end
 % 
 %        current_sim = current_sim + 1;
 %        if mod(current_sim, 20) == 0
 %            fprintf('已完成: %d / %d 次仿真\n', current_sim, total_sims);
 %        end
 



% 绘制原点辅助线
xline(0, 'k--', 'LineWidth', 1.5);
yline(0, 'k--', 'LineWidth', 1.5);

% 设置标题和坐标轴
title(sprintf('Vehicle Phase Plane (F_{xf} = %.0f N, \\delta = %d^\\circ, V = %.1f m/s, \\mu = %.2f)', Fxf_val, delta_deg, V, mu), 'FontSize', 14);
xlabel('Sideslip Angle \beta (rad)', 'FontSize', 12, 'FontWeight', 'bold');
ylabel('Yaw Rate r (rad/s)', 'FontSize', 12, 'FontWeight', 'bold');

% 限制坐标轴显示范围，让图面更聚焦
xlim([-0.5, 0.5]); 
ylim([-1.5, 1.5]);

% 优化坐标轴显示风格
set(gca, 'FontSize', 12, 'LineWidth', 1.2);

fprintf('绘制完成！\n');

2026.3.16 版本三：满屏撒点法         

对版本二进行了求解器修正和曲线平滑
仍需解决存在莫名断电，曲线不连续等问题。

% Simulink 单工况相平面高密度绘制脚本
clear; clc; close all;

%% 1. 车辆模型参数 (Vehicle Parameters)
m  = 1600;      % 整车质量 (kg)
Iz = 2800;      % 横摆转动惯量 (kg·m^2)
a  = 1.17;      % 质心到前轴距离 (m)
b  = 1.69;      % 质心到后轴距离 (m)
L  = a + b;     % 轴距
Cf = 148600;    % 前轮侧偏刚度 (N/rad)
Cr = 97600;     % 后轮侧偏刚度 (N/rad) 
g = 9.8;
h = 0.685;
mu = 0.85;

%% 2. 定义全局工况参数 (单工况)
mdl = 'three_model';      % Simulink 模型文件名
load_system(mdl);         % 预加载模型，加快循环速度

V = 10;                   % 车速 (m/s)
sim_time = '3';           % 仿真时长 (秒)

% [修改点1]：直接指定你要的具体工况：前轮纵向力 2000N，前轮转角 0度
Fxf_val = 2000;           
delta_deg = 0;            
delta_rad_val = delta_deg * pi / 180; % 转换为弧度

%% 3. 设定初始状态网格 (加密)
% [修改点2]：将网格点数从 7 提高到 15 (或者你可以改成 20，但注意仿真时间会变长)
% 15 * 15 = 225 个初始点，画出来的线簇会非常饱满
beta_grid = linspace(-0.6, 0.6, 20);     % 侧偏角初始值范围 (稍微拓宽了一点方便看全貌)
gamma_grid = linspace(-1.8, 1.8, 20);    % 横摆角速度初始值范围 

%% 4. 开始跑 Simulink 仿真并绘图
% 创建一个较大的单图窗口
figure('Name', 'Simulink Single Phase Plane', 'Position', [100, 100, 800, 600]);
hold on; box on; grid on;

total_sims = length(beta_grid) * length(gamma_grid);
current_sim = 0;
fprintf('开始仿真，总共需要运行 %d 次...\n', total_sims);

%遍历所有初始条件网格
for b0 = beta_grid
    for r0 = gamma_grid
        % 将当前循环的参数推送到 MATLAB 基础工作区
        assignin('base', 'beta_init', b0);
        assignin('base', 'gamma_init', r0);
        assignin('base', 'F_xf_phase', Fxf_val);
        assignin('base', 'delta_phase', delta_rad_val);
        assignin('base', 'V_phase', V);
        assignin('base', 'mu', mu);
       
        % 运行 Simulink 仿真
simout = sim(mdl, 'StopTime', sim_time, 'ReturnWorkspaceOutputs', 'on');
% 提取仿真结果
beta_traj = simout.beta_out.Data;
r_traj = simout.gamma_out.Data;
% 绘制当前起点的相轨迹

plot(beta_traj, r_traj, 'b-', 'LineWidth', 0.5);

% 画个黑点标记起始位置
%plot(beta_traj(1), r_traj(1), 'k.', 'MarkerSize', 5);

% 打印进度 (可选，不想看可以注释掉)
current_sim = current_sim + 1;
if mod(current_sim, 20) == 0
fprintf('已完成: %d / %d 次仿真\n', current_sim, total_sims);
end
    end
end
 % for b0 = beta_grid
 %    for r0 = gamma_grid
 %        % [注意这里：把变量名统一为你在 Simulink 里实际填写的名字]
 %        assignin('base', 'beta_init', b0);
 %        assignin('base', 'gamma_init', r0);       % 修改为 r_init
 %        assignin('base', 'Fxf_phase', Fxf_val);     % 修改为 Fxf
 %        assignin('base', 'delta_phase', delta_rad_val);
 %        assignin('base', 'V', V);             % 修改为 V
 %        assignin('base', 'mu', mu);
 % 
 %        try
 %            % 运行 Simulink 仿真 (加入超时和错误捕捉)
 %            simout = sim(mdl, 'StopTime', sim_time, 'ReturnWorkspaceOutputs', 'on');
 % 
 %            % 提取仿真结果
 %            beta_traj = simout.beta_out.Data; 
 %            r_traj = simout.gamma_out.Data;       % 确保对应 Simulink To Workspace 里的 r_out
 % 
 %            % 绘制当前起点的相轨迹
 %            plot(beta_traj, r_traj, 'b-', 'LineWidth', 0.8);
 %            plot(beta_traj(1), r_traj(1), 'k.', 'MarkerSize', 5);
 % 
 %        catch ME
 %            % 如果仿真发散崩溃，打印警告并继续下一个点，不会中断整个程序
 %            warning('初始点 (beta=%.2f, r=%.2f) 仿真发散或失败。跳过此点。', b0, r0);
 %            % 可选：画个红叉标记发散点
 %            plot(b0, r0, 'rx', 'MarkerSize', 8, 'LineWidth', 1.5);
 %        end
 % 
 %        current_sim = current_sim + 1;
 %        if mod(current_sim, 20) == 0
 %            fprintf('已完成: %d / %d 次仿真\n', current_sim, total_sims);
 %        end
 



% 绘制原点辅助线
xline(0, 'k--', 'LineWidth', 1.5);
yline(0, 'k--', 'LineWidth', 1.5);

% 设置标题和坐标轴
title(sprintf('Vehicle Phase Plane (F_{xf} = %.0f N, \\delta = %d^\\circ, V = %.1f m/s, \\mu = %.2f)', Fxf_val, delta_deg, V, mu), 'FontSize', 14);
xlabel('Sideslip Angle \beta (rad)', 'FontSize', 12, 'FontWeight', 'bold');
ylabel('Yaw Rate r (rad/s)', 'FontSize', 12, 'FontWeight', 'bold');

% 限制坐标轴显示范围，让图面更聚焦
xlim([-0.5, 0.5]); 
ylim([-1.5, 1.5]);

% 优化坐标轴显示风格
set(gca, 'FontSize', 12, 'LineWidth', 1.2);

fprintf('绘制完成！\n');

版本四：边缘发车方法：

% Simulink 单工况相平面高密度绘制脚本 (边缘发车版)
clear; clc; close all;

%% 1. 车辆模型参数 (Vehicle Parameters)
m  = 1600;      % 整车质量 (kg)
Iz = 2800;      % 横摆转动惯量 (kg·m^2)
a  = 1.17;      % 质心到前轴距离 (m)
b  = 1.69;      % 质心到后轴距离 (m)
L  = a + b;     % 轴距
Cf = 148600;    % 前轮侧偏刚度 (N/rad)
Cr = 97600;     % 后轮侧偏刚度 (N/rad) 
g = 9.8;
h = 0.685;
mu = 0.85;

%% 2. 定义全局工况参数 (单工况)
mdl = 'three_model';      % Simulink 模型文件名
load_system(mdl);         % 预加载模型，加快循环速度
V = 10;                   % 车速 (m/s)

% [关键修改 1]：增加仿真时长。因为点都在边缘，跑回中心需要更多时间
sim_time = '10';          % 仿真时长从 3 秒增加到 10 秒

Fxf_val = 2000;           
delta_deg = 0;            
delta_rad_val = delta_deg * pi / 180; % 转换为弧度

%% 3. 设定初始状态网格 (边缘发车策略)
% [关键修改 2]：不再使用全覆盖网格，而是只在相平面的四个外边缘撒点
b_max = 0.6;  b_min = -0.6;  % 侧偏角边界
r_max = 1.8;  r_min = -1.8;  % 横摆角速度边界
N_edge = 125;                 % 每条边上的起点数量 (20*4 = 80 个点，跑起来比原来 400 个点快很多，但图更漂亮)

% 构造四条边的初始点坐标 [beta, r]
edge_top    = [linspace(b_min, b_max, N_edge)', repmat(r_max, N_edge, 1)];
edge_bottom = [linspace(b_min, b_max, N_edge)', repmat(r_min, N_edge, 1)];
edge_left   = [repmat(b_min, N_edge, 1), linspace(r_min, r_max, N_edge)'];
edge_right  = [repmat(b_max, N_edge, 1), linspace(r_min, r_max, N_edge)'];

% 合并为一个总的初始点矩阵
initial_points = [edge_top; edge_bottom; edge_left; edge_right];
total_sims = size(initial_points, 1);

%% 4. 开始跑 Simulink 仿真并绘图
% 创建一个较大的单图窗口
figure('Name', 'Simulink Single Phase Plane', 'Position', [100, 100, 800, 600]);
hold on; box on; grid on;

current_sim = 0;
fprintf('开始仿真，总共需要运行 %d 次...\n', total_sims);

% [关键修改 3]：将双重 for 循环改为单层遍历矩阵的循环
for i = 1:total_sims
    b0 = initial_points(i, 1);
    r0 = initial_points(i, 2);
    
    % 将当前循环的参数推送到 MATLAB 基础工作区
    assignin('base', 'beta_init', b0);
    assignin('base', 'gamma_init', r0);
    assignin('base', 'F_xf_phase', Fxf_val);
    assignin('base', 'delta_phase', delta_rad_val);
    assignin('base', 'V_phase', V);
    assignin('base', 'mu', mu);
   
    try
        % 运行 Simulink 仿真 (这里我帮你把 try-catch 加回去了，防患于未然)
        simout = sim(mdl, 'StopTime', sim_time, 'ReturnWorkspaceOutputs', 'on');
        
        % 提取仿真结果
        beta_traj = simout.beta_out.Data;
        r_traj = simout.gamma_out.Data;
        
        % 绘制当前起点的相轨迹
        plot(beta_traj, r_traj, 'b-', 'LineWidth', 0.5);
        
    catch ME
        warning('边缘初始点 (beta=%.2f, r=%.2f) 仿真发散或失败。跳过。', b0, r0);
    end
    
    % 打印进度
    current_sim = current_sim + 1;
    if mod(current_sim, 20) == 0
        fprintf('已完成: %d / %d 次仿真\n', current_sim, total_sims);
    end
end
 
%% 5. 图像美化与收尾
% 绘制原点辅助线
xline(0, 'k--', 'LineWidth', 1.5);
yline(0, 'k--', 'LineWidth', 1.5);

% 设置标题和坐标轴
title(sprintf('Vehicle Phase Plane (F_{xf} = %.0f N, \\delta = %d^\\circ, V = %.1f m/s, \\mu = %.2f)', Fxf_val, delta_deg, V, mu), 'FontSize', 14);
xlabel('Sideslip Angle \beta (rad)', 'FontSize', 12, 'FontWeight', 'bold');
ylabel('Yaw Rate r (rad/s)', 'FontSize', 12, 'FontWeight', 'bold');

% 限制坐标轴显示范围，让图面更聚焦
xlim([-0.5, 0.5]); 
ylim([-1.5, 1.5]);

% 优化坐标轴显示风格
set(gca, 'FontSize', 12, 'LineWidth', 1.2);
fprintf('绘制完成！\n');

版本五：问题出在了轮胎模型上：
 

1. 致命的数学 Typo（导致轨迹断裂、杂糅的元凶）

请仔细看你代码里计算 term3 的这一行：

Matlab

term3 = -(Ca^2 / (27 * xi_safe^2 * mu^2 * Fz_safe^2)) * (zi^3); 

这里少乘了一个 $C_\alpha$！分子应该是 Ca^3，而不是 Ca^2。

这可不是一个无足轻重的小笔误。标准 Fiala 轮胎模型的完整公式是：

$F_y = -C_\alpha \tan\alpha + \frac{C_\alpha^2}{3 \xi \mu F_z} |\tan\alpha| \tan\alpha - \frac{C_\alpha^3}{27 \xi^2 \mu^2 F_z^2} \tan^3\alpha$

你代码里的错误会导致什么物理灾难？

在即将进入完全滑动区（abs(zi) 刚好等于 tan_alpha_sli）的临界点：

• 你的 term1 算出来是 $-3 \xi \mu F_z$

• 你的 term2 算出来是 $+3 \xi \mu F_z$

• 这两项刚好完全抵消成了 0！

• 因为你 term3 写成了 Ca^2，导致它算出来变成了一个极其微小的数值（大约是 $-\xi \mu F_z / C_\alpha$）。

也就是说，在侧偏角增大的过程中，轮胎力本来应该平滑地上升到 $-\xi \mu F_z$，但因为这个 Typo，轮胎力在到达边界的前一瞬间，突然掉到了接近 0！ 然后下一秒进入 else 逻辑，瞬间又突变回了最大摩擦力 $-\xi \mu F_z$。

这种高达几千牛顿的悬崖式突变，在相平面上就是一颗“炸弹”，求解器经过这里时瞬间“骨折”，所以你会看到乱七八糟的杂糅和尖锐的拐角。

2. 底层物理缺陷：分段函数与“硬饱和”

即使你把 Ca^2 改成了 Ca^3，修复了突变 Bug，标准的 Fiala 模型依然画不出文献里那种像丝绸一样的完美曲线。

因为 Fiala 模型是一个分段函数 (if/else)。

当侧偏角越过临界点后，侧向力直接变成了一个常数：-xi_safe * mu * Fz_safe * sign(alpha)。

这意味着在滑动区，侧向力对侧偏角的导数（变化率）瞬间变成了绝对的 0。

在微分方程和相平面拓扑学中，这种导数的不连续性（$C^1$ 不连续）必然导致相轨迹出现肉眼可见的“硬折角”。

---

🚀 如何修改？

方案：使用双曲正切拟合（复刻文献完美图的终极秘籍）

文献中画相平面，为了保证向量场绝对光滑，99% 都会用 $\tanh$（双曲正切）函数来代替分段模型。它能在 $0$ 点保持线性刚度 $C_\alpha$，并在远端无限平滑地趋近于附着极限，没有任何 if/else 的硬拐角。

你可以直接用下面这几行极其优雅的代码，整体替换掉你原来的 if/else 计算：(但这个是近视了一个轮胎模型并不是原来的fiala轮胎）

Matlab

function Fyi = improved_fiala_tire(alpha, Fz, Fx, mu, Ca)
    % 1. 安全保护
    Fz_safe = max(abs(Fz), 0.01);
    
    % 2. 纵向力修正系数 xi
    mu_Fz = mu * Fz_safe;
    inside_sqrt = max((mu_Fz)^2 - Fx^2, 0); 
    xi = sqrt(inside_sqrt) / mu_Fz;
    xi_safe = max(xi, 1e-5);
    
    % 3. 极限侧向力
    Fy_max = xi_safe * mu_Fz;
    
    % 4. 使用 tanh 函数进行全局无限平滑拟合 (完美替代 if/else 分段)
    % 公式: Fy = -Fy_max * tanh( (Ca * tan(alpha)) / Fy_max )
    Fyi = -Fy_max * tanh( (Ca * tan(alpha)) / Fy_max );
end

平滑后的fiala轮胎模型：
 

function Fyi = improved_fiala_tire(alpha, Fz, Fx, mu, Ca)
% #########################################################################
% 修复版：改进 Fiala 轮胎模型 (原汁原味，修复了导致相平面崩溃的 Bug)
% #########################################################################

    % 1. 垂直载荷保护
    Fz_safe = max(abs(Fz), 0.01);
    
    % 2. 计算摩擦圆极值与修正因子
    mu_Fz = mu * Fz_safe;
    inside_sqrt = max((mu_Fz)^2 - Fx^2, 0); 
    xi = sqrt(inside_sqrt) / mu_Fz;
    xi_safe = max(xi, 1e-5);
    
    % 3. ★ 核心保护 1：防止 tan(alpha) 走向无穷大 (Inf)
    % 在相平面大打滑工况下，限制 alpha 绝对值不超过 85 度 (约 1.48 rad)
    % 这样既保留了真实的 tan 几何关系，又绝对不会让求解器崩溃
    alpha_limit = min(max(alpha, -1.48), 1.48);
    zi = tan(alpha_limit);
    
    % 4. 计算滑动临界点的正切值
    tan_alpha_sli = (3 * xi_safe * mu_Fz) / Ca;
    
    % 5. 分段逻辑计算侧向力 Fyi
    if abs(zi) < tan_alpha_sli
        % --- 弹性区/半滑动区 ---
        term1 = -Ca * zi;
        term2 = (Ca^2 / (3 * xi_safe * mu_Fz)) * abs(zi) * zi;
        % ★ 核心保护 2：修复了 Ca^3 的 Typo！
        % 这一改动恢复了公式在边界处的平滑导数，消除了相平面的生硬折线
        term3 = -(Ca^3 / (27 * xi_safe^2 * mu_Fz^2)) * (zi^3); 
        Fyi = term1 + term2 + term3;
    else
        % --- 完全滑动区 ---
        % 在滑动区直接输出极限力，保持原来的逻辑不变
        Fyi = -xi_safe * mu_Fz * sign(alpha_limit);
    end
    
end

完整相平面代码：

% Simulink 单工况相平面高密度绘制脚本 (纯粹边缘发车版)
clear; clc; close all;

%% 1. 车辆模型参数 (Vehicle Parameters)
m  = 1600;      % 整车质量 (kg)
Iz = 2800;      % 横摆转动惯量 (kg·m^2)
a  = 1.17;      % 质心到前轴距离 (m)
b  = 1.69;      % 质心到后轴距离 (m)
L  = a + b;     % 轴距
Cf = 148600;    % 前轮侧偏刚度 (N/rad)
Cr = 97600;     % 后轮侧偏刚度 (N/rad) 
g = 9.8;
h = 0.685;
mu = 0.85;

%% 2. 定义全局工况参数 (单工况)
mdl = 'three_model';      % Simulink 模型文件名
load_system(mdl);         % 预加载模型，加快循环速度
V = 10;                   % 车速 (m/s)

% 仿真时长设为 20 秒，确保边缘出发的长线能跑到终点
sim_time = '20';          

Fxf_val = 2000;           
delta_deg = 0;            
delta_rad_val = delta_deg * pi / 180; % 转换为弧度

%% 3. 设定初始状态网格 (纯粹的边缘发车策略)
b_max = 0.6;  b_min = -0.6;  % 侧偏角边界
r_max = 1.8;  r_min = -1.8;  % 横摆角速度边界

% 纯边缘发车：为了保证画面饱满，将边缘点数加密至 150 个
N_edge = 100;                 
edge_top    = [linspace(b_min, b_max, N_edge)', repmat(r_max, N_edge, 1)];
edge_bottom = [linspace(b_min, b_max, N_edge)', repmat(r_min, N_edge, 1)];
edge_left   = [repmat(b_min, N_edge, 1), linspace(r_min, r_max, N_edge)'];
edge_right  = [repmat(b_max, N_edge, 1), linspace(r_min, r_max, N_edge)'];

% 合并总点集并去重 (去除四个角点重复，不再添加内部点)
initial_points = [edge_top; edge_bottom; edge_left; edge_right];
initial_points = unique(initial_points, 'rows');
total_sims = size(initial_points, 1);

%% 4. 开始跑 Simulink 仿真并绘图
% 创建一个较大的单图窗口，设为纯白背景适合放进论文
figure('Name', 'Simulink Pure Edge Phase Plane', 'Position', [100, 100, 800, 600], 'Color', 'w');
hold on; box on; grid on;

current_sim = 0;
fprintf('开始仿真，总共需要运行 %d 次...\n', total_sims);

for i = 1:total_sims
    b0 = initial_points(i, 1);
    r0 = initial_points(i, 2);
    
    % 将当前循环的参数推送到 MATLAB 基础工作区
    assignin('base', 'beta_init', b0);
    assignin('base', 'gamma_init', r0);
    assignin('base', 'F_xf_phase', Fxf_val);
    assignin('base', 'delta_phase', delta_rad_val);
    assignin('base', 'V_phase', V);
    assignin('base', 'mu', mu);
   
    try
        % 强制限制最大步长 MaxStep 为 0.005，保证线条平滑
        simout = sim(mdl, 'StopTime', sim_time, 'MaxStep', '0.005', 'ReturnWorkspaceOutputs', 'on');
        
        % 提取仿真结果
        beta_traj = simout.beta_out.Data;
        r_traj = simout.gamma_out.Data;
        
        % 文献级虚实线及颜色判定
        if abs(beta_traj(end)) > 0.5 || abs(r_traj(end)) > 1.6
            % 发散失控轨迹：浅蓝色 + 虚线
            plot(beta_traj, r_traj, '--', 'Color', [0.3010 0.7450 0.9330], 'LineWidth', 0.8);
        else
            % 收敛稳定轨迹：深蓝色 + 实线
            plot(beta_traj, r_traj, '-', 'Color', [0 0.4470 0.7410], 'LineWidth', 1.0);
        end
        
    catch ME
        % 屏蔽单点报错，防止整个图画废
        warning('初始点 (beta=%.2f, r=%.2f) 仿真发散或失败。跳过。', b0, r0);
    end
    
    % 打印进度
    current_sim = current_sim + 1;
    if mod(current_sim, 20) == 0
        fprintf('已完成: %d / %d 次仿真\n', current_sim, total_sims);
    end
end
 
%% 5. 图像美化与收尾
% 绘制原点辅助线
xline(0, 'k--', 'LineWidth', 1.5);
yline(0, 'k--', 'LineWidth', 1.5);

% 设置标题和坐标轴
title(sprintf('Vehicle Phase Plane (F_{xf} = %.0f N, \\delta = %d^\\circ, V = %.1f m/s)', Fxf_val, delta_deg, V), 'FontSize', 14);
xlabel('Sideslip Angle \beta (rad)', 'FontSize', 12, 'FontWeight', 'bold');
ylabel('Yaw Rate r (rad/s)', 'FontSize', 12, 'FontWeight', 'bold');

% 限制坐标轴显示范围，切除外面乱飞的线条，只展示核心区域
xlim([-0.5, 0.5]); 
ylim([-1.5, 1.5]);

% 优化坐标轴显示风格
set(gca, 'FontSize', 12, 'LineWidth', 1.2);
fprintf('纯边缘发车绘制完成！\n');

最终结果：
 

Fxf_val = 2000;

delta_deg = 0;

Fxf_val = 0;

delta_deg = 0;

Fxf_val = -2000;

delta_deg = 0;

Fxf_val =2000;

delta_deg = 5;

Fxf_val = 0;

delta_deg = 5;

Fxf_val = -2000;

delta_deg =5;

Fxf_val = -2000;

delta_deg = 15;

Fxf_val = 0;

delta_deg = 15;

Fxf_val = -2000;

delta_deg = 15;