Dynamic Predictive Model with Parametric Bias can be expressed by the following equation.

where t is the current time step, s is the state of the manipulated object and robot, u is the control command to the robot body, p is the parametric bias (PB), and h_dpmpb is the function representing the time series change in the state of the manipulated object and robot due to the control command. We use the state of the cloth and the state of the robot for s in this study of cloth manipulation. For the cloth state, we use z_t, which is the compressed value of the current image I_t by using AutoEncoder (Hinton and Salakhutdinov, 2006). Although the robot state should be different for each robot, we use f_t and l_t for the musculoskeletal humanoid used in this study (f_t and l_t represent the current muscle tension and muscle length). Thus, stT=(ztTftTltT). Also, we set uT=(θref,Tkref,T) (θ^ref represents the target joint angle and k^ref represents the target body stiffness value). The details of the target joint angle and body stiffness are described in Section 2.2. The parametric bias is a value that can embed implicit differences in dynamics, which are common for the same material and different for different materials. By collecting data using various cloth materials, information on the dynamics of the material of the manipulated object is embedded in p. Note that each sensor value is used as input to the network after normalization using all the obtained data.

In this study, the DPMPB consists of 10 layers: 4 fully-connected layers, 2 LSTM layers (Hochreiter and Schmidhuber, 1997), and 4 fully-connected layers, in order. The number of units is set to {N_u + N_s + N_p, 300, 100, 30, 30 (number of units in LSTM), 30 (number of units in LSTM), 30, 100, 300, N_s} (where N_{{u, s, p}} is the number of dimensions of {u, s, p}). The activation function is hyperbolic tangent and the update rule is Adam (Kingma and Ba, 2015). Regarding the AutoEncoder, the image is compressed by applying convolutional layers with kernel size 3 and stride 2 five times to a 128 × 96 binary image (the cloth part is extracted in color), then reducing the dimensionality to 256 and 3 units, in order, using the fully-connected layers, and finally restoring the image by the fully-connected layers and the deconvolutional layers. For all layers exempting the last layer, batch normalization (Ioffe and Szegedy, 2015) is applied, and the activation function is ReLU (Nair and Hinton, 2010). The dimension of p is set to 2 in this study, which should be sufficiently smaller than the number of cloth materials used for training. The execution period of Equation (1) is set to 5 Hz.

Regarding the musculoskeletal humanoid with variable stiffness mechanism used in this study, we will describe the details of θ^ref and k^ref. In this study, we use the following relationship of joint angle θ, muscle tension f, and muscle length l (Kawaharazuka et al., 2018).

h_bodyimage is trained using the actual robot sensor information. When using this trained network for control, the target joint angle θ^ref and the target muscle tension f^ref are determined and the corresponding target muscle length lref=hbodyimage(θref,fref) is calculated. However, since this value is the muscle length to be measured, not the target muscle length, l^send, which takes into account the muscle elongation in muscle stiffness control (Shirai et al., 2011), is sent to the actual robot (Kawaharazuka et al., 2018) in practice.

Here, we consider how f^ref is given. Normally, f^ref should be the value required to realize θ^ref. On the other hand, in Kawaharazuka et al. (2018), a constant value f^const is first given to f^ref for all muscles to achieve θ^ref to a certain degree, and then the current muscle tension f is given as f^ref to achieve θ^ref more accurately. The body stiffness can be changed to some extent depending on the value of f^const sent at the beginning. In addition, since the relationship between s and u is eventually acquired by learning in this study, it is not necessary to realize θ^ref precisely. Therefore, we use f^const as the body stiffness value k^ref included in u and operate the robot with lref=hbodyimage(θref,fconst). By acquiring training data while changing f^const, we can perform dynamic cloth manipulation considering the body stiffness value.

In the following, we show the results of experiments to see how much the operational stiffness of the arm changes with f^const. Figure 3 shows the operational stiffness ellipsoid of the left arm in the sagittal plane when θ^ref is set as the state with the elbow of the left arm of Musashi-W bent at 90 degrees and f^const is changed to {10, 30, 50, 70} [N]. The graph shows the displacement of the hand when a force of 1 N is applied from all directions. It can be seen that the size of the stiffness ellipsoid changes greatly with the change in f^const. Therefore, we use f^const as the body stiffness value k^ref in this study.

By setting u randomly, or operating the robot with GUI or VR device, we collect the data of s and u. For a single, coherent manipulation trial k, performed with the same cloth, the data D_k = {(s₁, u₁), (s₂, u₂), ⋯ , (s_Tk, u_Tk)} (1 ≤ k ≤ K, where K is the total number of trials and T_k is the number of time steps for the trial k). Then, we obtain the data D_train = {(D₁, p₁), (D₂, p₂), ⋯ , (D_K, p_K)} for training. p_k is the parametric bias for the trial k, a variable that has a common value during that one trial and a different value for different trials. We use this data D_train to train the DPMPB. In the usual training, only the network weight W is updated, but in this study, W and p_k are updated simultaneously as below,

where L is the loss function (mean squared error, in this study) and β is the learning rate. p_k will be embedded with the difference in dynamics in each trial, i.e., the dynamics of the manipulated cloth in this study. p_k is trained with an initial value of 0, and it is not necessary to directly provide parameters related to the dynamics of the cloth. Therefore, parametric bias can be applied to cases where the dynamics parameters are not known, as in the handling of actual cloth.

When manipulating a new cloth, if we do not know its correct dynamics, we will not be able to perform the manipulation correctly. Therefore, we need to obtain data on how the shape of the cloth changes when the cloth is manipulated, and estimate the dynamics of the cloth, i.e., p in this study, based on this data. First, we obtain the data D_new (D_train for one trial with K = 1) as in Section 2.3 when the cloth is manipulated by random commands, GUI, or the controller described later (Section 2.5). Using this data, we update only p while W is fixed, unlike the training process of Section 2.3. We do not execute Equation (3) but update only p by Equation (4) using the current p as the initial value. Since p is of a very low dimension compared to W, the overfitting problem can be prevented. In other words, by updating only the terms related to the cloth dynamics, we make the dynamics of DPMPB consistent with D_new. In this case, the update rule is MomentumSGD.

We describe the dynamic cloth manipulation control using DPMPB. First, we obtain the target image I^ref, which is the target state of the cloth, and compress it to z^ref by AutoEncoder. f^ref and l^ref are set to 0, and s^ref is generated by combining them. Next, we determine the number of expansions of DPMPB Nseqcontrol and set the initial value useqinit as useqopt, which is Nseqcontrol steps of u to be optimized (the abbreviation of u[t,t+Nseqcontrol-1]opt). We calculate the loss function as follows, and update useqopt from this value while W and p are fixed.

where h_loss is the loss function (to be explained later), sseqref is a time series of the target cloth states with s^ref arranged in Nseqcontrol steps, sseqpred is s[t+1,t+Nseqcontrol]pred, which is the time series of s predicted when useqopt is applied to the current state s_t, ||·||₂ is L2 norm, and γ is the learning rate. Equation (6) can be computed by the backpropagation method of neural networks, and Equation (7) represents the gradient descent method. In other words, the time series control command is updated to make the predicted time series of s closer to the target value. Note that W is obtained in Section 2.3, and p is obtained in Section 2.4. Here, γ can be constant, but in this study, we update useqopt using multiple values of γ for faster convergence. We use the value obtained by dividing [0, γ_max] into Nbatchcontrol equal logarithmic intervals (γ_max is the maximum value of γ) as γ, update useqopt using each γ, and adopt useqopt with the smallest L among them, repeating the process Nitercontrol times. An appropriate γ is selected for each iteration and Equation (7) is performed. Also, for the initial value useqinit of useqopt, we use the value optimized in the previous step useqprev (the abbreviation of u[t,t+Nseqcontrol-1]prev) by shifting it one step to the left and duplicating the last term, u{t+1,⋯,t+Nseqcontrol-1,t+Nseqcontrol-1}prev. This allows us to achieve faster convergence, taking into account the previous optimization results. The final obtained target value utopt at the current time step of useqopt is sent to the actual robot.

Here, we consider dynamic cloth manipulation tasks such as spreading out a bed sheet or a picnic sheet in the air. Let s^ref be the target value that contains the image z^ref of the unfolded state of the cloth, and we consider trying to achieve the target value. In this case, if the cloth does not spread well at the first attempt, it would be not possible to make minor adjustments, and it is necessary to generate a series of motions to try to spread the cloth from the beginning. In other words, even if the loss is lowered to some extent by the first attempt, it is necessary to go through a state where the loss is increased by large motions in order to lower the loss further. Therefore, if we just define sseqref as a vector of the same s^ref, we end up with useqopt that hardly moves after the first attempt. The same problem occurs when s^ref is defined as the state of spreading a cloth in the air since that state is not always realized. Thus, it is necessary to realize the target state by periodically spreading the cloth over and over again. Based on this feature, in this study, we determine the number of periodic steps of the periodic motion, Nperiodiccontrol, and change h_loss for each period. In this study, h_loss is expressed as follows.

where {z,f}seq{ref,pred} is the value of {z, f} extracted from sseq{ref,pred}, and w_loss is the weight constant. m_t (∈{0,1}Nseqcontrol) is a vector with 1 occurring at every Nperiodiccontrol step and 0, otherwise. It is shifted to the left at each time step, and 0 or 1 is inserted from the right depending on Nperiodiccontrol (e.g., if Nseqcontrol=Nperiodiccontrol=4, then (0, 1, 0, 0) → (1, 0, 0, 0) → (0, 0, 0, 1) in this order). This makes it possible to bring the cloth state z closer to the target value every Nperiodiccontrol step and enables dynamic cloth manipulation that handles periodic states which can only be realized momentarily. The second term on the right-hand side of Equation (8) is a term to minimize muscle tension as much as possible.

In this study, we set Nseqcontrol=8, Nbatchcontrol=30, γ_max = 1.0, Nitercontrol=3, w_loss = 0.001, and Nperiodiccontrol=8.

The robot and the types of cloth used in this study are shown in Figure 4. First, we conduct a cloth manipulation experiment using a simple 2-DOF robot in simulation with Mujoco (Todorov et al., 2012). Here, the cloth is modeled as a collection of 3 ×6 point masses, and the damping C_damp between those point masses and the weight C_mass of the entire cloth can be changed. Next, we use the actual robot Musashi-W, which is a musculoskeletal dual arm robot MusashiDarm (Kawaharazuka et al., 2019a) with a mechanum wheeled base and a z-axis slider. The head is equipped with Astra S camera (Orbbec 3D Technology International, Inc.). In this study, we use two kinds of cloth: soft type (5 mm thick) and hard type (2 mm thick) of polyethylene foam (TAKASHIMA color sheet, WAKI Factory, Inc.). Polyethylene foam was used instead of actual cloth in this study due to the joint speed limitation of Musashi-W. For the soft type, we prepare one or two sheets, and for the hard type, we prepare one, two, or three sheets (denoted as soft-1, soft-2, hard-1, hard-2, hard-3).

We move the arms of the simulated robot and Musashi-W only in the sagittal plane, with two degrees of freedom in the shoulder and elbow pitch axes. Although the control command of the simulated robot is two dimensional without hardware body stiffness, the control command of Musashi-W is three dimensional including the stiffness (hardware stiffness was difficult to reproduce in simulation). The image is compressed by AutoEncoder after color extraction, binarization, closing opening, and resizing (the front part of the cloth is red and the back part is a different color such as blue or black).

The cloth characteristics are changed to C_damp = {0.03, 0.05, 0.07} and C_mass = {0.05, 0.10, 0.15} while randomly commanding θ^ref (Random) for about 50 s, and while human operation by GUI (joint angle is given by GUI) for about 50 s. Using a total of 4,500 data points, we trained DPMPB with the number of LSTM expansions set to 30, the number of batches being 300, and the number of epochs being 300 (these parameters are empirically set). Principle Component Analysis (PCA) is applied to the parametric bias for each cloth obtained in the training, and its arrangement in the two-dimensional plane is shown in Figure 5. Note that even if p is two-dimensional, the axes of the principal components, PC1 and PC2, become more defined by applying PCA. It can be seen that the PBs are neatly arranged mainly on the axis of PC1 according to the magnitude of C_damp. On the other hand, for the axis of PC2, some PBs are not arranged according to the magnitude of C_mass, and we can see that the dynamics change due to C_mass is complex. Also, the contribution ratio in PCA is 0.65 for PC1 and 0.35 for PC2, indicating that the influence of C_mass on the dynamics of the cloth is small.

Here, with each cloth of (C_damp, C_mass) = {(0.05, 0.05), (0.07, 0.10), (0.03, 0.15)}, the Random motion and the online material estimation in Section 2.4 are executed for about 40 s. The trajectory of the parametric bias (traj) is shown in Figure 5. Note that the initial value of PB is 0 (the origin of PB is not necessarily at the origin of the figure, since PCA is applied). It can be seen that the current PB gradually approaches the respective PB values obtained during training in the dynamics of the current cloth. In particular, the accuracy of the online estimation is high for the axis of C_damp. On the other hand, for the axis of C_mass, though the material can be correctly recognized to some extent, the accuracy is less than for the axis of C_damp.

Next, we performed the control in Section 2.5 by setting the cloth characteristics to (C_damp, C_mass) = {(0.03, 0.05), (0.07, 0.15)} and by setting the current PB value to the PB obtained during the training of the object to be manipulated. Two images, Target-1 and Target-2 of Figure 6, are given as the target images of the cloth. We denote the case of control in this study as Control and the random trial as Random. For each of the two types of cloth, the Control and Random trials are performed for 50 s each, using the two target images. The percentage (y-axis) of ||zref-z||2 lower than a certain threshold (x-axis) is shown in Figure 6. Figure 6 indicates that the more the graph expands in the upper-left direction (the larger the y-axis value becomes when the x-axis value is low), the more the target image is realized. Here, in order to verify the reproducibility, all experiments including learning and control are conducted five times, and the mean and variance of these trials are shown in the graphs. Note that the mean of ||zref-z||2 in the initial condition shown in the left figure of Figure 4 is 1.55 for Target-1 and 1.36 for Target-2. In all cases, Control outperforms Random and succeeds in realizing the target image accurately. Since Target-1 is more difficult than Target-2, especially when (C_damp, C_mass) = (0.03, 0.05), the difference is clearly shown. Since the larger C_damp causes slower deformation of the cloth and makes it easier to lift up, the target image is relatively correctly realized when is set (C_damp, C_mass) to (0.07, 0.15), compared to when it is set to (0.03, 0.05). The reproducibility of performance with respect to learning and control is also high, especially in areas with a low threshold (x-axis), where the variance is small.

Next, assuming that the current cloth properties are (C_damp, C_mass) = (0.03, 0.10) and the current PB values are (C_damp, C_mass) = (0.07, 0.15) obtained during training, we simultaneously perform the online material estimation in Section 2.4 and control in Section 2.5. Note that the target image is set as Target-1 of Figure 6. The transition of PB (integrated-traj) is shown in Figure 5, and the transition of ||zref-z||2 is shown in Figure 7. It can be seen that the current PB gradually approaches the value obtained during training with (C_damp, C_mass) = (0.03, 0.10). Also, as PB becomes more accurate, ||zref-z||2 periodically shows smaller values. In other words, the online material estimation and the learning control can be performed together, and the control becomes more accurate as PB becomes more accurate.

Finally, we show that ||zref-z||2 is close to the actual distance between the raw images since the control and its evaluation are basically performed with the error of the latent variable z in this study. We use Symmetric Chamfer Distance (Borgefors, 1988), which represents the similarity between binary images, to calculate the distance between the current raw image I and its target value I^ref. The relationship between the logarithm of Chamfer Distance and ||zref-z||2 when the target image is set to Target-2 and the cloth is randomly moved for 90 s is shown in Figure 8. The two values are well correlated, with a correlation coefficient of 0.811.

Motion data was collected by randomly commanding θ^ref and k^ref (Random) for about 80 s, and while human operation with GUI (joint angle is given by GUI, joint stiffness is randomly given) for about 80 s, for the cloths of soft-1, soft-2, hard-1, and hard-3, respectively. The random command k^ref is set between [10, 70] [N]. Using a total of 3,000 data points, we trained DPMPB with the number of LSTM expansions being 20, the number of batches being 1,000, and the number of epochs being 600 (these parameters are empirically set). Principle Component Analysis (PCA) is applied to the parametric bias for each cloth obtained in the training, and its arrangement in the two-dimensional plane is shown in Figure 9. It can be seen that soft-1, soft-2, hard-1, and hard-3 are placed at the four corners, along with the thickness and stiffness of the cloth.

We execute the control in Section 2.5 and the online material estimation in Section 2.4 for about 30 s with each of the hard-1, hard-2, and hard-3 cloths. The trajectory of parametric bias is also shown in Figure 9. Note that the initial value of PB is 0 (the origin of PB is not necessarily at the origin of the figure, since PCA is applied). The updated PBs aligned in a straight line according to the thickness of the cloth. The updated PB of hard-3 is close to the PB of hard-3 obtained during training, while the updated PB of hard-1 deviated slightly from the PB of hard-1 obtained during training and was close to the PB of hard-3. The updated PB of hard-2, which is the data not used in the training, is located midway between hard-1 and hard-3. Therefore, we consider that the space of PB is self-organized according to the cloth thickness and cloth stiffness.

First, we performed the control in Section 2.5 with the value of PB set to the PB obtained during the training of the object to be manipulated. From the initial state of the cloth, as shown in the left figure of Figure 10, we use Section 2.5 to bring the cloth closer to the state of being spread out in the air as shown in the right figure of Figure 10. For hard-1 and soft-1, examples of transitions of the target joint angle θ^ref, the measured joint angle θ, and the muscle tension f^const representing the target body stiffness are shown in Figure 11. Note that θ_s−p and θ_e−p represent the pitch joint angles of the shoulder and elbow, respectively. The area surrounded by the red frame is the main part of the spreading motion, where the upraised shoulder and elbow are lowered down significantly (the larger the joint angle is, the lower the arm is). Here, if we look at f^const, we can see that in many cases, the stiffness is increased in the first half and decreased in the second half. This increases the hand speed and causes the cloth to spread out and float in midair, realizing the state shown in the right figure of Figure 10. The absolute values of the hand velocities when Equation (7) is not performed for k^ref and f^const is fixed (= 10 N) (w/o stiffness), and when the control in Section 2.5 is performed (w/ stiffness), are shown in Figure 12. The figure shows the mean and variance of the maximum hand velocities during five 20-s experiments from the state shown in the left figure of Figure 10. It can be seen that the speed is about 12% faster when the change in the body stiffness is used. Also, looking at the difference between soft-1 and hard-1 in Figure 11, in hard-1, θ_e−p and θ_s−p move almost simultaneously and the process from raising to lowering the hand is quick. On the other hand, in soft-1, θ_s−p starts to move slightly later than θ_e−p, and the process from raising to lowering the hand is often slow. The reason is that the soft cloth can maintain flight time even if the hand is lowered slowly at the end, while the hard cloth cannot maintain flight time unless the hand is lowered immediately. In this way, we can see that the way of manipulating the object varies depending on the difference in PB.

Next, we performed the cloth manipulation under several conditions and compared them. The trials with correct PBs (the PB of soft-1 is used for the manipulation of soft-1 and the PB of hard-1 is used for the manipulation of hard-1) are called Correct, the trials with wrong PBs (the PB of hard-1 is used for the manipulation of soft-1 and the PB of soft-1 is used for the manipulation of hard-1) are called Wrong, and random trials such as in Section 3.3 are called Random. The case where the PB is correct but Equation (7) is not performed for k^ref and f^const is fixed (10N) is called Correct w/o stiffness. Under these four conditions, the hard-1 and soft-1 cloths are manipulated for 25 s five times from the state shown in the left graph of Figure 10. The average of the percentage (y-axis) with ||zref-z||2 lower than a certain threshold (x-axis) is shown in Figure 13. The means and variances of the proportions of ||zref-z||2<0.5 taken from Figure 13 are shown in the upper graphs of Figure 14, the means and variances of the minimum value of ||zref-z||2 are shown in the middle graphs of Figure 14, and the means and variances of the maximum value of ||zref-z||2 are shown in the lower graphs of Figure 14. Figure 13 indicates that the more the graph expands in the upper-left direction (the larger the y-axis value becomes when the x-axis value is low), the more the target image is realized. Note that the mean of ||zref-z||2 for the initial condition shown in the left figure of Figure 10 is 0.86. For both hard-1 and soft-1, Correct is the best, and Random is the worst. Additionally, looking at the upper graphs of Figure 14, we can see the characteristics quantitatively. Although hard-1 is more difficult and, thus, has a lower realization rate in general, the results are consistent with the trend that Wrong and Correct w/o stiffness have lower realization rates than Correct. Looking at the middle graphs of Figure 14, the tendency of the minimum value of ||zref-z||2 is also consistent with the upper graphs of Figure 14, and the lower the minimum value is, the higher the percentage of realization is. On the other hand, for the lower graphs of Figure 14, Correct is the highest, and Random is the lowest. This may indicate that the more accurate the method can realize the target image, the more likely it is to produce a state that is different from the target image since the state must be changed significantly in order to realize the target image.

We performed a series of motion experiments incorporating dynamic cloth manipulation. Musashi-W picked up hard-1 and laid it on the table like a table cloth by the proposed method, and placed a basket of sweets on it. The scene is shown in Figure 15. First, the robot recognizes the point cloud of the cloth and grasps the cloth between its thumb and index by visual feedback. Next, the robot spreads the cloth on the desk by dynamic cloth manipulation, which has been previously trained on a similar setup. Finally, the robot recognizes the point cloud of the basket and grasps the basket between both hands by visual feedback, and successfully places it on the table.