Probability-Based Manipulability Score Using Forward-Kinematics Voxelization

2.1 Calculation of PBMS

The end-effector position of a manipulator can be computed using forward kinematics. When computing forward kinematics with randomly sampled joint angles, positions that can be reached through a greater number of joint angle combinations tend to produce denser point clouds, indicating a higher probability of achieving diverse manipulator configurations and orientations. In a voxelized workspace, voxels with denser point clouds are more likely to reflect regions where diverse configurations and orientations can be achieved, suggesting higher potential performance.

When a sufficiently large number of joint configurations are randomly sampled, this distribution approximates the result of uniformly sampling each joint over its range. In this case, the total number of possible configurations grows exponentially with the number of joints, forming a set of points in the form of n^D, where n is the number of sampled angles per joint and D is the number of joints. These points represent possible end-effector positions and are interpreted probabilistically. In this context, a sufficiently large sample size refers to the point at which further sampling no longer produces meaningful changes in the resulting distribution or PBMS value, indicating convergence of the Monte Carlo estimation^[20].

Adding one joint increases the total number of configurations by a factor of n. If we assume that this increase is evenly distributed across all voxels, it becomes equivalent to multiplying the number of points in each voxel by n. Therefore, an n-fold increase in the number of points within a voxel can be interpreted as the performance gain resulting from one additional degree of freedom. This assumption enables the use of a logarithmic transformation to quantify relative performance differences between manipulators at each voxel in terms of degrees of freedom.

Let N = n^D be the total number of sampled configurations. The base of the logarithm, a, is determined from N and the maximum score M defined by the user. It is chosen so that the increase in sample count associated with one additional degree of freedom, N^(1/D), maps directly to the corresponding score increment M/D. This increment represents the score contribution of one degree of freedom. The maximum number of points observed in a voxel is denoted as M_c. The bias term k shifts the score range so that the voxel with M_c attains the user-defined maximum score M in (1). This shifting ensures that all scores are normalized to a common baseline, enabling fair comparison across manipulators with different degrees of freedom. Based on these definitions, the number of points C within a specific voxel is counted, and the corresponding score S is computed as shown in (1).

When PBMS is used to compare the potential performance of manipulators with different degrees of freedom, a reference manipulator is first selected to define the baseline sampling density. A manipulator with more degrees of freedom inherently possesses a larger configuration space and can achieve a greater variety of end-effector positions and orientations. However, if both manipulators are evaluated with the same number of sampled configurations, the higher-DOF manipulator cannot fully express this structural advantage because its larger configuration space becomes under-sampled. As a result, the performance differences arising from the degrees of freedom would not be properly reflected in the PBMS.

To address this issue, the number of samples for the test manipulator must be scaled according to its degrees of freedom so that the sampling density per joint dimension remains consistent with the reference manipulator. This adjustment is performed using (2), where N_ref and D_ref denote the sample count and degrees of freedom of the reference manipulator, while N and D correspond to the test manipulator. By scaling N according to (2), the PBMS computed from (1) appropriately captures the performance differences that originate from the manipulators’ kinematic redundancies rather than from unequal sampling densities.

To ensure that the sample size is sufficiently large, sampling is performed iteratively in small batches, and convergence is monitored. In this evaluation, only the scores within the ISO cube, where most tasks are typically performed, are considered. A summary of this procedure is described in Algorithm 1.

Algorithm 1

PBMS computation Prcedure

Although the PBMS value in (1) is computed solely from the position-based count C, the sampling procedure also evaluates whether the range of feasible end-effector orientations has been sufficiently explored. A denser positional point cloud implicitly reflects a wider variety of achievable orientations, which is consistent with the notion of higher potential performance within a voxel. Therefore, orientation information is not used directly in the final PBMS score, but instead serves as a convergence criterion to prevent biased or insufficient angular coverage during sampling.

To evaluate convergence in orientation, it was assumed that each position voxel contains a corresponding orientation voxel to consider reachable orientation. Each orientation voxel is defined in a 3D space structured along the roll, pitch, and yaw axes, representing orientations for the corresponding position. At each iteration, the orientation is obtained from forward kinematics and mapped into this orientation voxel. Then, the number of reachable orientation voxel, which has any corresponding combination of roll, pitch, yaw angle for the given position, is counted. This number ranges from 1 to the cube of the number of divisions along each rotational axis. For instance, if each axis is divided into 10 intervals, the maximum number of reachable orientation voxels becomes 1,000. These bounds correspond to two extremes: a singular configuration such as a fully stretched manipulator with only one reachable orientation, and configurations where all orientations are accessible at a given position.

Convergence was evaluated for each voxel using a change ratio, defined as the relative change in results before and after each iteration. The change ratios corresponding to position and orientation, denoted as e_p and e_o, were calculated as shown in (3).

Here, C represents the number of points accumulated in each position voxel. O represents the numbers of reachable orientation voxels of orientation voxels for a given position voxel. C and O represent the updated values after the current iteration, while $$ C^{\prime }$$ and $$ O^{\prime }$$ are their corresponding values from the previous step. The indices i, j, k refer to the spatial dimensions x, y, z respectively. At each iteration, e_p and e_o, defined as the maximum change ratio across all position voxels, were calculated. The simulation was considered to be converged when e_p and e_o are below a certain threshold for five consecutive iterations.

3.1 Comparison of PBMS

A simulation was conducted to compare three manipulators using the proposed PBMS. All manipulators share an identical workspace and were evaluated within a common ISO cube to ensure fair comparison. These manipulators were implemented on a modular testbed named the MAT (Modularized Arm Testbed), which allows reconfiguration of joint types and link arrangements to facilitate comparative analysis. The first manipulator, a 4-DOF model, was used as the reference structure. The second manipulator has the same 4-DOF but features a different joint arrangement to investigate the effect of kinematic configuration on manipulability. The third manipulator is a 3-DOF model, included to analyze the influence of kinematic redundancy on workspace distribution and PBMS values. The geometric configurations of the MAT used in this study are illustrated in [Fig. 1]. From left to right, the models correspond to the reference manipulator, the 4-DOF test manipulator, and the 3-DOF test manipulator. Each actuator is highlighted in blue to indicate the locations of the active joints, and the link dimensions and joint spacings shown in the figure were used consistently throughout the forward-kinematics simulations and workspace analyses. Joint limits were not applied in order to focus on the intrinsic kinematic characteristics of each configuration.

[Fig. 1]

CAD models of the MAT (Modularized Arm Testbed) manipulators used in this study. From left to right: reference manipulator, 4-DOF test manipulator, and 3-DOF test manipulator. Blue components represent actuators, and the annotated dimensions (in mm) indicate the joint-to-joint spacing applied in the PBMS simulations.

To visualize the spatial distribution of scores, the computed PBMS values were mapped to each voxel in the workspace, as illustrated in [Fig. 2]. After PBMS computation, the workspace was color-coded voxel by voxel according to the relative score of each region. To facilitate interpretation, three transverse planes and one sagittal plane were used to slice the 3D space, and the voxel colors corresponding to each plane were projected onto 2D plots. In this simulation, the transverse planes were defined at z = 0 mm, 120 mm, and 240 mm corresponding to the shoulder, half-elbow and elbow levels, respectively, while the sagittal plane was defined at x = 0 mm. Each 2D plot in [Fig. 3] shows the score distribution of voxels on these planes for each manipulator. The shoulder plane corresponds to the position of the second joint, where the axes of the first, second, and third joints intersect in the reference structure. The elbow plane represents the end face of the ISO cube, approximately corresponding to the distance from the base to the elbow joint, while the half-elbow plane is defined at half of this distance. The corresponding heights of the shoulder, half-elbow, and elbow planes were set to 0 mm, 120 mm, and 240 mm, respectively, measured from the origin of the base frame, which is defined as the origin with coordinates, x = 0, y = 0 and z = 0. The scores for three manipulators were computed per voxel using equation (1). Since the 3-DOF manipulator had different degrees of freedom, the number of samples was adjusted to account for a one-degree-of-freedom difference. The maximum score of the reference manipulator was set to 100. Accordingly, the score difference corresponding to the effect of one degree of freedom was calculated as 25.0 points, derived by dividing the user-defined maximum score M by the number of degrees of freedom D.

[Fig. 2]

Illustration of voxel color distribution and slicing plane used for visualization of PBMS. Three transverse planes correspond to the shoulder, half elbow and elbow plane, respectively.

[Fig. 3]

Relative PBMS scores for each manipulator with respect to reference structure: (a) 4-DOF manipulator with different joint arrangement, and (b) 3-DOF manipulator. The heights of the transverse planes at the shoulder, half elbow and elbow were set to 0 mm, 120 mm, and 240 mm, respectively, measured from the base frame.

The evaluation space was defined as a bounding cube with sides of 1,100 mm, divided into 40 intervals along each axis to produce position voxels of 27.5 mm. However, scoring was performed only within the ISO cube, the largest inscribed cube fully enclosed within the manipulator’s workspace. In this case, the ISO cube had a side length corresponding to 18 intervals, or 495 mm. For the orientation space, each rotational axis was divided into 10 intervals, resulting in orientation voxels with a resolution of 36 degrees per axis. [Fig. 3] presents the relative PBMS scores for each manipulator compared to the reference structure. Two comparison cases are shown: (a) a 4-DOF manipulator with a different joint arrangement and (b) a 3-DOF manipulator. Each subplot visualizes the voxel-wise score differences using a color map across the same representative planes. The delta score indicates how much higher or lower each manipulator scores relative to the reference at each voxel. The color map is normalized based on the absolute score difference between each test manipulator and the reference. Purple regions represent areas where the relative score difference is negative, while orange regions indicate positive score differences. Voxels with values close to zero are shown in gray, representing regions where the test and reference manipulators exhibit nearly identical PBMS values. In this color scheme, higher color intensity corresponds to larger deviations in magnitude from the reference score, regardless of sign. As shown in [Fig. 3(a)], the 4-DOF test manipulator exhibited score distributions that were comparable to or slightly lower than the ideal value in the vicinity of the base frame. At the shoulder plane, regions within approximately 124 mm from the base frame appeared in gray, indicating that the scores were close to the expected value under the condition of equal degrees of freedom and therefore did not represent a significant loss. At the half-elbow plane, regions within about 80 mm from the base frame showed a similar pattern of purple or gray coloring, suggesting only a small deviation from the reference manipulator. In contrast, at the elbow plane, the overall distribution was visualized in orange or gray, implying slightly higher scores than the reference and indicating that performance differences became more noticeable toward the outer workspace.

The sagittal-plane plots also reveal that this pattern is distributed in a nearly spherical manner. The purple regions are mainly observed around the intersection point of the first through third joint axes in the shoulder area, implying that the score degradation of the test manipulator occurs near this region of the reference structure. In contrast, regions closer to the outer boundary of the workspace, particularly near the faces of the ISO cube, exhibit higher score distributions, indicating that the test manipulator performs relatively better toward the limits of its reachable workspace.

For the 3-DOF test manipulator, as illustrated in [Fig. 3(b)], the overall score distribution appeared predominantly in purple across all planes, indicating that most voxels exhibited lower PBMS values than the reference. Because this color corresponds to the expected decrease of approximately 25 points, the results suggest that the 3-DOF manipulator is generally affected by the degree-of-freedom difference throughout the workspace. The color distribution was mainly organized along the z-axis, rather than showing clear variations between transverse planes. Voxels located close to the z-axis appeared more grayish, indicating regions where the score difference approached near zero. This occurs because, when the reference manipulator aligns its first and third joint axes parallel to the z-axis, a local singularity arises, effectively reducing one degree of freedom. Consequently, the 3-DOF manipulator exhibited similar performance to the reference in this configuration.

These results suggest that, in the comparison between the two 4-DOF manipulators, the reference structure is more suitable for tasks near the base where the joint axes intersect, while the 4-DOF test manipulator performs better near the edges of the ISO cube, particularly around singularity-like configurations. The 3-DOF manipulator, on the other hand, shows relatively favorable score near the z-axis, where the singularity of the reference structure occurs. Nevertheless, both simulations indicate that the overall score variations remain small within the workspace when interpreted in light of the expected reduction due to the configuration and degree-of-freedom difference.

3.2 Convergence Criteria

[Fig. 4] presents the change ratio e_p and e_o computed at each iteration during the random sampling of end-effector positions and orientations, as performed using the equation (1) and (3) described in Algorithm 1. In each iteration, 5,000,000 random samples were used, and the threshold for determining convergence was set to 1%. Because the manipulators used in this study have three and four degrees of freedom, convergence was evaluated based solely on the position change ratio e_p. In the figure, the blue, green, and red lines correspond to the reference manipulator, the 4-DOF test manipulator, and the 3-DOF test manipulator, respectively. A total of 27 iterations were performed until convergence was achieved, resulting in 135,000,000 samples being used in the simulation. The simulation was conducted on an NVIDIA GeForce RTX 3060 GPU and completed in approximately 2 minutes.

[Fig. 4]

Change ratio ep and eo across iterations.

4.1 Experimental setup

The motion capture experiment was conducted using the modular testbed MAT to compare the tendencies between the simulation and the actual measurements, as show in [Fig. 5]. The optical motion-capture system consisted of VICON T160 cameras (16 MP resolution, 120 fps at full resolution) operating under Vicon Nexus software. The system employed Vicon’s proprietary gray-scale centroid-fitting algorithm for precise 2D tracking and real-time 3D reconstruction. Marker positions were captured with a typical spatial accuracy of approximately ±0.5 mm within the calibrated workspace. The MAT testbed was fabricated via 3D printing using PLA reinforced with carbon fiber (PLA-CF) to achieve a lightweight yet rigid structure. Actuation was provided by DYNAMIXEL PRO servos (H54-200-S500-R and H54-100-S500-R) to ensure accurate and repeatable motion across configurations.

[Fig. 5]

Manufactured model of MAT and experimental setup for motion capture. The marker at end-effector is tracked through the VICON cameras.

Because of the intrinsic limitations of the motion-capture process, it was not feasible to achieve the same large-scale sampling as in the simulation. Therefore, the experiment focused on comparing the overall tendencies under equivalent conditions, even though the end-effector could not reach the interior of the ISO cube. To explore a wide range of joint configurations within a limited time, rough joint limits of [-90°, ~90°] were assigned to all joints, and a relatively high joint velocity of approximately 15–25 deg/s was applied. Each configuration was measured for 300 seconds at 100 Hz, yielding a total of 30,000 samples of end-effector positions per configuration.

To ensure sufficient variation among joint combinations, each joint was actuated using a triangular-wave motion pattern with different periods. The periods were chosen such that their least common multiple exceeded the total measurement duration, preventing repetition of joint combinations within a single trial, as summarized in [Table 1]. These settings allowed the experiment to qualitatively replicate the sampling tendencies observed in the simulation, while maintaining comparable kinematic coverage within the reachable workspace.

[Table 1]

Experimental setup for each joint

4.2 Experimental Results and Analysis

To compare the tendencies between the simulation and the physical experiment, the Probability-Based Manipulability Score (PBMS) was computed within voxels of identical size. [Fig. 6] shows the PBMS results obtained from the motion-capture data of the MAT testbed. The motion trajectories recorded during the experiment formed a thin spherical shell centered around the base frame, representing the reachable surface swept by the end-effector. As in the simulation, the score distributions were analyzed on two transverse planes, the shoulder and the elbow heights, and on an additional upper plane (360 mm height) between the elbow level and the workspace boundary, where sufficient data were available.

[Fig. 6]

Relative PBMS scores for each manipulator with respect to reference structure from motion capture data: (a) 4-DOF manipulator with different joint arrangement, and (b) 3-DOF manipulator. The heights of the transverse planes at the shoulder, elbow and upper were set to 0 mm, 240 mm, and 360 mm, respectively, measured from the base frame.

The color map used for visualizing PBMS differences followed the same convention as in the simulation: orange indicates higher scores, purple indicates lower scores relative to the reference manipulator, and gray represents smaller deviations from the expected value. Although the motion covered a broad range of reachable postures, the plots exhibit asymmetric and sparse distributions because of the limited sample count relative to the dense voxel grid and the restricted joint-limit range.

For the 4-DOF test manipulator, the simulation results in [Fig. 3] showed that the score distribution tended to appear purple near the origin, gradually shifting toward orange near the ISO-cube boundaries, and finally turning gray again near the cube’s corner regions. The experimental results in [Fig. 6] demonstrate a similar trend: as the end-effector approaches the outer faces of the workspace, orange-colored regions become more prevalent, while purple areas emerge near the corners corresponding to the workspace boundaries.

For the 3-DOF test manipulator, the simulation indicated that the reference structure’s z-axis alignment produced a local singularity, yielding relatively advantageous scores in the surrounding region. This effect appeared as an hourglass-shaped distribution in the sagittal plane, with predominantly light-purple regions elsewhere. The experimental results in [Fig. 6] show a consistent tendency: orange-colored regions increasingly appear near the z-axis as the height increases toward the workspace boundary, which is particularly evident in the sagittal-plane plots where orange zones cluster along the central vertical axis.

[Fig. 7] presents snapshots of the experiment captured every 4 seconds, illustrating the motion of the MAT throughout the motion-capture trial. The motion began from a fully stretched configuration along the z-axis and proceeded with cyclic joint rotations following the periods and angular velocities listed in [Table 1]. These snapshots confirm that the MAT swept through a wide variety of trajectories within the reachable workspace during the experiment. Overall, the experimental PBMS distributions exhibited partially consistent spatial tendencies with the simulation results, particularly within the regions frequently reached by the end-effector.

[Fig. 7]

Sequential snapshots of the MAT obtained every 4 s using the motion-capture system. Starting from a fully stretched configuration along the z-axis, each joint rotated with the predefined periods and angular velocities listed in [Table 1].

5.1 Experimental Limitations and Insights

Because only 30,000 end-effector samples were available while the workspace was discretized into 40³ = 64,000 voxels, many voxels were left unsampled, particularly those located near the workspace boundaries. Consequently, quantitative comparison was feasible only in frequently visited regions. Furthermore, due to the joint-limit constraints in the experimental setup, none of the configurations reached the ISO cube; thus, the comparison focused on the overall spatial tendencies within the reachable workspace rather than complete volumetric coverage.

Despite these limitations, the experiment successfully reproduced the primary spatial tendencies observed in the simulation. The consistent PBMS distributions within the densely visited regions indicate that the proposed metric remains robust even under sparse sampling and restricted motion ranges. This result suggests that PBMS can serve as a reliable indicator for real-world manipulators where exhaustive sampling is impractical.

The experimental sparsity also reflects an inherent property of probability-based evaluations: local probability distributions can approximate global manipulability trends when the sampling diversity is maintained. Therefore, PBMS provides a scalable framework for both high-fidelity simulation and physical validation without requiring complete coverage of the configuration space. Future work will aim to enhance the experimental setup by expanding joint motion ranges, integrating end-effector orientation tracking, and increasing the sampling rate to improve volumetric coverage. The modular design of the MAT also allows further experiments with varied link geometries and actuation schemes, which will help generalize the applicability of PBMS to different manipulator types and redundancy levels.

5.2 Sparse Sampling Robustness Analysis

To determine how reliably the PBMS captures workspace-level tendencies under limited sampling, as observed in Section 5.1, additional simulations were conducted using significantly fewer samples than the fully converged case (135M samples). This analysis also identifies the approximate sampling threshold at which meaningful spatial trends begin to degrade. Five sampling conditions were tested: 5M, 500K, 50K, 30K (matching the MAT experiment), and 5K samples. For each case, the resulting PBMS distribution was compared against the converged reference using three criteria: RMSE (Root Mean Squared Error), Spearman’s rank correlation coefficient^[21], and Kendall’s rank correlation coefficient^[22].

RMSE was used to quantify the magnitude of numerical deviation from the fully converged PBMS distribution. Because PBMS is fundamentally a voxel-wise probability-based metric, RMSE provides a direct measure of how much each voxel’s estimated score drifts as the sample size becomes sparse.

Spearman’s rank correlation coefficient was used to evaluate whether the global monotonic trend of voxel-wise PBMS scores was preserved under reduced sampling density. As expressed in (4), Spearman’s coefficient computes the correlation between the rank-transformed vectors of the reference PBMS distribution X and the sparse-sampled distribution Y, where R(X) and R(Y) denote their respective rank-transformed values. The term $$ \mathrm{c}\mathrm{o}\mathrm{v}\left(\cdot \right)$$ is the covariance between the rank vectors, while σ_R(X) and σ_R(Y) represent their standard deviations. A value close to 1 indicates that the relative ordering of voxel scores is largely maintained.

Kendall’s rank correlation coefficient was additionally evaluated to quantify the stability of pairwise voxel ordering. As shown in (5), Kendall’s τ is defined using the number of concordant pairs C and discordant pairs D among all $$ \frac{1}{2}n\left(n-1\right)$$ possible voxel pairs. Because τ directly measures the proportion of voxel pairs whose ordering changes between the two PBMS distributions, it serves as a sensitive indicator of local ranking consistency. Together, Spearman’s and Kendall’s coefficients provide complementary insights into the trend preservation and ordering robustness of PBMS under sparse sampling. The results of these analyses are summarized in [Table 2].

[Table 2]

Robustness evaluation of PBMS under sparse sampling conditions. RMSE quantifies numerical deviation from the fully converged distribution (135M samples), while Spearman’s and Kendall’s rank correlation coefficients measure global trend preservation and voxel-wise ordering consistency, respectively.

Spearman’s rank correlation coefficient remains above 0.9 for sampling sizes down to 500K, indicating that the overall monotonic trend of voxel-wise PBMS scores is largely maintained. Kendall’s coefficient at 500K is 0.779, implying that approximately 11% of voxel pairs exhibit discordant ordering. Because Kendall’s metric is more sensitive to local ranking inconsistencies, this value suggests that most ranking changes occur among voxels with similar scores and therefore do not significantly distort the global spatial tendency.

Both correlation metrics exhibit a sharp decline at around 50K samples. This indicates that, below this threshold, ranking inconsistencies arise across wide regions of the workspace rather than being confined to locally similar voxels, marking the onset of structural breakdown in the PBMS trend. At 5K samples, Spearman’s coefficient drops to 0.333, implying that almost no meaningful spatial tendency remains.

Overall, these results indicate that meaningful trend preservation empirically requires a sampling size on the order of roughly ten times the number of voxels within the ISO cube. RMSE increases gradually as the sampling size decreases, reflecting greater deviation from the converged distribution. Moreover, smaller sample counts reduce the logarithmic base a in (1), lowering the effective score resolution and amplifying numerical discrepancies in the PBMS values.

5.3 Extension to Real-Robot Evaluation

To further validate the applicability of PBMS to real robotic systems, a simulation was conducted using the full-size humanoid robot TOCABI^[23] and a corresponding 7-DOF test manipulator. [Table 3] summarizes the Denavit–Hartenberg parameters of both manipulators. The test manipulator was derived by removing the shoulder yaw joint (the 4th joint) from TOCABI’s arm to investigate the effect of one fewer degree of freedom. To compensate for payload differences, the link lengths were slightly modified compared to the reference structure.

[Table 3]

Denavit–Hartenberg (DH) parameters of two single arms: (a) test manipulator, (b) reference manipulator (TOCABI arm)

[Fig. 8] presents the spatial distribution of PBMS for both manipulators. The color map was based on an expected reduction of 12.5 points, corresponding to the difference of one degree of freedom: orange represents higher scores, purple indicates lower scores, and white denotes nearly equivalent values. Transverse planes correspond to the shoulder, chest, and waist levels, positioned at 66 mm, –104 mm, and –274 mm, respectively, relative to the base frame.

[Fig. 8]

Relative PBMS scores for test manipulator with respect to humanoid robot TOCABI. The heights of the transverse planes at the shoulder, chest, and waist height were set to 66 mm, -104 mm, and -274 mm respectively, measure from the base frame.

As shown in [Fig. 8], the test manipulator exhibited slightly higher or comparable scores within approximately 50 mm of the base across all planes, represented by orange or white regions. In contrast, purple regions appeared between 50 mm and 250 mm at shoulder height, indicating reduced manipulability near the torso. At farther distances (250–400 mm), orange regions became dominant, suggesting improved performance in the arm-extended postures. The sagittal-plane visualization revealed that these regions formed a nearly spherical distribution. Notably, purple zones concentrated near the intersection of joints 2 through 4 in the reference manipulator, implying that the score degradation primarily occurs near the shoulder joint.

These findings demonstrate that the proposed PBMS can effectively distinguish the spatial advantages and disadvantages between manipulators with different kinematic structures, even when applied to a full-scale humanoid robot. This confirms the metric’s potential for practical evaluation and design optimization of redundant manipulators in real-world robotic systems.