1use anyhow::Result;
2use serde::{Deserialize, Serialize};
3
4use crate::config::SimulationConfig;
5use crate::metrics::{RunSummary, derived_metrics, summarize};
6use crate::model::step_state;
7use crate::monitor::{
8 FigureMetadata, StabilitySummary, admissible_status, build_figure_metadata, lyapunov_value,
9 normalized_authority_utilization, reduced_response_target, stability_summary,
10};
11use crate::scenarios::sample_control;
12use crate::state::{
13 DerivedMetricRecord, EventLatch, EventRecord, LimbBufferRecord, StepDiagnostics, SystemState,
14 TimeSeriesRecord,
15};
16
17#[derive(Debug, Clone, Serialize, Deserialize)]
18pub struct SimulationResult {
19 pub config: SimulationConfig,
20 pub time_series: Vec<TimeSeriesRecord>,
21 pub limb_buffers: Vec<LimbBufferRecord>,
22 pub events: Vec<EventRecord>,
23 pub summary: RunSummary,
24 pub derived_metrics: Vec<DerivedMetricRecord>,
25 pub stability_summary: StabilitySummary,
26 pub figure_metadata: FigureMetadata,
27}
28
29pub fn simulate(config: SimulationConfig) -> Result<SimulationResult> {
30 config.validate()?;
31
32 let mut state = SystemState::new(
33 config.model.pulse_energy_initial_j,
34 config.model.thermal_initial_k,
35 0.0,
36 0.0,
37 config.model.local_buffer_initial_j,
38 );
39 let mut time_series = Vec::new();
40 let mut limb_buffers = Vec::new();
41 let mut events = Vec::new();
42 let mut latch = EventLatch::default();
43 let mut previous_segment = "initial".to_string();
44 let mut previous_lyapunov = None;
45
46 push_records(
47 &config,
48 &mut time_series,
49 &mut limb_buffers,
50 &state,
51 &StepDiagnostics::zero(),
52 &mut previous_lyapunov,
53 );
54
55 while state.time_s < config.solver.duration_s - 1.0e-12 {
56 let remaining = config.solver.duration_s - state.time_s;
57 let dt_s = remaining.min(config.solver.dt_s);
58 let input = sample_control(&config.scenario, config.seed, state.time_s);
59 let outcome = step_state(&config.model, &state, &input, dt_s);
60
61 if outcome.diagnostics.active_segment != previous_segment {
62 events.push(EventRecord {
63 time_s: outcome.next_state.time_s,
64 event_type: "segment_transition".to_string(),
65 severity: "info".to_string(),
66 message: format!(
67 "Transitioned from '{}' to '{}'",
68 previous_segment, outcome.diagnostics.active_segment
69 ),
70 value: outcome.diagnostics.command_fraction,
71 threshold: 0.0,
72 });
73 previous_segment = outcome.diagnostics.active_segment.clone();
74 }
75
76 update_constraint_events(
77 &config,
78 &outcome.next_state,
79 &outcome.diagnostics,
80 &mut latch,
81 &mut events,
82 );
83
84 state = outcome.next_state;
85 push_records(
86 &config,
87 &mut time_series,
88 &mut limb_buffers,
89 &state,
90 &outcome.diagnostics,
91 &mut previous_lyapunov,
92 );
93 }
94
95 let summary = summarize(&config, &time_series, &limb_buffers, &events);
96 let derived_metrics = derived_metrics(&summary);
97 let stability_summary = stability_summary(
98 &time_series
99 .iter()
100 .map(|record| record.lyapunov_v)
101 .collect::<Vec<_>>(),
102 &time_series
103 .iter()
104 .map(|record| record.lyapunov_dv_dt)
105 .collect::<Vec<_>>(),
106 &time_series
107 .iter()
108 .map(|record| record.time_s)
109 .collect::<Vec<_>>(),
110 );
111 let figure_metadata = build_figure_metadata(&config);
112
113 Ok(SimulationResult {
114 config,
115 time_series,
116 limb_buffers,
117 events,
118 summary,
119 derived_metrics,
120 stability_summary,
121 figure_metadata,
122 })
123}
124
125fn push_records(
126 config: &SimulationConfig,
127 time_series: &mut Vec<TimeSeriesRecord>,
128 limb_buffers: &mut Vec<LimbBufferRecord>,
129 state: &SystemState,
130 diagnostics: &StepDiagnostics,
131 previous_lyapunov: &mut Option<(f64, f64)>,
132) {
133 let reduced_response_target_y_m =
134 reduced_response_target(diagnostics.command_fraction, &config.model);
135 let reduced_response_target_rate_mps = 0.0;
136 let reduced_response_error_m = state.y_m - reduced_response_target_y_m;
137 let reduced_response_error_rate_mps = state.v_mps - reduced_response_target_rate_mps;
138 let lyapunov_v = lyapunov_value(
139 config.model.mechanical_mass_kg,
140 diagnostics.stiffness,
141 reduced_response_error_m,
142 reduced_response_error_rate_mps,
143 );
144 let lyapunov_dv_dt = previous_lyapunov
145 .map(|(previous_time_s, previous_v)| {
146 let dt_s = (state.time_s - previous_time_s).max(1.0e-9);
147 (lyapunov_v - previous_v) / dt_s
148 })
149 .unwrap_or(0.0);
150 *previous_lyapunov = Some((state.time_s, lyapunov_v));
151
152 let authority_utilization = normalized_authority_utilization(
153 &config.model,
154 diagnostics.gain,
155 diagnostics.command_fraction,
156 diagnostics.delivered_ratio,
157 );
158 let admissible = admissible_status(
159 &config.model,
160 state.ep_j,
161 state.temperature_k,
162 state.y_m,
163 state.v_mps,
164 );
165
166 time_series.push(TimeSeriesRecord {
167 time_s: state.time_s,
168 ep_j: state.ep_j,
169 ep_gj: state.ep_j / 1.0e9,
170 temperature_k: state.temperature_k,
171 temperature_c: state.temperature_k - 273.15,
172 y_m: state.y_m,
173 v_mps: state.v_mps,
174 command_fraction: diagnostics.command_fraction,
175 active_segment: diagnostics.active_segment.clone(),
176 requested_actuator_power_w: diagnostics.requested_actuator_power_w,
177 delivered_actuator_power_w: diagnostics.delivered_actuator_power_w,
178 central_transfer_power_w: diagnostics.central_transfer_power_w,
179 commanded_recharge_power_w: diagnostics.commanded_recharge_power_w,
180 parasitic_loss_w: diagnostics.parasitic_loss_w,
181 heat_generation_w: diagnostics.heat_generation_w,
182 heat_rejection_w: diagnostics.heat_rejection_w,
183 gain: diagnostics.gain,
184 damping: diagnostics.damping,
185 stiffness: diagnostics.stiffness,
186 mechanical_force_n: diagnostics.mechanical_force_n,
187 acceleration_mps2: diagnostics.acceleration_mps2,
188 delivered_ratio: diagnostics.delivered_ratio,
189 saturation_fraction: diagnostics.saturation_fraction,
190 ep_dot_j_per_s: diagnostics.ep_dot_j_per_s,
191 temperature_dot_k_per_s: diagnostics.temperature_dot_k_per_s,
192 mechanical_power_w: diagnostics.mechanical_power_w,
193 authority_utilization,
194 reduced_response_target_y_m,
195 reduced_response_target_rate_mps,
196 reduced_response_error_m,
197 reduced_response_error_rate_mps,
198 lyapunov_v,
199 lyapunov_dv_dt,
200 raw_next_ep_j: diagnostics.raw_next_ep_j,
201 raw_next_temperature_k: diagnostics.raw_next_temperature_k,
202 raw_next_y_m: diagnostics.raw_next_y_m,
203 raw_next_v_mps: diagnostics.raw_next_v_mps,
204 ep_clamped: diagnostics.ep_clamped,
205 temperature_clamped: diagnostics.temperature_clamped,
206 y_clamped: diagnostics.y_clamped,
207 ydot_clamped: diagnostics.v_clamped,
208 near_admissible_boundary: admissible.near_boundary,
209 outside_admissible_region: admissible.outside_region,
210 admissible_margin_fraction: admissible.margin_fraction,
211 });
212 for (index, limb_flow) in diagnostics.limb_flows.iter().enumerate() {
213 limb_buffers.push(LimbBufferRecord {
214 time_s: state.time_s,
215 limb: limb_flow.limb.clone(),
216 buffer_energy_j: state.local_buffers_j[index],
217 buffer_energy_mj: state.local_buffers_j[index] / 1.0e6,
218 requested_power_w: limb_flow.requested_power_w,
219 transfer_power_w: limb_flow.transfer_power_w,
220 delivered_power_w: limb_flow.delivered_power_w,
221 saturation: limb_flow.saturation,
222 });
223 }
224}
225
226fn update_constraint_events(
227 config: &SimulationConfig,
228 state: &SystemState,
229 diagnostics: &StepDiagnostics,
230 latch: &mut EventLatch,
231 events: &mut Vec<EventRecord>,
232) {
233 let admissible = admissible_status(
234 &config.model,
235 state.ep_j,
236 state.temperature_k,
237 state.y_m,
238 state.v_mps,
239 );
240 let low_energy = state.ep_j < config.model.low_energy_threshold_j;
241 let high_temperature = state.temperature_k >= config.model.thermal_limit_k;
242 let local_buffer_low = state
243 .local_buffers_j
244 .iter()
245 .any(|energy| *energy < config.model.local_buffer_low_threshold_j);
246 let saturated = diagnostics.saturation_fraction > 0.0;
247 let admissible_outside = admissible.outside_region;
248 let admissible_margin = admissible.near_boundary && !admissible.outside_region;
249
250 transition_event(
251 events,
252 state.time_s,
253 "energy_low",
254 "warning",
255 "Pulse-layer energy dropped below the configured threshold.",
256 state.ep_j,
257 config.model.low_energy_threshold_j,
258 latch.low_energy,
259 low_energy,
260 );
261 transition_event(
262 events,
263 state.time_s,
264 "thermal_high",
265 "failure",
266 "Aggregate thermal state exceeded the configured threshold.",
267 state.temperature_k,
268 config.model.thermal_limit_k,
269 latch.high_temperature,
270 high_temperature,
271 );
272 transition_event(
273 events,
274 state.time_s,
275 "local_buffer_low",
276 "warning",
277 "At least one limb-local energy buffer dropped below the configured threshold.",
278 state
279 .local_buffers_j
280 .iter()
281 .copied()
282 .fold(f64::INFINITY, f64::min),
283 config.model.local_buffer_low_threshold_j,
284 latch.local_buffer_low,
285 local_buffer_low,
286 );
287 transition_event(
288 events,
289 state.time_s,
290 "actuator_saturation",
291 "warning",
292 "Local transfer limits and buffer state reduced delivered actuator power below request.",
293 diagnostics.delivered_ratio,
294 1.0,
295 latch.saturated,
296 saturated,
297 );
298 transition_event(
299 events,
300 state.time_s,
301 "admissible_region_breach",
302 "failure",
303 "Reduced-order state exited the configured admissible operating region.",
304 admissible.margin_fraction,
305 0.0,
306 latch.admissible_outside,
307 admissible_outside,
308 );
309 transition_event(
310 events,
311 state.time_s,
312 "admissible_margin_warning",
313 "warning",
314 "Reduced-order state approached the admissible operating region boundary.",
315 admissible.margin_fraction,
316 0.05,
317 latch.admissible_margin,
318 admissible_margin,
319 );
320
321 latch.low_energy = low_energy;
322 latch.high_temperature = high_temperature;
323 latch.local_buffer_low = local_buffer_low;
324 latch.saturated = saturated;
325 latch.admissible_outside = admissible_outside;
326 latch.admissible_margin = admissible_margin;
327}
328
329fn transition_event(
330 events: &mut Vec<EventRecord>,
331 time_s: f64,
332 event_type: &str,
333 severity: &str,
334 message: &str,
335 value: f64,
336 threshold: f64,
337 previous: bool,
338 current: bool,
339) {
340 if !previous && current {
341 events.push(EventRecord {
342 time_s,
343 event_type: event_type.to_string(),
344 severity: severity.to_string(),
345 message: message.to_string(),
346 value,
347 threshold,
348 });
349 }
350 if previous && !current {
351 events.push(EventRecord {
352 time_s,
353 event_type: format!("{event_type}_recovered"),
354 severity: "info".to_string(),
355 message: format!("{message} Condition recovered."),
356 value,
357 threshold,
358 });
359 }
360}