nyx_space/md/opti/multipleshooting/
equidistant_heuristic.rs

1/*
2    Nyx, blazing fast astrodynamics
3    Copyright (C) 2018-onwards Christopher Rabotin <christopher.rabotin@gmail.com>
4
5    This program is free software: you can redistribute it and/or modify
6    it under the terms of the GNU Affero General Public License as published
7    by the Free Software Foundation, either version 3 of the License, or
8    (at your option) any later version.
9
10    This program is distributed in the hope that it will be useful,
11    but WITHOUT ANY WARRANTY; without even the implied warranty of
12    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
13    GNU Affero General Public License for more details.
14
15    You should have received a copy of the GNU Affero General Public License
16    along with this program.  If not, see <https://www.gnu.org/licenses/>.
17*/
18
19use super::ctrlnodes::Node;
20use super::multishoot::MultipleShooting;
21pub use super::CostFunction;
22use crate::errors::TargetingError;
23use crate::md::prelude::*;
24use crate::{Orbit, Spacecraft};
25use log::error;
26
27impl<'a> MultipleShooting<'a, Node, 3, 3> {
28    /// Builds a multiple shooting structure assuming that the optimal trajectory is a straight line
29    /// between the start and end points. The position of the nodes will be update at each iteration
30    /// of the outer loop.
31    /// NOTE: this may cause some nodes to be below the surface of a celestial object if in low orbit
32    pub fn equidistant_nodes(
33        x0: Spacecraft,
34        xf: Orbit,
35        node_count: usize,
36        prop: &'a Propagator<SpacecraftDynamics>,
37    ) -> Result<Self, TargetingError> {
38        if node_count < 3 {
39            error!("At least three nodes are needed for a multiple shooting optimization");
40            return Err(TargetingError::UnderdeterminedProblem);
41        }
42
43        // Compute the direction of the objective
44        let mut direction = xf.radius_km - x0.orbit.radius_km;
45        if direction.norm() < f64::EPSILON {
46            return Err(TargetingError::TargetsTooClose);
47        }
48        let distance_increment = direction.norm() / (node_count as f64);
49        let duration_increment = (xf.epoch - x0.epoch()) / (node_count as f64);
50        direction /= direction.norm();
51
52        // Build each node successively (includes xf)
53        let mut nodes = Vec::with_capacity(node_count + 1);
54        let mut prev_node_radius = x0.orbit.radius_km;
55        let mut prev_node_epoch = x0.epoch();
56
57        for _ in 0..node_count {
58            // Compute the position we want.
59            let this_node = prev_node_radius + distance_increment * direction;
60            let this_epoch = prev_node_epoch + duration_increment;
61            nodes.push(Node {
62                x: this_node[0],
63                y: this_node[1],
64                z: this_node[2],
65                vmag: 0.0,
66                frame: x0.orbit.frame,
67                epoch: this_epoch,
68            });
69            prev_node_radius = this_node;
70            prev_node_epoch = this_epoch;
71        }
72        Ok(Self {
73            prop,
74            targets: nodes,
75            x0,
76            xf,
77            current_iteration: 0,
78            max_iterations: 50,
79            improvement_threshold: 0.01,
80            variables: [
81                Vary::VelocityX.into(),
82                Vary::VelocityY.into(),
83                Vary::VelocityZ.into(),
84            ],
85            all_dvs: Vec::with_capacity(node_count),
86        })
87    }
88}