相平面理解直观的比喻碗里的玻璃球想象一个底部平滑的碗这就好比你固定的车速、路面和转角它们决定了碗的形状。稳态如果玻璃球正好在碗底它就不动了。这就像你认为的“状态应该是固定的”。相轨迹动态演变如果我把玻璃球拿到碗壁的某个高处这就是给了它一个初始的侧偏角和横摆角速度然后松手。玻璃球会在碗里来回滚动它的位置和速度时刻在变化最终慢慢停在碗底。相平面图就是你把玻璃球放在碗壁上成百上千个不同的位置松手然后用相机把球滚动的所有的轨迹全部重叠拍在同一张照片上。怎么去实现呢如何联合simulink模型绘制相平面第一步改造你的 Simulink 模型关键为了让 MATLAB 脚本能控制 Simulink你需要让 Simulink 模型里的一些固定数值变成工作区变量。设置工况输入设为变量在 Simulink 中你的前轮转角、车速、纵向力、附着系数不要直接写死成数字。用Constant常数模块但里面的值分别填上变量名delta_phase、V、Fxf、mu。设置初始状态核心步骤找到你计算质心侧偏角 $\beta$ 和横摆角速度 $r$ 的那两个积分模块Integrator, 图标是1/s。双击它们把里面的Initial condition初始条件从 0 改成变量名$\beta$ 的积分模块初始条件设为beta_init$r$ 的积分模块初始条件设为r_init设置输出提取To Workspace在 $\beta$ 和 $r$ 的信号线上分别接一个To Workspace模块。变量名分别设为beta_out和r_out。Save format保存格式强烈建议选择Timeseries或Array下面的脚本以 Timeseries 为例。不同工况设置将Fxf分别设置为-2kN、0kN和2kN,其中Fxf-2kN表示前轮制动工况的漂移Fxf0kN表示后轮驱动车辆的漂移Fxf2kN表示前轮驱动工况的漂移。用质心侧偏角β和横摆角速度r表征整个车辆系统,并将这两个状态量作为横纵坐标,建立相平面图,相平面中的每一个点代表一个质心侧偏角状态值和一个横摆角速度状态值的组合,当给定不同初始状态值时,系统就会随着时间的推移不断发生变化,对应的系统状态就会在相平面图上绘制出一条曲线,不同的初始条件会产生不同的曲线,多条曲线构成了相轨迹簇。第一组路面附着系数μ 0.85和车辆速度V 10m/s,δ 0deg、δ 5deg和δ 15deg三种不同前轮转角下的相图问题一要想使用工作变量区中输出的变量通常存在out里需要在写脚本的时候先把他提出来。%把速度从out里取出来velout.y%速度数值timeout.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的这一行Matlabterm3 -(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轮胎Matlabfunction 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;
Simulink仿真漂移机理分析（二）：相图分析
相平面理解直观的比喻碗里的玻璃球想象一个底部平滑的碗这就好比你固定的车速、路面和转角它们决定了碗的形状。稳态如果玻璃球正好在碗底它就不动了。这就像你认为的“状态应该是固定的”。相轨迹动态演变如果我把玻璃球拿到碗壁的某个高处这就是给了它一个初始的侧偏角和横摆角速度然后松手。玻璃球会在碗里来回滚动它的位置和速度时刻在变化最终慢慢停在碗底。相平面图就是你把玻璃球放在碗壁上成百上千个不同的位置松手然后用相机把球滚动的所有的轨迹全部重叠拍在同一张照片上。怎么去实现呢如何联合simulink模型绘制相平面第一步改造你的 Simulink 模型关键为了让 MATLAB 脚本能控制 Simulink你需要让 Simulink 模型里的一些固定数值变成工作区变量。设置工况输入设为变量在 Simulink 中你的前轮转角、车速、纵向力、附着系数不要直接写死成数字。用Constant常数模块但里面的值分别填上变量名delta_phase、V、Fxf、mu。设置初始状态核心步骤找到你计算质心侧偏角 $\beta$ 和横摆角速度 $r$ 的那两个积分模块Integrator, 图标是1/s。双击它们把里面的Initial condition初始条件从 0 改成变量名$\beta$ 的积分模块初始条件设为beta_init$r$ 的积分模块初始条件设为r_init设置输出提取To Workspace在 $\beta$ 和 $r$ 的信号线上分别接一个To Workspace模块。变量名分别设为beta_out和r_out。Save format保存格式强烈建议选择Timeseries或Array下面的脚本以 Timeseries 为例。不同工况设置将Fxf分别设置为-2kN、0kN和2kN,其中Fxf-2kN表示前轮制动工况的漂移Fxf0kN表示后轮驱动车辆的漂移Fxf2kN表示前轮驱动工况的漂移。用质心侧偏角β和横摆角速度r表征整个车辆系统,并将这两个状态量作为横纵坐标,建立相平面图,相平面中的每一个点代表一个质心侧偏角状态值和一个横摆角速度状态值的组合,当给定不同初始状态值时,系统就会随着时间的推移不断发生变化,对应的系统状态就会在相平面图上绘制出一条曲线,不同的初始条件会产生不同的曲线,多条曲线构成了相轨迹簇。第一组路面附着系数μ 0.85和车辆速度V 10m/s,δ 0deg、δ 5deg和δ 15deg三种不同前轮转角下的相图问题一要想使用工作变量区中输出的变量通常存在out里需要在写脚本的时候先把他提出来。%把速度从out里取出来velout.y%速度数值timeout.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的这一行Matlabterm3 -(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轮胎Matlabfunction 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;
