Network Inference From Local Measurements: Application to Coordination of Groups of Mobile Three-Dimensional Printers

Conventional construction methods are often associated with physical and economical inefficiencies [1–3], thereby prompting interests in alternative construction paradigms. One of the most promising alternatives, three-dimensional construction printing (3DCP), leverages additive manufacturing technology to assemble a structure. Traditionally, 3DCP relies on layer-by-layer printing, using gantry- or robotic arm-based systems [4–6]. In the first case, the print nozzle is precisely localized and controlled in a 3D Cartesian space defined with respect to the gantry. In the second case, the material is extruded from a print nozzle attached to the end-effector of a stationary robotic arm. Despite their accuracy, these methods restrict the size of objects that can be printed to the size of the printer [7,8]. Numerous efforts have sought to overcome this limitation in scalability and achieve large-scale 3DCP. For example, in Ref. [9], a mobile 3D printing mechanism—consisting of a robotic arm with a print nozzle mounted on a mobile platform—was proposed in order to expand the printing space and print objects of sizes larger than that of the printer.

The limitation in scalability motivates the use of multi-agent printing systems based on multiple mobile 3D printers. To the best of our knowledge, multi-agent printing systems have only been explored in Refs. [10–12]. Therein, authors proposed a multi-agent 3D printing system capable of printing a large, single-piece structure. The printing tasks were allocated to each printer by means of a control system relying on a full connectivity between the agents. However, whether in urban or remote settings (such as in space exploration), construction inevitably suffers from uncertainties and disturbances. The environment is constantly evolving as more layers are printed: continuous communication among agents is difficult, if not impossible, to guarantee. As such, decentralization of the collective behavior of printing agents is a necessity to build robust, versatile multi-agent 3D printing systems.

In the decentralized collective 3DCP paradigm, several autonomous printers are tasked with printing together a large structure, whose layout is initially available. Without continuous communication with other printers, the printers cannot apprehend the progress of the process throughout the spatial extent of the structure. Hence, the printers must perform some form of sensing that allows them to reconstruct the actual printed shape of the structure, such that they can infer the progress of the other agents and act accordingly.

Here, we propose an active sensing framework that relies on scalar measurements from on-board sensors to infer the shape of the structure during printing. More specifically, printers shall use the physical medium they are creating to actively sense the addition of layers by other printers and therefore track the progress of the task. While several quantities could be measured by the printer, for simplicity, we focus on temperature, which is a scalar variable whose time evolution throughout the structure is governed by the classical Fourier equations [13]. The latter, once discretized, are equivalent to consensus dynamics for first-order integrators [14].

We focus on 2D geometries, which are representative of walls or other plate-like structures that could be printed by teams of mobile printers, similar to previous efforts in additive manufacturing [15–17]. In this vein, we approximate the shape of the structure through a 2D lattice, over which heat propagates. Ultimately, the problem of identifying the shape of the structure becomes that of reconstructing the topology of the lattice, a fundamental problem in network theory. Several methods for topology reconstruction have been proposed in the literature, for both deterministic and stochastic networked systems, see, for example, Refs. [18–22] and Refs. [23–28], respectively. A related problem, borrowing several methodological tools from these studies, is the inference of the size of a system, which is a first step before the complete topology can be determined [29–33].

To this end, we assume that the structure behaves as a linear time-invariant (LTI) system [34] during the addition of each layer by any of the printers. This hypothesis is valid for slow printing speeds, such that measurements of temperature are taken in a period of time short enough for new layers not to significantly affect the shape of the structure. In this way, the shape of the system is embedded within the time-invariant state-matrix that describes the heat propagation over the 2D lattice. Thus, by reconstructing the state-matrix of the system, we can also infer its topology. The problem then becomes one of identifying the state-matrix of the LTI system in order to reconstruct the shape of the structure.

We utilize the so-called subspace identification method to retrieve the state-matrix and other matrices of the system up to a similarity transformation, from input and output time-series [35]; the successful application of equivalent techniques to mechanical systems is documented in Refs. [33] and [36]. The identified, transformed state-matrix does not directly give access to the topology of the system: we need to reconstruct the original system matrices. Should one possess a priori the output-matrix, they could obtain the similarity transformation and from that the original system matrices [22]. If all the system matrices are unknown, the problem of reconstructing the original system matrices can be formulated as an inverse eigenvalue problem [37], where one searches for the specific topology that allows for matching the eigenvalues of the identified state-matrix.

The alternating projections algorithm provides a potential approach to solve this problem [38], but it is plagued by the need of constructing a database of graph matrices, whose number grows rapidly with the order of the system. To overcome this issue, we formulate the inference problem as an inverse eigenvalue problem within a smaller search space of topologies that are proximal to a nominal topology defined by previous instances of shape reconstruction or by the printing plan, which is known by each agent. This framework is valid when the differences between the actual and nominal topologies are small (that is, the frequency at which printers perform the shape reconstruction is sufficiently high). A major hypothesis in the solution of the inverse eigenvalue problem is that the specific topology we are reconstructing is univocally determined by its spectrum (DS-graph), that is, it is the only topology (up to a permutation) in the minimization space that generates a state-matrix with the same eigenvalues of the identified state-matrix [39].

We present a set of necessary identifiability conditions that need to be satisfied in order to retrieve the correct shape of the structure. Through numerical simulations, we confirm that, provided that the hypotheses are verified, our approach can reconstruct the shape of the structure up to a permutation. Further, we investigate the validity of the hypothesis that the topology is a DS-graph, by counting the number of nonisomorphic cospectral graphs in 2D lattice topologies with increasing number of blocks [40]. These simulations are an original extension of the studies on square polyominoes (plane figures formed by the union of square blocks sharing edges) [41], addressing the issue of counting the number of cospectral polyominoes.

The rest of this paper is organized as follows. We formulate the mathematical problem in Sec. 2. Next, we present a review of the subspace identification framework in Sec. 3. Section 4 introduces our approach to infer the state-matrix of the system, and hence the topology of the network. In Sec. 5, we investigate the identifiability of the system, and we validate the performance of the reconstruction through simulations. We conclude in Sec. 6 with a summary of the proposed framework and an outline for future work.

We consider an instance of decentralized collective 3DCP, illustrated in Fig. 1, where mobile 3D printers collectively build a 2D structure, in the ξ–η plane. Each 3D printer generates heat in the structure due to the printing process, propagating across the printed structure [42]. The 3D printers cannot communicate with each other, so that they seek to reconstruct the shape of the structure from measures of the temperature through on-board sensors. This operation is repeated by each printer at regular time intervals, such that shapes from previous reconstructions are available. In between reconstructions, the shape of the structure changes as printers progress with their printing jobs.

where τ denotes the time, $T(ξ,η,τ)$ is the temperature distribution, μ is the thermal diffusivity of the printing material, ρ is the material density, c is the material specific heat capacity, $q˙v(ξ,η,τ)$ is the known volumetric heat rate from the heat sources (printers), and $∇2(·)$ is the Laplacian operator. As a first approximation, we assume that the structure has adiabatic boundaries.

The heat equation in (1) can be discretized in space and time using the finite difference method [43]. For simplicity, we assume that the spatial discretization step is equal in both the ξ- and η-directions, $Δs=Δξ=Δη$⁠. We then consider a temporal discretization of the equation with the backward Euler method, with uniform time-step $Δτ$⁠. We hypothesize that the heat that serves as input is known, while the temperature is measured at a subset of the points thanks to on-board sensors.

We assume that the printing speed is slow compared to the time scale of heat propagation, such that new layers do not significantly affect the shape of the structure while the 3D printer takes the measurements. Whether or not a printing speed is slow depends on the thermal diffusivity of the printed material. Since we only consider 2D problems, we are interested in the maximum areal printing speed, defined as the areal speed at which the printer is spanning the 2D environment. The maximum potential areal printing speed is estimated from the maximum flowrate based on the depth of the structure; in our estimates, we utilize a depth on the order of one-tenth of a meter, which is typical of wall-like structures in civil construction. The printing speed can be safely regarded as slow if the maximum potential areal printing speed is higher than the thermal diffusivity of the printed material. We assume that measurements are taken at a sufficiently high rate, which is a reasonable assumption for the print of civil structures that span an area on the order of tens of square meters.

Based on these arguments, the assumption of slow printing speed should be considered with caution for concrete structures, since thermal diffusivity (⁠$10−6−10−5 m2/s$⁠) could be on the same order or lower than the maximum possible areal printing speed (⁠$10−5−10−4 m2/s$⁠) [44,45]. On the other hand, the hypothesis is safely verified for metallic structures, whereby the bound provided by thermal diffusivity (⁠$10−3−10−2 m2/s$⁠) is much higher than the maximum possible areal printing speed (⁠$10−5 m2/s$⁠) [46]. For this maximum areal rate, the time to print a square meter of structure is on the order of hours, such that it is plausible that measurements are taken at a sufficiently higher rate. When the printed structure is not solid, the thermal diffusivity of the printed material should be corrected to account for the reduced heat propagation, introducing a smaller effective thermal diffusivity [47,48], which could hinder the validity of our assumption. Nevertheless, the real applicability of our hypothesis is likely to extend to a broader range of scenarios, since the maximum potential areal printing speed is an upper bound of the actual printing speed, that printers attain only when printing at their maximum flowrate.

Under the hypothesis of slow printing speeds, we treat the system as an LTI system. The shape of the structure is described by a network $G={V,E}$ of N nodes $V={v1,…,vN}, E$ being the edge set. We collate the temperature of each node in the set at time-step t in a state vector $x(t)∈ℝN$⁠, the m inputs from the heat generated by the 3D printers into an input vector $u(t)∈ℝm$⁠, and the temperature measurements on l nodes from the on-board sensors into an output vector $y(t)∈ℝl$⁠.

where $A∈ℝN×N$ is the state-matrix, $B∈ℝN×m$ is the input-matrix, and $C∈ℝl×N$ is the output-matrix.

Our approach aims at inferring the topology of the network from the input and output time-series ${u(0),…,u(M−1)}$ and ${y(0),…,y(M−1)}$⁠, respectively, where M is the length of the time-series that we measure. The inference of the topology of the network requires first to estimate the number of nodes in the topology and second to identify how these nodes are connected. To this end, we utilize basic knowledge about the system dynamics to reconstruct the topology by means of system identification techniques. The topology is assumed to be embedded in the unknown state-matrix A that describes the dynamic evolution over the network.

We formulate the topology reconstruction problem as follows: find the adjacency matrix $A$ of the 2D lattice topology generating a state-matrix $A$ from Eq. (4) such that the system in Eq. (2) provides the measured output ${y(0),…,y(M−1)}$ for an input ${u(0),…,u(M−1)}$⁠. As a first step to solve this problem, we utilize the input–output time-series to retrieve the state-matrix $AT$ up to a similarity transformation (⁠$AT=T−1AT$ where $T∈ℝN×N$ is a nonsingular transformation matrix) through the subspace identification method [35].

We rely on the subspace identification method to retrieve the system matrices up to a similarity transformation [35]. While sufficient to capture the dynamics of the system, the transformed state-matrix does not readily provide the topology of the network, for which we need the original state-matrix; we will address this issue in Sec. 4.

The subspace identification method relies on the geometrical properties of the LTI system expressed in Eq. (2) to obtain the matrices of the LTI system up to a similarity transformation. To correctly identify the system matrices, this approach requires the system to satisfy a set of identifiability conditions. First, the system must be minimal, that is, it must be both observable and controllable [34]. Observability refers to the ability of reconstructing the state of a system from measurements of its output. Controllability refers to the ability of the input to drive the state of the system to the zero state in a finite time interval. Second, the input signal $u(t)$ in Eq. (2a) must be persistently exciting [50]: the spectrum of the input must contain a sufficiently large number of harmonics.

The LTI in Eq. (2) is observable if and only if $rank(O)=N$⁠, and is controllable if and only if $rank(C)=N$ [34].

This expression provides a practical condition to verify whether the input signal is persistently exciting. Note that, since $Up∈ℝms×p, p≥ms$ is a necessary condition for the system to be persistently excited; given that s > N, then also p > N.

Note that if s < N or p < ms, it is not possible to estimate N correctly.

In this case, the condition s > N guarantees that $UR(1:(s−1)l,1:N̂)$ has rank equal to $N̂$⁠, such that $[UR(1:(s−1)l,1:N̂)]†UR(1:(s−1)l,1:N̂)=IN̂$ [35]. As a result, $AT$ is not rank-deficient through our estimation.

However, the transformed state-matrix $AT$ does not provide information regarding the structural properties of the network described in the system. In Sec. 4, we propose a framework to reconstruct the actual state-matrix and therefore the network topology by solving an inverse eigenvalue problem.

where $Λ∈ℝN̂×N̂$ is a diagonal matrix containing the sorted eigenvalues of the matrices $AT$ and A.

where $N̂$ is the estimated order of the network obtained from Eq. (15), $S$ is the minimization space that contains the family of all possible graphs of size $N̂$⁠, and $Λ=eig(AT)$ is computed from Eq. (19). To accurately infer $N̂$ and $AT$⁠, the LTI system must be identifiable: (1) the system must be minimal (that is, observability Eq. (5) and controllability matrices Eq. (6) are full rank [34]), and (2) the input-signal must be persistently exciting [51] as in Eq. (12).

We investigate the minimality of the system, and hence the reconstructability of the topology of the structure, under conditions of increasing complexity on the system input- and output-matrices B and C: (1) fully actuated and measured systems (⁠$B=C=IN$⁠); (2) partially actuated and fully measured systems (⁠$B≠IN$ and $C=IN$⁠); and (3) partially actuated and measured systems (⁠$B≠IN$ and $C≠IN$⁠).

In this case, all the nodes in the network are actuated and measured (⁠$B=C=IN$⁠). The system is minimal as it is observable (⁠$rank(O)=N$⁠) and controllable (⁠$rank(C)=N$⁠). The identifiability of the system depends only on the excitation condition, whereby the input signal must satisfy the persistent excitation condition in Eq. (12).

The topology is finally reconstructed by plugging Eq. (22) in Eq. (4) and inverting this relation to obtain the adjacency matrix. Due to numerical approximations, $Â$ may contain elements that are not exact integers: in practice we round each element to the closest integer.

The chosen assumptions on the input- and output-matrices are restrictive for collective 3DCP applications. For large scale structures, the printers cannot actuate and measure all the nodes along the spatial extent of the structure. Hence, it is necessary to loosen such restriction to ensure practicality of our framework for shape reconstruction.

Similar to the previous case, all the nodes in the network are measured (⁠$C=IN$⁠), so that the system satisfies the observability condition. We loosen the restriction on the input-matrix and assume only a subset of the nodes are actuated. Hence, the identifiability of the system is not guaranteed and depends on the choice of actuated nodes such that the system is controllable. In case of successful selection of the nodes to ensure minimality, we leverage knowledge of the output-matrix and reconstruct the topology using Eq. (22), similar to the previous case.

Despite loosening the restriction on the system inputs, we still rely on the unrealistic assumption on the output-matrix (⁠$C=IN$⁠). In reality, each printer has access to only a subset of the nodes, located in its local environment. In fact, we cannot assume any prior knowledge of the input- and output-matrices as the topology and size of the network are initially unknown. As such, it is important to investigate the inference and reconstruction of the shape of the structure for the general case of partially actuated and measured system (⁠$B≠IN$ and $C≠IN$⁠).

Here, we actuate and measure only a subset of the nodes. Hence, the identifiability of the system depends on the choice of actuated and measured nodes such that the system is controllable and observable. Further, in the previous two cases, we use knowledge about the output-matrix to obtain the estimated state-matrix. However, in this case the output-matrix is unknown and it is not possible to compute the transformation matrix T in order to obtain the estimated state-matrix $Â$⁠. As such, we make use of the relation in Eq. (20) that holds for similar matrices and solve an inverse eigenvalue problem to reconstruct the network topology.

Several challenges arise when solving the problem in Eq. (23). First, the problem is constrained to the integer domain since all entries of the adjacency matrix must be binary. Second, even in solving the problem as a mixed-integer nonlinear programing problem, the objective function is nonconvex due to the existence of cospectral nonisomorphic graphs, which share the same spectrum but do not have the same shape. Consequently, the optimization problem must be constrained beyond the spectral properties of the network in order for the estimation to converge to the desired topology.

where $Â22∈BδN̂×δN̂$ encodes the connections between the nodes of the unknown additive perturbation, and $Â12∈BÑ×δN̂$ considers the connections between the two sets of nodes. The minimization only considers variations in matrices $Â12$ and $Â22$⁠, whose number of elements is much smaller than that of $Â$⁠. Thus, we reformulate the integer optimization problem presented in Eq. (23) over a smaller minimization space.

It is important to note that the minimization searches for a 2D lattice topology that matches the eigenvalue set $Λ$ obtained from the identification. To guarantee the convergence to the exact topology of the system, the topology must be distinguishable by its spectrum across the minimization space [39]. A graph is determined by its spectrum (DS-graph) with respect to a generalized adjacency matrix (the actual adjacency matrix or the combinatorial Laplacian) if it does not have a nonisomorphic cospectral counterpart with respect to that matrix [39]. While 2D lattice topologies are not DS-graphs in general, we hypothesize that the topology we need to reconstruct is a DS-graph over the reduced minimization space. As such, we expect the estimated topology to match the topology of the network up to a permutation $P∈BN̂×N̂$$(A=PÂPT)$⁠.

In this section, we conduct a series of numerical simulations to elucidate the feasibility of our approach for shape reconstruction, as well as the plausibility of our hypotheses. First, we investigate the conditions under which the system is minimal. Second, we examine the plausibility of our hypothesis that the topology we seek to reconstruct is a DS-graph over the space of perturbed topologies. To this end, we seek to enumerate the number of 2D lattice topologies (polyominoes) for an increasing number of unit blocks, showing that only a small fractions of these topologies are not DS-graphs. Finally, we show over a few different examples that our approach can accurately reconstruct complex printed shapes.

We start by investigating the influence of the number of actuated and sensed nodes, in terms of fraction of the total number of nodes, on the estimated size of network $N̂$⁠. Specifically, we generated rectangular topologies embodied in 2D lattice representation of different orders: N = 30, 50, and 80. During each iteration, we also randomize uniformly the number of nodes on each side of the rectangular topologies, while maintaining constant the overall number of nodes. We first considered a fully observable system, $C=I$⁠, with a persistently exciting input signal generated using the randi function on matlab. Further, we selected a fraction of the nodes m/N as actuated nodes. The nodes were chosen uniformly at random and the simulations were performed for different selection of m/N. For each combination of N and m/N, simulations were repeated for 50 trials. We found that $m/N≥0.2$ guarantees that the system is controllable (Fig. 2).

Similarly, in order to investigate the minimum number of measured nodes to identify the system, we considered a fully controllable system (⁠$B=I$⁠) with partial observation. We repeated the simulations for different choices of measured nodes l/N. Similar to the case of the partially controllable system, Fig. 3 shows that for $l/N≥0.3$ the system is always observable.

In Sec. 4, we hypothesized that the 2D lattice topology we seek to reconstruct is a DS-graph over the space of perturbed topologies. In other words, by matching the eigenvalues of the state-matrix, we could uniquely determine the network topology among the entire set of perturbed topologies, up to a permutation. While for special graphs some results on cospectral graphs are available [39], there is no general test for finding cospectral graphs from properties of the graph matrices or topology.

To support the validity of our hypothesis, we provide evidence through combinatorial simulations. Specifically, we conducted simulations to count the cardinality of the set of all possible topologies that can be formed for a given number of square blocks $NB$⁠. We then grouped isomorphic topologies in subsets such that any two graphs selected from independent subsets were nonisomorphic. Finally, we compared the spectrum of all subsets to check for the existence of nonisomorphic cospectral graphs.

In order to generate the set of nonisomorphic graphs, one can rely on a brute force approach, searching for all possible graphs and grouping each graph with its isomorphic counterparts. However, this approach is computationally very demanding. We relied on an exhaustive approach leveraging the symmetries in 2D lattice topologies to generate all their isomorphic counterparts. This procedure has been proposed in Ref. [52] for counting hexagonal and triangular polyominoes [41]. The symmetries of 2D lattice topologies include four rotations (0, 90, 180, and 270 deg), and for each of these rotations a reflection about the horizontal and vertical axes passing through a central block (Fig. 4).

Consider a set of nonisomorphic (nonidentical) graphs consisting of $NB−1$ blocks. For each graph, we add a block to an adjacent grid and performed these transformations about a reference block of the graph (denoted by a red dot in Fig. 4, in the second row from the bottom and middle column) to obtain its isomorphisms, keeping record of the explored adjacent grids in these transformations. In each step, we examine a new unexplored adjacent grid until all the adjacent grids are explored. As such, we divide the set of all possible graphs into subsets, each containing a list of isomorphic graphs. We repeat these steps for all possible graph configurations, while discarding graphs already generated from the transformations. By doing so, we substantially decrease the computational time in generating the subsets.

We generated the families of 2D lattice topologies for $NB$ from three to ten; we summarize the results in Table 1 where NI-graphs denotes the number of nonisomorphic graphs, CS-graphs (⁠$A$⁠) the number of cospectral graphs with respect to the adjacency matrix, and CS-graphs (⁠$L$⁠) the number of cospectral graphs with respect to the Laplacian matrix. The counting results on the number of nonisomorphic graphs are in agreement with those presented in [41,53]. Up to ten blocks, we did not find any nonisomorphic cospectral topology with respect to the adjacency or the Laplacian matrix. We cannot exclude cospectral nonisomorphic graphs for large values of $NB$⁠. These results support the claim that only a small fraction of 2D lattice topologies are cospectral nonisomorphic graphs, such that our hypothesis that the topology we want to reconstruct is a DS-graph over the space of perturbed topologies is viable.

We demonstrate the feasibility of our shape reconstruction framework through a series of different examples. These examples include both symmetric and asymmetric structures with nominal symmetric and asymmetric topologies with a different number of layers affected by uncertainty, associated with common types of printing inaccuracies. In these simulations, we selected an appropriate rescaling of the variables such that μ = 1. In addition, we hypothesized that the inputs were applied and outputs were measured on the same subset of nodes, covering 30% of the nominal size $Ñ$ and chosen uniformly at random. The input time-series were chosen to be binary input signals generated using idinput in matlab, validated to be persistently exciting as per Eq. (12).

The minimization was performed with a brute-force algorithm. First, we generated the entire set of perturbed topologies about the known nominal topology, based on the difference between the number of nodes in the nominal and actual topologies. Specifically, we started from the nominal topology and explored all the ways of appending nodes until the required node count was reached. Next, we evaluated the objective function in Eq. (20) for each perturbed topology. We then selected the estimated topology by choosing the perturbed topology that minimizes the cost function in Eq. (20). We demonstrate two examples of topologies consisting of 52 and 58 nodes in Fig. 5. In both cases, we were successful in exactly inferring the actual shape of the printed structure starting from its nominal topology and available measurements.

The examples we presented are representative of common inaccuracies in layer-by-layer printing [54]. Figure 5(a) demonstrates the case when the printed structure has distorted layers due to induced deformations and residual stresses. This is common during sequential deposition of layers during the printing process, which introduces deformations that either exceed dimensional tolerances or physically prevent additional material from being deposited [54]. On the other hand, Fig. 5(b) demonstrates deformations due to insufficient substrate support where the deposited layer is not able to serve as a solid foundation for successive layers deposited above it [54].

While the proposed framework has a low computational cost in terms of computational time, we acknowledge computational limitations that arose in populating the Hankel matrix in Eq. (13) imposed by the conditions on the length of segments s and number of partitions p. As the structure grows larger with the deposition of more layers (N increases), more inputs are required to ensure the system is controllable (m increases, see Fig. 2) and longer time-series are needed to retrieve the size of the system in Eq. (14) (s increases, since s > N). Hence, the number of partitions increases (p increases, since $p≥ms$⁠). The size of the Hankel matrices exceeded the admissible memory allocation in matlab of 20 gigabytes for large topologies (N > 100).

In this work, we consider the problem of collective mobile 3D construction printing, where a group of 3D printers collectively prints a structure. We aim at taking a first step toward the decentralization of these systems by proposing an active sensing framework that leverages modification in the environment to reconstruct the shape of the printed structure. By doing so, we facilitate communication among the agents to allow individual agents to track and assess the evolution of the printed structure and eventually take corrective actions.

The shape of the printed structure is approximated by a 2D lattice, which we seek to reconstruct. We rely on input and output time-series of the heat propagation across the structure and previous shape reconstructions or knowledge of the printing plan to infer the actual shape of the structure. By modeling heat propagation as a linear time-invariant system, we employ system identification techniques to reconstruct the lattice topology. First, we leverage the subspace identification method [55] to retrieve the system matrices up to a similarity transformation, provided that the system meets a series of identifiability conditions: the system must be minimal and persistently excited. We then solve an inverse eigenvalue problem to reconstruct the 2D lattice topology, over a reduced minimization space that leverages the shape reconstructed in the previous measurement instance or the a priori knowledge of the printing plan by the agent to define a series of potential perturbed topologies.

We investigate the reconstruction problem under the following conditions imposed on the system input- and output-matrices B and C: (1) fully actuated and measured systems (⁠$B=C=IN$⁠); (2) partially actuated and fully measured systems (⁠$B≠IN$ and $C=IN$⁠); and (3) partially actuated and measured systems (⁠$B≠IN$ and $C≠IN$⁠). We further investigate the inference problem under partial measurement and actuation. Specifically, we demonstrate the minimality of the system by varying the number of driver and observer nodes as a fraction of the total number of nodes. We show that it is sufficient to actuate and measure $30%$ of the nodes in the network to ensure a minimal representation condition, necessary for successful identification.

In the proposed framework, we hypothesize that the 2D lattice topology we seek to reconstruct is a DS-graph over the space of perturbed topologies: it is sufficient to have knowledge of the eigenvalues of the graph with respect to a graph matrix in order to infer its eigenvectors. We provide evidence through combinatorial simulations to support this hypothesis. For topologies consisting of $NB$ square blocks, we generate the set of all possible 2D lattice topologies, leveraging their symmetries, similar to what was done in Ref. [52]. We find that, up to $NB=10$⁠, there are no nonisomorphic cospectral topologies with respect to the adjacency or the Laplacian matrix. We further demonstrate the feasibility of our shape reconstruction framework through simulated examples.

Although our framework accurately reconstructs the exact topology of the system, it is not free of limitations that we plan to address in future work. First, while the proposed methodology improves on previous techniques in terms of scalability, it is still computationally unfeasible for large networks. To address this issue, one may consider including a time component to the print plan known a priori. By predicting the trajectory of the printing nozzle, the minimization space can be further reduced such that the perturbation are dictated by the printing direction. Large networks could also cause the Hankel matrix in Eq. (13) to become ill-conditioned, such that the subsequent RQ and SVD factorizations could be numerically unstable. In this case, algorithms for computing the factorizations should be carefully selected to ensure relevant results even for ill-conditioned matrices [56,57]. Second, the proposed framework assumes that the discrepancy between the nominal and actual topology is small. This assumption may be restrictive if there are significant errors in the printing. We plan to extend our framework by incorporating partial connectivity between the 3D printers, based on physical proximity, allowing some of the agents to directly share information and identify larger errors in the printing. Third, the hypothesis of slow printing speed should be relaxed to apply our framework to any printable material, for which the time-scale of conduction may be similar to the one of printing. A potential approach to relax this constraint is considering the topology as a time-varying network, in which the time-scale of the addition of nodes is governed by the printing speed. Fourth, environmental factors may cause issues in the implementation of the current shape reconstruction algorithm. Not only may the 3D printed structure be subject to other forms of heat exchange, such as convection and irradiation, but also measurements from sensors may be noisy. Addressing this limitation requires the extension of our framework to stochastic noisy systems, in line with previous efforts in the literature [22]. Last, while our investigation provides an insight on conditions that guarantee the minimal system representation, one is often restricted on the choice of actuation and sensing nodes. An interesting approach to solve this problem could be designing the initial printing plan to ensure a minimal system representation.

Overall, our work demonstrates the possibility of actively utilizing a physical medium to enable indirect transfer of information between agents, a framework that could be extended beyond collective 3DCP to design and understand collective systems. One example is in modeling the interactions in fish schools utilizing the physical medium to maintain a collective swimming formation [58]. One can think of hydrodynamic interactions between the fish as perturbations introduced in the environment that dictate the formation maintained by the school [59], working in conjunction or at odds with visual cues. Surprisingly, some species [60] can school in the absence of vision, suggesting that active modification of the environment may be of critical importance to the school to maintain formation.

• National Science Foundation, Directorate for Engineering (Award No. CMMI 1932187; Funder ID: 10.13039/100000084).

The datasets generated and supporting the findings of this article are obtainable from the corresponding author upon reasonable request.