3.1.1. Learnable fuzzy gain-scheduled PID (LF-PID) formulation

Consider a robotic system with $n$ controllable joints. For each joint $i$ , we use a PID controller whose gains are scheduled online by a learnable fuzzy system:

(1) \begin{equation} u_i(t) = K_{p,i}(t)\, e_i(t) + K_{i,i}(t) \int _0^t e_i(\tau ) d\tau + K_{d,i}(t)\, \dot {e}_i(t) \end{equation}

where $e_i(t) = q_{ref,i}(t) - q_i(t)$ is the tracking error.

We parameterize gain scheduling as a bounded correction around a meta-initialized base gain:

(2) \begin{equation} K_{\, ,i}(t) = \bar {K}_{\, ,i} + \Delta K_{\, ,i}(\mathbf{x}_i(t)), \quad \, \in \{p,i,d\}. \end{equation}

The fuzzy system input is a per-joint, scale-normalized error state

(3) \begin{equation} \mathbf{x}_i(t) = \big [s_{e,i} e_i(t),\; s_{\dot {e},i} \dot {e}_i(t),\; s_{I,i} \!\int _0^t e_i(\tau )d\tau \big ]^\top ,\quad \mathbf{x}_i \in \mathbb{R}^3, \end{equation}

where $(s_{e,i}, s_{\dot {e},i}, s_{I,i})$ are learned per-joint scaling factors.

Shared fuzzy membership partitions: To preserve common error semantics across platforms and joints, we use shared membership functions for each input dimension (e.g., Gaussian partitions with shared centers $c_{m,d}$ and widths $\sigma _{m,d}$ for $d \in \{e,\dot {e},I\}$ ).

Per-joint Takagi–Sugeno consequents: We employ a differentiable Takagi–Sugeno (TS) fuzzy model. For each rule $r$ , the activation weight is computed as a product of membership values, $w_r(\mathbf{x}_i)$ , and each gain correction is an activation-weighted affine function:

(4) \begin{equation} \Delta K_{\, ,i}(\mathbf{x}_i) = \frac {\sum _{r} w_r(\mathbf{x}_i)\,\big (\mathbf{a}_{\, ,i,r}^\top \mathbf{x}_i + b_{\, ,i,r}\big )}{\sum _{r} w_r(\mathbf{x}_i) + \epsilon }, \end{equation}

where $(\mathbf{a}_{\, ,i,r}, b_{\, ,i,r})$ are learned per-joint consequent parameters and $\epsilon$ ensures numerical stability. In practice, we enforce bounded and nonnegative gains via standard parameterizations (e.g., softplus for base gains and tanh-bounded corrections), which improves training stability during online adaptation.

Complete parameter set definition: Given a robot with $n$ joints and a fuzzy system with $R$ rules (derived from $M_e \times M_{\dot {e}} \times M_I$ membership functions per input dimension), the complete LF-PID parameter vector $\theta$ is partitioned into three subsets:

1. 1. Base gains $\bar {\boldsymbol{{K}}} = [\bar {K}_{p,1},\ldots ,\bar {K}_{p,n}, \bar {K}_{i,1},\ldots ,\bar {K}_{i,n}, \bar {K}_{d,1},\ldots ,\bar {K}_{d,n}]^\top \in \mathbb{R}^{3n}$ – meta-initialized nominal gains per joint

2. 2. Input scaling factors $\boldsymbol{{s}} = [s_{e,1}, s_{\dot {e},1}, s_{I,1}, \ldots , s_{e,n}, s_{\dot {e},n}, s_{I,n}]^\top \in \mathbb{R}^{3n}$ – per-joint normalization of fuzzy inputs

3. 3. TS consequent parameters $\boldsymbol{{c}} = \{(\mathbf{a}_{\, ,i,r}, b_{\, ,i,r})\}$ for $\, \in \{p,i,d\}$ , $i=1,\ldots ,n$ , $r=1,\ldots ,R$ – total dimension $d_c = 3nR(3+1) = 12nR$ (3 gains $\times$ n joints $\times$ R rules $\times$ (3 affine coefficients + 1 bias))

The complete parameter vector dimension is

(5) \begin{equation} D_{\text{total}} = 3n + 3n + 12nR = 6n + 12nR = 6n(1 + 2R) \end{equation}

For our experiments with $M_e=M_{\dot {e}}=M_I=3$ membership functions per input, we have $R = 3^3 = 27$ rules, yielding $D_{\text{total}} = 6n(1+54) = 330n$ parameters. For Franka Panda (n = 9), this gives 2,970 parameters; for Laikago (n = 12), 3,960 parameters.

Stage-specific parameter subsets:

• • Meta-learning stage: Predicts initialization for all three subsets $\theta _{\text{meta}} = [\bar {\boldsymbol{{K}}}, \boldsymbol{{s}}, \boldsymbol{{c}}]$ , dimension $D_{\text{total}} = 330n$

• • Zero-shot deployment (Meta-LF-PID): Uses $\theta _{\text{meta}}$ directly without adaptation

• • RL adaptation stage: Refines only input scaling factors $\boldsymbol{{s}}$ (3n parameters) while keeping base gains $\bar {\boldsymbol{{K}}}$ and consequents $\boldsymbol{{c}}$ fixed at meta-initialized values. This compact adaptation space enables fast convergence (1M steps, 10 minutes) and prevents catastrophic forgetting of cross-platform knowledge encoded in consequents.

3.1.2. Optimization objective

Our goal is to find LF-PID parameters $\theta$ (base gains, per-joint scaling factors, and TS consequents) that minimize the tracking error across different robot platforms and operating conditions:

(6) \begin{equation} \theta ^* = \arg \min _{\theta } \mathbb{E}_{r \sim \mathcal{R}} \left [ \mathcal{L}_r(\theta ) \right ] \end{equation}

where $\mathcal{R}$ is a distribution over robot platforms and $\mathcal{L}_r(\theta )$ is the tracking error for robot $r$ with parameters  $\theta$ .

3.2.1. Stage 1: Meta-learning for LF-PID initialization

The meta-learning stage learns a mapping $f_{\phi }\;:\; \mathcal{F} \rightarrow \Theta$ from robot feature space $\mathcal{F}$ to LF-PID parameter space $\Theta$ (base gains, per-joint scaling factors, and TS consequent parameters).

Robot feature extraction: For each robot platform, we extract 10-dimensional physical features:

(7) \begin{equation} \mathbf{f} = [n_{dof}, m_{total}, \mathbf{I}, \mathbf{L}, \mathbf{r}_{com}, \boldsymbol{\mu }] \in \mathbb{R}^{10} \end{equation}

where $n_{dof}$ is degrees of freedom, $m_{total}$ is total mass, $\mathbf{I}$ are inertia tensor components, $\mathbf{L}$ are link lengths, $\mathbf{r}_{com}$ are center-of-mass positions, and $\boldsymbol{\mu }$ are friction coefficients. Features are normalized to facilitate cross-platform learning.

Figure 1.

Meta-learning network architecture for LF-PID initialization. The network maps 10D robot features $\mathbf{f}$ to the complete LF-PID parameter vector $\hat {\boldsymbol{\theta }} = [\hat {\bar {\boldsymbol{{K}}}}, \hat {\boldsymbol{{s}}}, \hat {\boldsymbol{{c}}}] \in [0,1]^D$ ( $D=330n$ ) via three structured output heads: base gains $\bar {\boldsymbol{{K}}} \in \mathbb{R}^{3n}$ , input scaling factors $\boldsymbol{{s}} \in \mathbb{R}^{3n}$ , and TS consequent parameters $\boldsymbol{{c}} \in \mathbb{R}^{d_c}$ . Sigmoid activation ensures bounded outputs; denormalization maps to physical ranges.

Network architecture: We employ a hierarchical feedforward neural network with two encoder layers (256D each) and one hidden layer (128D), followed by a parameter head that outputs the LF-PID initialization vector, as illustrated in Figure 1. The forward propagation is defined as:

(8) \begin{align} \mathbf{h}_1 &= \text{ReLU}(\text{LayerNorm}(W_1 \mathbf{f} + b_1)) \in \mathbb{R}^{256} \\[-8pt] \nonumber \end{align}

(9) \begin{align} \mathbf{h}_2 &= \text{ReLU}(\text{LayerNorm}(W_2 \mathbf{h}_1 + b_2)) \in \mathbb{R}^{256} \\[-8pt] \nonumber \end{align}

(10) \begin{align} \mathbf{h}_{hidden} &= \text{ReLU}(W_3 \mathbf{h}_2 + b_3) \in \mathbb{R}^{128} \\[-8pt] \nonumber \end{align}

(11) \begin{align} \hat {\boldsymbol{\theta }} &= \sigma (W_{\theta } \mathbf{h}_{hidden} + b_{\theta }) \in [0,1]^{D} \\[8pt] \nonumber\end{align}

where LayerNorm ensures training stability across diverse robot morphologies, $\sigma$ is the sigmoid activation ensuring bounded outputs in $[0,1]$ , and $D = D_{\text{total}} = 330n$ is the dimensionality of the complete LF-PID parameter vector (base gains, input scaling factors, and TS consequent parameters). For a variable-DOF architecture, the output layer dimension is dynamically configured: $D_{\text{Franka}} = 2970$ for 9-DOF, $D_{\text{Laikago}} = 3960$ for 12-DOF. The predicted parameters are then denormalized to the actual control range via element-wise affine transformations calibrated to typical PID gain ranges and consequent magnitudes.

Overfitting prevention: Despite high output dimensionality, the network exhibits strong generalization (train–validation NMAE gap $\lt$ 0.5%, Appendix E) through L2 regularization ( $10^{-5}$ ), LayerNorm, and early stopping. The effective degrees of freedom are constrained by low-dimensional input features (10D) and learned manifold structure, enabling sample-efficient training on 232 samples.

The network is trained on a diverse dataset of robot configurations and their optimized LF-PID parameter vectors using weighted mean squared error loss:

(12) \begin{equation} \mathcal{L}_{meta} = \frac {1}{N}\sum _{v=1}^{N} w_v \left \|\boldsymbol{\theta }_v^* - \hat {\boldsymbol{\theta }}_v\right \|_2^2 \end{equation}

where $\boldsymbol{\theta }_v^*$ denotes optimized LF-PID parameters obtained via hybrid optimization (DE + Nelder–Mead) on each virtual robot, $\hat {\boldsymbol{\theta }}_v$ are predicted parameters, and $w_v$ are weights inversely proportional to the optimization error of each sample, prioritizing high-quality controllable configurations. The dataset consists of $N=232$ high-quality virtual robot variants (filtered from 303 generated samples) through physics-based data augmentation.

3.2.2. Stage 2: Reinforcement learning for online adaptation

The RL stage fine-tunes a compact subset of LF-PID parameters online (e.g., a small subset of TS consequents or scaling factors) to handle model uncertainties and disturbances.

State space: The state observation at time $t$ includes:

(13) \begin{equation} \mathbf{s}_t = [q_t, \dot {q}_t, e_t, \theta _t, q_{ref,t}, \dot {q}_{ref,t}] \end{equation}

where $q_t \in \mathbb{R}^n$ are joint positions, $\dot {q}_t \in \mathbb{R}^n$ are velocities, $e_t \in \mathbb{R}^n$ are tracking errors, $\theta _t$ are current LF-PID parameters being adapted, and $q_{ref,t}, \dot {q}_{ref,t}$ are reference trajectory information.

Action space (Reproducible configuration): The RL agent outputs relative adjustments to the per-joint input scaling factors $\boldsymbol{{s}} \in \mathbb{R}^{3n}$ only, while keeping base gains $\bar {\boldsymbol{{K}}}$ and TS consequents $\boldsymbol{{c}}$ fixed at meta-initialized values:

(14) \begin{equation} \mathbf{a}_t = [\Delta s_{e,1}, \Delta s_{\dot {e},1}, \Delta s_{I,1}, \ldots , \Delta s_{e,n}, \Delta s_{\dot {e},n}, \Delta s_{I,n}]^\top \in [\!-0.2, 0.2]^{3n} \end{equation}

The scaling factors are updated via multiplicative adjustment:

(15) \begin{equation} s_{\text{type},i}^{t+1} = s_{\text{type},i}^{t} \cdot (1 + \Delta s_{\text{type},i}^t), \quad \text{type} \in \{e, \dot {e}, I\}, \quad i=1,\ldots ,n \end{equation}

with clipping to ensure $s_{\text{type},i} \in [0.1, 10.0]$ to prevent numerical instability. This yields action dimension $d_a = 3n$ (27 for Franka 9-DOF, 36 for Laikago 12-DOF), consistent with the action space definition $\mathbf{a}_t \in [\!-0.2, 0.2]^{3n}$ .

Rationale for compact adaptation: Adapting only input scaling factors (rather than all 330n parameters) provides three benefits: (1) Sample efficiency – enables convergence in 1 M steps versus $\gt$ 10M for full parameter space, (2) stability preservation – avoids catastrophic forgetting of cross-platform fuzzy rules learned during meta-training, and (3) interpretability – scaling adjustments directly reflect how each joint’s error magnitude should be weighted in the fuzzy inference, facilitating deployment diagnostics. An ablation study comparing scales-only vs consequent-subset vs combined adaptation (Appendix D) demonstrates that scales-only achieves the best balance of performance gain (16.6%), sample efficiency, and training stability (0/5 divergence rate vs 2–3/5 for higher-dimensional spaces).

Reward function: We design a reward function balancing tracking accuracy, smoothness, and stability:

(16) \begin{equation} r_t = -\alpha _1 \frac {\|e_t\|}{\sqrt {n}} - \alpha _2 \frac {\|\dot {q}_t\|}{\sqrt {n}} - \alpha _3 \|\mathbf{a}_t\| \end{equation}

where $\alpha _1=10.0$ , $\alpha _2=0.1$ , $\alpha _3=0.1$ are weighting coefficients. Normalization by $\sqrt {n}$ ensures consistent scaling across different DOF platforms. The reward is clipped to $[\!-100, 10]$ to prevent numerical instability. Complete reward function design and coefficient selection rationale are detailed in Appendix A.

Training algorithm: We employ proximal policy optimization (PPO) [Reference Schulman, Wolski, Dhariwal, Radford and Klimov9] with the following hyperparameters: learning rate $1 \times 10^{-4}$ , discount factor $\gamma =0.99$ , GAE parameter $\lambda =0.95$ , entropy coefficient $0.02$ , and 8 parallel environments. Training proceeds for 1,000,000 timesteps, requiring approximately 10 minutes on a standard CPU.

3.3.1. Augmentation procedure

Starting from $K$ base robot platforms, we generate $M$ virtual variants for each base robot by perturbing physical parameters within constrained ranges:

Algorithm 1.

Physics-based data augmentation

Algorithm 1 Long description

An algorithm for physics-based data augmentation. The algorithm takes a base robot URDF and perturbation ranges as input and returns an augmented dataset. It involves sampling perturbations for mass, length, inertia, friction, and damping parameters. For each base robot, it samples these perturbations, generates a virtual robot with perturbed parameters, optimizes the LF-PID for the virtual robot, and adds the robot and its optimized parameters to the augmented dataset. The process is repeated for a specified number of iterations.

Design rationale and physical validity criteria: The perturbation ranges are carefully chosen to:

1. 1. Maintain physical plausibility: $\pm 10\%$ mass variation reflects realistic manufacturing tolerances (ISO 9283 [28]); $\pm 15\%$ inertia variation captures load/payload changes and rotational inertia uncertainty; $\pm 20\%$ friction and $\pm 30\%$ damping ranges match typical joint wear and lubrication variability [Reference Cho and Kim29]

2. 2. Ensure dynamic feasibility: Each generated variant passes four validation checks before inclusion:

   • • Positive-definiteness: Mass matrix $M(q) \succ 0$ verified at 100 random configurations

   • • Joint limit compliance: All perturbed link lengths satisfy workspace reachability constraints

   • • Simulation stability: 1,000-step forward simulation without numerical divergence (energy drift $\lt$ 1%)

   • • Controllability: Optimized LF-PID achieves tracking error $\mathcal{L}_v(\theta ^*_v) \lt 30^\circ$ (controllability threshold). This threshold is derived from a stability-based criterion: for sinusoidal reference trajectories with amplitude $A \in [0.2, 0.8]$ rad, a 30 $^\circ$ (0.52 rad) RMS error corresponds to a normalized tracking fidelity below 35%, indicating that the controller fails to maintain meaningful trajectory following. At this error level, the closed-loop system operates near the boundary of the linear stability region of the PD inner loop, where position-control mode may exhibit limit-cycle behavior. Sensitivity to this threshold choice is analyzed in Table II.

3. 3. Cover diverse dynamics: Parameter perturbations span natural frequency variations of 0.8–1.3 $\times$ baseline, ensuring heterogeneous training distribution

Quality filtering transparency: From 303 initially generated variants, 71 samples (23.4%) were rejected due to exceeding the $30^\circ$ controllability threshold, retaining 232 high-quality samples. To assess potential selection bias, we analyzed the distribution of physical parameters pre- and post-filtering:

• • Parameter distribution shift: Mean mass scaling $\alpha _{\text{mass}}$ changed from 1.00 $\pm$ 0.058 (pre-filter) to 0.998 $\pm$ 0.056 (post-filter); inertia scaling from 1.00 $\pm$ 0.087 to 0.995 $\pm$ 0.083. Kolmogorov–Smirnov tests show no significant distribution shift ( $p \gt 0.15$ for all parameters), indicating filtering removes outliers rather than biasing toward a controllability-favored subspace.

• • Feature diversity preservation: Post-filtering dataset retains 94.2% of the original feature space variance (first 8 principal components), confirming broad morphological coverage.

• • Error distribution: Filtered samples exhibit tracking errors $13.9\pm 8.4^\circ$ (mean $\pm$ std), spanning easy ( $5^\circ$ ) to challenging ( $29^\circ$ ) configurations, providing gradient for meta-learning.

Table II.

Sensitivity of controllability threshold on data quality and downstream meta-learning performance (Franka 7-DOF). Retained = number of virtual samples passing the filter out of 303 total. Downstream MAE is the meta-learning prediction error (degrees) on the held-out test set.

Table II Long description

A table with five rows and five columns. The columns are labeled Threshold, Retained, Mean error, NMAE, and Downstream MAE. The rows are labeled with different threshold values: 10 degrees, 20 degrees, 30 degrees, 40 degrees, and No filter. The values in the table are as follows: Row 1: Threshold, 10 degrees; Retained, 98; Mean error, 6.2 plus or minus 2.1; NMAE, 52.8; Downstream MAE, 9.14. Row 2: Threshold, 20 degrees; Retained, 178; Mean error, 9.8 plus or minus 4.7; NMAE, 48.3; Downstream MAE, 8.02. Row 3: Threshold, 30 degrees; Retained, 232; Mean error, 13.9 plus or minus 8.4; NMAE, 47.1; Downstream MAE, 7.51. Row 4: Threshold, 40 degrees; Retained, 268; Mean error, 16.5 plus or minus 10.2; NMAE, 49.6; Downstream MAE, 7.89. Row 5: Threshold, No filter; Retained, 303; Mean error, 19.1 plus or minus 12.8; NMAE, 55.4; Downstream MAE, 8.73.

Threshold sensitivity analysis. To justify the $30^\circ$ controllability threshold, we perform a sensitivity study by varying the threshold from $10^\circ$ to $40^\circ$ (with a no-filter baseline) and evaluating the impact on dataset composition and downstream meta-learning performance. Table II summarizes the results.

The $30^\circ$ threshold achieves the optimal balance: it retains sufficient diversity (232 samples, 76.6% of total) while excluding pathologically uncontrollable configurations. Stricter thresholds ( $10^\circ$ , $20^\circ$ ) aggressively prune challenging but learnable samples, reducing dataset diversity and degrading downstream performance. Conversely, relaxing to $40^\circ$ or removing the filter introduces noisy, uncontrollable samples that corrupt the meta-learning landscape. The resulting U-shaped curve in downstream MAE confirms that $30^\circ$ sits at the minimum of the quality–diversity trade-off.

Reproducibility and detailed protocols: All validation procedures (mass matrix positive-definiteness, energy drift tests, KS distribution comparisons, PCA variance retention) are fully documented with computational parameters, random seeds, and API calls in Appendix B to enable independent verification.

Algorithm 2.

Hybrid LF-PID optimization for virtual robots

Algorithm 2 Long description

The algorithm describes a two-stage process for optimizing LF-PID parameters in virtual robots. The first stage involves a global search using differential evolution. It initializes a population of parameters uniformly within specified bounds. For each generation, it selects random indices for mutation and crossover operations. The trial parameters are evaluated using a PyBullet simulation, and the best parameters are selected for the next generation. This process repeats for 15 generations. The second stage involves local refinement using the Nelder-Mead method. It initializes a simplex around the best parameters from the global search. The simplex vertices are sorted by objective value, and operations like reflection, expansion, contraction, or shrinkage are applied. The simplex is updated based on Nelder-Mead rules, and the best vertex in the final simplex is returned as the optimal parameters.

3.3.2. Hybrid optimization strategy for LF-PID parameters

For each virtual robot, we employ a two-stage hybrid optimization strategy that combines global search and local refinement to find optimized LF-PID parameters efficiently and accurately.

Optimization objective: For each virtual robot $r_v$ , we aim to minimize the root mean square joint tracking error over a test trajectory:

(17) \begin{equation} \mathcal{L}_v(\theta ) = \sqrt {\frac {1}{T} \sum _{t=1}^{T} \sum _{i=1}^{n} \left ( q_{ref,i}(t) - q_i(t; \, \theta ) \right )^2} \end{equation}

where $T$ is the trajectory length (typically $T=2000$ steps at 240 Hz control frequency), $q_{ref,i}(t)$ is the reference trajectory for joint $i$ , and $q_i(t; \, \theta )$ is the actual trajectory achieved with LF-PID parameters $\theta$ .

Stage 1 – Global search via differential evolution: We employ differential evolution (DE) [Reference Storn and Price25], a population-based stochastic optimizer, for robust global exploration of the LF-PID parameter space:

(18) \begin{equation} \theta ^*_{global} = \arg \min _{\theta \in [\theta _{min}, \theta _{max}]} \mathcal{L}_v(\theta ) \end{equation}

with population size $N_{pop}=8$ and $N_{iter}=15$ iterations. DE is particularly effective for non-convex, noisy objective functions common in robotic control, as it maintains a population of candidate solutions and evolves them through mutation, crossover, and selection operations, effectively avoiding local minima.

Stage 2 – Local refinement via Nelder–Mead: To achieve high-precision convergence from the globally optimal region identified by DE, we apply local refinement using the Nelder–Mead simplex method [Reference Nelder and Mead26]:

(19) \begin{equation} \theta ^*_v = \arg \min _{\theta } \mathcal{L}_v(\theta ), \quad \text{from } \theta ^*_{global} \end{equation}

with $N_{iter}^{polish}=20$ polishing iterations.

Table III.

Pilot comparison of optimization algorithms on 20 representative virtual robots. All methods use identical simulation budget (2000-step trajectory evaluation). “Obj.” denotes the final objective $\mathcal{L}_v(\theta ^*)$ in degrees. “Evals” denotes the number of objective function evaluations (each requiring one forward simulation).

Table III Long description

A table comparing optimization algorithms on virtual robots. The table has four rows and four columns. The columns are labeled Method, Obj. (°), Evals, Time (s), and Notes. Row 1: DE + NM (ours), 13.9 ± 4.2, 128, 42 ± 12, Robust; no gradient needed. Row 2: CMA-ES, 13.2 ± 3.8, 210, 68 ± 18, Better mean; higher cost. Row 3: L-BFGS-B, 18.7 ± 9.1, 85, 28 ± 8, Fast; frequent local minima. Row 4: PSO (N = 20), 14.5 ± 5.0, 400, 125 ± 30, Slowest; marginal gain.

This hybrid approach leverages complementary strengths: DE provides robust global search while Nelder–Mead offers rapid local convergence. In our implementation using scipy.optimize.differential_evolution, this is achieved via the polish=True parameter, which automatically applies local optimization after the DE phase.

Rationale and efficiency: Pure DE requires many iterations (typically $\gt$ 200) for high-precision convergence, while pure local methods risk converging to poor local optima. The hybrid approach achieves both global robustness and local precision efficiently, completing optimization for each virtual robot in 30–60 s on a single CPU core. This optimization process provides ground-truth LF-PID parameters for meta-learning training. We filter out samples with optimization error $\mathcal{L}_v(\theta ^*_v) \gt 30^\circ$ to ensure data quality, retaining 232 high-quality samples from 303 initially generated variants. When parallelized across 23 cores, the entire 303-sample augmentation completes in approximately 5 minutes. Detailed parameters for data augmentation and LF-PID optimization are provided in Appendix A.

Choice of optimizer – comparison with alternatives: We chose DE + Nelder-Mead (NM) for its simplicity and favorable properties in our specific setting, but acknowledge that more advanced optimizers exist. We conducted a pilot comparison on 20 representative virtual robots (10 Franka-type, 10 KUKA-type) to evaluate three alternatives: CMA-ES [Reference Hansen and Ostermeier39], L-BFGS-B [Reference Byrd, Lu, Nocedal and Zhu40], and PSO [Reference Trelea8] (Table III).

Key findings: (1) CMA-ES achieves marginally better objective values ( $-$ 5%) but requires 64% more function evaluations; since each evaluation involves a full PyBullet simulation, this translates to proportionally higher wall time. (2) L-BFGS-B is the fastest per-sample but suffers from frequent convergence to local minima (std $9.1^\circ$ vs $4.2^\circ$ ), yielding 34% worse mean quality – problematic because poor-quality training labels directly degrade meta-learning. (3) PSO offers no advantage over DE at 3 $\times$ computational cost. Given that (i) DE + NM serves as an offline data-generation tool (not a deployment-time component), (ii) the entire 232-sample generation takes only $\sim$ 5 minutes with parallelization, and (iii) CMA-ES’s marginal quality gain does not justify doubling offline preprocessing time, the DE + NM choice represents a practical and effective trade-off. We note that CMA-ES would be a reasonable alternative if higher-quality training labels are desired and offline computation budget is less constrained.

Regarding L-BFGS-B for high-dimensional spaces (330 $n$ ): While L-BFGS-B is well-suited for smooth, high-dimensional optimization, our LF-PID objective $\mathcal{L}_v(\theta )$ is evaluated through black-box simulation (PyBullet forward dynamics) and does not provide analytical gradients. Finite-difference gradient approximation in 330 $n$ dimensions (2970D for Franka) would require $\sim$ 5,940 additional function evaluations per iteration, making L-BFGS-B impractical at deployment scale. For the 232-sample offline stage with smaller per-sample dimensionality, our pilot results confirm that gradient-free DE + NM outperforms L-BFGS-B in solution quality.