基于Davis阻力模型的列车运行仿真及优化分析

问题描述: 假设一辆列车在水平轨道上运行，从站台A运行至站台B，其间距为5144.7m，运行的速度上限为100km/h，列车质量为176.3t，列车旋转部件惯性的旋转质量因数ρ = 1.08，列车电机的最大牵引力为310KN，机械制动部件的最大制动力为 760KN。列车受到的阻力满足Davis阻力方程f = 2.0895 + 0.0098v + 0.006v^2 ，该公式中的速度单位为m/s，阻力单位为KN。请你通过建模方法编写matlab程序以获得列车运行过程的速度-距离曲线、牵引制动力-距离曲线、时间-距离曲线与能量消耗-距离曲线。注意:使用欧拉法求解微分方程，列车在A,B处速度为零，且牵引/制动力应为变力，并且需要考虑旋转质量因数。需要获取列车以最短时间到达站台B、在最短运行时间上分别增加10s、20s、50s、150s、300s到达站台B总共六组曲线。

解题思路:

对列车运行过程进行建模，包括受力分析和运动学分析。
使用欧拉法求解微分方程，得到速度-距离曲线、牵引制动力-距离曲线、时间-距离曲线与能量消耗-距离曲线。
考虑不同时间到达站台B的情况，分别绘制六组曲线。

具体步骤:

列车受力分析

列车受到的力包括：牵引力、制动力、重力和阻力。

其中，重力可以表示为：

$F_g = mg$

阻力可以根据给定的Davis阻力方程表示为：

$f = 2.0895 + 0.0098v + 0.006v^2 $

牵引力和制动力需要根据列车的运动状态进行判断，即当列车速度小于100km/h时，牵引力可以取310KN，制动力取0；当列车速度大于100km/h时，牵引力取0，制动力可以取760KN。

运动学分析

根据牛顿第二定律，列车的加速度可以表示为：

$a = \frac{F_{\text{net}}}{m}$

其中，$F_{\text{net}}$为列车受到的净力，$m$为列车质量。

根据旋转惯量的定义，可以得到旋转惯量的表达式：

$I = \rho \cdot m \cdot r^2$

其中，$\rho$为旋转质量因数，$r$为旋转部件到转轴的距离。

列车的运动状态可以表示为：

$\begin{aligned} v &= \frac{ds}{dt} \ a &= \frac{dv}{dt} \end{aligned}$

欧拉法求解微分方程

根据欧拉法的公式，可以得到列车运动状态的更新公式：

$\begin{aligned} v_{i+1} &= v_i + a_i \cdot \Delta t \ s_{i+1} &= s_i + v_i \cdot \Delta t \end{aligned}$

其中，$v_i$和$s_i$分别表示第$i$个时间步长时的速度和位移，$a_i$表示第$i$个时间步长时的加速度，$\Delta t$表示时间步长。

绘制曲线

根据上述分析，可以得到列车运行过程的速度-距离曲线、牵引制动力-距离曲线、时间-距离曲线和能量消耗-距离曲线。根据题目要求，需要考虑不同时间到达站台B的情况，分别绘制六组曲线。

Matlab程序如下:

% 参数定义
m = 176.3*1000; % 列车质量（kg）
rho = 1.08; % 旋转质量因数
max_traction_force = 310*1000; % 最大牵引力（N）
max_brake_force = 760*1000; % 最大制动力（N）
g = 9.8; % 重力加速度（m/s^2）
distance = 5144.7; % 站台间距（m）
max_speed = 100*1000/3600; % 速度上限（m/s）
time_step = 0.1; % 时间步长（s）

% Davis阻力方程
function f = davis_resistance(v)
f = 2.0895 + 0.0098*v + 0.006*v^2;
end

% 初始化变量
s = 0; % 初始距离（m）
v = 0; % 初始速度（m/s）
t = 0; % 初始时间（s）
energy_consumption = 0; % 初始能量消耗（J）

% 存储数据
s_data = [];
v_data = [];
t_data = [];
force_data = [];
energy_consumption_data = [];

% 运行时间
run_times = [0, 10, 20, 50, 150, 300];

for i = 1:length(run_times)
    % 初始化变量
s = 0;
v = 0;
t = 0;
energy_consumption = 0;
    
    % 存储数据
s_data = [];
v_data = [];
t_data = [];
force_data = [];
energy_consumption_data = [];
    
    % 循环运行至到达站台B
    while s < distance
        % 计算阻力
        resistance_force = davis_resistance(v);
        
        % 计算净力
        if v < max_speed
            traction_force = max_traction_force;
            brake_force = 0;
        else
            traction_force = 0;
            brake_force = max_brake_force;
        end
        net_force = traction_force - brake_force - resistance_force - m*g;
        
        % 计算加速度
        a = net_force / m;
        
        % 更新速度和距离
        v = v + a*time_step;
        s = s + v*time_step;
        
        % 计算能量消耗
        energy_consumption = energy_consumption + traction_force*v*time_step;
        
        % 更新时间
        t = t + time_step;
        
        % 存储数据
s_data = [s_data, s];
v_data = [v_data, v];
t_data = [t_data, t];
force_data = [force_data, traction_force - brake_force];
energy_consumption_data = [energy_consumption_data, energy_consumption];
    end
    
    % 绘制曲线
    figure
    subplot(2,2,1)
    plot(s_data, v_data)
    title('速度-距离曲线')
    xlabel('距离（m）')
    ylabel('速度（m/s）')
    
    subplot(2,2,2)
    plot(s_data, force_data)
    title('牵引制动力-距离曲线')
    xlabel('距离（m）')
    ylabel('牵引制动力（N）')
    
    subplot(2,2,3)
    plot(s_data, t_data)
    title('时间-距离曲线')
    xlabel('距离（m）')
    ylabel('时间（s）')
    
    subplot(2,2,4)
    plot(s_data, energy_consumption_data)
    title('能量消耗-距离曲线')
    xlabel('距离（m）')
    ylabel('能量消耗（J）')
    
    % 设置图例
    sgtitle(['运行时间：', num2str(run_times(i)), 's'])
    
    % 保存图片
saveas(gcf, ['run_time_', num2str(run_times(i)), '.jpg'])
e