# -*- coding: utf-8 -*- """Untitled162.ipynb Automatically generated by Colab. Original file is located at https://colab.research.google.com/drive/16l8lAe93PQPgmHhyZqpL8ga6aqNdsJfc """ #!/usr/bin/env python3 """Reproducible computations for the forest optimal-control article. This script reproduces the numerical material described in the manuscript "Optimal Control of a Forest Resource Model under Population Growth and Mining Activities". It generates: * Tables 4-6: baseline objective and terminal-state comparisons; * Table 3A: three-grid convergence check; * Tables 7-8A: horizon, initial-state, forest-capacity, and weight sensitivity; * Table 9: four-start direct single-shooting cross-check; * Figures 1-2 in PNG and PDF formats; * a machine-readable run manifest and optional verification report. The model variables are dimensionless indices on the illustrative scenario scale specified in the manuscript. The code does not use or imply empirical calibration to a geographical region. Examples -------- Run every experiment and verify the published values:: python Supplementary_Code_Forest_Optimal_Control_Revised.py \ --mode all --output-dir Supplementary_Outputs \ --seed 20260731 --verify Run only the inexpensive baseline calculations:: python Supplementary_Code_Forest_Optimal_Control_Revised.py \ --mode baseline --output-dir Supplementary_Outputs --verify Use ``--help`` for all command-line options. """ from __future__ import annotations import argparse import csv import hashlib import json import platform import sys from dataclasses import asdict, dataclass, replace from datetime import datetime, timezone from importlib import metadata from pathlib import Path from typing import Any, Callable, Iterable, Sequence import matplotlib matplotlib.use("Agg") import matplotlib.pyplot as plt import numpy as np from scipy.optimize import minimize SCRIPT_VERSION = "1.0.2" DEFAULT_SEED = 20260731 STATE_NAMES = ("F", "N", "P", "M") CONTROL_NAMES = ("u1", "u2") def trapezoid(values: np.ndarray, *, dx: float, axis: int = 0) -> np.ndarray: """Compatibility wrapper for NumPy versions before ``np.trapezoid``.""" method = getattr(np, "trapezoid", None) if method is not None: return method(values, dx=dx, axis=axis) return np.trapz(values, dx=dx, axis=axis) # pragma: no cover @dataclass(frozen=True) class ModelParameters: """Illustrative model coefficients and objective weights.""" L: float = 50.0 r: float = 1.0 r0: float = 2.0 gamma1: float = 0.01 gamma2: float = 0.01 s: float = 1.0 K: float = 50.0 s0: float = 2.0 pi: float = 0.01 lam: float = 1.0 theta0: float = 4.0 Lambda: float = 5.0 alpha: float = 0.3 mu0: float = 0.2 thetaE: float = 1.0 muR: float = 1.0 AF: float = 5.0 AP: float = 10.0 AM: float = 10.0 B1: float = 15.0 B2: float = 15.0 @dataclass class SweepSettings: """Forward-backward sweep settings.""" horizon: float = 40.0 dt: float = 0.01 initial_state: tuple[float, float, float, float] = (25.0, 25.0, 6.0, 25.0) initial_active_control: float = 0.5 relaxation: float = 0.5 tolerance: float = 1.0e-7 max_iterations: int = 3000 @dataclass class SweepResult: """State, control, and objective data returned by a sweep calculation.""" time: np.ndarray state: np.ndarray control: np.ndarray iterations: int converged: bool final_control_change: float objective: float settings: SweepSettings active_controls: tuple[bool, bool] def state_rhs( state: np.ndarray, control: np.ndarray, parameters: ModelParameters ) -> np.ndarray: """Right-hand side of equations (1)-(4).""" F, N, P, M = state u1, u2 = control p = parameters return np.array( [ F * ( p.r - p.r0 * F / p.L - p.gamma1 * M - p.gamma2 * N ), N * ( p.s - p.s0 * N / p.K + p.pi * p.gamma2 * F ), p.lam * N - (p.theta0 + p.thetaE * u1) * P, p.Lambda + p.alpha * P - (p.mu0 + p.muR * u2) * M, ], dtype=float, ) def adjoint_rhs( adjoint: np.ndarray, state: np.ndarray, control: np.ndarray, parameters: ModelParameters, ) -> np.ndarray: """Right-hand side of the adjoint equations (17)-(20).""" F, N, P, M = state eta1, eta2, eta3, eta4 = adjoint u1, u2 = control p = parameters return np.array( [ p.AF - eta1 * ( p.r - 2.0 * p.r0 * F / p.L - p.gamma1 * M - p.gamma2 * N ) - eta2 * p.pi * p.gamma2 * N, eta1 * p.gamma2 * F - eta2 * ( p.s - 2.0 * p.s0 * N / p.K + p.pi * p.gamma2 * F ) - eta3 * p.lam, eta3 * (p.theta0 + p.thetaE * u1) - eta4 * p.alpha - p.AP, eta1 * p.gamma1 * F + eta4 * (p.mu0 + p.muR * u2) - p.AM, ], dtype=float, ) def running_cost( state: np.ndarray, control: np.ndarray, parameters: ModelParameters ) -> np.ndarray: """Running integrand of objective functional (12).""" p = parameters F = state[..., 0] P = state[..., 2] M = state[..., 3] u1 = control[..., 0] u2 = control[..., 1] return ( p.AP * P + p.AM * M - p.AF * F + 0.5 * p.B1 * u1**2 + 0.5 * p.B2 * u2**2 ) def integrate_state( controls: np.ndarray, parameters: ModelParameters, settings: SweepSettings, ) -> np.ndarray: """Integrate the state system forward using fourth-order Runge-Kutta.""" dt = settings.dt intervals = controls.shape[0] - 1 states = np.empty((intervals + 1, 4), dtype=float) states[0] = np.asarray(settings.initial_state, dtype=float) for index in range(intervals): current = states[index] midpoint_control = 0.5 * (controls[index] + controls[index + 1]) k1 = state_rhs(current, controls[index], parameters) k2 = state_rhs(current + 0.5 * dt * k1, midpoint_control, parameters) k3 = state_rhs(current + 0.5 * dt * k2, midpoint_control, parameters) k4 = state_rhs(current + dt * k3, controls[index + 1], parameters) states[index + 1] = current + dt * (k1 + 2.0 * k2 + 2.0 * k3 + k4) / 6.0 return states def integrate_adjoint( states: np.ndarray, controls: np.ndarray, parameters: ModelParameters, settings: SweepSettings, ) -> np.ndarray: """Integrate the adjoint system backward from eta(T)=0 using RK4.""" dt = settings.dt intervals = controls.shape[0] - 1 adjoints = np.empty((intervals + 1, 4), dtype=float) adjoints[-1] = 0.0 backward_step = -dt for index in range(intervals, 0, -1): current = adjoints[index] midpoint_state = 0.5 * (states[index] + states[index - 1]) midpoint_control = 0.5 * (controls[index] + controls[index - 1]) k1 = adjoint_rhs(current, states[index], controls[index], parameters) k2 = adjoint_rhs( current + 0.5 * backward_step * k1, midpoint_state, midpoint_control, parameters, ) k3 = adjoint_rhs( current + 0.5 * backward_step * k2, midpoint_state, midpoint_control, parameters, ) k4 = adjoint_rhs( current + backward_step * k3, states[index - 1], controls[index - 1], parameters, ) adjoints[index - 1] = current + backward_step * ( k1 + 2.0 * k2 + 2.0 * k3 + k4 ) / 6.0 return adjoints def solve_forward_backward_sweep( parameters: ModelParameters, settings: SweepSettings, active_controls: tuple[bool, bool] = (True, True), ) -> SweepResult: """Compute a stationary candidate by the documented relaxed sweep.""" intervals_float = settings.horizon / settings.dt intervals = int(round(intervals_float)) if not np.isclose(intervals_float, intervals, rtol=0.0, atol=1.0e-12): raise ValueError("horizon/dt must be an integer") time = np.linspace(0.0, settings.horizon, intervals + 1) controls = np.full( (intervals + 1, 2), settings.initial_active_control, dtype=float ) for control_index, active in enumerate(active_controls): if not active: controls[:, control_index] = 0.0 if not any(active_controls): states = integrate_state(controls, parameters, settings) objective = float( trapezoid(running_cost(states, controls, parameters), dx=settings.dt) ) return SweepResult( time=time, state=states, control=controls, iterations=0, converged=True, final_control_change=0.0, objective=objective, settings=settings, active_controls=active_controls, ) converged = False final_change = float("inf") iteration = 0 states = np.empty((intervals + 1, 4), dtype=float) for iteration in range(1, settings.max_iterations + 1): states = integrate_state(controls, parameters, settings) adjoints = integrate_adjoint(states, controls, parameters, settings) candidate = np.column_stack( ( adjoints[:, 2] * parameters.thetaE * states[:, 2] / parameters.B1, adjoints[:, 3] * parameters.muR * states[:, 3] / parameters.B2, ) ) candidate = np.clip(candidate, 0.0, 1.0) for control_index, active in enumerate(active_controls): if not active: candidate[:, control_index] = 0.0 updated = ( settings.relaxation * candidate + (1.0 - settings.relaxation) * controls ) final_change = float(np.max(np.abs(updated - controls))) controls = updated if final_change < settings.tolerance: converged = True break states = integrate_state(controls, parameters, settings) objective = float( trapezoid(running_cost(states, controls, parameters), dx=settings.dt) ) return SweepResult( time=time, state=states, control=controls, iterations=iteration, converged=converged, final_control_change=final_change, objective=objective, settings=settings, active_controls=active_controls, ) class SweepCache: """Avoid repeated calculations shared by tables and figures.""" def __init__(self) -> None: self._results: dict[tuple[Any, ...], SweepResult] = {} @staticmethod def _key( parameters: ModelParameters, settings: SweepSettings, active_controls: tuple[bool, bool], ) -> tuple[Any, ...]: return ( tuple(asdict(parameters).values()), settings.horizon, settings.dt, tuple(settings.initial_state), settings.initial_active_control, settings.relaxation, settings.tolerance, settings.max_iterations, active_controls, ) def solve( self, parameters: ModelParameters, settings: SweepSettings, active_controls: tuple[bool, bool], ) -> SweepResult: key = self._key(parameters, settings, active_controls) if key not in self._results: self._results[key] = solve_forward_backward_sweep( parameters, settings, active_controls ) return self._results[key] def scalar(value: Any) -> Any: """Convert NumPy scalar objects to JSON/CSV-friendly Python values.""" if isinstance(value, np.generic): return value.item() return value def write_csv(path: Path, rows: Sequence[dict[str, Any]]) -> None: """Write a sequence of dictionaries with stable column order.""" if not rows: return path.parent.mkdir(parents=True, exist_ok=True) fields = list(rows[0].keys()) with path.open("w", newline="", encoding="utf-8") as handle: writer = csv.DictWriter(handle, fieldnames=fields) writer.writeheader() for row in rows: writer.writerow({key: scalar(value) for key, value in row.items()}) def write_json(path: Path, value: Any) -> None: """Write indented UTF-8 JSON.""" path.parent.mkdir(parents=True, exist_ok=True) with path.open("w", encoding="utf-8") as handle: json.dump(value, handle, indent=2, ensure_ascii=False, allow_nan=False) handle.write("\n") def trajectory_rows(result: SweepResult, scenario: str) -> list[dict[str, Any]]: """Return one machine-readable row per grid node.""" return [ { "scenario": scenario, "time": time, "F": state[0], "N": state[1], "P": state[2], "M": state[3], "u1": control[0], "u2": control[1], } for time, state, control in zip( result.time, result.state, result.control, strict=True ) ] def mean_controls(result: SweepResult) -> np.ndarray: """Composite-trapezoidal time means of the controls.""" return trapezoid( result.control, dx=result.settings.dt, axis=0, ) / result.settings.horizon def objective_independent_components(result: SweepResult) -> dict[str, float]: """Unweighted time averages used in Tables 7A and 8A.""" horizon = result.settings.horizon dt = result.settings.dt states = trapezoid(result.state, dx=dt, axis=0) / horizon controls_squared = trapezoid(result.control**2, dx=dt, axis=0) / horizon return { "I_F_over_T": float(states[0]), "I_P_over_T": float(states[2]), "I_M_over_T": float(states[3]), "I_u1_squared_over_T": float(controls_squared[0]), "I_u2_squared_over_T": float(controls_squared[1]), } def run_baseline( output_dir: Path, cache: SweepCache, parameters: ModelParameters, settings: SweepSettings, ) -> dict[str, SweepResult]: """Run four baseline scenarios and write Tables 4-6.""" definitions = ( ("Without control", (False, False)), ("u1", (True, False)), ("u2", (False, True)), ("u1 + u2", (True, True)), ) results: dict[str, SweepResult] = {} for label, active in definitions: print(f"[baseline] {label}", flush=True) result = cache.solve(parameters, settings, active) if not result.converged: raise RuntimeError(f"Sweep did not converge for {label}") results[label] = result no_control = results["Without control"] table4_rows = [] table5_rows = [] table6_rows = [] trajectory_data: list[dict[str, Any]] = [] for label, _active in definitions: result = results[label] reduction_j = ( 100.0 * (no_control.objective - result.objective) / no_control.objective if label != "Without control" else 0.0 ) table4_rows.append( { "scenario": label, "J": result.objective, "reduction_relative_to_no_control_percent": reduction_j, "iterations": result.iterations, "converged": result.converged, "final_control_change": result.final_control_change, } ) table5_rows.append( { "scenario": label, "F_terminal": result.state[-1, 0], "N_terminal": result.state[-1, 1], "P_terminal": result.state[-1, 2], "M_terminal": result.state[-1, 3], } ) trajectory_data.extend(trajectory_rows(result, label)) if label != "Without control": table6_rows.append( { "scenario": label, "increase_F_terminal_percent": 100.0 * (result.state[-1, 0] - no_control.state[-1, 0]) / no_control.state[-1, 0], "reduction_P_terminal_percent": 100.0 * (no_control.state[-1, 2] - result.state[-1, 2]) / no_control.state[-1, 2], "reduction_M_terminal_percent": 100.0 * (no_control.state[-1, 3] - result.state[-1, 3]) / no_control.state[-1, 3], "reduction_J_percent": reduction_j, } ) write_csv(output_dir / "table4_objective_values.csv", table4_rows) write_csv(output_dir / "table5_terminal_state_values.csv", table5_rows) write_csv(output_dir / "table6_effectiveness.csv", table6_rows) write_csv(output_dir / "baseline_trajectories.csv", trajectory_data) return results def run_grid_convergence( output_dir: Path, cache: SweepCache, parameters: ModelParameters, settings: SweepSettings, ) -> list[dict[str, Any]]: """Run the Table 3A three-grid convergence calculation.""" rows: list[dict[str, Any]] = [] for dt in (0.02, 0.01, 0.005): print(f"[grid] dt={dt:.3f}", flush=True) grid_settings = replace(settings, dt=dt) result = cache.solve(parameters, grid_settings, (True, True)) if not result.converged: raise RuntimeError(f"Sweep did not converge for dt={dt}") rows.append( { "dt": dt, "J": result.objective, "F_terminal": result.state[-1, 0], "P_terminal": result.state[-1, 2], "M_terminal": result.state[-1, 3], "iterations": result.iterations, } ) write_csv(output_dir / "table3A_grid_convergence.csv", rows) return rows def run_sensitivity( output_dir: Path, cache: SweepCache, baseline_parameters: ModelParameters, baseline_settings: SweepSettings, ) -> dict[str, list[dict[str, Any]]]: """Run Tables 7-8A exactly as specified in the manuscript.""" table7_definitions: list[ tuple[str, dict[str, float], float, tuple[float, float, float, float]] ] = [ ("B1 = B2 = 7.5", {"B1": 7.5, "B2": 7.5}, 40.0, (25, 25, 6, 25)), ("B1 = B2 = 15 (baseline)", {}, 40.0, (25, 25, 6, 25)), ("B1 = B2 = 30", {"B1": 30.0, "B2": 30.0}, 40.0, (25, 25, 6, 25)), ("T = 20", {}, 20.0, (25, 25, 6, 25)), ("T = 60", {}, 60.0, (25, 25, 6, 25)), ("X0 = (20,20,4,20)", {}, 40.0, (20, 20, 4, 20)), ("X0 = (30,30,8,30)", {}, 40.0, (30, 30, 8, 30)), ("L = 20", {"L": 20.0}, 40.0, (25, 25, 6, 25)), ("L = 40", {"L": 40.0}, 40.0, (25, 25, 6, 25)), ("L = 60", {"L": 60.0}, 40.0, (25, 25, 6, 25)), ] table7_rows: list[dict[str, Any]] = [] table7a_rows: list[dict[str, Any]] = [] for label, changes, horizon, initial_state in table7_definitions: print(f"[sensitivity-7] {label}", flush=True) parameters = replace(baseline_parameters, **changes) settings = replace( baseline_settings, horizon=horizon, initial_state=tuple(float(value) for value in initial_state), ) result = cache.solve(parameters, settings, (True, True)) if not result.converged: raise RuntimeError(f"Sweep did not converge for {label}") controls = mean_controls(result) table7_rows.append( { "test": label, "horizon": horizon, "F_terminal": result.state[-1, 0], "P_terminal": result.state[-1, 2], "M_terminal": result.state[-1, 3], "mean_u1": controls[0], "mean_u2": controls[1], "iterations": result.iterations, } ) table7a_rows.append( { "test": label, **objective_independent_components(result), } ) table8_definitions: list[tuple[str, dict[str, float]]] = [ ("AF = 2.5", {"AF": 2.5}), ("AF = 10", {"AF": 10.0}), ("AP = 5", {"AP": 5.0}), ("AP = 20", {"AP": 20.0}), ("AM = 5", {"AM": 5.0}), ("AM = 20", {"AM": 20.0}), ("B1 = 7.5", {"B1": 7.5}), ("B1 = 30", {"B1": 30.0}), ("B2 = 7.5", {"B2": 7.5}), ("B2 = 30", {"B2": 30.0}), ("B2 = 60", {"B2": 60.0}), ("B2 = 150", {"B2": 150.0}), ] table8_rows: list[dict[str, Any]] = [] table8a_rows: list[dict[str, Any]] = [] for label, changes in table8_definitions: print(f"[sensitivity-8] {label}", flush=True) parameters = replace(baseline_parameters, **changes) result = cache.solve(parameters, baseline_settings, (True, True)) if not result.converged: raise RuntimeError(f"Sweep did not converge for {label}") controls = mean_controls(result) saturation_fraction = float( trapezoid( (result.control[:, 1] >= 0.99).astype(float), dx=result.settings.dt, ) / result.settings.horizon ) table8_rows.append( { "test": label, "F_terminal": result.state[-1, 0], "P_terminal": result.state[-1, 2], "M_terminal": result.state[-1, 3], "mean_u1": controls[0], "mean_u2": controls[1], "S2": saturation_fraction, "iterations": result.iterations, } ) table8a_rows.append( { "test": label, **objective_independent_components(result), } ) write_csv(output_dir / "table7_sensitivity.csv", table7_rows) write_csv(output_dir / "table7A_objective_independent_components.csv", table7a_rows) write_csv(output_dir / "table8_weight_sensitivity.csv", table8_rows) write_csv(output_dir / "table8A_objective_independent_components.csv", table8a_rows) return { "table7": table7_rows, "table7A": table7a_rows, "table8": table8_rows, "table8A": table8a_rows, } def make_figures(output_dir: Path, results: dict[str, SweepResult]) -> None: """Generate manuscript-style state and control profile figures.""" labels = { "Without control": "Without control", "u1": "u1", "u2": "u2", "u1 + u2": "u1 + u2", } colors = { "Without control": "#20B13D", "u1": "#202020", "u2": "#D62F89", "u1 + u2": "#1464F4", } linestyles = { "Without control": "--", "u1": "--", "u2": "--", "u1 + u2": "--", } markers = { "Without control": "h", "u1": "s", "u2": "4", "u1 + u2": "p", } horizon = next(iter(results.values())).settings.horizon # Figure 1 – State profiles figure, axes = plt.subplots(2, 2, figsize=(10, 7.5), constrained_layout=True) y_labels = [ "Forest-resource condition index", "Human-population index", "Population-pressure index", "Mining-activity index", ] for state_index, axis in enumerate(axes.flat): for scenario, result in results.items(): axis.plot( result.time, result.state[:, state_index], label=labels[scenario], color=colors[scenario], linestyle=linestyles[scenario], linewidth=1.7, marker=markers[scenario], markevery=max(1, len(result.time) // 16), markersize=4.5, ) axis.set_xlabel(f"Years\n({chr(97 + state_index)})") axis.set_ylabel(y_labels[state_index]) axis.set_xlim(0.0, horizon) axis.legend(fontsize=8) figure.suptitle("State-index profiles under the four baseline scenarios") figure.savefig(output_dir / "figure1_state_profiles.png", dpi=300) figure.savefig(output_dir / "figure1_state_profiles.pdf") plt.close(figure) # Figure 2 – Control profiles figure = plt.figure(figsize=(10, 7.5), constrained_layout=True) grid = figure.add_gridspec(2, 4) axis_u1 = figure.add_subplot(grid[0, 0:2]) axis_u2 = figure.add_subplot(grid[0, 2:4]) # Panel ketiga ditempatkan di tengah di bawah dua panel pertama. axis_both = figure.add_subplot(grid[1, 1:3]) u1_result = results["u1"] axis_u1.plot(u1_result.time, u1_result.control[:, 0], color="#d62728") axis_u1.set_ylabel("Control profile") u2_result = results["u2"] axis_u2.plot(u2_result.time, u2_result.control[:, 1], color="#1f77b4") axis_u2.set_ylabel("Control profile") both_result = results["u1 + u2"] axis_both.plot( both_result.time, both_result.control[:, 0], label="u1", color="#d62728" ) axis_both.plot( both_result.time, both_result.control[:, 1], label="u2", color="#1f77b4" ) axis_both.set_ylabel("Control profile") axis_both.legend() for panel_index, axis in enumerate((axis_u1, axis_u2, axis_both)): axis.set_xlabel(f"Years\n({chr(97 + panel_index)})") axis.set_xlim(0.0, horizon) axis.set_ylim(-0.02, 1.05) figure.savefig(output_dir / "figure2_control_profiles.png", dpi=300) figure.savefig(output_dir / "figure2_control_profiles.pdf") plt.close(figure) def state_jacobian( state: np.ndarray, control: np.ndarray, parameters: ModelParameters ) -> np.ndarray: """Jacobian of the state dynamics with respect to the state.""" F, N, P, M = state u1, u2 = control p = parameters return np.array( [ [ p.r - 2.0 * p.r0 * F / p.L - p.gamma1 * M - p.gamma2 * N, -p.gamma2 * F, 0.0, -p.gamma1 * F, ], [p.pi * p.gamma2 * N, p.s - 2.0 * p.s0 * N / p.K + p.pi * p.gamma2 * F, 0.0, 0.0], [0.0, p.lam, -(p.theta0 + p.thetaE * u1), 0.0], [0.0, 0.0, p.alpha, -(p.mu0 + p.muR * u2)], ], dtype=float, ) def control_jacobian( state: np.ndarray, parameters: ModelParameters ) -> np.ndarray: """Jacobian of the state dynamics with respect to (u1,u2).""" _F, _N, P, M = state p = parameters return np.array( [ [0.0, 0.0], [0.0, 0.0], [-p.thetaE * P, 0.0], [0.0, -p.muR * M], ], dtype=float, ) def direct_augmented_rhs( augmented: np.ndarray, control: np.ndarray, segment_index: int, segments: int, parameters: ModelParameters, ) -> np.ndarray: """State, objective, and variational equations for direct shooting.""" variables = 2 * segments state = augmented[:4] state_sensitivity_start = 5 state_sensitivity_end = state_sensitivity_start + 4 * variables state_sensitivity = augmented[ state_sensitivity_start:state_sensitivity_end ].reshape(4, variables) objective_sensitivity = augmented[state_sensitivity_end:] selector = np.zeros((2, variables), dtype=float) selector[0, 2 * segment_index] = 1.0 selector[1, 2 * segment_index + 1] = 1.0 dynamics = state_rhs(state, control, parameters) cost = float(running_cost(state, control, parameters)) derivative_state_sensitivity = ( state_jacobian(state, control, parameters) @ state_sensitivity + control_jacobian(state, parameters) @ selector ) cost_state_gradient = np.array( [-parameters.AF, 0.0, parameters.AP, parameters.AM], dtype=float ) cost_control_gradient = np.array( [parameters.B1 * control[0], parameters.B2 * control[1]], dtype=float ) derivative_objective_sensitivity = ( cost_state_gradient @ state_sensitivity + cost_control_gradient @ selector ) return np.concatenate( ( dynamics, np.array([cost]), derivative_state_sensitivity.ravel(), derivative_objective_sensitivity, ) ) def evaluate_direct_shooting( decision: np.ndarray, parameters: ModelParameters, initial_state: Sequence[float], horizon: float = 40.0, segments: int = 40, substep: float = 0.1, ) -> tuple[float, np.ndarray, np.ndarray]: """Evaluate direct shooting objective and variational-equation gradient.""" if decision.size != 2 * segments: raise ValueError("decision vector has the wrong length") segment_length = horizon / segments substeps_float = segment_length / substep substeps = int(round(substeps_float)) if not np.isclose(substeps_float, substeps, rtol=0.0, atol=1.0e-12): raise ValueError("segment length/substep must be an integer") variables = 2 * segments augmented = np.zeros(5 + 4 * variables + variables, dtype=float) augmented[:4] = np.asarray(initial_state, dtype=float) controls = decision.reshape(segments, 2) for segment_index, control in enumerate(controls): for _ in range(substeps): k1 = direct_augmented_rhs( augmented, control, segment_index, segments, parameters ) k2 = direct_augmented_rhs( augmented + 0.5 * substep * k1, control, segment_index, segments, parameters, ) k3 = direct_augmented_rhs( augmented + 0.5 * substep * k2, control, segment_index, segments, parameters, ) k4 = direct_augmented_rhs( augmented + substep * k3, control, segment_index, segments, parameters, ) augmented += substep * (k1 + 2.0 * k2 + 2.0 * k3 + k4) / 6.0 objective = float(augmented[4]) gradient = augmented[5 + 4 * variables :].copy() terminal_state = augmented[:4].copy() return objective, gradient, terminal_state class DirectObjective: """Cache direct objective and gradient calls made by SLSQP.""" def __init__( self, parameters: ModelParameters, initial_state: Sequence[float], horizon: float, segments: int, substep: float, ) -> None: self.parameters = parameters self.initial_state = tuple(initial_state) self.horizon = horizon self.segments = segments self.substep = substep self._decision: np.ndarray | None = None self._result: tuple[float, np.ndarray, np.ndarray] | None = None def evaluate(self, decision: np.ndarray) -> tuple[float, np.ndarray, np.ndarray]: if self._decision is None or not np.array_equal(decision, self._decision): self._decision = np.asarray(decision, dtype=float).copy() self._result = evaluate_direct_shooting( self._decision, self.parameters, self.initial_state, self.horizon, self.segments, self.substep, ) assert self._result is not None return self._result def fun(self, decision: np.ndarray) -> float: return self.evaluate(decision)[0] def jac(self, decision: np.ndarray) -> np.ndarray: return self.evaluate(decision)[1] def interval_average_controls(result: SweepResult, segments: int = 40) -> np.ndarray: """Average the fine-grid sweep controls over equal direct-shooting segments.""" segment_length = result.settings.horizon / segments averaged = np.empty((segments, 2), dtype=float) for segment in range(segments): left = segment * segment_length right = (segment + 1) * segment_length mask = (result.time >= left - 1.0e-12) & (result.time <= right + 1.0e-12) local_time = result.time[mask] local_control = result.control[mask] if local_time.size < 2: midpoint = left + 0.5 * segment_length averaged[segment] = np.column_stack( [ np.interp(midpoint, result.time, result.control[:, index]) for index in range(2) ] )[0] else: averaged[segment] = trapezoid( local_control, dx=result.settings.dt, axis=0, ) / segment_length return np.clip(averaged, 0.0, 1.0) def run_direct_crosscheck( output_dir: Path, cache: SweepCache, parameters: ModelParameters, settings: SweepSettings, seed: int, ) -> list[dict[str, Any]]: """Run the four-start direct single-shooting check reported in Table 9.""" sweep = cache.solve(parameters, settings, (True, True)) segments = 40 substep = 0.1 rng = np.random.default_rng(seed) starts = [ ("sweep-downsampled start", interval_average_controls(sweep, segments).ravel()), ("zero start", np.zeros(2 * segments, dtype=float)), ("one start", np.ones(2 * segments, dtype=float)), (f"pseudo-random start (seed {seed})", rng.random(2 * segments)), ] rows: list[dict[str, Any]] = [ { "method": "Forward-backward sweep (RK4)", "objective_J": sweep.objective, "iterations": sweep.iterations, "difference_from_sweep_percent": 0.0, "success": True, "status": 0, "function_evaluations": "", "gradient_evaluations": "", } ] bounds = [(0.0, 1.0)] * (2 * segments) for start_label, initial_decision in starts: print(f"[direct] {start_label}", flush=True) objective = DirectObjective( parameters, settings.initial_state, settings.horizon, segments, substep, ) result = minimize( objective.fun, initial_decision, method="SLSQP", jac=objective.jac, bounds=bounds, options={"ftol": 1.0e-8, "maxiter": 100, "disp": False}, ) terminal_state = objective.evaluate(result.x)[2] difference = 100.0 * (float(result.fun) - sweep.objective) / sweep.objective rows.append( { "method": f"Single shooting: {start_label}", "objective_J": float(result.fun), "iterations": int(result.nit), "difference_from_sweep_percent": difference, "success": bool(result.success), "status": int(result.status), "function_evaluations": int(result.nfev), "gradient_evaluations": int(result.njev), } ) controls = result.x.reshape(segments, 2) control_rows = [ { "start": start_label, "segment": index + 1, "time_start": index * settings.horizon / segments, "time_end": (index + 1) * settings.horizon / segments, "u1": controls[index, 0], "u2": controls[index, 1], } for index in range(segments) ] safe_label = ( start_label.replace(" ", "_") .replace("(", "") .replace(")", "") .replace("-", "_") ) write_csv(output_dir / f"direct_controls_{safe_label}.csv", control_rows) write_json( output_dir / f"direct_result_{safe_label}.json", { "start": start_label, "seed": seed, "success": bool(result.success), "status": int(result.status), "message": str(result.message), "iterations": int(result.nit), "objective": float(result.fun), "terminal_state": { name: float(value) for name, value in zip(STATE_NAMES, terminal_state, strict=True) }, }, ) write_csv(output_dir / "table9_direct_single_shooting.csv", rows) return rows def add_check( checks: list[dict[str, Any]], name: str, observed: float, expected: float, tolerance: float, ) -> None: """Append a scalar numerical verification check.""" absolute_error = abs(observed - expected) checks.append( { "check": name, "observed": observed, "expected": expected, "absolute_error": absolute_error, "tolerance": tolerance, "passed": absolute_error <= tolerance, } ) def verify_results( output_dir: Path, baseline: dict[str, SweepResult] | None, grid_rows: list[dict[str, Any]] | None, sensitivity: dict[str, list[dict[str, Any]]] | None, direct_rows: list[dict[str, Any]] | None, ) -> bool: """Compare computed values with the numerical values printed in the article.""" checks: list[dict[str, Any]] = [] if baseline is not None: expected_objectives = { "Without control": 13585.205, "u1": 12657.053, "u2": 1768.848, "u1 + u2": 1427.949, } expected_terminal = { "Without control": (10.1491, 25.0254, 6.2563, 34.3814), "u1": (10.6136, 25.0265, 5.6692, 32.5897), "u2": (17.2946, 25.0433, 6.2608, 6.4763), "u1 + u2": (17.3628, 25.0434, 5.6751, 6.2302), } for scenario, expected in expected_objectives.items(): add_check( checks, f"Table 4 J: {scenario}", baseline[scenario].objective, expected, 0.0015, ) for index, state_name in enumerate(STATE_NAMES): add_check( checks, f"Table 5 {state_name}(T): {scenario}", float(baseline[scenario].state[-1, index]), expected_terminal[scenario][index], 0.00015, ) if grid_rows is not None: expected_grid = { 0.02: (1427.958084, 17.362732, 5.675394, 6.230808, 25), 0.01: (1427.948807, 17.362759, 5.675138, 6.230188, 25), 0.005: (1427.946486, 17.362762, 5.675074, 6.230104, 25), } for row in grid_rows: expected = expected_grid[float(row["dt"])] for field, expected_value in zip( ("J", "F_terminal", "P_terminal", "M_terminal"), expected[:4], strict=True, ): add_check( checks, f"Table 3A {field}: dt={row['dt']}", float(row[field]), float(expected_value), 0.000015, ) add_check( checks, f"Table 3A iterations: dt={row['dt']}", float(row["iterations"]), float(expected[4]), 0.0, ) if sensitivity is not None: key_expected = { "L = 20": (6.948, 5.670, 6.233, 0.876, 0.996), "B1 = B2 = 15 (baseline)": (17.363, 5.675, 6.230, 0.886, 0.996), } rows7 = {str(row["test"]): row for row in sensitivity["table7"]} for label, expected in key_expected.items(): row = rows7[label] for field, expected_value in zip( ("F_terminal", "P_terminal", "M_terminal", "mean_u1", "mean_u2"), expected, strict=True, ): add_check( checks, f"Table 7 {field}: {label}", float(row[field]), expected_value, 0.0015, ) rows8 = {str(row["test"]): row for row in sensitivity["table8"]} key_expected8 = { "B2 = 60": (17.124, 5.674, 8.335, 0.894, 0.930, 0.070), "B2 = 150": (16.577, 5.673, 10.775, 0.939, 0.658, 0.020), } for label, expected in key_expected8.items(): row = rows8[label] for field, expected_value in zip( ("F_terminal", "P_terminal", "M_terminal", "mean_u1", "mean_u2", "S2"), expected, strict=True, ): add_check( checks, f"Table 8 {field}: {label}", float(row[field]), expected_value, 0.0015, ) if direct_rows is not None: for row in direct_rows[1:]: add_check( checks, f"Table 9 J: {row['method']}", float(row["objective_J"]), 1429.007, 0.002, ) passed = all(bool(check["passed"]) for check in checks) report = { "passed": passed, "checks_passed": sum(bool(check["passed"]) for check in checks), "checks_total": len(checks), "checks": checks, "note": ( "SLSQP iteration counts may differ slightly across SciPy versions; " "verification uses the converged objective rather than iteration count." ), } write_json(output_dir / "verification_report.json", report) write_csv(output_dir / "verification_checks.csv", checks) return passed def dependency_version(package: str) -> str: """Return installed dependency version without failing manifest generation.""" try: return metadata.version(package) except metadata.PackageNotFoundError: return "not installed" def source_identity() -> tuple[str, str | None]: """Return the source name and hash, including pasted-notebook execution. Python defines ``__file__`` when this program is run from a saved file, but not when its contents are pasted directly into a Jupyter or Colab cell. In that case there is no source file to hash, so the manifest records an explicit notebook-cell identity and a null hash. """ source = globals().get("__file__") if source is not None: script_path = Path(str(source)).resolve() if script_path.is_file(): return ( script_path.name, hashlib.sha256(script_path.read_bytes()).hexdigest(), ) return "notebook-cell", None def create_manifest( output_dir: Path, arguments: argparse.Namespace, parameters: ModelParameters, settings: SweepSettings, generated_files: Iterable[Path], ) -> None: """Record parameters, environment, seed, and hashes of generated outputs.""" script_name, script_sha256 = source_identity() files = [] for path in sorted(generated_files): if path.is_file() and path.name != "run_manifest.json": files.append( { "name": path.name, "sha256": hashlib.sha256(path.read_bytes()).hexdigest(), "size_bytes": path.stat().st_size, } ) manifest = { "script": script_name, "script_version": SCRIPT_VERSION, "script_sha256": script_sha256, "run_timestamp_utc": datetime.now(timezone.utc).isoformat(), "command_line": sys.argv, "mode": arguments.mode, "seed": arguments.seed, "verify_requested": bool(arguments.verify), "skip_direct": bool(arguments.skip_direct), "skip_figures": bool(arguments.skip_figures), "parameters": asdict(parameters), "sweep_settings": asdict(settings), "direct_settings": { "segments": 40, "piecewise_constant_control_pairs": 40, "rk4_substep_Years": 0.1, "method": "SLSQP", "gradient": "variational equations", "ftol": 1.0e-8, "max_iterations": 100, "starts": ["sweep-downsampled", "zeros", "ones", "pseudo-random"], }, "environment": { "python": platform.python_version(), "platform": platform.platform(), "numpy": dependency_version("numpy"), "scipy": dependency_version("scipy"), "matplotlib": dependency_version("matplotlib"), }, "generated_files": files, } write_json(output_dir / "run_manifest.json", manifest) def build_parser() -> argparse.ArgumentParser: """Build the documented command-line interface.""" parser = argparse.ArgumentParser( description=( "Reproduce the baseline, sensitivity, grid-convergence, direct " "single-shooting, and figure calculations reported in the article." ) ) parser.add_argument( "--mode", choices=("all", "baseline", "grid", "sensitivity", "direct", "figures"), default="all", help="experiment group to run (default: all)", ) parser.add_argument( "--output-dir", type=Path, default=Path("Supplementary_Outputs"), help="directory for CSV, JSON, PNG, and PDF outputs", ) parser.add_argument( "--seed", type=int, default=DEFAULT_SEED, help=f"pseudo-random direct-shooting seed (default: {DEFAULT_SEED})", ) parser.add_argument( "--verify", action="store_true", help="compare computed values with the values printed in the article", ) parser.add_argument( "--skip-direct", action="store_true", help="with --mode all, skip the more expensive four-start direct check", ) parser.add_argument( "--skip-figures", action="store_true", help="with --mode all, skip PNG/PDF figure generation", ) parser.add_argument( "--version", action="version", version=f"%(prog)s {SCRIPT_VERSION}", ) return parser def running_in_notebook() -> bool: """Return True when execution is hosted by Jupyter or Google Colab.""" return "ipykernel" in sys.modules or "google.colab" in sys.modules def parse_runtime_arguments( parser: argparse.ArgumentParser, argv: Sequence[str] | None, ) -> argparse.Namespace: """Parse CLI arguments while tolerating Jupyter's injected kernel flags. Jupyter and Google Colab launch notebooks with arguments such as ``-f /path/to/kernel.json``. Those arguments belong to the notebook kernel, not to this program. They are ignored only in a detected notebook session; ordinary command-line execution remains strict so that misspelled options still produce an error. """ if argv is not None: return parser.parse_args(argv) if running_in_notebook(): arguments, notebook_arguments = parser.parse_known_args() if notebook_arguments: print( "[notebook] Ignoring kernel arguments: " + " ".join(notebook_arguments), flush=True, ) return arguments return parser.parse_args() def main(argv: Sequence[str] | None = None) -> int: """Command-line entry point.""" parser = build_parser() arguments = parse_runtime_arguments(parser, argv) output_dir = arguments.output_dir.resolve() output_dir.mkdir(parents=True, exist_ok=True) parameters = ModelParameters() settings = SweepSettings() cache = SweepCache() baseline: dict[str, SweepResult] | None = None grid_rows: list[dict[str, Any]] | None = None sensitivity: dict[str, list[dict[str, Any]]] | None = None direct_rows: list[dict[str, Any]] | None = None modes_requiring_baseline = {"all", "baseline", "direct", "figures"} if arguments.mode in modes_requiring_baseline: baseline = run_baseline(output_dir, cache, parameters, settings) if arguments.mode in {"all", "grid"}: grid_rows = run_grid_convergence(output_dir, cache, parameters, settings) if arguments.mode in {"all", "sensitivity"}: sensitivity = run_sensitivity(output_dir, cache, parameters, settings) if arguments.mode in {"all", "direct"} and not arguments.skip_direct: direct_rows = run_direct_crosscheck( output_dir, cache, parameters, settings, arguments.seed ) if arguments.mode in {"all", "figures"} and not arguments.skip_figures: assert baseline is not None print("[figures] state and control profiles", flush=True) make_figures(output_dir, baseline) verification_passed = True if arguments.verify: print("[verify] comparing generated values with the article", flush=True) verification_passed = verify_results( output_dir, baseline, grid_rows, sensitivity, direct_rows ) generated_files = list(output_dir.iterdir()) create_manifest(output_dir, arguments, parameters, settings, generated_files) print(f"Outputs written to: {output_dir}", flush=True) if arguments.verify: status = "PASSED" if verification_passed else "FAILED" print(f"Verification: {status}", flush=True) return 0 if verification_passed else 1 if __name__ == "__main__": exit_code = main() if not running_in_notebook(): raise SystemExit(exit_code)