汽车纵向动力学补充内容(状态空间、常微分方程数值解)

常微分方程的数值解
问题描述： 问题描述：dy/dt=f(t, y), y(t0)=y0 对常微分方程的求解是基于一阶方程进行的， 一阶方程进行的 MATLAB 对常微分方程的求解是基于一阶方程进行的，通 常采用Runge-Kutta方法，所对应的MATLAB 命令为ode 常采用Runge-Kutta方法，所对应的MATLAB 命令为ode Runge 方法 的缩写），例如ode23 ），例如ode23、 （Odinary Differential Equation 的缩写），例如ode23、 ode45、ode23s、ode23tb、ode15s、ode113等 ode45、ode23s、ode23tb、ode15s、ode113等，分别用于求 解不同类型的微分方程，如刚性方程和非刚性方程等。 解不同类型的微分方程，如刚性方程和非刚性方程等。 ode23、ode45一般用于求解非时变常微分方程 一般用于求解非时变常微分方程（ ode23、ode45一般用于求解非时变常微分方程（包括非线 性方程） 性方程） ode23s、ode45tb、ode15s、ode113一般用于求解非线 ฀ ode23s、ode45tb、ode15s、ode113一般用于求解非线 性、时变常微分方程
主程序： 主程序：main.m
function dy=WheelDyn(t,x) global Mw Rd Iw Uw0 Omega0 Tb0 g Ki Kd Smin Smax S0 mu_h mu_g dy=zeros(3,1); % x(1)=Omega, x(2)=Uw, x(3)=Tb; Sb=(x(2)-x(1)*Rd)/x(2); if(Sb<=S0) mu=mu_h*Sb/S0; else mu=(mu_h-mu_g*S0)/(1-S0)-(mu_h-mu_g)*Sb/(1-S0); end Fxb=Mw*g*mu; dy(1)=(Fxb*Rd-x(3))/Iw; %dw dy(2)=-Fxb/Mw; % dUw if(Sb>Smax) dy(3)=-Kd; % dTb 减压 elseif(Sb<Smin) dy(3)=Ki; % dTb 增压 else 子程序： 子程序：WheelDyn.m dy(3)=0; % dTb 保压 end
-kd (λ≥λ0 )
µh
uw − ω rb s= uw
1 Fxb = µ Fz = µ ⋅ mg 4
µg= s sb ( sb ≤ s0 ) 0 µ = µh − µ g s0 − µh − µ g s ( s > s ) b b 0 1 − s0 1 − s0
基于1/4车辆模型的 基于 车辆模型的ABS系统设计与仿真 车辆模型的 系统设计与仿真
dω Iw = Fxb rd − Tb dt
系统动态 模型
duw Fxb mw = Fxb dt (λ < λ0 ,增压) dTb ki = dt −kd (λ ≥ λ0 ,减压)
ω r d
Fz
Tb
选择状态变量x 选择状态变量 1 = ω, x2 = uw , x3 = Tb dx1 = (Fxb *rd-x3) /Iw (1) dx2 = Fxb /mw w dx3 = ki (λ< λ0 ) (3) (2)
子程序： 子程序：VanderPol.m
function dy= VanderPol (t, y) dy(1) = y(2)*(y(2)^2-1 dy(2) = -(y(1)^2-1)* y(2)-y(1); dy=[dy(1);dy(2)];
四阶龙格——库塔数值方法 库塔数值方法 四阶龙格
在区间[x 内预估多个点上的斜率值K 在区间 i,xi+1]内预估多个点上的斜率值 1、 内预估多个点上的斜率值 K2、……Km，并用他们的加权平均数作为平 均斜率的近似值， 均斜率的近似值，能构造出具有很高精度的高 阶计算公式。经数学推导、求解， 阶计算公式。经数学推导、求解，可以得出四 阶龙格－库塔公式， 阶龙格－库塔公式，也就是在工程中应用广泛 的经典龙格－库塔算法： 的经典龙格－库塔算法： yi＋1 = yi+h*( K1+ 2*K2 +2*K3+ K4)/6 ＋ K1 =f(xi, yi) K2=f(xi+h/2, yi +h* K1 /2) K3=f(xi+h/2, yi +h*K2/2) K4 =f(xi+h, yi +h*K3)
Runge-Kutta编程 编程 方法求解
子程序： 子程序：VanderPol.m
function dy= VanderPol (t, y) dy(1) = y(2)*(y(2)^2-1 dy(2) = -(y(1)^2-1)* y(2)-y(1); dy=[dy(1);dy(2)];
主程序： 主程序：main.m
子程序： 子程序：Psub.m
derPol方程 方程: Van derPol方程: y"+(y2-1)y'+y=0
令 y1’=y2 ; y2’= y1”=-(y12-1)*y2-y1
t0=0; tn=20; y0=[0; 0.25]; h=0.001; t = t0 : h : tn; n = length (t); j = 1; for i = 1 : n t1 = t0 + h; K1 = VanderPol(t0, y0); K2 = VanderPol (t0 + h/2, y0 + h*K1/2); K3 = VanderPol (t0 + h/2, y0 + h*K2/2); K4 = VanderPol (t0 + h, y0 + h*K3); y1 = y0 + (h/6)*(K1 + 2*K2 + 2*K3 + K4); yy1(j)=y1(1); yy2(j)=y1(2); t0=t1; y0=y1; j=j+1; end t=0:h:tN; subplot(121), plot(t, yy1, t, yy2); grid subplot(122), plot(yy1), yy2); grid
1 Fxb = µ Fz = µ ⋅ mg 4
clc; clear all; global Mw Rd Iw Uw0 Omega0 Tb0 g Ki Kd Smin Smax S0 mu_h mu_g Mw=300; Rd=0.25; Iw=12; Uw0=30; Omega0=120; Tb0=600; g=9.81; Ki=4500; Kd=5000; Smin=0.18; Smax=0.22;S0=0.2; mu_h=0.8; mu_g=0.6; % x(1)=Omega, x(2)=Uw, x(3)=Tb; x0=[Omega0;Uw0;Tb0]; options=odeset(‘Events’,@events); %开启事件判断功能 开启事件判断功能 [t,x]=ode45(@WheelDyn,[0:0.01:5],x0,options); Sb=(x(:,2)-x(:,1)*Rd)./x(:,2); figure(1); plot(t,x(:,2),'r'); hold on; plot(t,x(:,1)*Rd); axis([0 5 0 30]); figure(2); plot(t,Sb); axis([0 5 0 1.0]);
events事件子程序：events.m 事件子程序： 事件子程序 function [value,isterminal,direction] = events(t,x) value = [x(1);x(2)]; % Detect height = 0 isterminal = [1;1]; % Stop the integration direction = [0;0]; % Negative direction only
Matlab环境 环境ODE命令方法求解 环境 命令方法求解
[t, y]=ode45(‘funname’, tspan, y0, options, p1,p2….); Example: 主程序： 主程序：main.m
t0 = 0; tn= 20; tol= 1e-6; y0 = [0; 0.25]; [t, y]=ode45(‘VanderPol’, [t0, tn], y0, tol);
y(t0)=y0， 对一阶常微分方程 dy/dt=f(t, y), y(t0)=y0， 四阶Runge Kutta方法求解流程 Runge四阶Runge-Kutta方法求解流程
t0=?; y0=?; h=?; N=? % 输入初始条件、计算步长和迭代次数 输入初始条件、 tt=[t0]; yy=[y0]; % 输出 1, y1 输出t for i = 1 : N t1 = t0 + h; K1 = f (t0，y0) = Psub(t0, y0); K2 = f (t0 + h/2, y0 + hK1/2) = Psub (t0 + h/2, y0 + h*K1/2); K3 = f (t0 + h/2, y0 + hK2/2) = Psub (t0 + h/2, y0 + h*K2/2); K4 = f (t0 + h, y0 + hK3) = Psub (t0 + h, y0 + h*K3); y1 = y0 + (h/6)*(K1 + 2K2 + 2K3 + K4); tt=[tt,t1]; yy=[yy, y1]; t0 =t1; y0 =y1; end 主程序： 主程序：main.m Plot(tt,yy) %输出数据或图形 输出数据或图形 function dx=Psub(t,y) dx=[]; dx(1)=f1(t,y); dx(2)=f2(t,y); ... dx(n)=fn(t,y); dx=[dx(1); dx(2); ...dx(n)];