简介这是一份基于MATLAB的自行车运动学模型仿真源码面向车辆动力学学习者、MATLAB使用者和智能车/自动驾驶方向初学者用于实现自行车左转、右转的动力学建模、仿真与可视化。资源共4个文件均为m脚本包含主程序Runme.m与功能函数整包仅4KB结构紧凑适合逐行研读和二次开发。目前已有292人学习下载。源码完整覆盖了自行车动力学基础、数学模型建立、状态空间表示、Simulink/Stateflow仿真环境搭建、转向控制器设计、可视化输出与源码调试优化等环节读者可通过修改轮距、前叉角度等关键参数直观观察转向稳定性差异并结合轨迹图、速度图深入理解动态行为。同时可掌握MATLAB函数封装、数据处理、运行调试等编码技巧为今后开展自动驾驶稳定性分析、骑行安全研究及控制策略优化提供可复用的实验基础。1. 自行车运动学模型的MATLAB仿真从左右转开始最靠谱做车辆路径规划和移动机器人避障的人几乎都绕不过自行车运动学模型但要亲手用MATLAB把它跑起来最容易卡住的其实不是方程本身而是坐标系、正负号和转向输入的配合。经常有人把delta设为正值结果轨迹往右偏于是怀疑公式抄错。实际上只要统一用右手系约定x向右、y向上、航向角逆时针为正那么delta大于0就是左转小于0就是右转。这篇文章围绕一个能直接运行的源码展开先给出状态方程和最小函数再做左右转的仿真脚本最后讲参数边界和验证技巧适合用MATLAB做毕设和原型验证的人。2. 自行车运动学模型的状态方程与MATLAB坐标约定2.1 左转右转的正负号怎么定先定坐标系。MATLAB的plot命令默认以屏幕横向为x轴、纵向为y轴且x向右、y向上因此theta从x轴正方向开始逆时针方向为正。在这个约定下前轮转角delta取正值表示前轮向左偏折车辆沿着左侧圆弧行驶也就是左转delta取负值表示前轮向右偏折车辆右转。很多现成源码把左右转方向写反大多是因为没有统一“正方向”。运动学模型的标准形式是dx/dt v * cos(theta) dy/dt v * sin(theta) dtheta/dt v * tan(delta) / L其中v是后轴中心的平动速度L是前后轮轴距delta是前轮转角。这个模型把车辆看成刚性连接前后轮的单车也叫单轨模型。需要强调x、y方向的推进只取决于车辆当前的航向角thetadelta并不直接出现在dx、dy里而是通过dtheta改变航向角再影响下一时刻的x、y方向。这就是自行车运动学模型和普通质点模型的本质区别。2.2 参考点选在后轴中心而不是前轴这里的常见陷阱是“参考点”。上面的方程对应的是后轴中心轨迹车辆速度方向与车身纵轴一致所以速度分解干净利落。如果参考点换到前轴中心方程就会变成dx/dt v * cos(theta delta) dy/dt v * sin(theta delta) dtheta/dt v * sin(delta) / L看着更“对称”但注意此时前轮速度方向与车身夹角为delta轮胎被要求无侧偏得到的转弯半径和转向角关系不再单调直观。在绝大多数路径规划和自动驾驶的横向控制文献里默认使用的是后轴中心模型。原因有三一是低速泊车和窄路掉头场景下后轴中心能更好描述车体“尾巴”的运动二是转向角delta到横摆角速度的映射在正负方向上完全对称调试左右转时逻辑简单三是后续做Stanley、纯追踪等横向控制器时误差量定义在后轴中心最自然。所以本文所有代码均采用后轴中心版本。2.3 模型参数表与最小转弯半径这里给出常用的一组参数方便做第一个仿真参数符号示例值单位说明轴距L1.2m前后轮触地点的距离自行车L较短车速v2.0m/s后轴中心的平动速度可设为常数最大前轮转角delta_max30deg前轮摆角物理限位需要转为rad初始航向角theta00rad车头方向与x轴夹角仿真步长dt0.01s欧拉积分的时间间隔根据dtheta v*tan(delta)/L最小转弯半径与delta_max的关系是R_min L / tan(delta_max)。在MATLAB里计算L 1.2; delta_max deg2rad(30); R_min L / tan(delta_max)运行结果是2.0784 m意味着这组参数下车辆以满舵30度持续行驶会画出一个半径约2.08 m的圆。这个数值后面用来验证左转/右转轨迹是否正确。注意delta趋近90度时tan(delta)发散运动学模型失效所以delta_max一般不超过60度。物理上自行车的转向角也很少超过45度。3. 用MATLAB写最小仿真源码欧拉积分驱动左右转3.1 把状态递推封装成一个独立函数先把运动学模型写成一个可复用的MATLAB函数文件bicycle_step.m。输入当前状态s[x;y;theta]加上速度v、前轮转角delta、轴距L和步长dt输出下一时刻状态s_new。function s_new bicycle_step(s, v, delta, L, dt) % bicycle_step 按自行车运动学模型向前推进一步 % 输入: % s - 3x1状态向量 [x; y; theta]theta单位rad % v - 后轴中心速度单位m/s % delta - 前轮转角单位rad正值左转负值右转 % L - 轴距单位m % dt - 积分步长单位s % 输出: % s_new - 更新后的状态向量 s_new zeros(3,1); s_new(1) s(1) v * cos(s(3)) * dt; s_new(2) s(2) v * sin(s(3)) * dt; s_new(3) s(3) (v / L) * tan(delta) * dt; % 将航向角归一化到 [-pi, pi) s_new(3) atan2(sin(s_new(3)), cos(s_new(3))); end代码只有四行核心运算但这里有三个容易写错的地方。第一theta是弧度制调用前必须通过deg2rad把30度转成0.5236 rad否则tan和cos输入的是“带度数的数”轨迹会完全错乱。第二航向角归一化没有用wrapToPi因为wrapToPi属于Mapping Toolbox不保证所有MATLAB环境都有atan2(sin,cos)是纯基础函数任何版本都能用。第三s_new(3)必须在更新x、y之后做角度归一化因为当前步的航向角已经参与了x、y推进归一化只影响下一步的cos/sin取值。3.2 主脚本先左转30度再右转30度有了状态递推函数主脚本只需要定义转向角随时间变化的序列。下面这段源码前2秒令delta30度左转2到4秒回正直行4到6秒令delta-30度右转6到8秒再回正。% simulate_turn.m L 1.2; % 轴距单位m v 2.0; % 车速单位m/s dt 0.01; % 仿真步长单位s T 8; % 总时长单位s t 0:dt:T; delta zeros(size(t)); delta(t 2) deg2rad(30); delta(t 2 t 4) 0; delta(t 4 t 6) deg2rad(-30); delta(t 6 t 8) 0; s zeros(3, numel(t)); s(:,1) [0; 0; 0]; % 初始位置(0,0)航向角0 for i 1:numel(t)-1 s(:, i1) bicycle_step(s(:,i), v, delta(i), L, dt); end figure; plot(s(1,:), s(2,:), b-, LineWidth, 1.6); axis equal; grid on; xlabel(x (m)); ylabel(y (m)); title(Bicycle kinematics: left turn then right turn);delta数组的长度是numel(t)循环只取delta(i)去推进第i步所以最后一个时间点的delta不会被用到这是正常现象。axis equal会让x、y方向比例一致否则画出来的圆弧会被拉伸成椭圆左右转的对称性就看不出来。左转段使用delta0右转段使用delta0如果用自己的代码时发现轨迹方向相反只需要把所有delta取负号即可。3.3 从轨迹上验证左右转是否成功运行脚本后轨迹看起来是一条先逆时针弯、后顺时针弯的S形曲线。可以用数值来确认左转段持续2秒横摆角速度omega vtan(delta)/L 2tan(30°)/1.2 ≈ 0.962 rad/s2秒后航向角理论变化量是1.924 rad右转段用同样的-30度变化量是-1.924 rad。因此第6秒车辆航向角应回到接近0左右转的圆弧部分中心角相同半径相同关于x轴大致对称。在MATLAB命令窗口里可以检查theta_at_2s s(3, find(t2)) theta_at_6s s(3, find(t6))如果第二个值接近第一个值的相反数说明左右转的符号逻辑没问题。注意这里因为中途有2秒直行航向角在直行段保持不变所以theta_at_2s约等于1.91 radtheta_at_6s约等于-0.01 rad。如果输出结果里theta_at_6s仍然接近1.91说明转向输入没有在4秒后生效通常是数组索引或时间判断写错了。4. 左右转仿真的参数调节与可视化增强4.1 用斜坡转向代替阶跃左右转更接近真实车辆第3章的delta序列是阶跃信号从前轮转角0直接跳到30度在实际车辆上不可能瞬间完成。MATLAB仿真虽然不会报错但轨迹在切换瞬间会出现明显“尖角”不符合运动规律。常见做法是给前轮转角加一个过渡斜坡例如用1秒时间完成0到30度的转向。将第3章源码里的delta生成部分替换为下面这段即可。t_ramp 1; delta zeros(size(t)); % 0~1s 从0升到30度 delta(t 1) deg2rad(30) * (t(t 1) / t_ramp); % 1~2s 保持30度 delta(t 1 t 2) deg2rad(30); % 2~3s 从30度回正 delta(t 2 t 3) deg2rad(30) * (1 - (t(t 2 t 3) - 2) / t_ramp); % 3~4s 直行 delta(t 3 t 4) 0; % 4~5s 从0转到-30度右转 delta(t 4 t 5) -deg2rad(30) * ((t(t 4 t 5) - 4) / t_ramp); % 5~6s 保持-30度 delta(t 5 t 6) -deg2rad(30); % 6~7s 从-30度回正 delta(t 6 t 7) -deg2rad(30) * (1 - (t(t 6 t 7) - 6) / t_ramp); % 7~8s 直行 delta(t 7 t 8) 0;t_ramp是转向过渡时间单位s。t_ramp越大转向动作越慢轨迹的圆弧半径在开始阶段并不是常数而是从半径无穷大逐渐收敛到L/tan(delta_max)t_ramp越小越接近阶跃响应。用这段斜坡替换阶跃后左转段从0到30度需要1秒完成相当于前轮“拧”过去而不是“跳”过去。如果要描述液压转向或机械转向的执行机构延迟还可以再加一阶惯性环节但斜坡输入已经能覆盖大部分运动学仿真需求。4.2 画上航向箭头让左右转方向一目了然仿真轨迹本身只能看出弯曲方向看不出车头朝向。MATLAB的quiver函数可以在轨迹上按固定间隔画小箭头每个箭头指向该时刻theta的方向。idx 1:20:numel(t); figure; plot(s(1,:), s(2,:), b, LineWidth, 1.4); hold on; quiver(s(1,idx), s(2,idx), cos(s(3,idx)), sin(s(3,idx)), 0.3, Color, [0.85 0.33 0.10]); axis equal; grid on;quiver调用中的第5个参数0.3是箭头缩放因子数值太大箭头会互相重叠太小则看不清方向。搭配grid on可以看到左转段箭头逐渐向逆时针方向偏转右转段箭头逐渐向顺时针方向偏转。如果还想看到车辆轮廓可以在每个采样时刻用两个点分别画出前轮和后轮位置。后轮中心就是s(:,idx)前轮中心可由后轮中心沿theta方向平移L得到x_front s(1,idx) L * cos(s(3,idx)); y_front s(2,idx) L * sin(s(3,idx));用plot把前后轮连成一条短线就能看到一个简化自行车在轨迹上“走”对演示左转右转非常有帮助。4.3 步长、车速与转角的约束关系欧拉积分是显式格式每一步的截断误差会累积。在自行车运动学模型里关键量是每个步长内航向角的增量dtheta (v*tan(delta)/L)*dt。下面表格给出不同车速下的建议步长上限。车速 v (m/s)delta_max30° 时横摆角速度 (rad/s)建议最大dt (s)2.00.9620.055.02.4050.028.03.8480.01粗略约束是每步航向角变化不超过0.1 rad约5.7度所以dt 0.1 / omega。以v8 m/s为例omega3.848 rad/sdt需小于0.026 s10 ms步长是安全的选择。如果只跑原型验证10ms已经够用如果要做多车轨迹对比或嵌入式移植步长需要进一步缩小或者改用ode45做自适应积分。还要注意当delta接近90度时tan值快速上升同样的dt误差会更大因此仿真代码里最好限制delta_max不超过60度否则数值结果会失去物理意义。5. 左右转仿真的3个验证技巧欧拉误差、镜像对称和航向角归一化5.1 用ode45交叉验证欧拉积分把同一个运动学模型交给MATLAB内置的ode45重新积分再和欧拉结果比较。先定义方程函数function ds bicycle_ode(t, s, v, L, delta_fun) d delta_fun(t); ds [v*cos(s(3)); v*sin(s(3)); v/L*tan(d)]; end调用时传入匿名函数delta_fun返回该时刻的转角delta_fun (t) deg2rad(30); [t45, s45] ode45((t,s) bicycle_ode(t,s,1.2,2.0,delta_fun), [0 2], [0 0 0]); s_euler s(1,1:201); % 欧拉法前2秒dt0.01 max(abs(s_euler - s45(:,1)))因为欧拉误差与dt同阶dt0.01时终点误差通常在毫米量级。若偏差过大先查是theta单位没转换还是L符号有误。5.2 把左转轨迹做镜像检查左右转对称分别用delta30度和-30度跑相同时间右转轨迹的y坐标取负后应与左转轨迹完全重合。在脚本末尾加上断言assert(max(abs(s_left(2,:) s_right(2,:))) 1e-6, 左右转不对称);如果运行时不报错说明正负号约定一致如果报错大概率是delta或theta的方向混了。5.3 航向角归一化与unwrap配合使用代码里用atan2把theta限制在[-pi,pi)长时间左转会让航向角从pi跳变到-pi画theta-time曲线时会出现锯齿。分析累计转角时需要先把角度展开theta_unwrapped unwrap(s(3,:));unwrap会检测相邻角度差超过pi的跳变并自动加减2pi。把theta_unwrapped和s(3,:)一起画能看到原始角度在边界跳变展开后的角度单调递增或递减对判断车辆转了多少圈很有用。本文还有配套的精品资源点击获取