Neural Inverse Source Problem

Authors: youngsun  wi
jayjun  lee
miquel  oller
nima  fazeli
Venue: robot-learning_org_CoRL_2024_Conference
Published: 2024-06-07 04:02:34

Abstract: Reconstructing unknown external source functions is an important perception capability for a large range of robotics domains including manipulation, aerial, and underwater robotics. In this work, we propose a Physics-Informed Neural Network (PINN) based approach for solving the inverse source problems in robotics, jointly identifying unknown source functions and the complete state of a system given partial and noisy observations.  Our approach demonstrates several advantages over prior works (Finite Element Methods (FEM) and data-driven approaches): it offers flexibility in integrating diverse constraints and boundary conditions; eliminates the need for complex discretizations (e.g., meshing); easily accommodates gradients from real measurements; and does not limit performance based on the diversity and quality of training data. We validate our method across three simulation and real-world scenarios involving up to 4th order partial differential equations (PDEs), constraints such as Signorini and Dirichlet, and various regression losses including Chamfer distance and L2 norm.

Main paper: 
 introduction: 
 We are interested in differential equations g with an unknown external source f : g(x, Φ, ∇ x Φ, ∇ 2x Φ, ...) = f (x), Φ : x -→ Φ(x).(EQUATIONOur objective is to determine the external source function f (x) ∈ R q and the full mapping from spatial/temporal coordinates x ∈ R s to the quantity of interest Φ(x) ∈ R r from known governing equations g and partial/noisy observations. This problem, known as an inverse source problem [2, 3, 4] , is of significant importance in various scientific and engineering domains such as signal processing [5] , fluid dynamics [6, 7] , optics [8] , and more.To address these inverse source problems in robotics, a common strategy is to use finite element analysis integrated with optimization techniques [9, 10] or physics engine [11, 12] . While these approaches focus on computation efficiencies, they suffer two notable limitations. First, they rely on complex discretization and/or meshing, which affects precision and introduces significant complexity. Second, these approaches require unique treatment for different equations and constraints, which significantly limits their generality. Recent progress in representation learning addresses these issues via learning priors over the simulated source functions, where the data implicitly contains all the constraints and boundary conditions. For example, prior works in manipulation simulate extrinsic contact, a source function exerting forces on deformable [13, 14, 15, 16] or rigid objects [17, 18] . Another example involves simulating sensor deformations as a source function transducing electric signals [19] . However, the performance of these approaches relies heavily on the diversity and quality of training data and does not guarantee adherence to the physics, especially with out-of-distribution inputs. Latent force models [20, 21] present an alternative approach using Gaussian processes to identify unknowns in differential equations; however, they are restricted to linear models, and it remains unclear how complex constraints could be effectively integrated.In this paper, we propose a method for solving the inverse source problem based on Physics-Informed Neural Networks [1] that enables the simultaneous inference of both the source function and the full state given partial and noisy observations. Our method offers several advantages over prior work in robotics: compatibility with various nth-order differential equations (including ordinary and partial differential equations), flexibility in integrating various constraints, initial values, and boundary conditions, elimination of the need for non-trivial discretizations like meshing, ease of accommodating gradients from real measurements, and does not limit performance based on the diversity training data. We open-source our codes and a real-world datasets at URL. 
related works: 
 Physics Informed Machine Learning: Physics-Informed Neural Networks (PINNs) have shown promising results in solving differential equations, including ODEs, PDEs, and IDEs. PINN [1] was initially limited to solving 1st or 2nd order differential equations [22, 23, 24] . Recently, [22, 25] succeeded in solving the forward problem of 4th-order spatio-temporal differential equations by incorporating residual layers [22, 23] . Our network architecture builds upon [25, 22] 's work for its capability in solving complex higher-order differential equations. In contrast, our goal is to recover the source and unobserved system states from partial observations and handle constraints.Inverse Problems in PINN: Inverse problems in PINNs primarily focus on system parameter or coefficient identification in sparse and noisy data regimes [26, 27, 28, 29, 1, 30, 31, 32, 33] . Parameter identification is fundamentally a different problem from source function identification; it is about the internal properties and inherent dynamics of the system, whereas we are interested in external influences that drive the system, affecting the transient and steady-state behaviors. Recently, a thermonuclear study [30] showed promising results in identifying the forcing function of a simple first-order spatio-temporal differential equation without any constraints. To the best of our knowledge, no existing work in PINN has focused on inverse source problems with constraints, validation on high-order differential equations, and applications to robotics.Inverse Source Problem in Robotics: External source identifications have seen considerable attention in the robotics from both first-principles and learning perspectives. In addition to the work cited in the Introduction, we refer the reader to Appendix. A.1. 
methodology: 
 3.1. Problem Formulation: Our goal is to learn a mapping Ψ : x -→ (Φ(x), f (x)) that satisfies:EQUATIONwhere x is a spatial/temporal coordinates we have partial/noisy access to, Φ(x) is the fully/partialy observable quantity of interest(e.g., deformation field, electric field, etc.), f (x) is a Lipschitz continuous source function, and Ω n , Ω m are subsets of a bounded domain Ω. When we approximate the mapping with a neural network as Ψ θ (x) and f θ (x), this problem becomes solving for the network parameters θ, satisfying the differential equations g(x, Φ,∇ x Φ θ , ∇ 2 x Φ θ , . . .) -f θ (x) = 0 and M constraints {C m (x, Φ θ , ∇ x Φ θ , ∇ 2 x Φ θ , . . . , f θ (x)} M m=1 .We cast this problem into a loss:EQUATIONwhere the first term represents residual losses (L r ) of the differential equation (DE), and the second term represents the constraint loss. The constraint loss includes various forms of regression losses (L reg ) enforcing desired output at specific coordinates (e.g., Chamfer Distance, Mean Square Error), boundary conditions (e.g., Dirichlet, Neumann, Signorini), and equality and inequality constraints.We use a modified multi-layer perceptron (MLP), which has been shown to excel at learning PINNs [22, 23] . The modified MLP's l-th layer has 2 steps: The first step is a fully connected layer with weight W (l) , bias b (l) , and a tanh activation function acting on the previous layer output z (l-1) . The second step weights the output z ′(l) with c and d from separate MLPs with a tanh activation function, as in Fig. 1 . Here, we propose the last fully connected layer unique to each dimension of y and f . Please read Appendix. E.1 and E.4 for details.EQUATIONIncorporating Inductive Bias into Network Design: Many practical robotics systems have low-frequency signals, and building a network with bandwidth matching the underlying signal is essential for robustness to signal noise and successful training. We found that removing the input mappings helps the model be robust to measurement noise (Appendix-E.2), whereas prior works in AI4science [25, 22] highlight the necessity of the input mappings, particularly when focusing on high-frequency signals.Ensuring inputs and outputs have reasonable scales is crucial for successful training [25, 34, 35] . We normalize inputs to be within a unit cube [-1, 1] s and each dimension of the outputs to have similar scales. For instance, if we scale inputs in [m] unit by some value k, we scale Young's modulus [N/m 2 ] by 1 k 2 and depth measurements in [m] unit by k. 3.5. Loss Weighting: We found that an effective training strategy is to quickly regress to the partial measurements and then gradually solve for the forcing function. One way to control the convergence speed of each loss term is by updating the loss weights every hundred epochs using each loss term's gradients [25] ; however, this approach is expensive in time and space, especially with highly nonlinear and high-order residual losses. Instead, we recommend fixed loss weightings satisfying 50λ r ∥∇ θ L r ∥ ≤ λ reg ∥∇ θ L reg ∥, where ∇ θ L reg and ∇ θ L r are the gradients of the regression loss and the residual loss. Please see Appendix-E.3 for details. 
experiments: 
 We present three inverse source problem examples in both simulation and real-world that vary in several key aspects: the order of the differential equations, the degree of observability (ranging from partial to full measurements), and the levels of noise in the data. Additionally, we demonstrate system parameter identification using the same framework in the membrane example.In this example, we validate our method on a well-known non-linear DE system and demonstrate the robustness of our method to different levels of measurement noise. The double pendulum in Fig. 3 is defined over a single variable, time x, with the quantities of interest being the cart's displacement q, and the angle of link1 ψ and link2 ϕ. The black-box external source function is the force f applied to the cart. The system comprises three differential equations without constraints:EQUATIONEQUATIONh 2 q cos ψ + h 4 ψ + h 5 φ cos(ψ -ϕ) -h 7 sin ψ + h 5 φ2 sin(ψ -ϕ) = 0 (8)EQUATIONwhere h 1 through h 8 are scalars determined from Appendix. B.1.We construct three datasets without and with observation noise drawn from N (0, 0.01) and N (0, 0.03) as visualized in Fig. 3 . Each dataset comprises 24 examples of inverted pendulum balancing with various initial conditions (see Appendix. B.2). Each example runs from 0 to 4 seconds at 50 fps, resulting in a total of 200 timestamps (u m , ψ m , ϕ m ) ts=0:200 with the initial condition.Training: When we use the shorthand notations e 1 θ , e 2 θ , and e 3 θ for the left-hand sides of Eqs. 7, 8, and 9, respectively, the loss function becomes:arg min θ λ 1 |Ω| Ω (e 1 θ -f θ ) 2 + (e 2 θ ) 2 + (e 3 θ ) 2 + λ 2 |Ω m | Ωm (u θ -u m ) 2 + (ψ θ -ψ m ) 2 + (ϕ θ -ϕ m ) 2 .Here, the time domain Ω is {x|0 ≤ x ≤ 4}, where the domain Ω m ⊂ Ω is at the observations. The parameters λ 1 and λ 2 are weights for the PDE residual loss and the regression loss, respectively.In this example, we demonstrate how our method effectively handles a highly ill-posed inverse problem given partial and noisy observations. This example deals with a softbody system interacting with a rigid terrain with unknown geometry. The goal is to predict 1) a full deformed softbody geometry and 2) contact pressure, given wrist wrench w m ∈ R 3 and partial/noisy segmented pointcloud P m .A softbody with linear elasticity interacting with a rigid environment is known as a Signorini Problem. This problem is defined in 3D coordinates u of an undeformed softbody x ∈ R 3 . The quantity of interest u(x) ∈ R 3 is the deformation field, and the forcing function f (x) ∈ R 3 is contact pressure applied perpendicular to the surface. Following [36] , our problem is:EQUATIONEQUATIONwhere σ(x) is the stress tensor, n is the normal vector at the query, ρ is the density, g is the gravitational constant, D is the elasticity tensor, and ε(x) is the strain tensor. Domain Ω is the entire volume of interest with the boundary ω b . The domain Ω s is the four sides of the cubed sponge, where we assume no contact for simplification. Eq. 11 represent the force equilibrium at the contact and the infinitesimal volume at x, respectively. f (x) is always greater than 0 because normal tractions can only be compressive. Detailed formulations can be found in Appendix-B.3.We select 20 examples from [13] 's real-world dataset, where each data sample consists of a single deformed sponge's partial and noisy measurement (P m , w m ). The sponge is a 46 mm cube with a Young's modulus E = 1.1 × 10 4 and a Poisson's ratio ν = 0.1. The sponge was attached to the Pandas Franka Emika robot with a force-torque sensor (ATI-Gamma) mounted at the wrist. This dataset includes ground truth contact locations for evaluation, obtained by directly observing the contact location through a transparent acrylic plate.Training: Our model takes the wrist wrench and partial point clouds as input, and it predicts the full deformation and contact pressure. If we define the entire 3D query space as Ω = [-0.43, 0.43] 3 , all six surfaces as Ω b , and four sides as Ω s , the loss function is defined as: where CD is a single directional Chamfer Distance from P m to estimated deformed surface pointcloud P θ = {x + u(x)|x ∈ Ω b } and w θ is a wrist wrench estimation calculated from the estimated contact pressure as illustrated in Appendix. C.2. The third term represents an inductive bias that there are no contacts at the sides, and the fourth term represents a Signorini constraint, where σ n = n • σ(x)n. The last two terms are residual losses from Eq. 11.EQUATIONHere, we emphasizes the capacity of our method to solve high-order and complex differential equations' inverse source problems. We show our model's ability to perform system identification with the same framework and robustness to partial and noisy observations, unlike the prior works [9, 10] .The classical theory of Von Karman [37] provides a 4th-order PDE model of a thin elastic plate system subjected to large deflections under external forces. This problem is defined in 2D in-plane coordinates α and β. The quantity of interest is the 3D deformation [u, w] corresponding to the undeformed in-plane coordinates of the membrane. The forcing function f (x) ∈ R represents the contact pressure applied perpendicular to the surface normal. The governing equation and the boundary conditions of the membrane sensor are given by:EQUATIONwhereD = Et 3 /12(1 -ν 2) is a constant of flexural rigidity consisting of Young's modulus E and Poisson's ratio ν, t is the plate thickness, the biharmonic operator∆ 2 is ∂ 4 ∂α 4 + ∂ 4 ∂β 4 + 2 ∂ 4 ∂α 2 ∂β 2 , σ αβ is the shear stress, u ∈ R 2 representsthe in-plane displacement, w denotes the out-of-plane deflection, and p is the traverse load per area from air pressure and gravity. The membrane is subject to a Dirichlet boundary condition in the domain Ω b . Please refer to Appendix B.4 for further details.Dataset: Our dataset D = {P m , p m } consists of partial and noisy point clouds P i and internal air pressures from the pressure sensor p m . We recorded 24 data samples from 6 different geometries with 4 interactions per object, as shown in Fig. 2 . For each interaction, the Soft Bubble randomly rotates in yaw ∈ [-90 • , 90 • ] and lowers until the contact force reaches 8 ∼ 11 N. We use a Kuka arm equipped with an ATI-gamma F/T sensor mounted on the wrist to move the Soft Bubble. The reaction wrench measurement from the F/T sensors is used only for evaluation purposes of our contact pressure estimation. The Soft Bubble's membrane has a thickness of t = 0.45 mm and is an ellipse with major axis a = 0.06 m and minor axis b = 0.04 m.System Identification: Our network can also perform system identifications via optimizing for unknown parameters. When we use an inflated sensor pointcloud without contact, the sourcing term is all zero f θ = 0 and Young's modulus is the only unknown parameter. The loss function becomesEQUATIONwhere g E,θ is the left hand side of Eq. 14,Ω = {(α, β, w)|( α a ) 2 + ( β b ) 2 ≤ 1}, Ω b is the boundary of Ω, CD is one-directional Chamfer Distance from noisy partial pointcloud measurement P m to our predicted full pointcloud P θ = (α, β, 0) + (u θ , w θ )|(α, β) ∈ Ω}. L1 Error Full State Source Functionq ψ ϕ fNo noise 2.117e-3 1.587e-3 7.735e-4 3.825e-4 N (0, 0.01) 2.505e-3 2.808e-3 1.528e-3 1.640e-3 N (0, 0.03) 4.683e-3 3.945e-3 5.632e-31.262e-2 Training: Given a pointcloud under contact P m , the loss function for contact pressure estimation isEQUATIONThe only difference between Eq. 16 is non zero f θ and that the Young's modulus is not updating. 
double pendulum problem: 
 Metrics: We use the L1 error to assess the full states and the sourcing function estimations at unseen time queries with ×4 dense samplings.Result:The results in Tab. 1 and Fig. 3 demonstrate that our model can almost perfectly predict the ground truth for both the full states and the source function when there is no measurement noise. As noise is introduced, our model's accuracy gradually decreases. This trend is more pronounced in the source function estimation because accurately determining n-th order derivatives from noisy data becomes significantly more challenging with noise. Nonetheless, average force errors of 1.262e-2[N] from noise N (0, 0.03) are still a very small fraction of the scale of the ground truth (ranging from -0.3 to 0.3 [N]) as visualized in Fig. 3 .Metrics: We consider a point to be in contact when the estimated contact pressure exceeds 1,500 Pa, found to be optimal through grid search. The average Chamfer Distance between the ground truth and the estimated contact patch in our approach is 45 mm 2 , whereas the baseline achieves 32 mm 2 .Baseline: Our SOTA baseline [13] tackles this ill-posed problem by training a neural network on a dataset comprising a total of 3,000 sponges-environment interactions. One-third of the data includes sponge-box interactions, similar to their real-world dataset's sponges-tables interactions. They utilize Isaac Gym's softbody simulation [38] , based on finite element method (FEM), which is non- differentiable and can only address forward problems. The baseline [13] takes wrist wrench, partial point clouds, and a trial code as inputs, and produces signed distances and binary contact as outputs.Result: We demonstrate that our results are both qualitatively and quantitatively comparable to the SOTA baseline [13] (Fig. 4 ). This is impressive considering that our model had to infer from an infinite number of contact possibilities from partial and noisy measurements and physics equations, whereas the baseline [13] relies on strong priors obtained from the data it was trained on.System Identification: We used ν = 0.5 like typical rubber [9] and p = 103, 320P a from the pressure sensor measurement. The resulting Young's modulus is E = 341, 260P a.We evaluate the predicted contact pressure with 1) contact patch and 2) net contact force estimation. For contact patch evaluations, we use Intersection over Union (IoU) and bidirectional Chamfer Distance (CD). We classify a point in contact if the predicted contact pressure satisfies f θ (x) > max(2500, 2|Ω| Ω f θ d|A|), where 1 |Ω| Ω f θ d|A| is the average contact pressure across the entire bubble surface area [9] . Here, the ground truth is the intersection of the Soft Bubble pointcloud and the object mesh, using known transformations. For contact force evaluations, we use L2 norm between the estimated wrench at the wrist (Fig. 2 ) following Appendix. C.2 and the ground truth.Baselines: Our SOTA baseline [9] is a Finite-element (FE) based approach, highly specialized in the Soft Bubble's 3D contact force and contact location estimation. Unlike our approach, [9] requires correspondence tracking between mesh before and after contact, a special system identification process requiring a total of 205 data points, and a result refinement step for handling partial and noisy observations using a convex optimization (CVXPY). We utilized their finest mesh resolution with 749 vertices along with their best-identified model calibration results.Results: Fig. 5 shows examples of bubble-object interactions and the resulting contact pressure estimation from our model. Tab. 2 indicates that our method excels at contact patch predictions for all shapes except for Y, producing 15% higher IoU and 22.2% smaller chamfer distances when compared to the baseline. Fig. 6 shows visualizations of both our and the baseline's contact patch estimations, highlighting that our predicted contact masks are much smoother, less noisy, and more precise at following the shapes of the objects. Although the predicted contact pressures outperformed results on contact patch estimation, the predicted contact force produced a 1.083 N higher contact force L2 error, which is about 10.8% of the average contact force scale, 9.99 N. Table 2 : Soft Bubble's contact patch and contact force estimation on six different shapes on our approach and a FEM-based baseline [9] . Force is a L2 norm of the contact force error. 
future works and limitations: 
 While our approach has shown effective in solving inverse source problems in some real-world robotics applications, we have yet to explore real-time inference capacity and representing multiple DE solutions with single networks. Interesting future work is to parameterize multiple DE solutions, similar to methods used in [13, 39, 16] , and integrate the latent space inference for potential realtime applications. Additionally, while our method relies on the known governing equations of the system, an exciting direction involves replacing these known differential equations with surrogate models identified from data [40, 41, 42] . This could further enhance the applicability of our method towards unidentified systems. In the field of robotics, there has been a growing focus on identifying unknown sources that influence system behavior. Particularly in the manipulation domain, the system inputs can originate from two primary sources: the robot itself and the environment. Inputs from the robot are typically measurable through various onboard sensors, such as force/torque (F/T) sensors or visio-tactile sensors.In contrast, inputs from the environment are not directly measurable, making their identification a crucial area of research. A significant topic in this research is the inference or identification of contact formations as unknown system inputs. For example, studies such as [13, 14, 16] have focused on identifying contact locations of elastic tools (e.g., sponges, spatulas) using partial point clouds and wrench data. Other studies, like those by [17, 43, 44, 18] , have investigated extrinsic contact identification using visio-tactile measurements. In other domains of robotics, unknown system inputs often take the form of external disturbances. For instance, in underwater and aerial robotics, estimating current or wind velocity is essential for achieving optimal robot motion [45, 46] .Additionally, across various robotic applications, identifying the actual input to a sensor (such as force or deformation) from the transduced signals (like electric signals or depth images) is another form of an inverse source problem. For instance, [19] focused on estimating tactile sensor deformations from electric signals, while [9, 10] studied contact pressure estimation using depth camera outputs from sensors. It's important to note that these previous works utilize either data-driven approaches [13, 14, 16, 17, 43, 44, 18, 19] or traditional methods often based on discretizations [9, 10, 45, 46] . These methods have certain limitations, which are discussed in Sec. 1.The parameters used in Eqs. 7, 8, and 9 are as follows:EQUATIONEQUATION)where m 1 is the mass of the cart, m 2 is the mass of the first link, m 3 is the mass of the second link, l 1 is the length of the first link, l 2 is the length of the second link, and g is the gravitational constant (9.81 m/s 2 ). The cart's mass is m 0 = 0.05 kg, and the lengths of the two links are l 1 = l 2 = 0.5 m with masses m 1 = m 2 = 0.005 kg.We give angle offsets in ψ or ϕ range from [-0. Here, we show detailed formulation used for Eq. 13 implementation. First, 3x3 tensors σ(x) and ε(x) are EQUATIONEQUATIONEQUATIONLinear stress and strain tensor relationship σ(x) = D • ε(x) was in practice implemented following Voigt notation for computation efficiency as follows:EQUATIONMoreover, the divergence is formulated asEQUATIONFollowing [47, 48] , we get the following equations:EQUATIONEQUATIONEQUATIONEQUATIONThe objects we selected mimic real-world items, as shown in Fig. 7 . For example, our cylinder resembles cans from the YCB object set (e.g., Pringles can, tuna can). Our dataset was collected in real-world conditions, with 3D-printed shapes specifically designed for precise quantitative analysis of contact location estimations. Using real-world objects (e.g., YCB dataset) for evaluations poses a challenge because obtaining their precise world poses is often difficult. This difficulty stems from object pose tracking methods, such as AprilTags, which are prone to errors caused by camera distortions, depth noise, and occlusions. In contrast, our 3D-printed objects are bolted to the table with known world locations, allowing us to access their precise positions with less than 1 mm of error. In summary, when using general 3D shapes in real-world scenarios, our model will still make predictions, but there will be significant degradation in our evaluation metrics due to the lack of access to ground truth. 4 : Specifications of training and network hyper-parameters. s is the input dimension, q + r is the network output dimension, L is the number of layers, and hf is the hidden features. For all approaches, we can fit all data points in a single batch.We train our model on NVIDIA RTX A6000 GPU. For all experiments, we use Adam Optimizer with a learning rate 1e-3 with a decay rate 0.9 per every 5,000 iteration steps. Tab. 4 shows experiment-specific details.Given a set of query points (α, β), we perform a Delaunay triangulation [49] to generate 2D triangle meshes. Using deformation field predictions from our network, we displace the vertices and calculate the deformed triangle area as described in line 6 of Alg. 1. Then, we calculate the average contact pressure force exerted at each triangle following line 7 of Alg. 1. The total wrench is the summation of the forces from all triangles as described in line 8 of Alg. 1.For the FEM baseline method [9] , we use the same optimized values for the parameters: force penalty k force = 0.3322 N -1 , displacement penalty k disp = 537592 m -2 , and membrane-thickness normalized Young's modulus E = 856 Pa, as reported in [9] and its associated codebase.Input: queries -an array of query points (α, β) on deflated bubble ellipse (0.06, 0.04) Input: deformation -an array of 3D deformation [u, w] at (α, β) from PINN Input: pressure -an array of pressure value p at (α, β) from PINN 1: contact force ← 0 2: for ∆ in DelaunayTriangulation(queries) do We show that our approach can be applied even when the right-hand side of the equation consists of the source function and its derivatives:EQUATIONThe RLC circuit in Fig. 8 is described by the equation:EQUATION)where v is the voltage at the capacitor and f is the black box source function. For training the dataset, we simulated the RLC circuit with the ground truth source function 6e -3t , resulting in voltage measurements v m = 44 3 e -2t + 1 3 e -5t -9e -3t . We simulated the circuit for 2 seconds and generated 1k measurements. The source function estimation from training produces only a 5.450 × 10 -3 V absolute mean error, aligning well with the ground truth as visualized in Fig. 8 .In this section, we evaluate our approach and FEM against observation noise, as shown in Fig. 11 and Tab. 5. Specifically, we added Gaussian noise with zero mean and standard deviations of 0.0001, 0.001, 0.002, and 0.003 [m] to depth measurements. Our method demonstrates robustness up to a standard deviation of 0.002 [m], showing only a 12% drop in accuracy, whereas the FEM-based approach and the baseline [9] show performance decreases of 36% and 69%, respectively. This is an impressive result because the baseline [9] explicitly accounts for observation noise via a group lasso optimization problem, while our approach leverages the inherent noise regularization ability of neural networks without task-specific noise removal treatments. 11 : Comparison of contact pressure estimation between our method, FEM, and the baseline [9] . Each column shows the contact pressure distribution given raw depth measurements and with Gaussian noise having standard deviations of 0.0001, 0.001, 0.002, and 0.003 [m]. The baseline includes an internal 'regularization' process to correct point cloud measurement noise. The 'FEM' variant is our modification of [9] , which omits the 'regularization' phase and relies solely on FEM for solving the inverse problem. The light blue regions in the FEM and [9] results indicate contact locations found via thresholding.We propose to add a fully connected layer for each output dimension, a component absent in the architecture of prior works [25, 22] . We demonstrate that this additional layer brings significant improvements, particularly in source function estimation. Qualitatively, in Fig. 9 's contact pressure estimation results, our approach with the final layer shows sharper and more accurate contact patch estimations, closely matching the ground truth. In contrast, the contact pressures estimated by [25, 22] are widely dispersed and fail to capture the shapes of the objects in contact. Furthermore, quantitatively, for three randomly selected examples, we observe an increase in Intersection over Union (IoU) from 0.498 to 0.661 for the flower, from 0.504 to 0.780 for the cylinder, and from 0.298 to 0.667 for the 90-degree curve, respectively, as shown in Fig. 9 and Fig. 10 .In this section, we demonstrate how the sinusoidal input mapping [25] proposed in prior works makes a model susceptible to measurement noise, which is inevitable in real-world experiments. Using 6 : Ablation studies on input mapping given measurement noise N (0, 0.01). We report the L1 norms of the full state estimation and forcing function estimation. q and f are in meters, and ψ and ϕ are in radians.the same double pendulum system as described in Sec. 4.1 with measurement noise N (0, 1), we compare the results with and without the sinusoidal input mapping of [25] . We utilize four examples with initial offsets of -0.03 in ψ, ϕ, ψ, φ, respectively. Fig. 12 qualitatively demonstrates that our model without the input embedding exhibits superior robustness in force function estimation under input noise. Quantitatively, Tab. 6 indicates that using input mapping produces 47% more errors in force term estimation compared to not using it (ours).Moreover, we demonstrate how sinusoidal input mapping degrades the forcing term estimation results in contact location estimation performance. Our proposed architecture, without the input mapping, increased the IoU score from 0.538 to 0.661, 0.606 to 0.780, and 0.399 to 0.667, as visualized in Fig. 9 . This represents a performance increase of 22%, 28%, and 67%, respectively.Conventionally, as shown in [22, 25] , prior works have insisted on keeping the norm of each weighted loss's gradients equal to each other for a fixed number of epochs. However, these weighting methods, which have been classically used in AI4science, are not only 1) expensive in time and space, especially with high-order and non-linear differential equations, but also 2) ineffective for the Inverse Source Problem. This is because the solution to inverse source problems tends to remain near the randomly initialized forcing term while quickly minimizing the DE residual loss. We propose a new approach that involves using a fixed loss term that produces higher loss gradients (by a factor of 50) on regression losses that use partial or noisy observations at the beginning of the training. This approach facilitates the network in first converging to the left side of the DE using observations, so that the forcing term converges to the identified left-side value via the slowly converging DE residual loss. We show the effectiveness of our loss weighting in Fig. 13 and Tab. 7. The results, for instance, show that the loss weighting with λr λreg = 100 shows 90 % more accuracy with different loss weightings, ranging from λr λreg = 500 (second column) to 5 (sixth column). The first column of the first row shows the ground truth interaction between the membrane and the 3D-printed objects (pink), with the point cloud observation visualized in green. The second to sixth columns display the contact pressure estimation results, with the ground truth contact patch outlined in blue. The first column of the second row presents the ground truth contact estimation (blue) overlaid with a partially noisy point cloud observation (green). The color scale of the contact pressure distributions matches that of Fig. 4 and Fig. 5 . We observe that contact location detection for both cases significantly deteriorates qualitatively when λr λreg < 50. A quantitative analysis is provided in Tab. 7. 7 : Quantitative evaluations of results of different loss weighting in Fig. 13 , where the evaluation metric for the soft body is IoU and for the membrane is Chamfer Distance (CD).in IoU for the membrane example and 54 % less errors in softbody example compare to λr λreg = 1 suggested in the prior works [22, 25] .We evaluate our method's ability to handle high-order differential equations through the membrane equation's forward problem. This example involves a simplified membrane system described by the equation:EQUATIONwith Dirichlet boundary conditions. The main differences between Eq. 30 and Eq. 14 are that Eq. 30 does not account for membrane stress in the force balance equation, which holds true only if the deformation is small enough and in-plane deformation u = 0. Moreover, Eq. 30 assumes there is no contact, i.e., f = 0. This simplified system differs from the actual system, resulting in larger deformations. However, removing the stress term makes solution derivation much easier, which is useful for quantitative analysis on deformation predictions. Here, we solve a forward problem given non-updating system parameters (E, ν, p, D, t) as follows:EQUATION 
 
 Appendix:  


Decision: **Post-Rebuttal Guidance**
Guidance: ACCEPT. The author rebuttal was very strong, with precise clarifications on what exactly the contributions are, additional comparisons to baselines and ablations. There was 1 strong accept and 1 weak accept. The reviewer (**FvMB**) who provided a weak reject has not acknowledged/responded to the rebuttal nor updated the score - I had posted my thoughts to reconsider their score in light of the rebuttal.

I would like to convey this message to the authors, which was brought up by the reviewer (**dC7P**): 
"The authors should also discuss and cite other recent works in the robotics community that have used PINN, for planning [1], and control [2,3].

[1] NTFields: Neural Time Fields for Physics-Informed Robot Motion Planning [2] Neural Optimal Control Using Learned System Dynamics [3] Value Iteration in Continuous Actions, States and Time [4] Physics-Informed Neural Networks for Continuum Robots: Towards Fast Approximation of Static Cosserat Rod Theory"

**Pre-Rebuttal Summary**

Strengths
- Elegant solution
- Strong quantitative evaluation, including on the real world

Weaknesses
- Limited technical novelty
- The presentation of the paper could be more clear
- Questions on generalization to more complex real data (not 3d printed shapes for example)

Dear authors, looks like the reviews are mixed and have several clarifying questions - please take a look.

Reviews: review: title: title: Review of CoRL submission 423
paper_summary: paper_summary: This work proposes to solve inverse source problem of PDE for robotics applicaiton, solving source function and state mapping at the same time. The proposed method is based on PINN and many examples are given to show the results, including one real-robot example.
main_review: main_review: The motivation of the paper is sound and clear, and it is important to adopt neural networks to find the solutions of natural dynamics as underlying PDEs in robotics applications. The experiments are extensive and detailed, showing that the PINN-based method can give results for different robotic tasks. However, there are some weaknesses below.

- The core idea is simple and not novel: neural networks are used to parameterize the unknown PDE parts as well as the solution, and the losses are based on PDE residual errors and boundary conditions. The pipeline for training PINNs seems to be typical in the AI4Science area. Therefore, the contribution of introducing PINNs into the robotics area seems to be incremental.

- In different examples, the NN architecture is the same and the only difference lies in different losses and data, which further weakens the contribution and cannot provide insights into how PINNs will be designed differently for different robot tasks.

- The first example is a common non-linear control problem, and it is expected to compare the proposed method with other off-the-shelf control or RL or other learning-based methods, to fairly show the superiority of PINNs through this traditional robot control example.

- The organization of the paper can be improved. The example section is too long and includes some details about the method, experiment setup, and results, which look a bit messy and not coherent from the view of the whole paper.
questions: questions: See weakness above.
limitations: limitations: 2
review_summary: review_summary: Incremental contribution, lack of insight for PINNs in robotics, poor oganization.
title: title: Introduce a existing tool to robotics area
paper_summary: paper_summary: This paper proposes to use PINN to solve the inverse source problem in robotics
main_review: main_review: This work uses PINN to solve the inverse problem in the robotics area, e.g. recover the shape from tactile.

The paper demonstrates extensive experiments including simulation and real-world setups. 

PINN method is new for the inverse problems in robotics and it may lead to multiple following works. 

However, the technique in the paper is not very novel. There are similar PINN network structure designs [1]. 

I also have some questions about the comparison with FEM results. PINN is known for its slower convergence and low accuracy compared with the FEM method [2]. I am curious about the better performance of PINN in the paper. Is it due to the noisy observation data? Will PINN have better performance for fully observed data?

Also, the scalability is another problem. Is it possible to scale to a real-world dataset but not 3D-printed shapes?

[1] PirateNets: Physics-informed Deep Learning with Residual Adaptive Networks

[2] Can Physics-Informed Neural Networks beat the Finite Element Method?
questions: questions: 1. Could you explore more on FEM results? Is it due to the noisy data? Can PINN outperform FEM in simulated experiments without noise?
limitations: limitations: 3
review_summary: review_summary: The paper introduces a new tool to the inverse problem in robotics. It can have a good impact for the following methods.
title: title: Well written paper on solving inverse source problems with PINNs.
paper_summary: paper_summary: In this paper, the authors introduce a novel way to solve inverse source problems by exploiting Physics Informed Neural Networks with a modified Loss function. This structure allows to reconstruct the external source function based on partial and noise observations of the system. The proposed method is evaluated on simulations and real-world datasets. In comparison to existing approaches, it shows superior performance regarding the accuracy of the predictions.
main_review: main_review: The paper addresses an important problem in robotics which is the identification of the source function that acts on a PDE or ODE system. Even though the approach seems simple, the simulation and real-world experiments demonstrates the effectiveness of the proposed method. Instead, I would rather say that the simplicity is actually the beauty of the approach. In this way, I think that this is an original idea which can have some impact on the future research on inverse source problems.

Major comments:

- I think that this approach suffers from the main problem of PINNs which is that since the physics are enforced by the loss function evaluated at some points, a non-zero loss means that there are (minor) violations of the physics / constraints at the evaluation points and there is no guarantee what happens in between the points. 

- Furthermore, the literature review does not cover all relevant areas as it’s limited to NN approaches. For instance, there is a rich literature on latent force models that address the same issue, see, for instance

Alvarez, Mauricio A., David Luengo, and Neil D. Lawrence. "Linear latent force models using Gaussian processes." IEEE transactions on pattern analysis and machine intelligence 35.11 (2013): 2693-2705.

- Nevertheless, the experimental verification makes a strong case that the presented method can generate reasonable results in the shown use cases

Minor
47: typo “limited to solvinng"

UPDATE AFTER REBUTTAL: Thank you for addressing all my concerns. I updated my score and recommend to accept the paper.
questions: questions: -	Can you please extend the literature section including other approaches such as latent force models?
limitations: limitations: 3
review_summary: review_summary: Simple but effective approach to learn the source function of an ODE or PDE with PINNs with convincing numerical evaluations.