03_geo_raise_optim/
raise_optim.rs

1#![doc = include_str!("./README.md")]
2extern crate log;
3extern crate nyx_space as nyx;
4extern crate pretty_env_logger as pel;
5
6use anise::{
7    almanac::{Almanac, metaload::MetaFile},
8    constants::{
9        celestial_objects::{MOON, SUN},
10        frames::{EARTH_J2000, IAU_EARTH_FRAME},
11    },
12};
13use hifitime::{Epoch, TimeUnits, Unit};
14use log::info;
15use nyx::{
16    Spacecraft,
17    cosmic::{GuidanceMode, Mass, MetaAlmanac, Orbit, SRPData},
18    dynamics::{
19        GravityField, OrbitalDynamics, SolarPressure, SpacecraftDynamics,
20        guidance::{Ruggiero, Thruster},
21    },
22    io::gravity::GravityFieldData,
23    md::prelude::{Objective, OrbitalElement, StateParameter},
24    propagators::{ErrorControl, IntegratorOptions, Propagator},
25};
26use radiate::*;
27use std::{error::Error, sync::Arc};
28
29// Shared state struct for the fitness evaluation to avoid reading files thousands of times
30struct SharedState {
31    almanac: Arc<Almanac>,
32    harmonics: Arc<GravityField>,
33    srp_dyn: Arc<SolarPressure>,
34}
35
36impl SharedState {
37    fn new() -> Result<Self, Box<dyn Error>> {
38        let almanac = Arc::new(MetaAlmanac::latest().map_err(Box::new)?);
39
40        let mut jgm3_meta = MetaFile {
41            uri: "http://public-data.nyxspace.com/nyx/models/JGM3.cof.gz".to_string(),
42            crc32: Some(0xF446F027),
43        };
44        jgm3_meta.process(true)?;
45
46        let harmonics = GravityField::from_stor(
47            almanac.frame_info(IAU_EARTH_FRAME)?,
48            GravityFieldData::from_cof(&jgm3_meta.uri, 4, 4, true)?,
49        );
50        let srp_dyn = SolarPressure::default_flux(EARTH_J2000, almanac.clone())?;
51
52        Ok(Self {
53            almanac,
54            harmonics,
55            srp_dyn,
56        })
57    }
58}
59
60fn main() -> Result<(), Box<dyn Error>> {
61    pel::init();
62
63    /* let (prop_usage_kg, penalty) = evaluate_weights(
64        &[0.22033301, 0.8512096, 0.49421895],
65        60.0,
66        Arc::new(SharedState::new()?),
67    )
68    .unwrap();
69
70    println!("Best weight prop usage = {prop_usage_kg:.3} kg \t penalty = {penalty:.3}"); */
71
72    // Set up shared state (read large files only once!)
73    let shared_state = Arc::new(SharedState::new()?);
74
75    // Set up the genetic algorithm optimization
76    let codec = FloatCodec::vector(3, 0.1_f32..1.0_f32); // 3 weights for SMA, Ecc, Inc
77    let problem = EngineProblem {
78        objective: radiate::Objective::Multi(vec![Optimize::Minimize, Optimize::Minimize]), // NSGA2 Multi Objective
79        codec: Arc::new(codec),
80        fitness_fn: Some(Arc::new(move |weights: Vec<f32>| {
81            // Full 60 days propagation for evaluating the actual performance, but running fast due to shared state
82            let (prop_usage, penalty) =
83                evaluate_weights(&weights, 60.0, shared_state.clone()).unwrap_or((1e6, 1e6));
84            Score::from(vec![prop_usage as f32, penalty as f32])
85        })),
86        raw_fitness_fn: None,
87    };
88
89    let mut engine = GeneticEngine::<FloatChromosome<f32>, Vec<f32>>::builder()
90        .population_size(20)
91        .parallel()
92        .multi_objective(vec![Optimize::Minimize, Optimize::Minimize])
93        .problem(problem)
94        .survivor_selector(NSGA2Selector::new())
95        .build();
96
97    // Wrap the engine with the UI
98    let final_generation = engine.run(|generation: &Generation<FloatChromosome<f32>, Vec<f32>>| {
99        let scores = generation.score().as_slice();
100        println!(
101            "[ {:?} ]: Best Score: Prop usage {:.3} kg, Penalty {:.3}",
102            generation.index(),
103            scores[0],
104            scores[1]
105        );
106        generation.index() >= 5
107    });
108
109    let best_weights = final_generation
110        .value()
111        .iter()
112        .map(|w| format!("W: = {w}"))
113        .collect::<Vec<String>>()
114        .join(", ");
115    let best_score = final_generation
116        .score()
117        .iter()
118        .enumerate()
119        .map(|(i, w)| format!("S[{i}]: = {w}"))
120        .collect::<Vec<String>>()
121        .join(", ");
122    println!("Optimization finished. Best weights: [{best_weights}] -> Best score: [{best_score}]");
123
124    // Evaluate these weights.
125    let best_weights: Vec<f32> = final_generation.value().to_vec();
126
127    let (prop_usage_kg, penalty) =
128        evaluate_weights(&best_weights, 60.0, Arc::new(SharedState::new()?)).unwrap();
129
130    println!("Best weight prop usage = {prop_usage_kg:.3} kg \t penalty = {penalty:.3}");
131
132    Ok(())
133}
134
135fn evaluate_weights(
136    weights: &[f32],
137    prop_time_days: f64,
138    state: Arc<SharedState>,
139) -> Result<(f64, f64), Box<dyn Error>> {
140    let ηthresholds: Vec<f64> = weights.iter().map(|w| *w as f64).collect();
141
142    let eme2k = state.almanac.frame_info(EARTH_J2000).unwrap();
143    let epoch = Epoch::from_gregorian_utc_hms(2024, 2, 29, 12, 13, 14);
144
145    let orbit = Orbit::keplerian(24505.9, 0.725, 7.05, 0.0, 0.0, 0.0, epoch, eme2k);
146
147    let sc = Spacecraft::builder()
148        .orbit(orbit)
149        .mass(Mass::from_dry_and_prop_masses(1000.0, 1000.0))
150        .srp(SRPData::from_area(3.0 * 6.0))
151        .thruster(Thruster {
152            isp_s: 4435.0,
153            thrust_N: 0.472,
154        })
155        .mode(GuidanceMode::Thrust)
156        .build();
157
158    let prop_time = prop_time_days * Unit::Day;
159
160    let objectives = &[
161        Objective::within_tolerance(
162            StateParameter::Element(OrbitalElement::SemiMajorAxis),
163            30_000.0,
164            20.0,
165        ),
166        Objective::within_tolerance(
167            StateParameter::Element(OrbitalElement::Eccentricity),
168            0.001,
169            5e-5,
170        ),
171        Objective::within_tolerance(
172            StateParameter::Element(OrbitalElement::Inclination),
173            0.05,
174            1e-2,
175        ),
176    ];
177
178    let ctrl = Ruggiero::from_ηthresholds(objectives, &ηthresholds, sc)?;
179
180    let mut orbital_dyn = OrbitalDynamics::point_masses(vec![MOON, SUN]);
181    orbital_dyn.accel_models.push(state.harmonics.clone());
182
183    let sc_dynamics = SpacecraftDynamics::from_model(orbital_dyn, state.srp_dyn.clone())
184        .with_guidance_law(ctrl.clone());
185
186    let (final_state, _traj) = Propagator::rk89(
187        sc_dynamics.clone(),
188        IntegratorOptions::builder()
189            .min_step(10.0_f64.seconds())
190            .tolerance(1e-8)
191            .error_ctrl(ErrorControl::RSSCartesianStep)
192            .build(),
193    )
194    .with(sc, state.almanac.clone())
195    .for_duration_with_traj(prop_time)?;
196
197    let prop_usage = sc.mass.prop_mass_kg - final_state.mass.prop_mass_kg;
198
199    let mut penalty = 0.0;
200    for obj in objectives {
201        let (achieved, error) = obj.assess(&final_state)?;
202        if !achieved {
203            penalty += error.abs();
204        }
205        info!("{obj} error: {error:.3}, achieved? {achieved}");
206    }
207
208    info!("{ηthresholds:?} -> {prop_usage:.3} kg\tpenalty = {penalty:.3}");
209
210    Ok((prop_usage, penalty * 1000.0))
211}
03_geo_raise_optim/raise_optim.rs

03_geo_raise_optim/
raise_optim.rs
