closedloop_multi_oe_ruggiero.rs 13.5 KB
Newer Older
1 2 3 4 5 6 7 8
extern crate hifitime;
extern crate nalgebra as na;
extern crate nyx_space as nyx;

use self::hifitime::{Epoch, SECONDS_PER_DAY};
use self::nyx::celestia::{bodies, Cosm, Geoid, State};
use self::nyx::dynamics::celestial::CelestialDynamics;
use self::nyx::dynamics::propulsion::{Propulsion, Thruster};
Christopher Rabotin's avatar
Christopher Rabotin committed
9
use self::nyx::dynamics::spacecraft::{Spacecraft, SpacecraftState};
10 11
use self::nyx::dynamics::thrustctrl::{Achieve, Ruggiero};
use self::nyx::dynamics::Dynamics;
12
use self::nyx::propagators::events::{EventKind, EventTrackers, OrbitalEvent, SCEvent};
13 14 15 16 17 18 19 20 21 22 23 24 25 26 27
use self::nyx::propagators::{PropOpts, Propagator, RK4Fixed};

/// NOTE: Herein shows the difference between the QLaw and Ruggiero (and other control laws).
/// The Ruggiero control law takes quite some longer to converge than the QLaw.

#[test]
fn qlaw_as_ruggiero_case_a() {
    // Source: AAS-2004-5089

    let mut cosm = Cosm::from_xb("./de438s");
    cosm.mut_gm_for_geoid_id(bodies::EARTH, 398_600.433);
    let earth = cosm.geoid_from_id(bodies::EARTH);

    let start_time = Epoch::from_gregorian_tai_at_midnight(2020, 1, 1);

28 29
    let orbit =
        State::<Geoid>::from_keplerian(7000.0, 0.01, 0.05, 0.0, 0.0, 1.0, start_time, earth);
30 31 32 33

    let prop_time = 39.91 * SECONDS_PER_DAY;

    // Define the dynamics
34
    let dynamics = CelestialDynamics::two_body(orbit);
35 36 37 38 39 40 41 42 43 44

    // Define the thruster
    let lowt = vec![Thruster {
        thrust: 1.0,
        isp: 3100.0,
    }];

    // Define the objectives
    let objectives = vec![
        Achieve::Sma {
45
            target: 42000.0,
46 47 48 49 50 51 52 53
            tol: 1.0,
        },
        Achieve::Ecc {
            target: 0.01,
            tol: 5e-5,
        },
    ];

54
    // Track these events
55 56 57 58
    let tracker = EventTrackers::from_events(vec![
        SCEvent::orbital(OrbitalEvent::new(EventKind::Sma(42000.0))),
        SCEvent::orbital(OrbitalEvent::new(EventKind::Ecc(0.01))),
    ]);
59

60
    let ruggiero = Ruggiero::new(objectives, orbit);
61 62 63 64

    let dry_mass = 1.0;
    let fuel_mass = 299.0;

65
    let prop_subsys = Propulsion::new(ruggiero, lowt, true);
66

67
    let mut sc = Spacecraft::with_prop(dynamics, prop_subsys, dry_mass, fuel_mass);
68 69
    println!("{:o}", orbit);

70
    let mut prop = Propagator::new::<RK4Fixed>(&mut sc, &PropOpts::with_fixed_step(10.0));
71
    prop.event_trackers = tracker;
72
    prop.until_time_elapsed(prop_time);
73
    println!("{}", prop.event_trackers);
74
    let final_state = prop.dynamics.celestial.state();
75
    let fuel_usage = fuel_mass - sc.fuel_mass;
76 77 78
    println!("{:o}", final_state);
    println!("fuel usage: {:.3} kg", fuel_usage);

79 80 81 82
    assert!(
        sc.prop.as_ref().unwrap().ctrl.achieved(&final_state),
        "objective not achieved"
    );
83

84
    assert!((fuel_usage - 93.449).abs() < 1.0);
85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100
}

#[test]
fn qlaw_as_ruggiero_case_b() {
    // Source: AAS-2004-5089
    let cosm = Cosm::from_xb("./de438s");
    let earth = cosm.geoid_from_id(bodies::EARTH);

    let start_time = Epoch::from_gregorian_tai_at_midnight(2020, 1, 1);

    let orbit =
        State::<Geoid>::from_keplerian(24505.9, 0.725, 7.05, 0.0, 0.0, 0.0, start_time, earth);

    let prop_time = 160.0 * SECONDS_PER_DAY;

    // Define the dynamics
101
    let dynamics = CelestialDynamics::two_body(orbit);
102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124

    // Define the thruster
    let lowt = vec![Thruster {
        thrust: 0.350,
        isp: 2000.0,
    }];

    // Define the objectives
    let objectives = vec![
        Achieve::Sma {
            target: 42165.0,
            tol: 20.0,
        },
        Achieve::Ecc {
            target: 0.001,
            tol: 5e-5,
        },
        Achieve::Inc {
            target: 0.05,
            tol: 5e-3,
        },
    ];

125
    let ruggiero = Ruggiero::new(objectives, orbit);
126 127 128 129

    let fuel_mass = 1999.9;
    let dry_mass = 0.1;

130
    let prop_subsys = Propulsion::new(ruggiero, lowt, true);
131

132
    let mut sc = Spacecraft::with_prop(dynamics, prop_subsys, dry_mass, fuel_mass);
133 134
    println!("{:o}", orbit);

135
    let mut prop = Propagator::new::<RK4Fixed>(&mut sc, &PropOpts::with_fixed_step(10.0));
136 137 138
    prop.until_time_elapsed(prop_time);

    let final_state = prop.dynamics.celestial.state();
139
    let fuel_usage = fuel_mass - sc.fuel_mass;
140 141 142
    println!("{:o}", final_state);
    println!("fuel usage: {:.3} kg", fuel_usage);

143 144 145 146
    assert!(
        sc.prop.as_ref().unwrap().ctrl.achieved(&final_state),
        "objective not achieved"
    );
147

148
    assert!((fuel_usage - 223.515).abs() < 1.0);
149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164
}

#[test]
fn qlaw_as_ruggiero_case_c() {
    // Source: AAS-2004-5089
    let cosm = Cosm::from_xb("./de438s");
    let earth = cosm.geoid_from_id(bodies::EARTH);

    let start_time = Epoch::from_gregorian_tai_at_midnight(2020, 1, 1);

    let orbit =
        State::<Geoid>::from_keplerian(9222.7, 0.2, 0.573, 0.0, 0.0, 0.0, start_time, earth);

    let prop_time = 3.0 * SECONDS_PER_DAY;

    // Define the dynamics
165
    let dynamics = CelestialDynamics::two_body(orbit);
166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184

    // Define the thruster
    let lowt = vec![Thruster {
        thrust: 9.3,
        isp: 3100.0,
    }];

    // Define the objectives
    let objectives = vec![
        Achieve::Sma {
            target: 30000.0,
            tol: 1.0,
        },
        Achieve::Ecc {
            target: 0.7,
            tol: 5e-5,
        },
    ];

185
    let ruggiero = Ruggiero::new(objectives, orbit);
186 187 188 189

    let fuel_mass = 299.9;
    let dry_mass = 0.1;

190
    let prop_subsys = Propulsion::new(ruggiero, lowt, true);
191

192
    let mut sc = Spacecraft::with_prop(dynamics, prop_subsys, dry_mass, fuel_mass);
193 194
    println!("{:o}", orbit);

195
    let mut prop = Propagator::new::<RK4Fixed>(&mut sc, &PropOpts::with_fixed_step(10.0));
196 197 198
    prop.until_time_elapsed(prop_time);

    let final_state = prop.dynamics.celestial.state();
199
    let fuel_usage = fuel_mass - sc.fuel_mass;
200 201 202
    println!("{:o}", final_state);
    println!("fuel usage: {:.3} kg", fuel_usage);

203 204 205 206
    assert!(
        sc.prop.as_ref().unwrap().ctrl.achieved(&final_state),
        "objective not achieved"
    );
207

208
    assert!((fuel_usage - 41.742).abs() < 1.0);
209 210 211 212 213
}

#[test]
#[ignore]
fn qlaw_as_ruggiero_case_d() {
214
    // Broken: https://gitlab.com/chrisrabotin/nyx/issues/103
215 216 217 218 219 220 221 222 223 224 225 226
    // Source: AAS-2004-5089
    let cosm = Cosm::from_xb("./de438s");
    let earth = cosm.geoid_from_id(bodies::EARTH);

    let start_time = Epoch::from_gregorian_tai_at_midnight(2020, 1, 1);

    let orbit =
        State::<Geoid>::from_keplerian(24505.9, 0.725, 0.06, 0.0, 0.0, 0.0, start_time, earth);

    let prop_time = 113.0 * SECONDS_PER_DAY;

    // Define the dynamics
227
    let dynamics = CelestialDynamics::two_body(orbit);
228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254

    // Define the thruster
    let lowt = vec![Thruster {
        thrust: 89e-3,
        isp: 1650.0,
    }];

    // Define the objectives
    let objectives = vec![
        Achieve::Sma {
            target: 26500.0,
            tol: 1.0,
        },
        Achieve::Inc {
            target: 116.0,
            tol: 5e-3,
        },
        Achieve::Ecc {
            target: 0.7,
            tol: 5e-5,
        },
        Achieve::Raan {
            target: 360.0 - 90.0,
            tol: 5e-3,
        },
    ];

255
    let ruggiero = Ruggiero::new(objectives, orbit);
256 257 258 259

    let fuel_mass = 67.0;
    let dry_mass = 300.0;

260
    let prop_subsys = Propulsion::new(ruggiero, lowt, true);
261

262
    let mut sc = Spacecraft::with_prop(dynamics, prop_subsys, dry_mass, fuel_mass);
263 264
    println!("{:o}", orbit);

265
    let mut prop = Propagator::new::<RK4Fixed>(&mut sc, &PropOpts::with_fixed_step(10.0));
266 267 268
    prop.until_time_elapsed(prop_time);

    let final_state = prop.dynamics.celestial.state();
269
    let fuel_usage = fuel_mass - sc.fuel_mass;
270 271 272
    println!("{:o}", final_state);
    println!("fuel usage: {:.3} kg", fuel_usage);

273 274 275 276
    assert!(
        sc.prop.as_ref().unwrap().ctrl.achieved(&final_state),
        "objective not achieved"
    );
277 278 279 280 281 282 283

    assert!((fuel_usage - 23.0).abs() < 1.0);
}

#[test]
#[ignore]
fn qlaw_as_ruggiero_case_e() {
284
    // Broken: https://gitlab.com/chrisrabotin/nyx/issues/103
285 286 287 288 289 290 291 292 293 294 295 296
    // Source: AAS-2004-5089
    let cosm = Cosm::from_xb("./de438s");
    let earth = cosm.geoid_from_id(bodies::EARTH);

    let start_time = Epoch::from_gregorian_tai_at_midnight(2020, 1, 1);

    let orbit =
        State::<Geoid>::from_keplerian(24505.9, 0.725, 0.06, 0.0, 0.0, 0.0, start_time, earth);

    let prop_time = 400.0 * SECONDS_PER_DAY;

    // Define the dynamics
297
    let dynamics = CelestialDynamics::two_body(orbit);
298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328

    // Define the thruster
    let lowt = vec![Thruster {
        thrust: 89e-3,
        isp: 1650.0,
    }];

    // Define the objectives
    let objectives = vec![
        Achieve::Sma {
            target: 26500.0,
            tol: 1.0,
        },
        Achieve::Ecc {
            target: 0.7,
            tol: 5e-5,
        },
        Achieve::Inc {
            target: 116.0,
            tol: 5e-3,
        },
        Achieve::Raan {
            target: 270.0,
            tol: 5e-3,
        },
        Achieve::Aop {
            target: 180.0,
            tol: 5e-3,
        },
    ];

329
    let ruggiero = Ruggiero::new(objectives, orbit);
330 331 332 333

    let fuel_mass = 1999.9;
    let dry_mass = 0.1;

334
    let prop_subsys = Propulsion::new(ruggiero, lowt, true);
335

336
    let mut sc = Spacecraft::with_prop(dynamics, prop_subsys, dry_mass, fuel_mass);
337 338
    println!("{:o}", orbit);

339
    let mut prop = Propagator::new::<RK4Fixed>(&mut sc, &PropOpts::with_fixed_step(10.0));
340 341 342
    prop.until_time_elapsed(prop_time);

    let final_state = prop.dynamics.celestial.state();
343
    let fuel_usage = fuel_mass - sc.fuel_mass;
344 345 346
    println!("{:o}", final_state);
    println!("fuel usage: {:.3} kg", fuel_usage);

347 348 349 350
    assert!(
        sc.prop.as_ref().unwrap().ctrl.achieved(&final_state),
        "objective not achieved"
    );
351 352 353 354 355 356 357 358 359 360 361 362

    assert!((fuel_usage - 23.0).abs() < 1.0);
}

#[test]
fn qlaw_as_ruggiero_case_f() {
    // Source: AAS-2004-5089
    /*
        NOTE: Due to how lifetime of variables work in Rust, we need to define all of the
        components of a spacecraft before defining the spacecraft itself.
    */

Christopher Rabotin's avatar
Christopher Rabotin committed
363 364 365 366
    // We'll export this trajectory as a POC. Adding the needed crates here
    extern crate csv;
    use std::sync::mpsc;
    use std::sync::mpsc::{Receiver, Sender};
367
    use std::thread;
Christopher Rabotin's avatar
Christopher Rabotin committed
368

369 370 371 372 373 374 375 376 377 378 379
    let cosm = Cosm::from_xb("./de438s");
    let earth = cosm.geoid_from_id(bodies::EARTH);

    let start_time = Epoch::from_gregorian_tai_at_midnight(2020, 1, 1);

    let orbit =
        State::<Geoid>::from_keplerian(15378.0, 0.01, 98.7, 0.0, 0.0, 0.0, start_time, earth);

    let prop_time = 30.0 * SECONDS_PER_DAY;

    // Define the dynamics
380
    let dynamics = CelestialDynamics::two_body(orbit);
381 382 383 384 385 386 387 388 389 390 391 392 393

    // Define the thruster
    let lowt = vec![Thruster {
        thrust: 89e-3,
        isp: 1650.0,
    }];

    // Define the objectives
    let objectives = vec![Achieve::Ecc {
        target: 0.15,
        tol: 1e-5,
    }];

394
    let ruggiero = Ruggiero::new(objectives, orbit);
395 396 397 398

    let fuel_mass = 67.0;
    let dry_mass = 300.0;

399
    let prop_subsys = Propulsion::new(ruggiero, lowt, true);
400

401
    let mut sc = Spacecraft::with_prop(dynamics, prop_subsys, dry_mass, fuel_mass);
402 403
    println!("{:o}", orbit);

Christopher Rabotin's avatar
Christopher Rabotin committed
404 405 406 407 408 409 410 411 412 413 414 415
    let (tx, rx): (Sender<SpacecraftState>, Receiver<SpacecraftState>) = mpsc::channel();

    // Set up the writing channel
    thread::spawn(move || {
        let mut wtr = csv::Writer::from_path("./rugg_case_f.csv").expect("could not create file");
        while let Ok(rx_state) = rx.recv() {
            // Serialize only the orbital state
            wtr.serialize(rx_state.orbit)
                .expect("could not serialize state");
        }
    });

416
    let mut prop = Propagator::new::<RK4Fixed>(&mut sc, &PropOpts::with_fixed_step(10.0));
Christopher Rabotin's avatar
Christopher Rabotin committed
417
    prop.tx_chan = Some(&tx);
418 419 420
    prop.until_time_elapsed(prop_time);

    let final_state = prop.dynamics.celestial.state();
421
    let fuel_usage = fuel_mass - sc.fuel_mass;
422 423 424
    println!("{:o}", final_state);
    println!("fuel usage: {:.3} kg", fuel_usage);

425 426 427 428
    assert!(
        sc.prop.as_ref().unwrap().ctrl.achieved(&final_state),
        "objective not achieved"
    );
429

430
    assert!((fuel_usage - 10.378).abs() < 1.0);
431 432 433 434 435 436 437 438 439 440 441 442 443
}

#[test]
fn ruggiero_iepc_2011_102() {
    // Source: IEPC 2011 102
    let cosm = Cosm::from_xb("./de438s");
    let earth = cosm.geoid_from_id(bodies::EARTH);

    let start_time = Epoch::from_gregorian_tai_at_midnight(2020, 1, 1);

    let orbit =
        State::<Geoid>::from_keplerian(24396.0, 0.7283, 7.0, 1.0, 1.0, 1.0, start_time, earth);

444
    let prop_time = 105.0 * SECONDS_PER_DAY;
445 446

    // Define the dynamics
447
    let dynamics = CelestialDynamics::two_body(orbit);
448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465 466 467 468 469 470

    // Define the thruster
    let lowt = vec![Thruster {
        thrust: 89e-3,
        isp: 1650.0,
    }];

    // Define the objectives
    let objectives = vec![
        Achieve::Sma {
            target: 42164.0,
            tol: 20.0,
        },
        Achieve::Inc {
            target: 0.001,
            tol: 5e-3,
        },
        Achieve::Ecc {
            target: 0.001,
            tol: 5e-5,
        },
    ];

471
    let ruggiero = Ruggiero::new(objectives, orbit);
472 473 474 475

    let fuel_mass = 67.0;
    let dry_mass = 300.0;

476
    let prop_subsys = Propulsion::new(ruggiero, lowt, true);
477

478
    let mut sc = Spacecraft::with_prop(dynamics, prop_subsys, dry_mass, fuel_mass);
479 480
    println!("{:o}", orbit);

481
    let mut prop = Propagator::new::<RK4Fixed>(&mut sc, &PropOpts::with_fixed_step(10.0));
482 483 484
    prop.until_time_elapsed(prop_time);

    let final_state = prop.dynamics.celestial.state();
485
    let fuel_usage = fuel_mass - sc.fuel_mass;
486 487 488
    println!("{:o}", final_state);
    println!("fuel usage: {:.3} kg", fuel_usage);

489 490 491 492
    assert!(
        sc.prop.as_ref().unwrap().ctrl.achieved(&final_state),
        "objective not achieved"
    );
493 494 495 496

    // WARNING: Paper claims this can be done with only 49kg of fuel.
    assert!((fuel_usage - 49.0).abs() < 1.0);
}