Thank you for your insightful and constructive comments! We discuss some of your questions and concerns below.

1. Rigid body dynamics & action space

While we can directly integrate the rigid body dynamics into the implicit-time integration scheme, the system will become quite cumbersome and therefore we decide to process the deformable and rigid system separately. Currently the robot manipulator is modelled by position-based control, i.e., the action space of rigid body end-effectors are defined as SE(3) movements within each timestep, and therefore the elastomers are attached to moving rigid boundaries controlled by the actions.

We are still working on further modelling of rigid body dynamics. For the rigid manipulated objects, we plan to support them by adding an stiff orthogonality potential energy on a deformable mesh as in ABD[1]. We also plan to include the penalty energies of joint constraints to model the dynamics behavior of articulated body system.

2. Constrainted simulation

We support the constraints of fixing vertex positions. In Folding and Forming, we fix one end of the paper on the table. In Separating and Following, we fix two manipulator in a gripper shape.

The technique to simulate the constraints is simple. To fix one vertex (or even one dimension of this vertex), we simply set the derivative and hessian matrix related to this dimension to zero before solving the equation. This will result in an unchanged position for this vertex while modeling it's force and influence to other vertices.

3. Real-world manipulation

We add more real-world manipulation cases identical to Folding tasks in simulation. We use a rigid manipulator to interact with a curved paper and try to fold on different creases.

Questions

Q1: The contact model uses quadratic energy instead of barrier. How does the system correctly identify that 2 triangles are penetrating? Based on the paper, it seems it is based on normal direction, however, this does not seem correct since the normal would flip if one of the contacting shell is flipped.

To model contact behavior, we apply quadratic penelty energies on a dynamic set of contact constraints. The reason of not using a barrier function is that we observe that the quadratic model converges much faster than the log-barrier model on complex scenes involving multiple layers of thin-shell material. With quadratic energies, we actually make non-penetration a soft constraint with a large stiffness coefficient. When a constraint is established between a vertex and a shell surface, we maintain a "side" property between them, whose value is one of the shell's two sides. This property remains unchanged until the vertex's projection falls on the shell's boundary. When dealing with contact between the vertex and the shell, all faces of the shell will be assigned with a normal direction pointing towards the corresponding side. In this way, we avoid the ambiguity problem of shell surface's normal direction, and also the penetration problem caused by the soft constraint, since a vertex slightly penetrating into a surface will be repulsed towards a consistent direction.

Q2: The paper describes the system as simulating tetrahedral and triangular elements, how is the rigid-body action space achieved?

It is addressed above in point 1.

Q3: Newton's method requires solving a large linear system. How is this implement on the GPU? Is it solved through sparse matrix algorithms, conjugate gradient methods, or some other method? What is the speed/performance of this simulator? What are the bottlenecks?

For linear equation solving, we use a python binding of the cuSPARSE library, which supports solving sparse matrices using Cholesky Factorization on GPU. The speed is around 2.74 FPS in one case and the bottlenecks are 1) linear equation solving at each iteration of Newton's method and 2) Hessian matrices building and SPD projecting. The details could be found in general response point 4.

[1] Lan, Lei, et al. "Affine body dynamics: Fast, stable & intersection-free simulation of stiff materials." arXiv preprint arXiv:2201.10022 (2022).

We sincerely appreciate your comments. Please feel free to let us know if you have further questions.

Best, Authors

We appreciate your positive and insightful comments! Below, we address your concerns in detail.

1. Labeling inverse tasks

We have separate those tasks in tables and figures in our main paper. Thank you for pointing them out!

2. Conclusion and future directions

We draw a conclusion that RL methods are of low data efficiency in most of the times and it would be good to seek a combination of sample-based methods (maybe RL) and gradient based method. More discussions about this are shown in General Response point 7.

Questions

Q1: Interested to know what is the speed of the simulation? Is multiple GPU parallelization supported by the simulation?

The speed is approximately 2.74 FPS in a test case. Currently, there are no obstacles preventing us from running parallel environments on different GPUs. However, we haven't pursued this option because the GD method does not derive any benefits from parallelization. More details can be found in General Response point 4.

We sincerely appreciate your comments. Please feel free to let us know if you have further questions.

Best, Authors

We appreciate the positive and insightful comments from you! We adress your concerns in details below.

1. Typo and figure quality

Thank you for pointing them out, we have made the revisions.

2. Real-to-sim disparity

While we are using tactile sensor and force feedbacks as observations in the real-world, we compute the disparity by the L2 loss between the simulated tactile markers (projected into 2D images given the camera parameters) and the images collected in real-world, and by the L2 loss between the simulated forces and the real-world force feedbacks. The loss of several experiments and 2D images of markers are ploted in the appendix.

3. Hybrid method details

This information has been incorporated into General Responses point 5. In essence, we apply a fixed proportion to the steps of CMA-ES and GD in the hybrid method. Additionally, we introduce specific constraints to shape the behavior and enhance the learning process of CMA-ES.

Questions

Q1: How well does the simulation environment run in practice, is it achieving something close to real-time, or how long does simulation take?

We test the speed under a specific setting, and it's 2.74 FPS for forward simulation. Details are included in General Response point 4.

Q2: When parallelizing the simulations on GPU, does that only benefit CMA-ES and RL, or is GD also gaining some advantage from the parallelized environments?

The feature of massive parallelism provided by Taichi is designed to parallelize the computing process within a single environment, resulting in faster speeds. Consequently, all algorithms benefit from this improvement.

It's worth noting that it's also possible to run multiple environments on different GPUs, which provides advantages for CMA-ES and RL (while GD does not gain advantages from this feature). In our testing, we ensure a fair comparison by using a single GPU to run a single experiment.

We additional explanations have convinced you of the merits of our work. Please do not hesitate to contact us if you have other concerns.

We appreciate your time! Thank you so much!

Best, Authors

Thank you for the extensive replies to all the reviewers, and for grouping the concerns into specific categories. Excellent clarification of the solver that was implemented. The detailed runtimes now make the usability very clear for practical users of the framework. Thank you to the authors for the extra time on the rebuttal!

We appreciate the positive and constructive comments from you! We have modified our paper according to your comments.

Success Indicators While we design continuous rewards without a distinct boundary between success and failure and primarily assess the success of a policy through visualization, we also present scores for human-designed policies, theoretical analyses, and success indicators in General Responses point 6.

Questions

Q1: In Sec 4.1, it seems that the observation includes a set of points. Does the order of points matter? How do the learning algorithms (SAC, PPO) used in the paper handle the observation if the input observations (points) are unordered?

The order of points matters in the observation. As we uniformly sample points in the full state, the sequence of selected nodal points remains constant within a given task.

Q2: In Sec 4.1, does the 6DoF stand for the pose change in SE(3), and does the 7DoF include one more DoF for gripper closeness (mimic behaviors for two finger tips)?

Yes, the 6DoF stand for the pose change in SE(3) (x, y, z, roll, pitch, yaw), and the 7DoF include one more DoF for gripper closeness.

Q3: In Sec 5.1, the authors mentioned that they reported "the maximum scores achieved within a specified number of timesteps". How do the authors count "timesteps" for CMA-ES and GD methods? What's the meaning of "episode" (x-axis) for GD and CMA-ES in Figure 3?

The "episode" means a rollout which is a full sequence of actions from the initial state to a final state.

The "timesteps" here has the same meaning of the previous "episode". As a result, we can simply count the number of full action sequences executed in GD, CMA-ES, and RL to count the episode. We will change them into a uniform one to avoid confusing. Thank you for pointing it out.

Q4: Does it mean that the benchmark currently focuses on "optimization" rather than "generalization"?

The experiment setting we currently test on indeed focuses on "optimization" rather than "generalization".

However, we still want to share vision about how to address generalization in General Response point 7. Concretely, we can switch the initial state and parameters in one task and evaluate close-loop policies. Another way is to use open-loop method to collect successful trajectories under different initial state, parameters or even different tasks, then we can use imitation learning with those data to obtain a generalist policy.

Q5: In Sec 5.2, for "our gradient-based method converges these parameters towards a uniform outcome.", do the authors mean that the method converges to similar solutions despite different initialization?

Yes. While the Card still have a respectively large standard deviation in the end, it is clearly shown in Sliding that they converges to similar solutions.

Please let us know if you have any further questions for our paper. We really appreciate your time! Thank you!

Best, Authors

7.Hybrid method In this work, we focus more on how to utilize the gradient in simulation to train an open-loop policy within limited time. While RL policy can generalize better to unseen state, open-loop methods employing differentiable physics prove to be more data-efficient. CMA-ES, GD, and our hybrid method all operate as open-loop trajectory optimization approaches. The concept of combining CMA-ES and GD with differentiable physics is proposed to alleviate the issue of local minima associated with the GD method in scenarios characterized by rich contact interactions, while design of RL policy considers no physics properties. More details of this hybrid method is shown in the general response.

Questions

Q1: The paper states that the goal is for this to be used as a benchmark. Is the idea to use this similar to gym and atari game setups, or more like the benchmarks used in the computer vision community with tasks on withheld data to counteract overfitting on a set of tasks?

It is similar to gym and atari game setups, where users can generate their own data and build their own experiment.

Q2: What solution quality / success rate do the various rewards correspond to?

As outlined in General Response point 6, we have established a set of scores to signify success in our benchmark. It is crucial to note, however, that we consider these tasks to be continuous rather than discrete. Our primary method of assessing success revolves around visualizations, and the success indicators are primarily employed for comparative purposes.

Q3: How fast does the simulator run, is it real-time, faster, slower?

It's about 2.74 FPS as discussed in the General Response point 4.

Q4: What is the limit of object complexity that can be simulated?

We employ 3D meshes to simulate objects, enabling the simulation of any object that can be represented using a 3D mesh. It's important to note, however, that simulations may experience decreased speed and increased memory requirements when dealing with objects with a high number of vertices.

We wish that our response has addressed your concerns, and turns your assessment to the positive side. If you have any more questions, please feel free to let us know during the rebuttal window.

Best, Authors

Thank you for your insightful and constructive comments! We have added additional experiments and modified our paper according to your comments.

1.Simulation details

This is addressed in General Response point 3. The detailed equations are too extensive to include in the main paper.

2.What is the level of success for each task

As demonstrated in the General Response point 6, we integrate the theoretically highest scores with human-intuitive policies to derive the success score for each task. The success rate is then computed based on the mean and standard deviation for the CMA-ES + GD method.

3.Actual definition of task objective Although we have detailed the reward function in Appendix B.2, we provide the exact formulas of the reward functions for further clarification.

Lifting: Let $B$ be the set of vertex positions on the block and $B^T$ be the target vertex positions, and $B[i]$ be the $i_{th}$ vertex in $B$. $reward=-\sum_{i}\{(B[i].x - B^T[i].x)^2 + (B[i].y - B^T[i].y)^2 + (B[i].z - B^T[i].z)^2\}$

Separating: Let $B$ be the set of vertex positions on the block and $C$ be the set of vertex positions on the paper. Since there are $256$ vertices in the shell and $144$ vertices in the block, we multiply the sume of block vertices with $256 / 144$. $reward=256 / 144 * \sum_{i}B[i].x - \sum_{i}C[i].x$

Following: $reward=-\sum_{i}B[i].x$

Folding-U: Let U be the set of rest angle (angle at rest position) on the upper crease and L be the set of rest angle on the lower crease. $reward=\sum_{i}U[i] - \sum_{i}L[i]$

Folding-L: $reward=-\sum_{i}U[i] + \sum_{i}L[i]$

Pick-Folding: Let M be the set of rest angle (angle at rest position) on the middle of the paper. $reward=\sum_{i}M[i]$

Forming: Let $S$ be the set of vertex positions on the shell and $S^T$ is the target vertex positions, and $S[i]$ be the $i_{th}$ vertex in $S$. $reward=-\sum_{i}\{(S[i].x - S^T[i].x)^2 + (S[i].y - S^T[i].y)^2 + (S[i].z - S^T[i].z)^2\}$

Sliding: Let $S$ be the set of vertex positions on the bottom shell. $reward=-\sum_{i}S[i].x$

Bouncing: Let $S$ be the set of vertex positions on the shell. $reward=\sum_{i}S[i].z$

Card: Let $C$ be the set of vertex positions on the card. $reward=-\sum_{i}C[i].x$

4.The presentation of the results

The symbol $\pm$ in the table denotes the standard deviation resulting from multiple runs. Notably, the GD method exhibits no standard deviation as it lacks randomness in its multiple executions. In the context of inverse design tasks, the optimization is solely focused on the parameters. Consequently, we evaluate a distinct set of methods for these tasks, which explains the absence of certain rows in the table.

We will ensure to articulate these details more explicitly in the main paper. Thank you for bringing this to our attention.

5.Conveying the goal

Our goal of the experiment section is to delve into task behavior to highlight the inherent challenges, and analyze the strengths and limitations of various methods in addressing these challenges. This analysis serves as the groundwork for the future development of innovative algorithms for thin-shell manipulation tasks. We've add this to our main paper, thank you for metioning it.

6.Real world experiment

• [Real-to-sim system identification] We have introduced more detailed experimental results for the real-to-sim system identification tasks [Appendix C]. These results are analyzed through the observation of the 2D marker pose differences, presented qualitatively, and the loss curve, presented quantitatively. Although tactile simulation is not the primary focus of this paper, we argue that it demonstrates satisfactory alignment with the real world. We have encountered challenges in accurately retrieving small-scale bending stiffness from the real world using tactile sensors. This difficulty arises from the sensitivity of results to the Poisson ratio of the tactile material. We deem that a better way to estimate small-scale bending stiffness is to first use vision method to obtain it's state and subsequently estimating it using the balance between gravitational and bending forces.

• [Sim-to-real manipulation tasks] In the real-world execution of the Separating task, our paper discusses that the bending stiffness of the thin shell significantly influences the final policy. Successfully deploying this policy to the real world serves as evidence of our effective approximation of the bending stiffness of the card paper. We have also incorporated additional real-world manipulation results. In the Folding task, we employ a rigid manipulator with 3 Degrees of Freedom (DoFs) and train a new policy in the simulator, which is subsequently deployed into the real world. The outcome illustrates our success in capturing the dynamics of the paper and accomplishing the task in the real-world setting.

4. Back-Propagation

This part has been briefly discussed in Appendix A.2. The equation to be solved in forward simulation is

where $\mathbf x^t$ and $\mathbf x^{t-1}$ are node positions of $t_{th}$ and $(t-1)_{th}$ timestep,

$\overline{\mathbf \theta}^t$ are undeformed angles for $t_{th}$ timestep, $\mathbf x_{fixed}^t$ are positions of fixed nodes in the beginning of $(t+1)_{th}$ timestep,

and $\mathbf\eta$ is the material parameters. Here, $\overline{\mathbf \theta}^t$ is actually a function of $\mathbf x^{t}$ and $\overline{\mathbf \theta}^{t - 1}$, and the friction force is a function of $\mathbf x_{i}$ and $\mathbf x_{i - 1}$. DiffPD [6] discusses the general idea of how to pass gradient in the implicit time integration:

With $\frac{\partial \mathbf x}{\partial \mathbf y}$ in hand, we can backpropagate the gradient from $x^{t+1}$ to $y^t$. Given $y^t=2x^t-x^{t-1}+h^2M^{-1}f_{ext}$, we can deliver the gradient to both $x^t$ and $x^{t-1}$.

Similarly, by substituting $y$ in the above equation with other variables like $\mathbf x_{fixed}^t$ or $\overline{\mathbf \theta}^t$, we can calculate the gradient. The remaining challenge is computing the term $\frac{\partial \nabla \mathbf g(\mathbf x)}{\partial \mathbf y}$ in the first line when changing $y$ into other variables. For these variables, this term will involve the Jacobian matrix for the force, as shown in Equation (1).

For $\overline{\mathbf \theta}^t$ and $\mathbf\eta$, their relationship to the forces is linear, allowing us to calculate the Jacobians directly. Regarding $\mathbf x_{fixed}^t$, the Jacobians are computed during forward simulation, where they are set to zeros. We can apply the same procedure in forward simulation and retain the Jacobians during the backpropagation process.

Concerning $\mathbf x^t$, it primarily influences the friction force since we utilize the pressure force, contact surface, and the original projected point from the last timestep. We observe that [7] calculates the derivative for friction similarly to our approach. The distinction lies in our choice not to pass the gradient through the surface change during backpropagation

Again, thank you for your insightful and constructive comments!

[1] Li, Minchen, et al. "Incremental potential contact: intersection-and inversion-free, large-deformation dynamics." ACM Trans. Graph. 39.4 (2020): 49.

[2] Smith, Breannan, Fernando De Goes, and Theodore Kim. "Stable neo-hookean flesh simulation." ACM Transactions on Graphics (TOG) 37.2 (2018): 1-15.

[3] Eftychios Sifakis and Jernej Barbic. Fem simulation of 3D deformable solids: A practitioner’s guide to theory, discretization and model reduction. ACM SIGGRAPH 2012 Courses, 2012. 4

[4] Grinspun, Eitan, et al. "Discrete shells." Proceedings of the 2003 ACM SIGGRAPH/Eurographics symposium on Computer animation. 2003.

[5] Tamstorf, Rasmus, and Eitan Grinspun. "Discrete bending forces and their Jacobians." Graphical models 75.6 (2013): 362-370.

[6] Du, Tao, et al. "Diffpd: Differentiable projective dynamics." ACM Transactions on Graphics (TOG) 41.2 (2021): 1-21.

[7] Huang, Zizhou, et al. "Differentiable solver for time-dependent deformation problems with contact." arXiv preprint arXiv:2205.13643 (2022).

3. Thin-Shell Energy

Following [4], we apply the stretching energy $U_s$ and the bending energy $U_b$ for thin-shell materials. We further divide the stretching energy into $U_s^e$ and $U_s^a$, representing stretching energy for edges and areas of triangles. The energy formula is as follows:

where $e$ represents all the edges, $A$ represents all the triangles, and $\theta_e$ represents the angle on edge $e$. $|\overline{e}|$, $|\overline{A}|$, $|\overline{h_e}|$, and $\overline{\theta_e}$ respectively represent the rest length of the edge, the rest area for the triangle, one-third of the average height of the two adjacent triangles of the edge, and the rest angle on this edge. $K_e$, $K_a$, and $K_b$ are the coefficients.

For $U_s^e$, it's simple to calculate the derivatives and Hessian matrix. The derivative of the length $|e|$ to the nodal point positions $x_0$ and $x_1$ is $\frac{\partial |e|}{\partial x_0}=(x_0 - x_1)/|e|$. The Hessian matrix will require the calculation of $\frac{\partial^2 |e|}{\partial a_0^2}$ and $\frac{\partial^2 |e|}{\partial a_0 \partial b_0}$ due to symmetry, where $x_0=(a_0,b_0,c_0)$.

For $U_s^a$, it becomes more complex to compute the derivatives and Hessian matrix. Assume the triangle is $(x_0, x_1, x_2)$ and $x_0=(a_0, b_0, c_0)$. Similar to the calculation of the contact energy, we use math tools to generate the symbolic formula of the derivative $\frac{\partial |A|}{\partial a_0}$ and hard code it with Taichi; then, we can calculate all the derivatives due to symmetry. For the Hessian matrix, we use math tools to generate formulas for $\frac{\partial^2 |A|}{\partial a_0^2}$, $\frac{\partial^2 |A|}{\partial a_0 \partial a_1}$, $\frac{\partial^2 |A|}{\partial a_0 \partial b_0}$, $\frac{\partial^2 |A|}{\partial a_0 \partial b_1}$; then, we can calculate the full matrix using symmetry.

For $U_b$, the Hessian matrix is too complex to compute, even when we try to use some math tools. As a result, we follow [5] to compute the derivatives and Hessian matrix for bending energy. This paper discusses how to compute the Hessian matrix of bending energy in a simplified form.

We provide additional details about the implementation of our simulation, specifically addressing the process and methods we used to compute the derivatives and Hessian matrix to address your concerns.

1. Contact Energy According to Section 3.2, the contact energies $U_r$ and $U_f$ are defined as

where the distance $d$ is a signed distance between vertex $\mathbf p$ and its projection to a triangle $(\mathbf p_1,\mathbf p_2,\mathbf p_3)$, whose normal direction is represented by $(\mathbf p_2-\mathbf p_1)\times(\mathbf p_3-\mathbf p_1)$. $k_r$ stands for a constant stiffness parameter, $\mu$ stands for a friction coefficient, and $\lambda$ stands for the contact force along normal direction at last timestep. The signed distance $d$ is then calculated as

in which the 3D vertices are represented as column vectors. Let $\mathbf x=[\mathbf p_1,\mathbf p_2,\mathbf p_3,\mathbf p]\in\mathbb R^{12}$ denote the involved variables, we implemented tool functions in Taichi to calculate

together with the gradients $\partial D/\partial\mathbf x,\partial C/\partial\mathbf x$ and Heissan matrices $\partial^2 D/\partial \mathbf x^2,\partial^2 C/\partial\mathbf x^2$. With the help of SymPy library for symbolic computation, we optimized the computation graph and hard coded it into Taichi snippets. The repulsive energy $U_r$ is then a composition function on $D(\mathbf x),C(\mathbf x)$, and therefore we use the chain rule to calculate its gradient and Hessian matrix.

For the friction energy $U_f(\mathbf u)$, we follow IPC [1] to calculate the $f_0(\epsilon_v,\|\mathbf u\|)$ function, its derivative $f_1(\epsilon_v,\|\mathbf u\|)$ and second derivative $f_2(\epsilon_v,\|\mathbf u\|)$. Once we get the $\partial U_f/\partial \|\mathbf u\|$ and $\partial^2 U_f/\partial \|\mathbf u\|^2$, we transform it back to the tangential space to get $\partial U_f/\partial\mathbf u$ and $\partial^2U_f/\partial\mathbf u^2$.

2. Elastic Energy

We follow Stable Neo-Hookean [2] to define the elastic energy as

in which $\mathbf F$ denote the deformation gradient of one tetrahedral element, $I_C=\text{tr}(\mathbf F^T\mathbf F)$ and $J=\det \mathbf F$. The motivation of replacing the original $\log J$ term in original Neo-Hookean model to $J-1$ is that $J-1$ is stable even for $J\le0$. Therefore, even if the whole tetrahedral element is inverted, this Stable Neo-Hookean model still works robustly and can provide gradients & Hessian matrices to recover the element (See the original paper of Stable Neo-Hookean [2]). This situation is common in our simulator, since the robot manipulator involves a moving boundary and can easily cause invertion when the moving speed is large. We observe that Stable Neo-Hookean perfectly addresses this problem.

The gradient $\partial U_e/\partial \mathbf F$ can be computed easily. To compute the Hessian matrix $\partial ^2U_e/\partial\mathbf F^2$, we first enumerate the differential $\delta F_{ij}$ for single entry and then calculate the correponding force differential $\delta(\partial U_e/\partial\mathbf F)$. We follow a well-known tutorial on FEM [3] to calculate the force differentials, and then concatenate them into the Hessian matrix $\partial^2 U_e/\partial \mathbf F^2$.

6.Success indicator

Defining the level of success in certain tasks proves challenging due to uncertainty regarding the optimal performance achievable within the current settings using the best policy. Currently, we gauge the success of a policy by visually inspecting the results. Nevertheless, we aim to enhance analytical insight by presenting the scores of human-designed heuristic policies alongside theoretically optimal scores. Subsequently, we can formulate a success indicator based on an analysis of these scores.

• [Human Designed] To elucidate the meaning of reward values, we present the scores of human-designed policies for each task in the table below. Notably, generated policies occasionally outperform human-designed counterparts, particularly in tasks requiring meticulous manipulation and where certain policies exhibit unexpected effectiveness (as observed in Separating).
• [Oracle] In addition to human-designed policies, we also present the theoretically maximum scores for each task, disregarding certain physical constraints. For Pick-Folding, we calculate this score by assuming the entire middle curve is fully folded. In the case of Separating, we assume the paper is affixed to the manipulator, while the block experiences no friction force with the paper. For Following, we assume both the manipulator and the block are affixed to the paper. In both Folding-L and Folding-U, we assume that the target curve is fully folded while the other one is completely flattened. For tasks involving Lifting and Forming, we assume a perfect match with the target. It is important to note that although certain Oracle scores, such as those for Pick-Folding and Separating, are exceptionally high, they are not physically feasible within the current experimental settings. We present this information solely to illustrate the approximate score range for our tasks.
• [Success Score] To make it further concrete, we design the score of success for each tasks. Concretely, we define the score of success as $min(HumanDesign, Oracle/2)$ for Pick-Folding, Separating Folllowing and $Oracle/2$ for Folding-U and Folding-L since human designed policies perform poorly in those two tasks. And the $-3e-4$ as the successful score of Forming and Lifting.
• [Success Rate] Finally, we calculate the success rate of our CMA-ES + GD method according to the mean and standard deviation.

Table A.

7. Conclusion and future directions Through the evaluation of various methods in our benchmark tasks of thin-shell manipulation, we have concluded that reinforcement learning algorithms exhibit limited data efficiency. This inefficiency arises from the sensitivity of thin-shell materials and the heightened complexity of their states. Simultaneously, gradient-based methods face challenges in reliability within contact-rich scenarios.While sample-based methods, such as CMA-ES, generally perform well, they suffer from reduced data efficiency in later episodes and intricate settings. Consequently, we propose a hybrid approach that combines the strengths of CMA-ES and gradient-based methods to achieve superior results.

In terms of future work, we recommend exploring closed-loop control under varied parameters, initializations, and tasks. Currently, our focus has primarily been on open-loop methods like CMA-ES and gradient-based methods. We posit that an avenue worth exploring involves leveraging the synergy between reinforcement learning and gradient methods [1] for improved data efficiency in closed-loop control. Additionally, an alternative strategy could involve collecting successful trajectories under diverse settings using open-loop methods and subsequently employing a neural network for imitation learning. This approach aims to address generalization challenges rather than optimization problems.

We hope our responses below convincingly address all reviewers’ concerns. We thank all reviewers’ time and efforts again!

[1] Mora, Miguel Angel Zamora, et al. "Pods: Policy optimization via differentiable simulation." International Conference on Machine Learning. PMLR, 2021.

4. Parallelization and Running Speed of Simulation

In ThinShellLab, both the forward and backward passes are parallelized in a single GPU. The major bottlenecks of the running time are 1) linear equation solving at each iteration of Newton's method and 2) Hessian matrices building and SPD projecting. For Hessian matrices, we use Taichi programming language to process the submatrices for independent elements in parallel. We also implement the QR iteration in Taichi to project each submatrix onto a manifold of SPD (symmetric positive definite) matrices. For linear equation solving, we use a python binding of the cuSPARSE library, which supports solving sparse matrices using Cholesky Factorization on GPU. The rest parts of simulation (e.g., line search and collision detection) are fast enough and therefore can be negligible. The collision detection is efficiently solved by space hashing using Taichi's sparse data structure.

To benchmark the running speed of simulation, we test our simulator on the Forming task, which contains a typical scene with 1698 DoFs involving both volumetric and thin-shell deformable bodies. The test is done on a single RTX 4090. In average, the forward / backward passes take 36.41 s / 4.67 s respectively for 100 timesteps, i.e., framerates of 2.74 fps / 21.41 fps. In comparison, DiffCloth reports their framerates at 0.73 fps / 2.72 fps using the same Newton-Cholesky method and 7.14 fps / 17.54 fps using the fast PD acceleration, both under a scene with 2469 DoFs. C-IPC reports a framerate at 0.024 fps for forward simulation of a thin-shell scene with 18.3K DoFs.

In practice, it takes about 3 hours to train a policy in one task.

While we run our experiments on a single GPU, it is also possible to run multiple environments on different GPUs in parallel. We didn't do this because gradient-based trajectory optimization method do not benefit from parallel environments.

5. Details and analysis of hybrid method We list the detailed design and analysis of the hybrid method as follows.

• Parameters We opt for a population size of $40$ in CMA-ES with zero-initialized trajectories. The initial variance is set to 1.0, corresponding to 0.0003 in each dimension of manipulator movement per timestep. The maximum movement within one timestep is either 0.001 or 0.002, depending on the task. Generally, we execute CMA-ES for 80% of the total episodes, where episodes are defined by the number of rollouts. For the Pick-Folding task, which demands more extensive training, we set the CMA-ES episode count to 1000 and supplement it with an additional 150 episodes of gradient descent.
• Design To enhance the effectiveness of the CMA-ES method, we filter out non-feasible states, such as instances where there is no contact between manipulators and thin shell or when the manipulator is experiencing excessive force. In these situations, we directly assign a low reward value and save the time for executing the rest actions.
• Analysis Some curves in Figure 3 exhibit large variance. Specifically, in Folding-U, Following, and Separating, the variance arises from the sensitivity of thin-shell manipulation results to minor changes, as previously discussed in the introduction section. In the case of Pick-Folding, although sensitivity contributes to the variance, the primary factor is occasional failure of the CMA-ES method. Failures result in rewards near 0, and successful runs yield rewards greater than 6, leading to substantial variance in the outcomes.