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};
25
26impl<'a> MultipleShooting<'a, Node, 3, 3> {
27 /// Builds a multiple shooting structure assuming that the optimal trajectory is a straight line
28 /// between the start and end points. The position of the nodes will be update at each iteration
29 /// of the outer loop.
30 /// NOTE: this may cause some nodes to be below the surface of a celestial object if in low orbit
31 pub fn equidistant_nodes(
32 x0: Spacecraft,
33 xf: Orbit,
34 node_count: usize,
35 prop: &'a Propagator<SpacecraftDynamics>,
36 ) -> Result<Self, TargetingError> {
37 if node_count < 3 {
38 error!("At least three nodes are needed for a multiple shooting optimization");
39 return Err(TargetingError::UnderdeterminedProblem);
40 }
41
42 // Compute the direction of the objective
43 let mut direction = xf.radius_km - x0.orbit.radius_km;
44 if direction.norm() < f64::EPSILON {
45 return Err(TargetingError::TargetsTooClose);
46 }
47 let distance_increment = direction.norm() / (node_count as f64);
48 let duration_increment = (xf.epoch - x0.epoch()) / (node_count as f64);
49 direction /= direction.norm();
50
51 // Build each node successively (includes xf)
52 let mut nodes = Vec::with_capacity(node_count + 1);
53 let mut prev_node_radius = x0.orbit.radius_km;
54 let mut prev_node_epoch = x0.epoch();
55
56 for _ in 0..node_count {
57 // Compute the position we want.
58 let this_node = prev_node_radius + distance_increment * direction;
59 let this_epoch = prev_node_epoch + duration_increment;
60 nodes.push(Node {
61 x: this_node[0],
62 y: this_node[1],
63 z: this_node[2],
64 vmag: 0.0,
65 frame: x0.orbit.frame,
66 epoch: this_epoch,
67 });
68 prev_node_radius = this_node;
69 prev_node_epoch = this_epoch;
70 }
71 Ok(Self {
72 prop,
73 targets: nodes,
74 x0,
75 xf,
76 current_iteration: 0,
77 max_iterations: 50,
78 improvement_threshold: 0.01,
79 variables: [
80 Vary::VelocityX.into(),
81 Vary::VelocityY.into(),
82 Vary::VelocityZ.into(),
83 ],
84 all_dvs: Vec::with_capacity(node_count),
85 })
86 }
87}