fastsim_core/drive_cycle/
maneuvers.rs

1use std::collections::HashSet;
2
3use crate::{imports::*, simdrive::roadload::StepInfo, simdrive::SimParams, vehicle::Vehicle};
4
5use super::manipulation_utils::trapz_distance_for_step;
6use super::{
7    manipulation_utils::{
8        accel_array_for_constant_jerk, accel_for_constant_jerk, calc_constant_jerk_trajectory,
9        trapz_step_distances, trapz_step_start_distance, CoastTrajectory, CycleCache, PassingInfo,
10        RendezvousTrajectory,
11    },
12    Cycle,
13};
14
15#[serde_api]
16#[derive(Clone, Debug, Deserialize, Serialize, PartialEq)]
17#[non_exhaustive]
18#[serde(deny_unknown_fields)]
19#[cfg_attr(feature = "pyo3", pyclass(module = "fastsim", subclass, eq))]
20pub struct Maneuver {
21    // Cycle Instances
22    /// Cycle to apply maneuver to
23    #[serde(default)]
24    pub cyc: Cycle,
25    /// Reference cycle
26    #[serde(default)]
27    pub cyc0: Cycle,
28
29    // Chassis Data
30    /// Constant mass to assume for maneuvers
31    #[serde(default)]
32    pub mass: si::Mass,
33    /// Constant air density to assume for manuvers
34    #[serde(default)]
35    pub air_density: si::MassDensity,
36    /// Constant aerodynamic drag coefficient to assume
37    #[serde(default)]
38    pub drag_coef: si::Ratio,
39    /// Constant Frontal Area of vehicle to assume
40    #[serde(default)]
41    pub frontal_area: si::Area,
42    /// Constant Wheel Rolling Resistance Coefficient to assume
43    #[serde(default)]
44    pub wheel_rr_coef: si::Ratio,
45    /// Wheel inertia per wheel
46    #[serde(default)]
47    pub wheel_inertia: si::MomentOfInertia,
48    /// Number of wheels
49    #[serde(default)]
50    pub num_wheels: u8,
51    /// Wheel radius
52    #[serde(default)]
53    pub wheel_radius: si::Length,
54
55    // Solver Settings
56    /// max number of iterations allowed in setting achieved speed when trace
57    /// cannot be achieved
58    #[serde(default)]
59    pub ach_speed_max_iter: u32,
60    /// tolerance in change in speed guess in setting achieved speed when trace
61    /// cannot be achieved
62    #[serde(default)]
63    pub ach_speed_tol: si::Ratio,
64    /// Newton method gain for setting achieved speed
65    #[serde(default)]
66    pub ach_speed_solver_gain: f64,
67
68    // Coasting Parameters
69    /// whether to allow coasting or not.
70    #[serde(default)]
71    pub coast_allow: bool,
72    /// for testing: triggers coasting when vehicle passes the given speed
73    #[serde(default)]
74    pub coast_start_speed: si::Velocity,
75    /// speed at which mechanical braking will initiate during coasting maneuvers
76    #[serde(default)]
77    pub coast_brake_start_speed: si::Velocity,
78    /// acceleration assumed during braking for coast maneuvers
79    /// NOTE: should be negative
80    #[serde(default)]
81    pub coast_brake_accel: si::Acceleration,
82    /// if true, accuracy will be favored over performance for grade per step
83    /// estimates Specifically, for performance, grade for a step will be
84    /// assumed to be the grade looked up at step start distance. For accuracy,
85    /// the actual elevations will be used. This distinciton only makes a
86    /// difference for CAV maneuvers.
87    #[serde(default)]
88    pub favor_grade_accuracy: bool,
89    /// if true, coasting vehicle can eclipse the shadow trace (i.e., reference
90    /// vehicle in front)
91    #[serde(default)]
92    pub coast_allow_passing: bool,
93    /// maximum allowable speed under coast
94    #[serde(default)]
95    pub coast_max_speed: si::Velocity,
96    /// "look-ahead" time for speed changes to be considered to feature coasting
97    /// to hit a given stopping distance mark
98    #[serde(default)]
99    pub coast_time_horizon_for_adjustment: si::Time,
100
101    // IDM - Intelligent Driver Model, Adaptive Cruise Control version
102    /// if true, initiates the IDM - Intelligent Driver Model, Adaptive Cruise
103    /// Control version
104    #[serde(default)]
105    pub idm_allow: bool,
106    /// IDM algorithm: a way to specify desired speed by course distance
107    /// traveled. Can simulate changing speed limits over a driving cycle.
108    /// optional list of (distance (m), desired speed (m/s)).
109    #[serde(default)]
110    pub idm_desired_speed_by_distance: Option<Vec<(si::Length, si::Velocity)>>,
111    /// IDM algorithm: desired speed (m/s). Only used if
112    /// idm_v_desired_in_m_per_s_by_distance_m is NOT set (i.e., is None)
113    #[serde(default)]
114    pub idm_desired_speed: si::Velocity,
115    /// IDM algorithm: headway time desired to vehicle in front (s)
116    #[serde(default)]
117    pub idm_headway: si::Time,
118    /// IDM algorithm: minimum desired gap between vehicle and lead vehicle (m)
119    #[serde(default)]
120    pub idm_minimum_gap: si::Length,
121    /// IDM algorithm: delta parameter
122    #[serde(default)]
123    pub idm_delta: f64,
124    /// IDM algorithm: acceleration parameter
125    #[serde(default)]
126    pub idm_acceleration: si::Acceleration,
127    /// IDM algorithm: deceleration parameter
128    #[serde(default)]
129    pub idm_deceleration: si::Acceleration,
130
131    // Internal Fields
132    pub i: usize,
133    pub coast_delay_index: Vec<i32>,
134    #[serde(default)]
135    pub impose_coast: Vec<bool>,
136    #[serde(default)]
137    pub idm_target_speed_m_per_s: Vec<f64>,
138
139    pub cyc0_cache: CycleCache,
140}
141
142#[pyo3_api]
143impl Maneuver {
144    #[pyo3(name = "create_from")]
145    #[staticmethod]
146    /// create a maneuver object based on the given cycle and vehicle.
147    fn create_from_py(cyc: &Cycle, veh: &Vehicle) -> PyResult<Self> {
148        Ok(Maneuver::from(cyc, veh))
149    }
150
151    #[pyo3(name = "apply_maneuvers")]
152    /// apply all maneuvers to the cycle and return it.
153    fn apply_maneuvers_py(&mut self) -> PyResult<Cycle> {
154        self.apply();
155        let cyc = self.cyc.clone();
156        Ok(cyc)
157    }
158
159    #[pyo3(name = "is_coasting")]
160    /// return a vector of signals indicating 1.0 for coast, otherwise 0.0
161    fn is_coasting_py(&self) -> PyResult<Vec<f64>> {
162        let mut result = Vec::with_capacity(self.impose_coast.len());
163        for ic in &self.impose_coast {
164            if *ic {
165                result.push(1.0);
166            } else {
167                result.push(0.0);
168            }
169        }
170        Ok(result)
171    }
172}
173
174impl SerdeAPI for Maneuver {
175    #[cfg(feature = "resources")]
176    const RESOURCES_SUBDIR: &'static str = "maneuvers";
177}
178
179impl Init for Maneuver {
180    fn init(&mut self) -> Result<(), Error> {
181        self.i = 1;
182        let n = self.cyc.speed.len();
183        self.coast_delay_index = vec![0; n];
184        self.impose_coast = vec![false; n];
185        self.idm_target_speed_m_per_s = vec![0.0; n];
186        self.cyc0_cache = self.cyc0.build_cache();
187        Ok(())
188    }
189}
190
191impl Default for Maneuver {
192    fn default() -> Self {
193        Self {
194            cyc: Cycle::default(),
195            cyc0: Cycle::default(),
196            mass: si::Mass::default(),
197            air_density: si::MassDensity::default(),
198            drag_coef: si::Ratio::default(),
199            frontal_area: si::Area::default(),
200            wheel_rr_coef: si::Ratio::default(),
201            wheel_inertia: si::MomentOfInertia::default(),
202            num_wheels: u8::default(),
203            wheel_radius: si::Length::default(),
204            ach_speed_max_iter: 3,
205            ach_speed_tol: 1.0e-3 * uc::R,
206            ach_speed_solver_gain: 0.9,
207            coast_allow: false,
208            coast_start_speed: 0.0 * uc::MPS,
209            coast_brake_start_speed: 20.0 * uc::MPH,
210            coast_brake_accel: -2.5 * uc::MPS2,
211            favor_grade_accuracy: true,
212            coast_allow_passing: false,
213            coast_max_speed: 40.0 * uc::MPS,
214            coast_time_horizon_for_adjustment: 20.0 * uc::S,
215            idm_allow: false,
216            idm_desired_speed_by_distance: None,
217            idm_desired_speed: 75.0 * uc::MPH,
218            idm_headway: 1.0 * uc::S,
219            idm_minimum_gap: 2.0 * uc::M,
220            idm_delta: 4.0,
221            idm_acceleration: 1.0 * uc::MPS2,
222            idm_deceleration: 1.5 * uc::MPS2,
223            i: 1,
224            coast_delay_index: vec![0, 0],
225            impose_coast: vec![false, false],
226            idm_target_speed_m_per_s: vec![0.0, 0.0],
227            cyc0_cache: CycleCache::default(),
228        }
229    }
230}
231
232impl Maneuver {
233    /// Create maneuver object from cycle and vehicle.
234    pub fn from(cyc: &Cycle, veh: &Vehicle) -> Self {
235        let mut c = cyc.clone();
236        c.init().unwrap();
237        let mut v = veh.clone();
238        v.init().unwrap();
239        let cyc0 = c.clone();
240        let cyc0_cache = cyc0.build_cache();
241        let default_mass = 1200.0 * uc::KG;
242        let mass = if let Ok(Some(m)) = veh.mass() {
243            m
244        } else {
245            default_mass
246        };
247        // TODO[mok]: what is the proper way to get air_density from the
248        //            vehicle? The below returns 0.0...
249        // *veh.state.air_density.get_fresh(|| format_dbg!()).unwrap();
250        let air_density = 1.2 * uc::KGPM3;
251
252        // NOTE: wheel radius should exist as Some(value) after v.init() above.
253        let wheel_radius = veh.chassis.wheel_radius.unwrap();
254        let params = SimParams::default();
255        Self {
256            cyc: c,
257            cyc0,
258            cyc0_cache,
259            mass,
260            air_density,
261            drag_coef: veh.chassis.drag_coef,
262            frontal_area: veh.chassis.frontal_area,
263            wheel_rr_coef: veh.chassis.wheel_rr_coef,
264            wheel_inertia: veh.chassis.wheel_inertia,
265            num_wheels: veh.chassis.num_wheels,
266            wheel_radius,
267            ach_speed_max_iter: params.ach_speed_max_iter,
268            ach_speed_tol: params.ach_speed_tol,
269            ach_speed_solver_gain: params.ach_speed_solver_gain,
270            ..Default::default()
271        }
272    }
273
274    /// Apply the eco-coast maneuver to the given cycle with
275    /// the given reference cycle.
276    /// - cyc: cycle to modify
277    /// - cyc0: reference cycle
278    pub fn apply(&mut self) {
279        self.cyc.init().unwrap();
280        self.cyc0.init().unwrap();
281        self.i = 1;
282        let cyc_len = self.cyc.time.len();
283        self.coast_delay_index = vec![0; cyc_len];
284        self.impose_coast = vec![false; cyc_len];
285        self.idm_target_speed_m_per_s = vec![0.0; cyc_len];
286        self.cyc0_cache = self.cyc0.build_cache();
287        self.walk(cyc_len);
288    }
289
290    /// Walk from step to step for the maneuver simulation.
291    fn walk(&mut self, cyc_len: usize) {
292        while self.i < cyc_len {
293            self.step();
294        }
295        // NOTE: force dist and elev to recalculate
296        // TODO: need to investigate re-deriving grade from interpolation of
297        // elevation by distance.
298        self.cyc.dist = vec![];
299        self.cyc.elev = vec![];
300        self.cyc.init().unwrap();
301    }
302
303    /// Step: compute a single time-step for the maneuver simulation.
304    fn step(&mut self) {
305        if self.idm_allow {
306            self.idm_target_speed_m_per_s[self.i] = match &self.idm_desired_speed_by_distance {
307                Some(vtgt_by_dist) => {
308                    let mut found_v_target = vtgt_by_dist[0].1;
309                    let mut current_d = si::Length::ZERO;
310                    for (idx, d) in self.cyc.dist.iter().enumerate() {
311                        if idx > self.i {
312                            break;
313                        }
314                        current_d += *d;
315                    }
316                    for (d, v_target) in vtgt_by_dist {
317                        if current_d >= *d {
318                            found_v_target = *v_target;
319                        } else {
320                            break;
321                        }
322                    }
323                    found_v_target.get::<si::meter_per_second>()
324                }
325                None => self.idm_desired_speed.get::<si::meter_per_second>(),
326            };
327            self.set_speed_for_target_gap_using_idm(self.i);
328        }
329        if self.coast_allow {
330            self.set_coast_speed(self.i);
331            self.cyc.grade[self.i] = self.lookup_grade_for_step(self.i, None);
332        }
333        self.i += 1;
334    }
335
336    /// Set gap
337    /// - i: non-negative integer, the step index
338    ///
339    /// RETURN: None
340    ///
341    /// EFFECTS:
342    /// - sets the next speed (m/s)
343    ///
344    /// EQUATION:
345    /// parameters:
346    ///     - v_desired: the desired speed (m/s)
347    ///     - delta: number, typical value is 4.0
348    ///     - a: max acceleration, (m/s2)
349    ///     - b: max deceleration, (m/s2)
350    /// s = d_lead - d
351    /// dv/dt = a * (1 - (v/v_desired)**delta - (s_desired(v,v-v_lead)/s)**2)
352    /// s_desired(v, dv) = s0 + max(0, v*dt_headway + (v * dv)/(2.0 * sqrt(a*b)))
353    /// REFERENCE:
354    /// Treiber, Martin and Kesting, Arne. 2013. "Chapter 11: Car-Following Models Based on Driving Strategies".
355    ///     Traffic Flow Dynamics: Data, Models and Simulation. Springer-Verlag. Springer, Berlin, Heidelberg.
356    ///     DOI: <https://doi.org/10.1007/978-3-642-32460-4>
357    pub fn set_speed_for_target_gap_using_idm(&mut self, i: usize) {
358        // PARAMETERS
359        let v_desired_m_per_s = if self.idm_target_speed_m_per_s[i] > 0.0 {
360            self.idm_target_speed_m_per_s[i]
361        } else {
362            let mut v = self.cyc0.speed[0];
363            for vi in &self.cyc0.speed {
364                if *vi > v {
365                    v = *vi;
366                }
367            }
368            v.get::<si::meter_per_second>()
369        };
370        // DERIVED VALUES
371        self.cyc.speed[i] = self.next_speed_by_idm(
372            i,
373            self.idm_acceleration.get::<si::meter_per_second_squared>(),
374            self.idm_deceleration.get::<si::meter_per_second_squared>(),
375            self.idm_headway.get::<si::second>(),
376            self.idm_minimum_gap.get::<si::meter>(),
377            v_desired_m_per_s,
378            self.idm_delta,
379        ) * uc::MPS;
380    }
381
382    /// Calculate the next speed by the Intelligent Driver Model
383    /// - i: int, the index
384    /// - a_m_per_s2: number, max acceleration (m/s2)
385    /// - b_m_per_s2: number, max deceleration (m/s2)
386    /// - dt_headway_s: number, the headway between us and the lead vehicle in seconds
387    /// - s0_m: number, the initial gap between us and the lead vehicle in meters
388    /// - v_desired_m_per_s: number, the desired speed in (m/s)
389    /// - delta: number, a shape parameter; typical value is 4.0
390    ///
391    /// RETURN: number, the next speed (m/s)
392    ///
393    /// REFERENCE:
394    /// Treiber, Martin and Kesting, Arne. 2013. "Chapter 11: Car-Following Models Based on Driving Strategies".
395    ///     Traffic Flow Dynamics: Data, Models and Simulation. Springer-Verlag. Springer, Berlin, Heidelberg.
396    ///     DOI: <https://doi.org/10.1007/978-3-642-32460-4>.
397    #[allow(clippy::too_many_arguments)]
398    pub fn next_speed_by_idm(
399        &mut self,
400        i: usize,
401        a_m_per_s2: f64,
402        b_m_per_s2: f64,
403        dt_headway_s: f64,
404        s0_m: f64,
405        v_desired_m_per_s: f64,
406        delta: f64,
407    ) -> f64 {
408        if v_desired_m_per_s <= 0.0 {
409            return 0.0;
410        }
411        let a_m_per_s2 = a_m_per_s2.abs();
412        let b_m_per_s2 = b_m_per_s2.abs();
413        let dt_headway_s = dt_headway_s.max(0.0);
414        // we assume the vehicles start out a "minimum gap" apart
415        let s0_m = s0_m.max(0.0);
416        // DERIVED VALUES
417        let sqrt_ab = (a_m_per_s2 * b_m_per_s2).powf(0.5);
418        let v0_m_per_s = self.cyc.speed[i - 1].get::<si::meter_per_second>();
419        let v0_lead_m_per_s = self.cyc0.speed[i - 1].get::<si::meter_per_second>();
420        let dv0_m_per_s = v0_m_per_s - v0_lead_m_per_s;
421        let d0_lead_m = self.cyc0_cache.trapz_distances_m[(i - 1).max(0)] + s0_m;
422        let d0_m = trapz_step_start_distance(&self.cyc, i).get::<si::meter>();
423        let s_m = (d0_lead_m - d0_m).max(0.01);
424        let dt = (self.cyc0.time[i] - self.cyc0.time[i - 1]).get::<si::second>();
425        // IDM EQUATIONS
426        let s_target_m = s0_m
427            + ((v0_m_per_s * dt_headway_s) + ((v0_m_per_s * dv0_m_per_s) / (2.0 * sqrt_ab)))
428                .max(0.0);
429        let accel_target_m_per_s2 = a_m_per_s2
430            * (1.0 - ((v0_m_per_s / v_desired_m_per_s).powf(delta)) - ((s_target_m / s_m).powi(2)));
431        (v0_m_per_s + (accel_target_m_per_s2 * dt)).max(0.0)
432    }
433
434    /// For situations where cyc can deviate from cyc0, this method
435    /// looks up and accurately interpolates what the average grade over
436    /// the step should be. The achieved value is used to predict the
437    /// distance traveled over the step.
438    ///
439    /// NOTE:
440    /// If not allowing coasting (i.e., sim_params.coast_allow == False)
441    /// and not allowing IDM/following (i.e., self.sim_params.idm_allow
442    /// == False) then returns self.cyc.grade\[i\]
443    pub fn lookup_grade_for_step(&self, i: usize, speed_ach: Option<si::Velocity>) -> si::Ratio {
444        if self.cyc0_cache.grade_all_zero {
445            return 0.0 * uc::R;
446        }
447        if !self.coast_allow && !self.idm_allow {
448            return self.cyc.grade[i];
449        }
450        match speed_ach {
451            Some(v1) => {
452                let dt = self.cyc.time[i] - self.cyc.time[i - 1];
453                self.cyc0.average_grade_over_range(
454                    trapz_step_start_distance(&self.cyc, i),
455                    0.5 * (v1 + self.cyc.speed[i - 1]) * dt,
456                    Some(&self.cyc0_cache),
457                )
458            }
459            None => self.cyc0.average_grade_over_range(
460                trapz_step_start_distance(&self.cyc, i),
461                trapz_distance_for_step(&self.cyc, i),
462                Some(&self.cyc0_cache),
463            ),
464        }
465    }
466
467    /// Determine whether the vehicle should go into a 'coasting' state.
468    /// Normal coasting logic is that the vehicle will coast if it is
469    /// within coasting distance of a stop:
470    /// - if distance to coast from start of step <= distance to next stop
471    /// - AND distance to coast from end of step (using reference speed) is
472    ///   > distance to next step
473    /// - AND vehicle was at or above the speed to start braking
474    /// - AND at least four time-steps away from where braking would start
475    ///
476    /// NOTE: for the case when coast-start speed is used, we only worry
477    /// about if the vehicle is above the coast-start speed. This is mainly
478    /// for testing.
479    pub fn should_impose_coast(&mut self, i: usize) -> bool {
480        let v0 = self.cyc.speed[i - 1];
481        if self.coast_start_speed > si::Velocity::ZERO {
482            return v0 >= self.coast_start_speed;
483        }
484        if v0 < self.coast_brake_start_speed {
485            return false;
486        }
487        // distance to stop by coasting from start of step (i - 1)
488        let dtsc0 = self.calc_distance_to_stop_coast_v2(i);
489        if dtsc0.is_none() {
490            return false;
491        }
492        let dtsc0 = dtsc0.unwrap();
493        // distance to next stop (m)
494        let d0 = trapz_step_start_distance(&self.cyc, i);
495        let dts0 = self
496            .cyc0
497            .calc_distance_to_next_stop_from(d0, Some(&self.cyc0_cache));
498        let dtb = -0.5 * v0 * v0 / self.coast_brake_accel;
499        dtsc0 >= dts0 && dts0 >= (4.0 * dtb)
500    }
501
502    /// Calculate the distance to stop via coasting.
503    /// - i: the current index
504    ///
505    /// RETURN: if Some, the distance, else None
506    ///
507    /// NOTES:
508    /// - if None, that means there is no solution to a coast-down distance.
509    ///   This can happen due to being too close to the given stop or
510    ///   perhaps due to coasting downhill (i.e., will not stop).
511    /// - if Some, the distance in meters that the vehicle would freely coast
512    ///   is unobstructed. We do account fro grade between the current point
513    ///   and the end-point.
514    pub fn calc_distance_to_stop_coast_v2(&mut self, i: usize) -> Option<si::Length> {
515        let not_found = -1.0;
516        let v0 = self.cyc.speed[i - 1].get::<si::meter_per_second>();
517        let v_brake = self.coast_brake_start_speed.get::<si::meter_per_second>();
518        let a_brake = self.coast_brake_accel.get::<si::meter_per_second_squared>();
519        let ds = &self.cyc0_cache.trapz_distances_m;
520        let gs: Vec<f64> = self
521            .cyc0
522            .grade
523            .iter()
524            .map(|g| g.get::<si::ratio>())
525            .collect();
526        assert!(
527            ds.len() == gs.len(),
528            "Assumed lengths of distances and grades must equal. ds.len()={}, gs.len()={}",
529            ds.len(),
530            gs.len(),
531        );
532        let d0 = trapz_step_start_distance(&self.cyc, i).get::<si::meter>();
533        let mut grade_by_distance = Vec::with_capacity(ds.len());
534        for idx in 0..ds.len() {
535            if ds[idx] >= d0 {
536                grade_by_distance.push(gs[idx]);
537            }
538        }
539        let veh_mass_kg = self.mass.get::<si::kilogram>();
540        let air_density_kg_per_m3 = self.air_density.get::<si::kilogram_per_cubic_meter>();
541        let cdfa_m2 =
542            self.drag_coef.get::<si::ratio>() * self.frontal_area.get::<si::square_meter>();
543        let rrc = self.wheel_rr_coef.get::<si::ratio>();
544        let gravity_m_per_s2 = uc::ACC_GRAV.get::<si::meter_per_second_squared>();
545        // distance traveled while stopping via friction-braking (i.e., distance to brake)
546        let dtb = -0.5 * v_brake * v_brake / a_brake;
547        if v0 <= v_brake {
548            let result = -0.5 * v0 * v0 / a_brake;
549            return Some(result * uc::M);
550        }
551        let grade_mult = 10000;
552        let unique_grades = {
553            let mut result = HashSet::new();
554            for gr in &grade_by_distance {
555                let gr_to_store = (*gr * (grade_mult as f64)).round() as i32;
556                result.insert(gr_to_store);
557            }
558            result
559        };
560        if unique_grades.len() == 1 {
561            // if there is only one grade, there may be a closed-form solution
562            let unique_grade = (*unique_grades.iter().nth(0).unwrap() as f64) / (grade_mult as f64);
563            let theta = unique_grade.atan();
564            let c1 = gravity_m_per_s2 * (theta.sin() + rrc * theta.cos());
565            let c2 = (air_density_kg_per_m3 * cdfa_m2) / (2.0 * veh_mass_kg);
566            let v02 = v0 * v0;
567            let vb2 = v_brake * v_brake;
568            let mut d = not_found;
569            let a1 = c1 + c2 * v02;
570            let b1 = c1 + c2 * vb2;
571            if c2 == 0.0 {
572                if c1 > 0.0 {
573                    d = (1.0 / (2.0 * c1)) * (v02 - vb2);
574                }
575            } else if a1 > 0.0 && b1 > 0.0 {
576                d = (1.0 / (2.0 * c2)) * (a1.ln() - b1.ln());
577            }
578            if d != not_found {
579                let result = d + dtb;
580                return Some(result * uc::M);
581            }
582        }
583        let ct = self.generate_coast_trajectory(i);
584        if ct.found_trajectory {
585            Some(ct.distance_to_stop_via_coast_m * uc::M)
586        } else {
587            None
588        }
589    }
590
591    /// Set the coasting speed.
592    pub fn set_coast_speed(&mut self, i: usize) {
593        let tol = 1e-6;
594        let v0 = self.cyc.speed[i - 1].get::<si::meter_per_second>();
595        if v0 > tol && !self.impose_coast[i] && self.should_impose_coast(i) {
596            let ct = self.generate_coast_trajectory(i);
597            if ct.found_trajectory {
598                let d = ct.distance_to_stop_via_coast_m;
599                if d < 0.0 {
600                    for idx in i..self.cyc0.speed.len() {
601                        self.impose_coast[idx] = false;
602                    }
603                } else {
604                    self.apply_coast_trajectory(&ct);
605                }
606                if !self.coast_allow_passing {
607                    self.prevent_collisions(i, None);
608                }
609            }
610        }
611        if !self.impose_coast[i] {
612            if !self.idm_allow {
613                // NOTE: transforming to i32 so we can carry negative indices.
614                let i_i32 = i32::try_from(i).ok();
615                let target_idx = i_i32.map(|v| v - self.coast_delay_index[i]);
616                let target_idx = match target_idx {
617                    Some(ti) => {
618                        if ti < 0 {
619                            Some(0)
620                        } else {
621                            usize::try_from(ti).ok()
622                        }
623                    }
624                    None => None,
625                };
626                if let Some(ti) = target_idx {
627                    self.cyc.speed[i] = self.cyc0.speed[ti.min(self.cyc0.speed.len() - 1)];
628                }
629            }
630            return;
631        }
632        let v1_traj = self.cyc.speed[i].get::<si::meter_per_second>();
633        let v_brake = self.coast_brake_start_speed.get::<si::meter_per_second>();
634        if v0 > v_brake {
635            if self.coast_allow_passing {
636                // NOTE: We could be coasting downhill so could in theory go
637                // to a higher speed. Since we can pass, allow vehicle to go
638                // up to max coasting speed (m/s). The solver will show us what
639                // we can actually achieve.
640                self.cyc.speed[i] = self.coast_max_speed;
641            } else {
642                self.cyc.speed[i] =
643                    v1_traj.min(self.coast_max_speed.get::<si::meter_per_second>()) * uc::MPS;
644            }
645        }
646        // Solve for the actual coasting speed
647        let coast_speed = self.solve_step(i);
648        if self.impose_coast[i - 1] && v1_traj <= v_brake {
649            // NOTE: if we've been coasting for at least one step already
650            // and the current trajectory takes us below v_brake, we should
651            // use that as it is the brake trajectory.
652            self.cyc.speed[i] = v1_traj * uc::MPS;
653        } else {
654            self.cyc.speed[i] = coast_speed;
655        }
656        let v_tol = tol * uc::MPS;
657        let dt = self.cyc.time[i] - self.cyc.time[i - 1];
658        let accel_proposed = (self.cyc.speed[i] - self.cyc.speed[i - 1]) / dt;
659        if self.cyc.speed[i] < v_tol {
660            for idx in i..self.cyc0.speed.len() {
661                self.impose_coast[idx] = false;
662            }
663            self.set_coast_delay(i);
664            self.cyc.speed[i] = si::Velocity::ZERO;
665            return;
666        }
667        if (self.cyc.speed[i] - v1_traj * uc::MPS).abs() > v_tol {
668            let mut adjusted_current_speed = false;
669            let brake_speed_start_tol = 0.1 * uc::MPS;
670            if self.cyc.speed[i] < (self.coast_brake_start_speed - brake_speed_start_tol) {
671                let (_, num_steps) =
672                    self.cyc
673                        .modify_with_braking_trajectory(self.coast_brake_accel, i, None);
674                for idx in i..self.cyc0.speed.len() {
675                    self.impose_coast[idx] = idx < (i + num_steps);
676                }
677                adjusted_current_speed = true;
678            } else {
679                // NOTE: will the below work when coasting downhill and picking up speed?
680                let (traj_found, traj_n, traj_jerk_m_per_s3, traj_accel_m_per_s2) = self
681                    .calc_next_rendezvous_trajectory(
682                        i,
683                        self.coast_brake_accel,
684                        accel_proposed.min(0.0 * uc::MPS2),
685                    );
686                if traj_found {
687                    // adjust cyc to perform the trajectory
688                    let final_speed = self.cyc.modify_by_const_jerk_trajectory(
689                        i,
690                        traj_n,
691                        traj_jerk_m_per_s3 * uc::MPS3,
692                        traj_accel_m_per_s2 * uc::MPS2,
693                    );
694                    for idx in i..self.cyc0.speed.len() {
695                        self.impose_coast[idx] = idx < (i + traj_n);
696                    }
697                    adjusted_current_speed = true;
698                    let i_for_brake = i + traj_n;
699                    if (final_speed - self.coast_brake_start_speed).abs() < brake_speed_start_tol {
700                        let (_, num_steps) = self.cyc.modify_with_braking_trajectory(
701                            self.coast_brake_accel,
702                            i_for_brake,
703                            None,
704                        );
705                        for idx in i_for_brake..self.cyc0.speed.len() {
706                            self.impose_coast[idx] = idx < i_for_brake + num_steps;
707                        }
708                        adjusted_current_speed = true;
709                    } else {
710                        // TODO: logging feature?
711                        //#[cfg(feature = "logging")]
712                        println!("## WARNING ##");
713                        println!("final_speed={:?} not close to coast_brake_start_speed={:?} for i={:?}; i_for_brake={:?}, traj_n={:?}",
714                            final_speed, self.coast_brake_start_speed, i, i_for_brake, traj_n);
715                        //log::warn!(
716                        //    "final_speed={}i not close to coast_brake_start_speed={} for i={}; i_for_brake={}, traj_n={}",
717                        //    final_speed,
718                        //    self.coast_brake_start_speed,
719                        //    i,
720                        //    i_for_brake,
721                        //    traj_n
722                        //);
723                    }
724                }
725            }
726            if adjusted_current_speed {
727                if !self.coast_allow_passing {
728                    self.prevent_collisions(i, None);
729                }
730                self.cyc.speed[i] = self.solve_step(i);
731            }
732        }
733    }
734
735    /// Solve for coast speed and set.
736    pub fn solve_step(&mut self, i: usize) -> si::Velocity {
737        let dt = self.cyc.time[i] - self.cyc.time[i - 1];
738        let step_info = StepInfo {
739            dt,
740            speed_prev: self.cyc.speed[i - 1],
741            cyc_speed: self.cyc.speed[i],
742            grade_curr: self.cyc.grade[i],
743            air_density: self.air_density,
744            mass: self.mass,
745            drag_coef: self.drag_coef,
746            frontal_area: self.frontal_area,
747            wheel_inertia: self.wheel_inertia,
748            num_wheels: self.num_wheels,
749            wheel_radius: self.wheel_radius,
750            wheel_rr_coef: self.wheel_rr_coef,
751            pwr_prop_fwd_max: 0.0 * uc::KW,
752        };
753        let coast_speed = step_info.solve_for_speed(
754            self.ach_speed_max_iter,
755            self.ach_speed_tol,
756            self.ach_speed_solver_gain,
757        );
758        let max_coast_speed = self.coast_max_speed.min(coast_speed);
759        let brake_start_speed = self.coast_brake_start_speed + 0.1 * uc::MPS;
760        if coast_speed > brake_start_speed {
761            max_coast_speed
762        } else {
763            // NOTE: We follow trace below the brake start speed as the
764            // brake trajectory has been added to the cycle
765            self.cyc.speed[i]
766                .min(max_coast_speed)
767                .max(si::Velocity::ZERO)
768        }
769    }
770
771    /// Calculate next rendezvous trajectory for eco-coasting
772    /// - i: the index into cyc for the end of start-of-step (i.e., the step
773    ///   that may be modified; should be i)
774    /// - min_accel: the minimum acceleration permitted
775    /// - max_accel: the maximum acceleration permitted
776    ///
777    /// RETURN: (Tuple
778    ///     found_rendezvous: Bool, if True the remainder of the data is valid; if False, no rendezvous found
779    ///     n: positive integer, the number of steps ahead to rendezvous at
780    ///     jerk_m__s3: number, the Jerk or first-derivative of acceleration (m/s3)
781    ///     accel_m__s2: number, the initial acceleration of the trajectory (m/s2)
782    /// )
783    /// If no rendezvous exists within the scope, the returned tuple has False for the first item.
784    /// Otherwise, returns the next closest rendezvous in time/space
785    pub fn calc_next_rendezvous_trajectory(
786        &self,
787        i: usize,
788        min_accel: si::Acceleration,
789        max_accel: si::Acceleration,
790    ) -> (bool, usize, f64, f64) {
791        let tol = 1e-6;
792        let min_accel_m_per_s2 = min_accel.get::<si::meter_per_second_squared>();
793        let max_accel_m_per_s2 = max_accel.get::<si::meter_per_second_squared>();
794        // v0 is where n=0; i.e., i - 1
795        let v0 = self.cyc.speed[i - 1].get::<si::meter_per_second>();
796        let brake_start_speed_m_per_s = self.coast_brake_start_speed.get::<si::meter_per_second>();
797        let brake_accel_m_per_s2 = self.coast_brake_accel.get::<si::meter_per_second_squared>();
798        let time_horizon_s = self.coast_time_horizon_for_adjustment.get::<si::second>();
799        // distance_horizon_m = 1_000.0;
800        let not_found_n = 0;
801        let not_found_jerk_m_per_s3 = 0.0;
802        let not_found_accel_m_per_s2 = 0.0;
803        let not_found = (
804            false,
805            not_found_n,
806            not_found_jerk_m_per_s3,
807            not_found_accel_m_per_s2,
808        );
809        if v0 < (brake_start_speed_m_per_s + tol) {
810            // don't process braking
811            return not_found;
812        }
813        let (min_accel_m_per_s2, max_accel_m_per_s2) = if min_accel_m_per_s2 > max_accel_m_per_s2 {
814            (max_accel_m_per_s2, min_accel_m_per_s2)
815        } else {
816            (min_accel_m_per_s2, max_accel_m_per_s2)
817        };
818        let num_samples = self.cyc.speed.len();
819        let d0 = trapz_step_start_distance(&self.cyc, i).get::<si::meter>();
820        // a_proposed = (v1 - v0) / dt
821        // distance to stop from start of time-step
822        let dts0 = self
823            .cyc0
824            .calc_distance_to_next_stop_from(d0 * uc::M, Some(&self.cyc0_cache))
825            .get::<si::meter>();
826        if dts0 < 0.0 {
827            return not_found;
828        }
829        let dt = (self.cyc0.time[i] - self.cyc0.time[i - 1]).get::<si::second>();
830        // distance to brake from the brake start speed (m/s)
831        let dtb =
832            -0.5 * brake_start_speed_m_per_s * brake_start_speed_m_per_s / brake_accel_m_per_s2;
833        // distance to brake initialization from start of time-step (m)
834        let dtbi0 = dts0 - dtb;
835        if dtbi0 < 0.0 {
836            return not_found;
837        }
838        // Now, check rendezvous trajectories
839        let mut step_idx = i;
840        let mut dt_plan = 0.0;
841        let mut r_best_found = false;
842        let mut r_best_n = 0;
843        let mut r_best_jerk_m_per_s3 = 0.0;
844        let mut r_best_accel_m_per_s2 = 0.0;
845        let mut r_best_accel_spread_m_per_s2 = 0.0;
846        while dt_plan <= time_horizon_s && step_idx < num_samples {
847            dt_plan += dt;
848            let step_ahead = step_idx - (i - 1);
849            if step_ahead == 1 {
850                // for brake init rendezvous
851                let accel = (brake_start_speed_m_per_s - v0) / dt;
852                let v1 = (v0 + accel * dt).max(0.0);
853                let dd_proposed = ((v0 + v1) / 2.0) * dt;
854                if (v1 - brake_start_speed_m_per_s).abs() < tol && (dtbi0 - dd_proposed).abs() < tol
855                {
856                    r_best_found = true;
857                    r_best_n = 1;
858                    r_best_jerk_m_per_s3 = 0.0;
859                    r_best_accel_m_per_s2 = accel;
860                    break;
861                }
862            } else {
863                // rendezvous trajectory for brake-start -- assumes fixed time-steps
864                if dtbi0 > 0.0 {
865                    let r_bi_traj = calc_constant_jerk_trajectory(
866                        step_ahead,
867                        0.0,
868                        v0,
869                        dtbi0,
870                        brake_start_speed_m_per_s,
871                        dt,
872                    );
873                    let r_bi_jerk_m_per_s3 = r_bi_traj.jerk_m_per_s3;
874                    let r_bi_accel_m_per_s2 = r_bi_traj.acceleration_m_per_s2;
875                    if r_bi_accel_m_per_s2 < max_accel_m_per_s2
876                        && min_accel_m_per_s2 < r_bi_accel_m_per_s2
877                        && r_bi_jerk_m_per_s3 >= 0.0
878                    {
879                        let as_bi = accel_array_for_constant_jerk(
880                            step_ahead,
881                            r_bi_accel_m_per_s2,
882                            r_bi_jerk_m_per_s3,
883                            dt,
884                        );
885                        let as_bi_min = as_bi.iter().cloned().reduce(f64::min).unwrap_or(0.0);
886                        let as_bi_max = as_bi.iter().cloned().reduce(f64::max).unwrap_or(0.0);
887                        let accel_spread = (as_bi_max - as_bi_min).abs();
888                        let flag = as_bi_max < (max_accel_m_per_s2 + 1e-6)
889                            && as_bi_min > (min_accel_m_per_s2 - 1e-6)
890                            && (!r_best_found || (accel_spread < r_best_accel_spread_m_per_s2));
891                        if flag {
892                            r_best_found = true;
893                            r_best_n = step_ahead;
894                            r_best_accel_m_per_s2 = r_bi_accel_m_per_s2;
895                            r_best_jerk_m_per_s3 = r_bi_jerk_m_per_s3;
896                            r_best_accel_spread_m_per_s2 = accel_spread;
897                        }
898                    }
899                }
900            }
901            step_idx += 1;
902        }
903        if r_best_found {
904            (
905                r_best_found,
906                r_best_n,
907                r_best_jerk_m_per_s3,
908                r_best_accel_m_per_s2,
909            )
910        } else {
911            not_found
912        }
913    }
914
915    /// Coast Delay allows us to represent coasting to a stop when the lead
916    /// vehicle has already moved on from that stop.  In this case, the coasting
917    /// vehicle need not dwell at this or any stop while it is lagging behind
918    /// the lead vehicle in distance. Instead, the vehicle comes to a stop and
919    /// resumes mimicing the lead-vehicle trace at the first time-step the
920    /// lead-vehicle moves past the stop-distance. This index is the "coast delay index".
921    ///
922    /// Arguments
923    /// ---------
924    /// - i: the step index
925    ///
926    /// NOTE: Resets the coast_delay_index to 0 and calculates and sets the next
927    /// appropriate coast_delay_index if appropriate
928    pub fn set_coast_delay(&mut self, i: usize) {
929        let speed_tol = 0.01; // m/s
930        let dist_tol = 0.1; // m
931        for idx in i..self.cyc0.speed.len() {
932            // clear all future coast delays
933            self.coast_delay_index[idx] = 0;
934        }
935        let mut coast_delay = None;
936        if !self.idm_allow && self.cyc.speed[i].get::<si::meter_per_second>() < speed_tol {
937            let d0 = trapz_step_start_distance(&self.cyc, i).get::<si::meter>();
938            let d0_lv = self.cyc0_cache.trapz_distances_m[i - 1];
939            let dtlv0 = d0_lv - d0;
940            if dtlv0.abs() > dist_tol {
941                let mut d_lv = 0.0;
942                let mut min_dtlv = None;
943                for (idx, (&dd, &v)) in trapz_step_distances(&self.cyc0)
944                    .iter()
945                    .zip(self.cyc0.speed.iter())
946                    .enumerate()
947                {
948                    let dd = dd.get::<si::meter>();
949                    let v = v.get::<si::meter_per_second>();
950                    d_lv += dd;
951                    let dtlv = (d_lv - d0).abs();
952                    if v < speed_tol && (min_dtlv.is_none() || dtlv <= min_dtlv.unwrap()) {
953                        if min_dtlv.is_none()
954                            || dtlv < min_dtlv.unwrap()
955                            || (d0 < d0_lv && min_dtlv.unwrap() == dtlv)
956                        {
957                            let i_i32 = i32::try_from(i).unwrap();
958                            let idx_i32 = i32::try_from(idx).unwrap();
959                            coast_delay = Some(i_i32 - idx_i32);
960                        }
961                        min_dtlv = Some(dtlv);
962                    }
963                    if min_dtlv.is_some() && dtlv > min_dtlv.unwrap() {
964                        break;
965                    }
966                }
967            }
968        }
969        if let Some(cd) = coast_delay {
970            if cd < 0 {
971                let mut new_cd = cd;
972                for idx in i..self.cyc0.speed.len() {
973                    self.coast_delay_index[idx] = new_cd;
974                    new_cd += 1;
975                    if new_cd == 0 {
976                        break;
977                    }
978                }
979            } else {
980                for idx in i..self.cyc0.speed.len() {
981                    self.coast_delay_index[idx] = cd;
982                }
983            }
984        }
985    }
986
987    /// Generate a coast trajectory without actually modifying the cycle.
988    /// This can be used to calculate the distance to stop via coast using
989    /// actual time-stepping and changing grade.
990    pub fn generate_coast_trajectory(&mut self, i: usize) -> CoastTrajectory {
991        let v0 = self.cyc.speed[i - 1].get::<si::meter_per_second>();
992        let v_brake = self.coast_brake_start_speed.get::<si::meter_per_second>();
993        let a_brake = {
994            let result = self.coast_brake_accel.get::<si::meter_per_second_squared>();
995            if result > 0.0 {
996                -result
997            } else {
998                result
999            }
1000        };
1001        let ds = &self.cyc0_cache.trapz_distances_m;
1002        let d0 = trapz_step_start_distance(&self.cyc, i).get::<si::meter>();
1003        let mut distances_m = Vec::with_capacity(ds.len());
1004        let mut grade_by_distance = Vec::with_capacity(ds.len());
1005        for (idx, d) in ds.iter().enumerate() {
1006            if *d >= d0 {
1007                distances_m.push(*d - d0);
1008                grade_by_distance.push(self.cyc0.grade[idx].get::<si::ratio>());
1009            }
1010        }
1011        if distances_m.is_empty() {
1012            return CoastTrajectory {
1013                found_trajectory: false,
1014                distance_to_stop_via_coast_m: 0.0,
1015                start_idx: 0,
1016                speed_m_per_s: None,
1017                distance_to_brake_m: None,
1018            };
1019        }
1020        if v0 <= v_brake {
1021            return CoastTrajectory {
1022                found_trajectory: true,
1023                distance_to_stop_via_coast_m: -0.5 * v0 * v0 / a_brake,
1024                start_idx: i,
1025                speed_m_per_s: None,
1026                distance_to_brake_m: None,
1027            };
1028        }
1029        // dtb = distance to brake: distance traveled during the friction
1030        //       braking part of the maneuver
1031        let dtb = -0.5 * v_brake * v_brake / a_brake;
1032        let mut d = 0.0;
1033        let d_max = distances_m.last().unwrap() - dtb;
1034        let mut unique_grades = HashSet::with_capacity(ds.len());
1035        let grade_mult = 10000.0;
1036        for g in grade_by_distance.iter() {
1037            let grade = (g * grade_mult).round() as i32;
1038            unique_grades.insert(grade);
1039        }
1040        let unique_grade = if unique_grades.len() == 1 {
1041            let ug = unique_grades.iter().nth(0).unwrap();
1042            let ug = (*ug as f64) / grade_mult;
1043            Some(ug)
1044        } else {
1045            None
1046        };
1047        let has_unique_grade = unique_grade.is_some();
1048        let max_iter = 180;
1049        let iters_per_step = if self.favor_grade_accuracy { 2 } else { 1 };
1050        let mut new_speeds_m_per_s = Vec::with_capacity(max_iter as usize);
1051        let mut v = v0;
1052        let mut iter = 0;
1053        let mut idx = i;
1054        // dts0 = distance to next stop from d0
1055        let dts0 = self
1056            .cyc0
1057            .calc_distance_to_next_stop_from(d0 * uc::M, Some(&self.cyc0_cache))
1058            .get::<si::meter>();
1059        while v > v_brake
1060            && v >= 0.0
1061            && d <= d_max
1062            && iter < max_iter
1063            && idx < self.cyc0.speed.len()
1064        {
1065            let dt_s = (self.cyc0.time[i] - self.cyc0.time[i - 1]).get::<si::second>();
1066            let mut gr = match unique_grade {
1067                Some(g) => g,
1068                None => self.cyc0_cache.interp_grade(d + d0),
1069            };
1070            let mut k = self.calc_dvdd(v, gr);
1071            let mut v_next = v * (1.0 + 0.5 * k * dt_s) / (1.0 - 0.5 * k * dt_s);
1072            let mut vavg = 0.5 * (v + v_next);
1073            let mut dd: f64;
1074            for _ in 0..iters_per_step {
1075                k = self.calc_dvdd(vavg, gr);
1076                v_next = v * (1.0 + 0.5 * k * dt_s) / (1.0 - 0.5 * k * dt_s);
1077                vavg = 0.5 * (v + v_next);
1078                dd = vavg * dt_s;
1079                if self.favor_grade_accuracy {
1080                    gr = match unique_grade {
1081                        Some(g) => g,
1082                        None => {
1083                            let dist = (d + d0) * uc::M;
1084                            let delta_dist = dd * uc::M;
1085                            self.cyc0
1086                                .average_grade_over_range(dist, delta_dist, Some(&self.cyc0_cache))
1087                                .get::<si::ratio>()
1088                        }
1089                    };
1090                }
1091            }
1092            if k >= 0.0 && has_unique_grade {
1093                // there is no solution for coast-down -- speed will never decrease
1094                return CoastTrajectory {
1095                    found_trajectory: false,
1096                    distance_to_stop_via_coast_m: 0.0,
1097                    start_idx: 0,
1098                    speed_m_per_s: None,
1099                    distance_to_brake_m: None,
1100                };
1101            }
1102            if v_next <= v_brake {
1103                break;
1104            }
1105            vavg = 0.5 * (v + v_next);
1106            dd = vavg * dt_s;
1107            // dtb = distance to break
1108            let dtb = -0.5 * v_next * v_next / a_brake;
1109            d += dd;
1110            new_speeds_m_per_s.push(v_next);
1111            v = v_next;
1112            if d + dtb > dts0 {
1113                break;
1114            }
1115            iter += 1;
1116            idx += 1;
1117        }
1118        if iter < max_iter && idx < self.cyc0.speed.len() {
1119            let dtb = -0.5 * v * v / a_brake;
1120            let dtb_target = (dts0 - d).max(0.5 * dtb).min(2.0 * dtb);
1121            // dtsc = distance to stop via coasting
1122            let dtsc = d + dtb_target;
1123            return CoastTrajectory {
1124                found_trajectory: true,
1125                distance_to_stop_via_coast_m: dtsc,
1126                start_idx: i,
1127                speed_m_per_s: Some(new_speeds_m_per_s),
1128                distance_to_brake_m: Some(dtb_target),
1129            };
1130        }
1131        CoastTrajectory {
1132            found_trajectory: false,
1133            distance_to_stop_via_coast_m: 0.0,
1134            start_idx: 0,
1135            speed_m_per_s: None,
1136            distance_to_brake_m: None,
1137        }
1138    }
1139
1140    /// Allply the given coasting trajectory to the drive cycle.
1141    fn apply_coast_trajectory(&mut self, coast_traj: &CoastTrajectory) {
1142        if !coast_traj.found_trajectory {
1143            return;
1144        }
1145        let num_speeds = match &coast_traj.speed_m_per_s {
1146            Some(speeds_m_per_s) => {
1147                for (di, &new_speed) in speeds_m_per_s.iter().enumerate() {
1148                    let idx = coast_traj.start_idx + di;
1149                    if idx >= self.cyc0.speed.len() {
1150                        break;
1151                    }
1152                    self.cyc.speed[idx] = new_speed * uc::MPS;
1153                }
1154                speeds_m_per_s.len()
1155            }
1156            None => 0,
1157        };
1158        let (_, n) = self.cyc.modify_with_braking_trajectory(
1159            self.coast_brake_accel,
1160            coast_traj.start_idx + num_speeds,
1161            coast_traj.distance_to_brake_m.map(|d| d * uc::M),
1162        );
1163        for di in 0..(self.cyc0.speed.len() - coast_traj.start_idx) {
1164            let idx = coast_traj.start_idx + di;
1165            self.impose_coast[idx] = di < num_speeds + n;
1166        }
1167    }
1168
1169    /// Calculates the derivative dv/dd (change in speed by change in distance)
1170    /// - speed_m_per_s: the speed at which to evaluate dv/dd (m/s)
1171    /// - grade: the road grade as a decimal fraction
1172    ///
1173    /// RETURN: number, the dv/dd for these conditions
1174    pub fn calc_dvdd(&self, speed_m_per_s: f64, grade: f64) -> f64 {
1175        let v = speed_m_per_s;
1176        if v <= 0.0 {
1177            return 0.0;
1178        }
1179        let (atan_grade_sin, atan_grade_cos) = if grade == 0.0 {
1180            (0.0, 1.0)
1181        } else {
1182            let atan_grade = grade.atan();
1183            (atan_grade.sin(), atan_grade.cos())
1184        };
1185        let g = uc::ACC_GRAV.get::<si::meter_per_second_squared>();
1186        let m = self.mass.get::<si::kilogram>();
1187        let rho_cdfa = self.air_density.get::<si::kilogram_per_cubic_meter>()
1188            * self.drag_coef.get::<si::ratio>()
1189            * self.frontal_area.get::<si::square_meter>();
1190        let rrc = self.wheel_rr_coef.get::<si::ratio>();
1191        -1.0 * ((g / v) * (atan_grade_sin + rrc * atan_grade_cos)
1192            + (0.5 * rho_cdfa * (1.0 / m) * v))
1193    }
1194
1195    /// Prevent collision between the vehicle in cyc and the one in cyc0.
1196    /// If a collision will take place, reworks the cyc such that a rendezvous occurs instead.
1197    ///
1198    /// # Arguments
1199    /// - i: int, index for consideration
1200    /// - passing_tol_m: None | float, tolerance for how far we have to go past the lead vehicle to be considered "passing"
1201    ///
1202    /// RETURN: Bool, True if cyc was modified
1203    fn prevent_collisions(&mut self, i: usize, passing_tol: Option<si::Length>) -> bool {
1204        let passing_tol_m = passing_tol.map(|pt| pt.get::<si::meter>()).unwrap_or(1.0);
1205        let pass_info = PassingInfo::from(&self.cyc, &self.cyc0, i, passing_tol);
1206        if !pass_info.passing_detected {
1207            return false;
1208        }
1209        let mut best = RendezvousTrajectory {
1210            found_trajectory: false,
1211            idx: 0,
1212            n: 0,
1213            full_brake_steps: 0,
1214            jerk_m_per_s3: 0.0,
1215            accel0_m_per_s2: 0.0,
1216            accel_spread: 0.0,
1217        };
1218        let a_brake_m_per_s2 = {
1219            let result = self.coast_brake_accel.get::<si::meter_per_second_squared>();
1220            if result > 0.0 {
1221                -result
1222            } else {
1223                result
1224            }
1225        };
1226        for full_brake_steps in 0..4 {
1227            for di in 0..(self.cyc.speed.len() - i) {
1228                let idx = i + di;
1229                if !self.impose_coast[idx] {
1230                    if idx == i {
1231                        break;
1232                    } else {
1233                        continue;
1234                    }
1235                }
1236                let n = pass_info.index - idx + 1 - full_brake_steps;
1237                if n < 2 {
1238                    break;
1239                }
1240                if (idx - 1 + full_brake_steps) >= self.cyc.speed.len() {
1241                    break;
1242                }
1243                let dt = pass_info.time_step_duration.get::<si::second>();
1244                let v_start_m_per_s = self.cyc.speed[idx - 1].get::<si::meter_per_second>();
1245                let dt_full_brake = (self.cyc.time[idx - 1 + full_brake_steps]
1246                    - self.cyc.time[idx - 1])
1247                    .get::<si::second>();
1248                let dv_full_brake = dt_full_brake * a_brake_m_per_s2;
1249                let v_start_jerk_m_per_s = (v_start_m_per_s + dv_full_brake).max(0.0);
1250                let dd_full_brake = 0.5 * (v_start_m_per_s + v_start_jerk_m_per_s) * dt_full_brake;
1251                let d_start_m =
1252                    trapz_step_start_distance(&self.cyc, idx).get::<si::meter>() + dd_full_brake;
1253                let pass_distance_m = pass_info.distance.get::<si::meter>();
1254                if pass_distance_m <= d_start_m {
1255                    continue;
1256                }
1257                let jerk_trajectory = calc_constant_jerk_trajectory(
1258                    n,
1259                    d_start_m,
1260                    v_start_jerk_m_per_s,
1261                    pass_info.distance.get::<si::meter>(),
1262                    pass_info.speed.get::<si::meter_per_second>(),
1263                    dt,
1264                );
1265                let mut accels_m_per_s2 = vec![];
1266                let mut trace_accels_m_per_s2 = vec![];
1267                for ni in 0..n {
1268                    if (ni + idx + full_brake_steps) >= self.cyc.time.len() {
1269                        break;
1270                    }
1271                    accels_m_per_s2.push(accel_for_constant_jerk(
1272                        ni,
1273                        jerk_trajectory.acceleration_m_per_s2,
1274                        jerk_trajectory.jerk_m_per_s3,
1275                        jerk_trajectory.step_duration_s,
1276                    ));
1277                    let index1 = ni + idx + full_brake_steps;
1278                    let index0 = index1 - 1;
1279                    let dvi = (self.cyc.speed[index1] - self.cyc.speed[index0])
1280                        .get::<si::meter_per_second>();
1281                    let dti = (self.cyc.time[index1] - self.cyc.time[index0]).get::<si::second>();
1282                    trace_accels_m_per_s2.push(dvi / dti);
1283                }
1284                let all_sub_coast = trace_accels_m_per_s2
1285                    .iter()
1286                    .copied()
1287                    .zip(accels_m_per_s2.iter().copied())
1288                    .fold(
1289                        true,
1290                        |all_sc_flag: bool, (trace_accel, accel): (f64, f64)| {
1291                            if !all_sc_flag {
1292                                return all_sc_flag;
1293                            }
1294                            trace_accel >= accel
1295                        },
1296                    );
1297                let (min_accel_m_per_s2, max_accel_m_per_s2) = {
1298                    if !accels_m_per_s2.is_empty() {
1299                        let mut a_min = accels_m_per_s2[0];
1300                        let mut a_max = accels_m_per_s2[0];
1301                        for a in &accels_m_per_s2 {
1302                            if *a < a_min {
1303                                a_min = *a;
1304                            }
1305                            if *a > a_max {
1306                                a_max = *a;
1307                            }
1308                        }
1309                        (a_min, a_max)
1310                    } else {
1311                        (0.0, 0.0)
1312                    }
1313                };
1314                let accept = all_sub_coast;
1315                let accel_spread = (max_accel_m_per_s2 - min_accel_m_per_s2).abs();
1316                if accept && (!best.found_trajectory || accel_spread < best.accel_spread) {
1317                    best = RendezvousTrajectory {
1318                        found_trajectory: true,
1319                        idx,
1320                        n,
1321                        full_brake_steps,
1322                        jerk_m_per_s3: jerk_trajectory.jerk_m_per_s3,
1323                        accel0_m_per_s2: jerk_trajectory.acceleration_m_per_s2,
1324                        accel_spread,
1325                    };
1326                }
1327            }
1328            if best.found_trajectory {
1329                break;
1330            }
1331        }
1332        if !best.found_trajectory {
1333            let new_passing_tol_m = if passing_tol_m < 10.0 {
1334                10.0
1335            } else {
1336                passing_tol_m + 5.0
1337            };
1338            if new_passing_tol_m > 60.0 {
1339                return false;
1340            }
1341            return self.prevent_collisions(i, Some(new_passing_tol_m * uc::M));
1342        }
1343        for fbs in 0..best.full_brake_steps {
1344            if (best.idx + fbs) >= self.cyc.time.len() {
1345                break;
1346            }
1347            let dt =
1348                (self.cyc.time[best.idx + fbs] - self.cyc.time[best.idx - 1]).get::<si::second>();
1349            let dv = a_brake_m_per_s2 * dt;
1350            let v_start = self.cyc.speed[best.idx - 1].get::<si::meter_per_second>();
1351            self.cyc.speed[best.idx + fbs] = (v_start + dv).max(0.0) * uc::MPS;
1352            self.impose_coast[best.idx + fbs] = true;
1353            self.coast_delay_index[best.idx + fbs] = 0;
1354        }
1355        self.cyc.modify_by_const_jerk_trajectory(
1356            best.idx + best.full_brake_steps,
1357            best.n,
1358            best.jerk_m_per_s3 * uc::MPS3,
1359            best.accel0_m_per_s2 * uc::MPS2,
1360        );
1361        for idx in (best.idx + best.n)..self.cyc0.speed.len() {
1362            self.impose_coast[idx] = false;
1363            self.coast_delay_index[idx] = 0;
1364        }
1365        true
1366    }
1367}
1368
1369#[cfg(test)]
1370mod tests {
1371    use crate::prelude::SimDrive;
1372
1373    use super::*;
1374
1375    #[test]
1376    fn test_that_coasting_works() {
1377        let udds = crate::drive_cycle::Cycle::from_resource("udds.csv", false).unwrap();
1378        let veh = crate::vehicle::Vehicle::from_resource("2012_Ford_Fusion.yaml", false).unwrap();
1379        let mut man = Maneuver::from(&udds, &veh);
1380        man.coast_allow = true;
1381        man.coast_start_speed = 20.0 * uc::MPS;
1382        man.coast_allow_passing = true;
1383        man.apply();
1384        let udds_mod = man.cyc;
1385        assert_eq!(udds_mod.time.len(), udds.time.len());
1386        assert_eq!(udds_mod.speed.len(), udds.speed.len());
1387        assert_eq!(udds_mod.dist.len(), udds.dist.len());
1388        let mut speeds_differ = false;
1389        for idx in 0..udds.time.len() {
1390            speeds_differ = udds_mod.speed[idx] != udds.speed[idx];
1391            if speeds_differ {
1392                break;
1393            }
1394        }
1395        assert!(speeds_differ);
1396    }
1397
1398    #[test]
1399    fn test_advanced_coasting() {
1400        let udds = crate::drive_cycle::Cycle::from_resource("udds.csv", false).unwrap();
1401        let veh = crate::vehicle::Vehicle::from_resource("2012_Ford_Fusion.yaml", false).unwrap();
1402        let mut man = Maneuver::from(&udds, &veh);
1403        man.coast_allow = true;
1404        man.coast_start_speed = 0.0 * uc::MPS;
1405        man.coast_brake_start_speed = 20.0 * uc::MPH;
1406        man.coast_brake_accel = -2.5 * uc::MPS2;
1407        man.favor_grade_accuracy = false;
1408        man.coast_allow_passing = true;
1409        man.coast_max_speed = 75.0 * uc::MPH;
1410        man.coast_time_horizon_for_adjustment = 120.0 * uc::S;
1411        man.apply();
1412        let udds_mod = man.cyc;
1413        let mut sd = SimDrive::new(veh, udds_mod, None);
1414        sd.walk().unwrap();
1415    }
1416
1417    #[test]
1418    fn test_cruise() {
1419        let udds = crate::drive_cycle::Cycle::from_resource("udds.csv", false).unwrap();
1420        let vavg = udds.average_speed(true);
1421        let veh = crate::vehicle::Vehicle::from_resource("2012_Ford_Fusion.yaml", false).unwrap();
1422        let mut man = Maneuver::from(&udds, &veh);
1423        man.idm_allow = true;
1424        man.idm_desired_speed = vavg;
1425        man.idm_headway = 1.0 * uc::S;
1426        man.idm_minimum_gap = 1.0 * uc::M;
1427        man.idm_delta = 4.0;
1428        man.idm_acceleration = 1.0 * uc::MPS2;
1429        man.idm_deceleration = 2.5 * uc::MPS2;
1430        man.coast_allow = false;
1431        man.apply();
1432        let udds_mod = man.cyc;
1433        let mut sd = SimDrive::new(veh, udds_mod, None);
1434        sd.walk().unwrap();
1435    }
1436
1437    #[test]
1438    fn test_cruise_and_coast() {
1439        let udds = crate::drive_cycle::Cycle::from_resource("udds.csv", false).unwrap();
1440        let vavg = udds.average_speed(true);
1441        let veh = crate::vehicle::Vehicle::from_resource("2012_Ford_Fusion.yaml", false).unwrap();
1442        let mut man = Maneuver::from(&udds, &veh);
1443        man.idm_allow = true;
1444        man.idm_desired_speed = vavg;
1445        man.idm_headway = 1.0 * uc::S;
1446        man.idm_minimum_gap = 1.0 * uc::M;
1447        man.idm_delta = 4.0;
1448        man.idm_acceleration = 1.0 * uc::MPS2;
1449        man.idm_deceleration = 2.5 * uc::MPS2;
1450        man.coast_allow = true;
1451        man.coast_brake_start_speed = 8.9408 * uc::MPS;
1452        man.coast_brake_accel = -2.5 * uc::MPS2;
1453        man.favor_grade_accuracy = true;
1454        man.coast_allow_passing = true;
1455        man.coast_time_horizon_for_adjustment = 120.0 * uc::S;
1456        man.apply();
1457        let udds_mod = man.cyc;
1458        let mut found_coast = false;
1459        for ic in man.impose_coast {
1460            if ic {
1461                found_coast = true;
1462                break;
1463            }
1464        }
1465        assert!(found_coast);
1466        let mut sd = SimDrive::new(veh, udds_mod, None);
1467        sd.walk().unwrap();
1468    }
1469}