## Abstract

Packing and routing separately are each challenging NP-hard problems. Therefore, solving the coupled packing and routing problem simultaneously will require disruptive methods to better address pressing-related challenges, such as system volume reduction, interconnect length reduction, ensuring non-intersection, and physics (thermal, hydraulic, or electromagnetic) considerations. Here we present a novel two-stage sequential design framework to perform simultaneous physics-based packing and routing optimization. Stage 1 generates interference-free initial layouts that are fed to stage 2 as starting points to perform continuous physics-based optimization. Three distinct strategies for stage 1 have been introduced recently, (1) the force-directed layout method (FDLM), (2) an extension of the shortest path algorithms (SPAs), and (3) a unique geometric topology (UGT) generation algorithm. In stage 2, a gradient-based topology optimization method is used to simultaneously optimize both component locations and interconnect routing paths. In addition to geometric considerations, this method supports optimization based on system behavior by including physics-based objectives and constraints. The proposed framework is demonstrated using three case studies. First, the layout generation methods developed for stage 1 are compared with respect to system performance metrics obtained from stage 2. Second, a multi-objective optimization problem using the epsilon-constraint method is solved to obtain Pareto optimal solutions. Third, an extension to multi-loop systems is demonstrated. In summary, the design automation framework integrates several elements together as a step toward a more comprehensive solution of 3D packing and routing problems with both geometric and physics considerations.

## 1 Introduction

Power electronics circuits and fluid cooling systems (and many other types of other engineering systems) are composed of components that exchange energy, as well as routing (interconnects) that facilitate energy transfer. One reason these systems are difficult to design is that they have many requirements, including, performance, cost, geometry, and volume restrictions. Identification of feasible designs can be exceptionally difficult in applications where space is limited, components and interconnects involve complicated geometries, and system performance depends on spatial relationships and multiple physics couplings. Current practice relies largely upon human expertise, design rules, modification of existing designs, and manual adjustments to solve these problems. This limits both the complexity of systems that can be designed involving non-trivial packing and routing decisions, as well as the realization of potentially improved functionality or performance. Many previous efforts have focused on creating design automation methods to address elements of the integrated packing and routing problem described here, but not the combined problem. In this paper, a two-stage design automation approach is presented that integrates several of these elements as a step toward a more comprehensive solution of 3D packing and routing problems with both geometric and physics considerations. Development of such methods for the formalization, automation, and prioritization of topological layouts for complex, multi-domain systems will offer benefits across several industry sectors.

Design automation methods for solving the optimal spatial packing and routing problem have been developed and studied previously in the context of many applications, such as vehicle assembly [1], electronic module layout design [2], 3D container loading [3], bin packing [4], computer animation [5], the layout of components in additive manufacturing [6], and automotive transmission design [7]. The packing problem in most cases is also considered as a layout placement problem [8], where component geometries can be arbitrary, with multiple types of design goals and spatial constraints satisfaction. For practical purposes, the minimization of layout cost functions is done under certain constraints imposed by design, fabrication, and operational requirements. Most layout algorithms are restricted to a certain class of systems and are as a whole intractable due to their combinatorial nature. Previous research algorithms used in 2D layout packing can be generally classified under four categories: genetic algorithms [9], heuristic methods [10], gradient descent algorithms [11], and simulated-annealing-based algorithms [12].

Furthermore, many efforts have addressed the interconnect routing problem specifically, where the component layout is fixed. Especially in the electrical engineering domain, many examples of 2D routing algorithms were developed for VLSI circuit layouts based on Manhattan rules and its variants [13]. Some other applications include aero-engine routing [14], ship pipe routing [15], chemical plant pipe routing [16], electrical wire routing in buildings [17], in developing CAD-based FPGA (field-programmable gate array) design tools [18], unmanned aerial vehicle navigation [19], and robotic path planning [20]. Optimization approaches have incorporated metrics such as packaging volume and mass properties [21], and have utilized solution methods such as simulated annealing [22,23], pattern search [24,25], genetic algorithms [26], ant colony optimization [27], and several other heuristic methods [28].

The main limitation of the above methods is that they address the component packing and interconnect routing problems separately, instead of in a combined manner that accounts for inherent coupling. In addition, most of the methods consider only geometric aspects of the problem, neglecting important physical system properties such as operating temperature, pressure drop, thermal loading, aerodynamic, and electromagnetic effects. Thus, existing methods may not extend well to the general coupled problem. In addition, the performance evaluation of the designs obtained from existing is left to human designers. The amount of time required for a designer to generate a feasible design and analyze its performance limits the ability of engineers to explore these complex design spaces within a constrained project timeline. Existing strategies can produce feasible designs, but they may not be optimal when considering all the system requirements, and the complexity of systems that can be considered is limited. In current practice, many aspects of layout and routing problems are solved manually, which severely limits design capabilities for systems involving complex packing and routing tasks (especially in cases with strong physics interactions). In this research expanded from our previous work in Ref. [29], coupled computational methods will be presented that have the potential to generate designs faster, with better system performance, and for higher-complexity systems than those designed with methods that require significant human input.

### 1.1 Objectives and Contributions.

The objective of the study presented here is to develop a novel computational design framework for optimally distributing system components and interconnect paths within a two-dimensional volume, subject to performance objectives that rely on physics coupling between continuum and lumped-parameter models. The primary contributions of the paper include: (1) a new two-stage sequential design framework for performing optimal packing and routing of 2D interconnected electro-thermal systems, (2) a rigorous comparison of three different automated feasible layout generation methods (stage 1) with respect to system performance metrics obtained in stage 2, and (3) a demonstration of multi-objective and multi-loop optimization thermal-fluidic system case studies. The design methodology presented here as follows:

1. Supports systems with different component geometries, port connectivity, and physics-based operating conditions, such as temperature, hydraulic pressure, and fluid flowrate.

2. Is applicable to small- to medium-scale systems (up to approximately 25 components). For a system-of-systems problem, each component can be considered as a subsystem. Extension to large-scale systems is a topic of ongoing work.

3. Uses a gradient-based topology optimization method based on differentiable geometric projections recently developed by the authors [30].

4. Supports geometric constraints to ensure non-interference between device-device, device-interconnect, and interconnect-interconnect combinations.

5. Determines the optimal packing and routing layouts for various system-level performance metrics, such as bounding box volume, total head loss, and maximum component temperature.

The remainder of this paper is organized as follows. Section 2 describes the two-stage framework and a 2D geometric layout design representation for a generic fluid-thermal system used in this research. Section 3 reviews the three feasible layout generation methods used for obtaining feasible layouts. The geometric projection layout method used for simultaneous optimization of packing and routing in stage 2 is presented in Sec. 4. Section 4.1 shows the various system-level objective function formulations for a simple hybrid unmanned aerial vehicle (UAV) example. A set of three different case studies are presented in Sec. 5 with the discussion of results in Sec. 6. Finally, the conclusion and future work topics are provided in Sec. 7.

## 2 Two-Stage Design Framework

Although discrete formulations have dominated earlier research efforts to solve the optimal spatial packing problem (OSPP), these approaches quickly become intractable, as both the packing problem and the routing problem are NP-hard. Rather than make incremental progress on established methods, we are proposing a novel strategy for the OSPP using a continuous representation applied to distinct geometric topologies. This enables the use of gradient-based methods to efficiently search the packing and routing design space. This methodology is centered on the use of simple geometric bars to approximately model both the component geometry and routing paths (see Fig. 1). Bars have favorable geometric properties for both the packing and routing problems, enabling simultaneous solution. This simple geometric representation with physics-based boundary conditions (heat source, heat sink, and forced convection) will be used as part of the two-stage framework. Specifically, the combined physics-based packing and routing problem can be subdivided into two important sequential stages as shown in Fig. 2: (1) generate spatially feasible geometric layouts without any interference between the components and interconnect routing (for example, pipes or ducts) and (2) to utilize these feasible layouts as initial designs for performing multi-physics optimization to meet system-level requirements such as packaging volume, operating temperature, fluid pressure, and cost. It is assumed here that the continuous optimization strategy cannot make topological changes, two classes of which are described below.

Fig. 1
Fig. 1
Close modal
Fig. 2
Fig. 2
Close modal

Geometric topology (GT) refers to the mathematical study of the properties that are preserved through continuous deformations, such as stretching and twisting, but not tearing or gluing. For instance, a circle is topologically equivalent to an ellipse, but a circle and a hollow ring have different topologies. Here, system architecture (SA) specifies what components comprise a system, and how ports on components are connected to specific ports on other components. For each system architecture, many geometric topologies may exist. For example, if an interconnect links ports P1 and P2, many options may exist for how this interconnect passes around various other interconnects and components in the system. The SA enumeration problem for electro-thermal systems was studied extensively in Ref. [31]. Figure 3 shows three layout options for a system with three components that have two ports each. All three have the same system architecture (the component connectivity matrix is the same), which is one aspect of system topology. Layouts A and B differ only in the interconnect path; they have the same geometric topology. Layout C has a distinct geometric topology; i.e., the red-colored interconnect cannot be continuously morphed on a plane to achieve the geometric topology exhibited in layouts A and B without either leaving the plane or requiring cuts or reconnections.

Fig. 3
Fig. 3
Close modal

## 3 Stage 1: Spatially Feasible Layout Generation Methods

Mathematical optimization methods often require initial designs to satisfy constraints (i.e., initial points that are feasible). The goal of stage 1 of this framework is to generate initial layout designs that are feasible with respect to geometric interference constraints. This feasibility is often required to perform associated physics-based analysis to assess performance. Stage 1 also provides an opportunity to investigate distinct geometric topologies (GTs), assuming that continuous optimization in stage 2 cannot make discrete switches between unique GTs. To develop a simultaneous optimization method to combine these two stages when using gradient-based methods, we need a continuous relaxation technique that can generate smooth gradients and is compatible with physics-based finite element analysis. Implementing this approach, however, will require a significant effort, and should be addressed as a topic for future work to compare it against the proposed sequential design framework. This task beyond the scope of this article.

Attaining feasible layouts in stage 1 can be achieved either by enforcing explicit interference constraints between objects or implicitly satisfying interference constraints via strategies such as the force-directed layout method (FDLM), shortest path-planning algorithm (SPA), or the unique geometric topology (UGT) enumeration method [32]. This section discusses these three (FDLM, SPA, and UGT) feasible layout generation methods in detail below.

### 3.1 Force-Directed Layout Method.

Force-directed graph drawing algorithms were originally developed by Tutte et al. [33]. They have been used for graph untangling in the context of several visualization applications, such as gene networks and computer graphics. The spring force-based graph layout algorithm is the most popularly used method [34,35] and is utilized here.

In the FDLM, a physical system with its components and interconnects is represented as a graph-based force model. The adjacent graph nodes have attractive forces between them and repulsive forces with remaining nodes. The entire graph is considered to be a system of particles (nodes) with spring forces acting between them. The initial random graph layout has high overall energy and is unstable because of the forces between the nodes. The end goal is to position these nodes so that each node has locally minimum energy. In other words, force-directed techniques seek an equilibrium state where the sum of the forces acting on each particle (node) is zero. At equilibrium state, all the springs are at their natural length l. This value of l in addition to the spring stiffness k (proportional to the spring force) ensures that the components and bars do not interfere. For electro-thermal system topologies, the nodes of the graph are the components centers, the interconnect bar centers, and their end points. The graph edges are formed by the springs in between these nodes. Figure 4 shows a 2D representation of the FDLM where all system elements are free to move only within a hard bounding box (HBB). Distributed inward boundary forces (shaded region) are applied from the walls to ensure the objects do not escape the HBB.

Fig. 4
Fig. 4
Close modal

Figure 5 shows examples of how the FDLM is applied to multi-component systems with and without interconnects, respectively. Figures 5(a) and 5(b) show both the initial and final layouts of a ten-component system without interconnects. The springs connecting the components are indicated as straight lines. Similarly, Figs. 5(c) and 5(d) show an example where component and interconnect routing untangling were both considered together. The spring network is not shown for better visibility of the topology. For multi-component systems with interconnects, the boundary forces have been removed to allow free untangling of the interconnect network. The advantage of the FDLM is that it generates feasible layouts without the need for explicit interference constraints. But the limitation with FDLM is that it does not support well complex interconnect routing with many bars. This is because only bar centers and endpoints are taken as nodes in the force-directed graph, and this makes untangling at other areas of the bars difficult where the magnitude of the spring forces is less.

Fig. 5
Fig. 5
Close modal

### 3.2 Shortest-Path Based Layout Algorithm.

Shortest path algorithms (SPAs) are grid-based methods that are widely used in applications such as robot path planning [36], obstacle avoidance [37], and the traveling salesman problem [38]. SPAs use grid-based techniques where obstacles occupy a group of cells within the grid known as closed cells, and the aim is to reach a target point from a source point avoiding the obstacles. To achieve the task the robot moves toward the target through the remaining empty cells on the grid. We exploit this property of the SPAs to perform interconnect routing for a multi-component system with a fixed architecture. It is initially assumed that the components are not overlapping and distant from each other. A SPA can be directly applied to find a feasible path between two ports on two different components treating components and defined interconnects as obstacles. This process is repeated sequentially for every connection to define all required routing paths. For an n-component system with two-port components, $(n2)$ combinations of routing sequences exist, and thus a maximum of $(n2)$ unique GTs. It is possible that unique sequences may generate the same GTs. Some sequences might fail if they do not have empty cells to form a feasible route. In that case, a finer grid can be used to perform more complex routings. Grid refinement, however, will increase the computational expense.

The A* algorithm, which is a well-known and efficient SPA [39], is utilized here to perform sequential interconnect routing. A* finds the shortest path between the source and the target nodes using a heuristic method. Figure 6 shows an example of the A* algorithm being used to generate a feasible layout (see Fig. 6(b)). The components are at fixed locations; restricted regions are indicated as black boxes. The interconnect routing sequence numbers are shown in Figs. 6(b) and 6(c). Figure 6(c) illustrates an infeasible final layout because the ports P1 and P2 shown here do not have a feasible connection path due to the coarse grid (60 × 40 elements) that limits the availability of empty cells.

Fig. 6
Fig. 6
Close modal

### 3.3 Unique Geometric Topology Enumeration.

Here, a new layout enumeration algorithm is presented that can systematically generate feasible unique geometric topologies (GTs) for given system architecture (SA). This exhaustive enumeration strategy serves as a benchmark for assessing geometric topology design space coverage of other strategies. Since a continuous mapping exists between any two layouts that share a GT, the complete set of GTs is an important set of starting design points for continuous layout optimization.

Consider a three-component system where each component has two ports, as shown in Fig. 7. Each of the three components is connected directly to the other two components in GT-A. This configuration is termed the directly connected GT. Figure 7 also shows three other unique GTs (namely, Indirect1, Indirect2, and Indirect3) for the same system architecture represented by GT-Direct. Note that in these GTs one interconnect circulates over a component to reach the other component. Such a configuration is termed an indirectly connected GT. Thus, a three-component and two-port interconnected system has a total of four unique GTs (one direct and three indirect). Similarly, for a four-component system with each component having two ports, we achieve a single unique direct GT and four unique indirect GTs, as shown in Fig. 8. By mathematical induction, for an n-component system with each having two ports each, 2D layouts involve n + 1 unique GTs. This method can also be applied to a multi-component system with each component having different shapes and a varying number of ports.

Fig. 7
Fig. 7
Close modal
Fig. 8
Fig. 8
Close modal

## 4 Stage 2: Physics-Based Layout Optimization

Stage 2 uses the feasible layouts obtained from stage 1 as initial designs for performing multi-physics optimization to meet system-level requirements, such as system volume, operating temperature, fluid pressure, cost, and reliability measures. This sequential process is based on the assumption that physics-based optimization methods are gradient-based and cannot navigate discrete changes in geometric topology. Here, the geometric projection method (GPM) of Norato et al. [40] is used to solve the continuous optimization problem. In the original projection method work, the new parameterization approach was used to perform continuum-based topology optimization of structures made of discrete elements while ensuring that the resulting design could be made from stock materials, such as beams with standard shapes and sizes. The geometric parameterization involves design variables that facilitate convenient derivation of lumped-parameter model sensitivities. This supports the use of physics-based analysis in evaluating designs for system-level performance. Similarly, the packing and routing problem contains discrete elements such as components with various shapes, interconnect networks, and other geometric features. Hence, for stage 2 optimization, we extended the GPM method successfully to solve the physics-based combined layout and routing optimization problem in our recent work [30]. The GPM was used to create routing designs that can be manufactured from standard circular cross-sectional pipes. Pairs of ports between components are connected via physical interconnects, and device topology is assumed to be given and unchanging. Each interconnect is represented in the GPM using one or more straight geometric segments. Increasing the number of segments in a connection supports consideration of curved or more complex interconnect geometries, but at the cost of increased computational expense. Thus, $x$, shown in the formulation below, represents the design optimization variables that include the port locations, start and end point locations of the interconnect segments, interconnect widths, and segment to segment nodal bend radii. The complete optimization problem formulation for stage 2 is given by
$minxf(x,T)$
(1a)
$s.t.:gphys(x,T)≤0$
(1b)
$gdd(x)≤0$
(1c)
$gsd(x)≤0$
(1d)
$gss(x)≤0$
(1e)
where
$K(x)T=P(x)$
(1f)
Here, $f(x,T)$ is the objective function and $g(x,T)$ are constraint functions. In general, these functions may depend both on design ($x$) and state ($T$) variables. The function f(·) can be any one of the candidate objectives such as bounding box area, head loss, device temperature, etc. The constraints $gphys(x,T)$ are constraints that depend both on design and on solutions to the physics models (i.e., the value of the state vector $T$). The interference constraints $gdd(x)$, $gsd(x)$, and $gss(x)$ prevent interference between two devices, one routing segment and one device, and two routing segments, respectively. These constraints are independent of any physics models, so they are all explicit functions of the design variables. The 1D lumped parameter model and the 2D finite element models were developed to simulate the physics response of the multi-component interconnect system. Temperature distribution was modeled on the continuum level using the finite element method. The strong form of the boundary value problem for steady-state heat conduction is given by
$∇⋅(κ∇T(x))+Q=0,xinΩ$
(2)
$T(x)=T*,xonΓT$
(3)
$n⋅(κ∇T(x))=q*,xonΓq$
(4)
where $κ$ is the matrix of thermal conduction coefficients, $T(x)$ is the temperature solution field, Q is heat flux per unit volume in the domain, and $n$ is the unit normal to the domain boundary. Temperature, T*, and heat flux, q*, boundary conditions are applied on the ΓT and Γq portions of the domain boundary, respectively. Detailed derivation of the finite element equations and implementation can be found in Ref. [41]. Here, we will skip to the final equation that solves for temperatures at the nodes of the finite element mesh, which is obtained by discretizing the boundary value problem in Eqs. (2)(4) using the finite element method.
$KT=P$
(5)
Equation (5) is solved for the temperature field vector $T$, where $K$ is the global thermal stiffness matrix assembled from element stiffness matrices, $Kel$, defined in Eq. (6), and $P$ is the global load vector assembled from element load vectors, $Pel$, defined in Eq. (7).
$Kel=∫ΩeBTκBdΩ−∫∂ΩhhNNTd∂Ωh$
(6)
$Pel=∫ΩeQNdΩ+∫∂ΩhhTenvNd∂Ωh$
(7)
Here, $N$ and $B$ are element shape function and shape function gradients, respectively. These equations also include convection boundary conditions on the ∂Ωh portion of the boundary. The temperature of the convecting fluid is Tenv, and the convection coefficient is h (assumed constant here). The sensitivity calculations for the objective functions and constraints are described in detail in recent work by the authors [30].

### 4.1 Stage 2 Illustration: Unmanned Aerial Vehicle Test Platform.

To illustrate stage 2, a notional power electronics cooling system for an unmanned aerial vehicle (UAV) as represented in Fig. 9 was optimized. This test platform has rich potential in multi-physics optimization for thermal, fluid, and mechanical packing. Both the system architecture and the geometric topology of the system are fixed during the optimization. The initial system layout is depicted in Fig. 10(a), and corresponding device properties are given in Table 1. The system consists of two battery packs, an AC/DC converter, and a heat exchanger. The battery packs and AC/DC converter add heat to the system, and the heat exchanger removes heat. A fixed-location inlet and outlet for the fluid loop are placed on the left edge. Boundary conditions for the thermal problem are shown in Fig. 11. The top and bottom edges of the domain have convection boundary conditions with a convection coefficient of h = 35.4 W/(m2K) and environment temperature of 0 °C. The convection boundaries on the top and bottom edges can be related to the surfaces of the UAV exposed to external cooling by the environment during flight. The right edge has a fixed temperature of 100 °C. For instance, this could be considered as the heat coming from the engine during flight, and the left edge is insulated. The maximum pipe diameter is allowed to be 0.03 m. There are two free points in each connection to allow for more complex interconnect routing paths. The optimization was solved using both head loss and bounding box objective functions.

Fig. 9
Fig. 9
Close modal
Fig. 10
Fig. 10
Close modal
Fig. 11
Fig. 11
Close modal
Table 1

Device properties

Device numberDescriptionQd (W/m2)Tmax (°C)
1Battery500030
2Battery500030
3AC/DC converter100070
4Heat exchanger−2000
Device numberDescriptionQd (W/m2)Tmax (°C)
1Battery500030
2Battery500030
3AC/DC converter100070
4Heat exchanger−2000
The optimization formulation for minimizing the system bounding box area is given by
$minxA(⋅)=(max(x¯)−min(x¯))(max(y¯)−min(y¯))$
(8a)
$subjectto:Tdi=1⋯4≤TdmaxMax.devicetemp.constraint$
(8b)
$Tf≤TfmaxMax.fluidtemp.constraint$
(8c)
$Hl≤HlcMax.headlossconstraint$
(8d)
$ggeo≤0Geometricconstraints$
(8e)
$where:K(x)T=P(x)Physics-basedmodelequations$
(8f)
where $x¯$ and $y¯$ are the set of x and y coordinates of device reference points and bar segment end points. The critical device temperatures, Tc, for all devices are shown in Table 1. Head loss constraints, $Hlc$, will be used for the bounding box objective function. Here, all geometric constraints from the generic optimization formulation have been lumped together into $ggeo(⋅)$. The function in Eq. (8a) represents a rectangular bounding box containing all free devices and interconnects that is aligned with the x and y axes. Similarly, the objective function for minimizing the pressure head loss objective function is given by
$minxf(⋅)=Hl$
(9a)
$subjectto:Tdi=1⋯4≤TdmaxMax.devicetemp.constraint$
(9b)
$Tf≤TfmaxMax.fluidtemp.constraint$
(9c)
$ggeo(⋅)≤0Geometricconstraints$
(9d)
$where:K(x)T=P(x)Physics−basedmodelequations$
(9e)

Fig. 12
Fig. 12
Close modal
Table 2

Power electronics cooling system optimization results

ObjectiveHl (m)Bounding box (m2)T1 (°C)T2 (°C)T3 (°C)
Bounding box1.500.31130.022.917.2
ObjectiveHl (m)Bounding box (m2)T1 (°C)T2 (°C)T3 (°C)
Bounding box1.500.31130.022.917.2

Note: Optimal objective function values are highlighted in bold.

## 5 Case Studies

### 5.1 Case Study 1: Comparison of Layout Generation Methods.

Sections 3 and 4 demonstrated methods for generating feasible layouts and the simultaneous physics-based optimization method. Here, we compare the FDLM, SPA, and UGT methods in terms of how their resulting layouts perform as start points in stage 2. A multi-start approach is used where all feasible layouts from stage 1 are used as initial designs for stage 2 optimization problems with bounding box area and pressure head loss objective functions separately to see how it scales over an increasing number of components n. For the case studies, we restrict the class of topologies to only n-component systems with each component having two ports each (say, nC-2P). The FDLM methods begin with both components and routing overlapped. The SPA and UGT methods have non-overlapping components initially.

The problem specifications are as follows: The multi-component interconnected system has n equally-sized square devices, with each device dissipating at least 3000 W/m2 of heat, and with each additional device individually dissipating 50 W/m2 more than its predecessor. For example, a system with ten devices will have the 10th device generating 3000 + (10 − 1)50 = 3450 W/m2 of heat. The devices are connected in a single loop. Each connection has one free intersection point. A short fixed pipe segment is attached to each port. This allows constraints between devices and all free pipe segments to be enforced. The same boundary conditions shown in the UAV illustration study are enforced here. The maximum temperature constraint on all devices is set to 30 °C. For the bounding box objective function, the head loss constraint is set to 5 m for all cases of n ∈ {3, 4, …, 9, 10}. The optimization problem terminates when one of the following convergence criteria is satisfied: (1) objective function step is within a prescribed tolerance, i.e., $|fn+1(x)−fn(x)|≤δf$ or (2) design variable vector change magnitude is within a prescribed tolerance, i.e., $‖xn+1−xn‖≤δx$. Here, δf and $δx$ are set to be 10−7 and 10−5, respectively. All computations in this work were performed using a workstation with an Intel Xeon E5-2660 CPU @2.00 GHz, 64 GB DDR4-2400 RAM, Windows 10 64-bit, and matlab 2019b.

Figure 13 shows the final optimal bounding box objective function values for all the initial feasible layouts. Each layout is represented on the plot using different bolded circles. The direct UGT is the best performer and the average values of the indirect UGTs are worse because of the longer interconnects that reduce the potential to shrink. The objective function values for the indirectly connected UGTs and the SPA layouts for a given nC-2P system are averaged to keep it consistent with the other solutions. Figure 14 shows the initial feasible and final optimal layouts of a six-component two-port system (6C-2P), corresponding to design points (a)–(d) in Fig. 13. Similarly, Fig. 15 shows the final pressure head loss objective function values for all the initial feasible layouts. Figure 16 shows the initial feasible and final design layouts of corresponding to points (a)–(d) in Fig. 15. UGT-Indirect is the best performer for the pressure head loss objective function because the longer interconnects it produces can help achieve smoother pipe bends. UGT-Direct does not significantly improve the head loss objective function as it contains shorter interconnects where a greater possibility of having sharper pipe bends exists. It should be noted that, for illustration purposes only, one representative indirect-UGT layout and one representative SPA layout are shown in Figs. 14 and 16, although the plots in Figs. 13 and 15 show averaged times for all indirect-UGTs and SPAs, respectively. In addition, the initial SPA layouts used in this case study are generated with some obstacles to perform routing so that the effect of interconnect complexity can be assessed during optimization. These obstacles or restricted regions are removed during optimization to make a fair comparison with other layout generation methods. Section 6 describes the important inferences from this case study in more detail. To achieve convergence given the prescribed tolerances, the average computation times for solving the bounding box minimization and the pressure head loss minimization problems were 2,345.5 s and 2,743.2 s, respectively.

Fig. 13
Fig. 13
Close modal
Fig. 14
Fig. 14
Close modal
Fig. 15
Fig. 15
Close modal
Fig. 16
Fig. 16
Close modal

Here, we also observed that the physics-based device temperature constraints for all feasible layouts were either active (equal to the allowable maximum temperature which is 30 °C), or well within the feasible range. For example, Table 3 shows the final device temperatures of the six components from each of the four feasible layouts shown in Fig. 14. In the final optimized layouts, devices near the boundaries had temperatures well below 30 °C, but those far from the boundaries were nearly active (or at least hotter) due to less convective cooling from ambient conditions.

Table 3

Final device temperatures of optimized layouts shown in Fig. 14 are satisfying the physics-based temperature constraints

Final layoutsT1 (°C)T2 (°C)T3 (°C)T4 (°C)T5 (°C)T6 (°C)
(a) Direct-UGT22.418.629.828.930.029.9
(b) FDLM23.220.727.722.623.4.029.3
(c) SPA2719.728.521.623.730.0
(d) Indirect-UGT29.728.926.423.024.229.0
Final layoutsT1 (°C)T2 (°C)T3 (°C)T4 (°C)T5 (°C)T6 (°C)
(a) Direct-UGT22.418.629.828.930.029.9
(b) FDLM23.220.727.722.623.4.029.3
(c) SPA2719.728.521.623.730.0
(d) Indirect-UGT29.728.926.423.024.229.0

### 5.2 Case Study 2: Bi-Objective Optimization Problem.

Multi-objective optimization studies are performed to analyze the trade-offs that exist between two or more conflicting performance metrics. The epsilon-constraint method tackles multi-objective problems by solving a series of single objective subproblems, where all but one objective are transformed into constraints. A Pareto front can be obtained efficiently for bi-objective problems. In this case study, the two competing objectives are bounding box area (f1) and pressure head loss (f2) functions. We treat the head loss objective function f2 as a constraint, varying within the range 2.2 ≤ ɛi ≤ 6.2 m over 40 equally spaced intervals and solve the corresponding bounding box minimization problems sequentially. The optimization is performed on a specific 5C-2P system, and its spatially feasible initial layout is shown in Fig. 17. We maintain the same physics-based boundary conditions, device sizes, heat dissipation rates, and temperature constraints as specified for case study in Sec. 5.1. The ɛ-constraint optimization formulation is as follows:
$minxf(⋅)=f1(x)$
(10a)
$subjectto:Tdi=1⋯4≤TdmaxMax.devicetemp.constraint$
(10b)
(10c)
$Tf≤TfmaxMax.fluidtemp.constraint$
(10d)
$ggeo(⋅)≤0Geometricconstraints$
(10e)
$where:K(x)T=P(x)Physics-basedmodelequations$
(10f)
$2.2≤εj≤6.2rangeofε-constraintvalues$
(10g)

The Pareto optimal solutions plotted in the objective function space are shown in Fig. 18. As the overall bi-objective design task is to minimize both the bounding box area and the pressure drop across interconnects, we desire points in the space close to the bottom-left corner (toward the origin).

Fig. 17
Fig. 17
Close modal
Fig. 18
Fig. 18
Close modal

Optimal solutions have a range of pressure drop values from 2.2 to 6.2 m, and it results in a range of bounding box values from 0.015 to 0.072 m2. The labels (a) through (d) that identify specific marked points in Fig. 18 correspond to the designs given in Figs. 17(a)17(d). These representative solutions (a) through (d) were chosen subjectively such that they are evenly distributed in the objective function space. Points (a) and (d) are the anchor points of the Pareto front. This case study indicates how sharp pipe bends cause more head loss across the system.

### 5.3 Case Study 3: Multi-Loop Optimization Example.

In this study, we demonstrate how the proposed framework also supports multi-loop optimization problems and can be used for more complex interconnected component layouts. A two-loop system optimization result is shown in Fig. 19, illustrating the efficacy of this framework. We maintain the same boundary conditions, device sizes, heat dissipation rates, and temperature constraints as defined for the case study in Sec. 5.1. Optimal layouts are obtained for the bounding box objective function (as shown in Fig 19(a), improving by 66.03% from 0.1443 m2 to 0.04901 m2), and for the pressure head loss objective function (as shown in Fig 19(b), improving by 40.01% from 5.98 m to 3.578 m). The two optimization studies required 1855.4 s and 2234.7 s of computation time, respectively, using the same hardware defined above. As a future task, this multi-loop system design capability will be expanded to apply to industry-relevant 2D packing and routing problems.

Fig. 19
Fig. 19
Close modal

## 6 Discussion

The primary findings from the above case studies (especially from Secs. 5.1 and 5.2) are as follows:

1. The directly connected topologies obtained from the UGTs have the best performance for the bounding box objective compared to all other kinds of topologies. This is because they can shrink very easily since there is no complex routing involved in such layouts. The indirectly connected GTs on average did not perform well for the bounding box objective function because their interconnects have longer lengths compared to the direct GTs. This restricts volume reduction significantly because there is a greater chance of the interconnects intercepting the components they are circulating.

2. The feasible layouts from the force-directed layout method (FDLM) performed well for the bounding box objective function after the direct GTs. The disadvantage of the FDLM is that it fails for complex interconnects. This turns out to be a hidden advantage in that the feasible layouts attained using FDLM are free from interconnect untangling and component overlapping. They are produced by a spring layout method that spread the system components uniformly until all the nodes reach minimal energy. This property may not be advantageous in other problem contexts.

3. The results for the SPAs with simpler topologies were good but worsened for more complex networks. This might be because larger-scale systems have interconnects that are very long, not allowing the components to pack closely.

4. The initial feasible layouts that had longer interconnects tend to have lower pressure (head) losses when pressure head was used as the objective. These layouts could better avoid sharp bends (the dominant factor in head loss), expanding within the bounding box as some of the interconnects may not be constrained between the components.

5. The direct GTs have higher head loss results. This is because the interconnects are constrained within the region between the components, and sharp bends might occur while avoiding intersections with the components.

6. The FDLM performs reasonably well for both objective functions. The interconnect untangling is done systematically, uniformly balancing all nodes in the FDLM. The feasible layouts that are produced are favorable designs and tend to be flexible in supporting both smooth bends and volume reduction.

Finally, it is important to note that the above observations only present a general notion of how the 2D layout generation methods perform in terms of the start points they provide. In the future, we plan to investigate a more comprehensive set of system types that can provide deeper insights for extending the design framework. These initial results indicate that the two-stage design framework helps obtain a variety of solutions, each having the potential to represent appropriate applications.

## 7 Conclusion

In this article, we have demonstrated a two-stage design automation approach that systematically enumerates and describes the set of feasible topological layouts of a fluid-thermal system and subsequently optimizes each layout via a gradient-based design optimization procedure that accounts for the physics-based performance of the system. Three distinct layout generation methods for generating feasible geometric layouts in stage 1 were described and demonstrated. A stage 2 physics-based optimization method was described. A notional power electronic cooling system for an unmanned aerial vehicle (UAV) was optimized in an illustrative study. The example demonstrates our technique for simultaneous optimization of device placement and interconnect routing, given a GT and initial layout from stage 1. Starting from the initial condition, both the devices and interconnects move through the domain in search of a locally optimal solution. The three layout generation methods have been compared over different system performance metrics. The UGT-Direct method was the best for the minimum bounding box objective because the direct interconnects between the components have relatively short lengths and reduced complexity, thus supporting enhanced volume reduction compared to other methods. The SPAs and UGT-Indirect methods generate layouts that, in general, have complex routing paths. Longer interconnect routing helps in minimizing the head loss in the pipes because lengthy paths can help reduce the occurrence of very sharp pipe bends in the final optimized layout. A bi-objective optimization problem with bounding box and head loss objective functions was solved in the second case study using the ɛ-constraint method to attain Pareto-optimal design solutions. The final case study demonstrated the solution of a multi-loop optimization problem using the proposed framework.

Future work items for FDLM include incorporating a continuous force field along the boundaries of the interconnect bars to prevent any interference. Several improvements can be made in stage 2, including, but not limited to, adding component rotation in the geometric projection method. This will provide more design freedom for optimization to attain better solutions. Furthermore, in future, we plan to incorporate the concept of irregularly shaped bounding boxes and components, multiple flow loops, pumps, fluid storage tanks, heat sinks, etc., to accommodate more complex system architectures. The electrical effects of the components were not considered in this work, and coupling with the already existing thermal and hydraulic effects is important for many practical systems. Multi-domain modeling and analysis is important to attain more realistic optimal solutions. Extension of the proposed framework to large-scale systems is part of ongoing work. Here, we plan to implement the UGT method for larger 2D systems by utilizing machine learning techniques. However, it should be noted that this work also aims to lay a foundation to aid design practitioners and engineers to gain deeper insights and to further develop the theory and methods required to tackle real-world 3D packing and routing problems, which present significant challenges beyond the 2D problems presented here.

## Acknowledgment

This material is based upon work supported by the National Science Foundation Engineering Research Center (NSFERC) for Power Optimization of Electro-Thermal Systems (POETS) with cooperative agreement EEC-1449548. The authors also sincerely thank Dr. Larry Zeidner and Dr. Neal Herring from the United Technologies Research Center (UTRC), as well as Dr. Joseph Zimmerman from CU Aerospace, for providing valuable feedback while serving as industry partner liaisons for this work.

## Conflict of Interest

There are no conflicts of interest.

## Data Availability Statement

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

## References

1.
Dong
,
H.
,
Guarneri
,
P.
, and
,
G.
,
2011
, “
Bi-Level Approach to Vehicle Component Layout With Shape Morphing
,”
ASME J. Mech. Des.
,
133
(
4
), p.
041008
. 10.1115/1.4003916
2.
Schafer
,
M.
, and
Lengauer
,
T.
,
1999
, “
Automated Layout Generation and Wiring Area Estimation for 3D Electronic Modules
,”
ASME J. Mech. Des.
,
123
(
3
), pp.
330
336
. 10.1115/1.1371478
3.
Yano
,
N.
,
Morinaga
,
T.
, and
Saito
,
T.
,
2008
, “
Packing Optimization for Cargo Containers
,”
2008 SICE Annual Conference
,
Tokyo, Japan
,
IEEE
, pp.
3479
3482
.
4.
Bansal
,
N.
,
Lodi
,
A.
, and
Sviridenko
,
M.
,
2005
, “
A Tale of Two Dimensional Bin Packing
,”
46th Annual IEEE Symposium on Foundations of Computer Science (FOCS’05)
,
Pittsburgh, PA
,
IEEE
, pp.
657
666
.
5.
Abdel-Malek
,
K. A.
,
Yeh
,
H. J.
, and
Maropis
,
N.
,
1998
, “
Determining Interference Between Pairs of Solids Defined Constructively in Computer Animation
,”
Engin. Comput.
,
14
(
1
), pp.
48
58
. 10.1007/BF01198974
6.
Panesar
,
A.
,
Brackett
,
D.
,
Ashcroft
,
I.
,
Wildman
,
R.
, and
Hague
,
R.
,
2015
, “
Design Framework for Multifunctional Additive Manufacturing: Placement and Routing of Three-Dimensional Printed Circuit Volumes
,”
ASME J. Mech. Des.
,
137
(
11
), p.
111414
. 10.1115/1.4030996
7.
Yin
,
S.
,
Cagan
,
J.
, and
Hodges
,
P.
,
2004
, “
Layout Optimization of Shapeable Components With Extended Pattern Search Applied to Transmission Design
,”
ASME J. Mech. Des.
,
126
(
1
), pp.
188
191
. 10.1115/1.1637663
8.
Cagan
,
J.
,
Degentesh
,
D.
, and
Yin
,
S.
,
1998
, “
A Simulated Annealing-Based Algorithm Using Hierarchical Models for General Three-Dimensional Component Layout
,”
Comput. Aided Des.
,
30
(
10
), pp.
781
790
. 10.1016/S0010-4485(98)00036-0
9.
Jain
,
S.
, and
Gea
,
H. C.
,
1998
, “
Two-Dimensional Packing Problems Using Genetic Algorithms
,”
Engin. Comput.
,
14
(
3
), pp.
206
213
. 10.1007/BF01215974
10.
López-Camacho
,
E.
,
Ochoa
,
G.
,
Terashima-Marín
,
H.
, and
Burke
,
E. K.
,
2013
, “
An Effective Heuristic for the Two-Dimensional Irregular Bin Packing Problem
,”
Ann. Oper. Res.
,
206
(
1
), pp.
241
264
. 10.1007/s10479-013-1341-4
11.
Sridhar
,
R.
,
Chandrasekaran
,
D.
,
Sriramya
,
C.
, and
Page
,
T.
,
2017
, “
Optimization of Heterogeneous Bin Packing Using Adaptive Genetic Algorithm
,”
IOP. Conf. Ser.: Mater. Sci. Eng.
,
183
(
03
), p.
012026
. 10.1088/1757-899X/183/1/012026
12.
Rao
,
R. L.
, and
Iyengar
,
S. S.
,
1994
, “
Bin-Packing by Simulated Annealing
,”
Comp. Math. Appl.
,
27
(
5
), pp.
71
82
.
13.
Koh
,
C.-K.
, and
,
P. H.
,
2000
, “
Manhattan or Non-Manhattan? A Study of Alternative Vlsi Routing Architectures
,”
Proceedings of the 10th Great Lakes Symposium on VLSI, GLSVLSI ’00, Association for Computing Machinery
,
Chicago, IL
,
Association for Computer Machinery
,
New York
, pp.
47
52
.
14.
Van der Velden
,
C.
,
Bil
,
C.
,
Yu
,
X.
, and
Smith
,
A.
,
2007
, “
An Intelligent System for Automatic Layout Routing in Aerospace Design
,”
Innovat. Syst. Softw. Engin.
,
3
(
2
), pp.
117
128
. 10.1007/s11334-007-0021-4
15.
Park
,
J.-H.
, and
Storch
,
R.
,
2002
, “
Pipe-routing Algorithm Development: Case Study of a Ship Engine Room Design
,”
Expert Syst. Appl. (UK)
,
23
(
3
), pp.
299
309
. 10.1016/S0957-4174(02)00049-0
16.
Guirardello
,
R.
, and
Swaney
,
R. E.
,
2005
, “
Optimization of Process Plant Layout with Pipe Routing
,”
Comput. Chem. Engin.
,
30
(
1
), pp.
99
114
. 10.1016/j.compchemeng.2005.08.009
17.
Liu
,
C.
,
2018
, “
Optimal Design of High-rise Building Wiring Based on Ant Colony Optimization
,”
Cluster Comput.
,
22
(
2
), pp.
1
8
.
18.
Betz
,
V.
, and
Rose
,
J.
,
1997
, “Vpr: a New Packing, Placement and Routing Tool for Fpga Research,”
Field-Programmable Logic and Applications
,
Luk
,
W.
,
Cheung
,
P. Y. K.
,
Glesner
,
M.
, eds.,
Springer Berlin Heidelberg
, pp.
213
222
.
19.
Tisdale
,
J.
,
Kim
,
Z.
, and
Hedrick
,
J. K.
,
2009
, “
Autonomous Uav Path Planning and Estimation
,”
IEEE Robot. Automat. Magaz.
,
16
(
2
), pp.
35
42
. 10.1109/MRA.2009.932529
20.
Jan
,
G. E.
,
Yin Chang
,
K.
, and
Parberry
,
I.
,
2008
, “
Optimal Path Planning for Mobile Robot Navigation
,”
IEEE/ASME Trans. Mechatron.
,
13
(
4
), pp.
451
460
. 10.1109/TMECH.2008.2000822
21.
Landon
,
M. D.
, and
Balling
,
R. J.
,
1994
, “
Optimal Packaging of Complex Parametric Solids According to Mass Property Criteria
,”
ASME J. Mech. Des.
,
116
(
2
), pp.
375
381
. 10.1115/1.2919389
22.
Szykman
,
S.
, and
Cagan
,
J.
,
1997
, “
Constrained Three-Dimensional Component Layout Using Simulated Annealing
,”
ASME J. Mech. Des.
,
119
(
1
), pp.
28
35
. 10.1115/1.2828785
23.
Szykman
,
S.
,
Cagan
,
J.
, and
Weisser
,
P.
,
1998
, “
An Integrated Approach to Optimal Three Dimensional Layout and Routing
,”
ASME J. Mech. Des.
,
120
(
3
), pp.
510
512
. 10.1115/1.2829180
24.
,
C.
,
Cagan
,
J.
, and
,
K.
,
2006
, “
Objective Function Effect Based Pattern Search–Theoretical Framework Inspired by 3D Component Layout
,”
ASME J. Mech. Des.
,
129
(
3
), pp.
243
254
. 10.1115/1.2406095
25.
Yin
,
S.
, and
Cagan
,
J.
,
2000
, “
An Extended Pattern Search Algorithm for Three-Dimensional Component Layout
,”
ASME J. Mech. Des.
,
122
(
1
), pp.
102
108
. 10.1115/1.533550
26.
Ren
,
T.
,
Zhu
,
Z.-L.
,
Dimirovski
,
G.
,
Gao
,
Z.-H.
,
Sun
,
X.-H.
, and
Yu
,
H.
,
2014
, “
A New Pipe Routing Method for Aero-engines Based on Genetic Algorithm
,”
Proc. Inst. Mech. Engin. Part G (J. Aerospace Eng.)
,
228
(
3
), pp.
424
434
. 10.1177/0954410012474134
27.
Qu
,
Y.
,
Jiang
,
D.
,
Gao
,
G.
, and
Huo
,
Y.
,
2016
, “
Pipe Routing Approach for Aircraft Engines Based on Ant Colony Optimization
,”
J. Aerospace Eng.
,
29
(
3
), p.
04015057
. 10.1061/(ASCE)AS.1943-5525.0000543
28.
Gulić
,
M.
, and
Jakobović
,
D.
,
2013
, “
Evolution of Vehicle Routing Problem Heuristics with Genetic Programming
,”
2013 36th International Convention on Information and Communication Technology
,
Opatija, Croatia
,
Electronics and Microelectronics (MIPRO)
, pp.
988
992
.
29.
,
S. R. T.
,
James
,
K. A.
, and
Allison
,
J. T.
,
2020
, “
A Novel Two-stage Design Framework for 2D Spatial Packing of Interconnected Components
,”
ASME 2020 International Design Engineering Technical Conferences, No. IDETC2020-22695
,
Online
,
Aug. 17–19 [Virtual]
.
30.
Jessee
,
A.
,
,
S. R. T.
,
Lohan
,
D. J.
,
Allison
,
J. T.
, and
James
,
K. A.
,
2020
, “
Simultaneous Packing and Routing Optimization Using Geometric Projection
,”
ASME J. Mech. Des.
,
142
(
11
), p.
111702
. 10.1115/1.4046809
31.
,
S. R. T.
,
Herber
,
D. R.
,
Pangborn
,
H. C.
,
Alleyne
,
A. G.
, and
Allison
,
J. T.
,
2019
, “
Optimal Flow Control and Single Split Architecture Exploration for Fluid-Based Thermal Management
,”
ASME J. Mech. Des.
,
141
(
8
), p.
083401
. 10.1115/1.4043203
32.
,
S. R. T.
,
Rodriguez
,
S. B.
,
James
,
K. A.
, and
Allison
,
J. T.
,
2020
, “
Automated Layout Generation Methods for 2D Spatial Packing
,”
ASME 2020 International Design Engineering Technical Conferences, No. IDETC2020-22627
,
Virtual
,
Aug. 17–19
,
ASME
.
33.
Tutte
,
W. T.
,
1963
, “
How to Draw a Graph
,”
Proc. Lond. Math. Soc.
,
s3–13
(
1
), pp.
743
767
. 10.1112/plms/s3-13.1.743
34.
,
T.
, and
Kawai
,
S.
,
1989
, “
An Algorithm for Drawing General Undirected Graphs
,”
Inform. Process. Lett.
,
31
(
1
), pp.
7
15
. 10.1016/0020-0190(89)90102-6
35.
Fruchterman
,
T. M. J.
, and
Reingold
,
E. M.
,
1991
, “
Graph Drawing by Force-Directed Placement
,”
Softw. Practice Exp.
,
21
(
11
), pp.
1129
1164
. 10.1002/spe.4380211102
36.
Chaari
,
I.
,
Koubaa
,
A.
,
Bennaceur
,
H.
,
Ammar
,
A.
,
Alajlan
,
M.
, and
Youssef
,
H.
,
2017
, “
Design and Performance Analysis of Global Path Planning Techniques for Autonomous Mobile Robots in Grid Environments
,”
,
14
(
2
), p.
1729881416663663
. 10.1177/1729881416663663
37.
Alexopoulos
,
C.
, and
Griffin
,
P. M.
,
1992
, “
Path Planning for a Mobile Robot
,”
IEEE Trans. Syst. Man Cyber.
,
22
(
2
), pp.
318
322
. 10.1109/21.148404
38.
Hoffman
,
K. L.
,
,
M.
, and
Rinaldi
,
G.
,
2013
,
Traveling Salesman Problem
,
Springer US
,
Boston, MA
, pp.
1573
1578
.
39.
Guruji
,
A. K.
,
Agarwal
,
H.
, and
Parsediya
,
D.
,
2016
, “
Time-Efficient a* Algorithm for Robot Path Planning
,”
Proc. Technol.
,
23
, pp.
144
149
.
3rd International Conference on Innovations in Automation and Mechatronics Engineering 2016, February, ICIAME 2016 05-06
.
40.
Norato
,
J.
,
Bell
,
B.
, and
Tortorelli
,
D.
,
2015
, “
A Geometry Projection Method for Continuum-Based Topology Optimization With Discrete Elements
,”
Comput. Methods Appl. Mech. Engin.
,
293
, pp.
306
327
. 10.1016/j.cma.2015.05.005
41.
Hughes
,
T. J. R.
,
2000
,
The Finite Element Method
,
Dover Publications
.