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}