Saeid BayatDepartment of Industrial and Enterprise Systems Engineering,
University of Illinois at Urbana-Champaign,
Urbana, IL, USA
email<EMAIL_ADDRESS>
Nastaran Shahmansouri Autodesk Research,
661 University Ave,
Toronto, Canada
email<EMAIL_ADDRESS>
Satya RT Peddada Department of Industrial and Enterprise Systems Engineering,
University of Illinois at Urbana-Champaign,
Urbana, IL, USA
email<EMAIL_ADDRESS>
Alexander Tessier Autodesk Research,
661 University Ave, Toronto, Canada
email<EMAIL_ADDRESS>
Adrian Butscher Autodesk Research,
661 University Ave,
Toronto, Canada
email<EMAIL_ADDRESS>
James T Allison Department of Industrial and Enterprise Systems Engineering,
University of Illinois at Urbana-Champaign,
Urbana, IL, USA
email<EMAIL_ADDRESS>
# Multi-split configuration design for fluid-based thermal management systems

(Version , )

###### Abstract

High power density systems require efficient cooling to maintain their thermal
performance. Despite this, as systems get larger and more complex, human
practice and insight may not suffice to determine the desired thermal
management system designs. To this end, a framework for automatic architecture
exploration is presented in this article for a class of single-phase, multi-
split cooling systems. For this class of systems, heat generation devices are
clustered based on their spatial information, and flow-split are added only
when required and at the location of heat devices. To generate different
architectures, candidate architectures are represented as graphs. From these
graphs, dynamic physics models are created automatically using a graph-based
thermal modeling framework. Then, an optimal fluid flow distribution problem
is solved by addressing temperature constraints in the presence of exogenous
heat loads to achieve optimal performance. The focus in this work is on the
design of general multi-split heat management systems. The architectures
discussed here can be used for various applications in the domain of
configuration design. The multi-split algorithm can produce configurations
where splitting can occur at any of the vertices. The results presented
include 3 categories of cases and are discussed in detail.

###### keywords:

Design Synthesis, Graph Generation, Optimization, Optimal Flow Control,
Thermal Management System

##  1 Introduction

Electrification of many systems in various domains has increased over the last
few decades [1, 2, 3]. These electronic devices are experiencing intensifying
miniaturization [4, 5, 6, 7, 8, 9], resulting in higher power densities.
Currently, microprocessors generate heat fluxes of over 102 W/cm2. Hot spots
on microelectronic chips can also generate heat fluxes of 1 kW/cm2 or more,
which result in excessive temperatures in local regions [4]. When devices are
operated at excessive temperatures, their performance and reliability are
negatively affected, resulting in their malfunction in the end [10, 11].
Consequently, the optimal design of thermal management systems is essential
for expediting design processes and achieving ambitious performance goals [12,
13].

Throughout the literature, the design of heat management systems has largely
focused on improving individual components [14, 15]. Yet, optimizing
individual components can result in a sub-optimal overall design.
Additionally, the design of a whole heat management system sometimes is only
relevant to a particular application [16]. Peddada et al. [14] took an
important step towards designing a class of heat management systems that are
not limited to a specific application. However, the systems considered in that
work are restricted to single-split system architectures, i.e., configurations
with a single source where thermal-devices of the system are located in
branches that only start from the source. There are, nevertheless, many
instances in which multi-split systems are required [17]. In a multi-split
architecture, branching can occur from the main source or any thermal-device
of the system. In this paper, we introduce a strategy for the generation and
the optimal design of general multi-split heat management systems using
graphs. Here, spatial considerations are also incorporated into the design
problem; the spatial data of the heat-devices of the system are used for
clustering and defining junctions where branching starts. After forming the
clusters, for each cluster, we find the Euclidean distance of cluster’s heat-
devices from the cluster’s centroid; the heat-device closest to the centroid
is defined as the cluster’s junction, refer to section 4.2 for more details.
The architectures discussed in the paper can be used across a wide range of
applications. Additionally, the design strategy presented here applies to both
single-split and multi-split architectures, supporting quantification of
trade-offs between multi-split system performance improvement and cost and
complexity increases.

Designs of these systems should meet transient response conditions [18] as
many electronic devices work in applications with time-varying workloads [4].
Accurate modeling of the design problem as an optimization problem, therefore,
requires treatment of system dynamics, including bounds and constraints on
dynamic behavior, as well as optimal control for active systems. Continuous
optimization can be applied to the combined design of physical and control
system design, a well-established problem in Control Co-Design (CCD) [19, 20,
21, 22], for a given system configuration. Distinct configurations can have
fundamentally different dynamics, design variables, and constraints, so each
configuration must be treated independently. Practical solution of such
problems requires automatic generation and solution of the optimization
problem for each configuration. Here, for every generated unique and feasible
multi-split system configuration a continuous optimization problem is
formulated and solved [23, 24, 14].

Thermal management system configuration design considered in this paper
belongs to an especially challenging class of optimization problems where
discrete decisions change the set of continuous decisions to be made. This
necessitates a nested discrete/continuous approach. Furthermore, the very
general nature of these problems prevent the use of established integer
programming methods that can be used for efficient solution of other problem
classes with spatial properties. At least, there are three possible approaches
for solving general physical systems configuration design problems. The first
approach, as described above and used in the studies presented here, entails
enumerating all unique and feasible configurations for a given design space,
solving the continuous optimization problem for each one, and then producing a
ranked set of design candidates. This is the only approach that can produce a
result that is a confirmed optimum. A second approach employs a population-
based optimization algorithm, such as a Genetic Algorithm (GA), or other
gradient-free search strategy to navigate the configuration space while still
solving the continuous optimization problem for each configuration. A third
approach is to utilize machine learning or artificial intelligence techniques
to restructure the system configuration design space [25, 26] such that it is
more tractable, while again nesting continuous optimization within this
search. Recent machine learning strategies have proven to be only marginally
better than random search [26]. General solution approaches for configuration
problems too large for enumeration and without special problem structure
remains elusive, and is truly an engineering design grand challenge.

The main contributions of the work presented in this article are:

  1. 1

Introduction of a new automated modeling strategy for multi-split thermal
management systems. This strategy has applications for configuration
generation across a wide range of systems.

  2. 2

Employing spatial information to cluster data and to define junctions.

  3. 3

Comparison of multi-split versus single-split configurations, providing
insight into the trade-offs between system performance and system complexity

  4. 4

In-depth study of multi-split system optimal results and analysis of the
system signals such as flow rates and flow temperatures at heat-device
locations that led to this optimal solution

This article continues as follows, in Section 2, we discuss the thermal
management system architectures studied in this work. The dynamic graph based
models of the thermal architectures are presented in Section 3. Section 4
describes the graph-based representation of multi-split architectures. The
formulation of the variable time horizon dynamic optimization problem is
explained in Section 5. Section 6 presents case studies using different
architectures and heat loads. Section 7 concludes with a summary of the design
methodology, guidelines for thermal management system design, and potential
future research topics.

##  2 System Description and Modeling

Figure 1: Class of problems considered in this paper. The systems include a
tank, a pump, valve(s), CPHXs in parallel and series, a LLHX, and a sink.

Figure 1 illustrates the class of problems considered in this article. The
purpose of this system is to manage the temperature of various heat generating
devices mounted on Cold Plate Heat Exchangers (CPHXs) through which a coolant
flows. The coolant is stored in a tank and is transferred to each branch by a
pump. Each branch has several valves and can divide the flow it receives into
its sub-branches. The coolant that passes through heat exchangers absorbs heat
and transfers it to the thermal sink through a Liquid-to-Liquid Heat Exchanger
(LLHX).

The class of architectures considered in this work produces multi-split
configurations as opposed to the single-split configurations generated in the
study by Peddada et al. [14]. A single-split assumption limits the
configuration search space; expanding the search space to include multi-split
configurations may enhance system performance. The multi-split enumeration
algorithm created in this work produces configurations where splits may be
made at the pump (source) or at any of the CPHX locations. Sections 3 and 4
describe the graph-based physical modeling of the thermal systems and the
representation and generation of multi-split architectures. Here the heat load
applied to each CPHX, as well as the inlet temperature and mass flow rate of
the thermal sink, have been assumed to be known exogenous inputs. In addition,
heat generating devices are assumed to have the same temperature as the wall
of the CPHX on which they are mounted. It is assumed that heat loss through
pipes is zero.

An optimal control problem is defined for each fluid-based thermal management
systems configuration generated using the multi-split algorithm. The control
problem seeks the optimal flow rate trajectory for each pipe that maximizes
system performance, while satisfying component temperature constraints. The
flow rates are controlled by valves. Dynamic system models incorporate the
thermal physics of advection, convection, and bi-directional advection. A
unique model is generated for each configuration design. We assume that system
performance is quantified by thermal endurance, i.e., the goal is to maximize
the time that the device is on while ensuring that all temperature bounds are
met. This is consistent with the approach in the study conducted by Peddada et
al. [14], supporting direct comparison with earlier single-split studies using
graph-based configuration representations.

The code workflow is illustrated in Fig. 2, providing an overview of the code
structure. In the subsequent sections, each part will be discussed in detail.
The variable Data represents the positions of CPHX in [$x,y,z$] coordinates.
The parameter numLevels determines the depth of the graph considered for
junction creation. For example, when set to 1, junctions can only be added to
branches connected to the tank. Multiple graphs are generated based on the
provided Data and numLevels, resulting in different configurations, such as
parallel, series, or a combination thereof. Among a set of $N$ graphs, the
variable configNum specifies a particular graph. The variable distrb denotes
the heat load at each node. Using the given Data, Heat Load, and Graph Config,
the Base Graph is generated, representing the connections between the tank,
junctions (denoted as $j_{2}$ and $j_{3}$ in Fig.2), and other CPHX nodes.
Subsequently, a Physics Graph is created by adding CPHX wall nodes, sink
nodes, and source power (heat load) to the Base Graph, along with the
underlying physics between these nodes, which include convection, advection,
and bidirectional advection. Based on the Physics Graph, an Optimal Open Loop
Control (OLOC) is defined, consisting of three components: system dynamics
($\mathbf{f}$), objective function ($J$), and path constraints ($G$). Here,
$t$ represents time, $\mathbf{\xi}$ denotes the state, and $\mathbf{u}$
represents the control signal. These components will be discussed in detail in
subsequent sections.

Figure 2: The workflow diagram that illustrates the steps involved in the code
execution. Starting with the given Data, Head Load, and Config Number, the
base graph and physics graph are generated. The OLOC problem is then defined
and solved to obtain the objective function value. The objective value
represents thermal endurance, with the incorporation of a penalty function for
control signal smoothness and convergence. Refer to Sec. 5 for a comprehensive
discussion

##  3 Dynamic Graph-Based Modeling

This paper uses the dynamic graph-based modeling framework discussed in
Peddada et al. [14]. Nodes represent components or fluids, and a temperature
is associated with each node. Edges represent power-flows between nodes. For
brevity, model details are omitted here, but the dynamics governing the main
components are presented in Fig. 3. The model includes advection, convection,
and bidirectional advection. Convection happens between CPHX wall nodes and
fluid nodes, advection occurs between fluid nodes, and bidirectional advection
occurs between the LLHX and sink. In Fig. 3, $c^{p}$ is the specific heat
capacitance, $\rho$ is the density of the fluid, $A_{s}$ is the convective
surface area, $h$ is the heat transfer coefficient, and $\dot{m}$ is the mass
flow rate of the fluid.

Figure 3: Notional example to illustrate the elements of thermal physics
included in this system model. Here, a node represents a temperature and an
edge represents a power flow.

###  3.1 Heat-transfer Model of Graph Nodes

The power-flow type in the tank is advection, which occurs between 1) the tank
and the fluid node on the CPHX, and 2) the tank and the primary side of the
LLHX. Advection and convection are present in the CPHX. Advection occurs
between 1) the CPHX fluid node and tank and 2) the CPHX fluid node and primary
side of the LLHX. Convection occurs between the CPHX fluid and the wall node.
The LLHX involves advection, convection, and bi-directional advection.
Advection occurs between 1) the primary side of the LLHX and tank, 2) the
primary side of the LLHX and the CPHX fluid node. Bi-directional advection
takes place between the secondary side of the LLHX and the sink. Convection
occurs between 1) the LLHX wall node and the primary side of LLHX and 2) the
LLHX wall node and the secondary side of LLHX. These are building blocks of
any larger system. Using these components, a model of the dynamics of any
complex system of the form illustrated in Figure 1 can be generated.

###  3.2 Graph-Based Model for Multi-split Architectures

Figure 4: Graph-based model for multi-split architectures studied in this
article. Figure 5: Graph illustrating the variables for a simple example
multi-split configuration. Types of power flow are illustrated, including the
distinction between independent and dependent flows.

Figure 4 shows a graph that represent the general class of multi-split
architectures studied in this paper. Here, $P^{s}$ represents the heat load. A
simple example problem of this class is shown in Fig. 5. Each of the two
junctions has 3 associated nodes. The first junction ($j_{2}$) distributes the
flow to 3 different branches, and the second junction ($j_{3}$) distributes
the flow to two branches, where 2 nodes in a branch are in series, and the
second branch has only 1 node. In this figure, states are denoted using the
symbol $\bm{\xi}$; here the states all correspond to node temperatures. For
example, $\xi_{1}$ is the temperature of node 1 and $\xi_{w-1}$ is the wall
temperature of the CPHX connected to node 1. The system involves both
independent and dependant flows, depicted in Fig. 5. Independent flows can be
controlled, whereas dependent flows can be calculated based on the independent
flows (see Fig. 5 for these equations). In the figure, the edges that carry
dependent power flows are shown with dashed-lines.

###  3.3 State Equations for the Graph-Based Model

The state equations for the thermal management systems are obtained by using a
methodology similar to the one introduced in references [14, 27]. The
methodology described there has been extended to accommodate multi-split
configurations. Equation (1) shows the state equation for this system. Here
$\mathbf{T}(t)$ is the vector of node temperatures (states), $T^{t}$
represents the sink temperature, which is considered to be a constant value
and is known, $\dot{m}_{p}(t)$ is the pump flow rate, $\dot{m}_{t}(t)$ is the
sink flow rate, $\dot{\mathbf{m}}_{f}(t)$ is composed of the sub-branch flow
rates, and $\mathbf{P}^{s}(t)$ is the vector of component heat loads. Here,
$\mathbf{C}$ is the diagonal matrix of heat capacitance, and $\mathbf{D}$
represents the connection of external sources to the system. Matrix
$\mathbf{A}$ is calculated using 3 matrices: 1) $\mathbf{\overline{M}}$, a
matrix derived from the graph incidence matrix representing the structural
mapping from power flows to the states, 2) a weighted incidence matrix
corresponding to convective power flow, and 3) $\mathbf{C}$. Matrix
$\mathbf{B}_{1}$ is obtained from $\mathbf{\overline{M}}$ and $\mathbf{C}$.
$\mathbf{B}_{2}$ is obtained from a weighted incidence matrix corresponding to
advective and bidirectional advective power flows. For a detailed explanation
of these matrices, please refer to reference [14].

$\displaystyle\dot{\bm{T}}(t)=\bm{A}\begin{bmatrix}\bm{T}(t)\\\
T^{t}(t)\end{bmatrix}+\bm{B}_{1}\left(diag\left(\bm{Z}\begin{bmatrix}\dot{m}_{p}(t)\\\
\dot{\bm{m}}_{f}(t)\\\
\dot{m}_{t}(t)\end{bmatrix}\right)\right)\bm{B}_{2}\begin{bmatrix}\bm{T}(t)\\\
T^{t}(t)\end{bmatrix}+\bm{C}^{-1}\bm{D}\bm{P}^{s}(t)$ (1)

##  4 Generating Multi-split Spatial Graphs

###  4.1 Graph Representation

The system configurations here can be represented by acyclic undirected
connected graphs where any two nodes are connected by only one simple path (
i.e. a tree). In the current representation, the tank is always the root node
and each CPHX is labeled with a number from 1 to N. Splittings occur at the
root or at any other node. In the corresponding symbolic representation,
branching is shown by parentheses and consequent CPHXs, i.e., nodes in one
branch are separated by a comma. Figure 6 shows the representation for two
configurations and their equivalent graphs.

(a) $0\,(1,2)(3)$ (b) $0\left(1,2(3,4)(5)\right)(7)\left(8(9,10)\right)(11)$

Figure 6: Two multi-split representations and their equivalent graphs.

The number of single-split configurations increases rapidly as the number of
heat-generating devices increases. Figure 7(a) shows the total number of
single-split configurations as a function of number of nodes in graphs.
Equation (2) quantifies the number of generated configurations for this
enumeration. When splits (or junctions) are added, the total number of
configurations increases more. Assuming there is only one-layer of junction
nodes, the number of configurations depends on the number of junctions as well
as the number of non-junction nodes. One can calculate the number of
configurations for a fixed number of non-junction nodes ($N$) and a variable
number of junctions, ranging from 1 to $J$, in a recursive manner. Table 1
summarizes the equations.

$\displaystyle G(N)=\sum_{k=0}^{N}{N\choose k}{N-1\choose k-1}(N-k)!$ (2)

Table 1: Number of graphs for a given number of junctions Junctions | Number of Graphs
---|---
J = 1 | $F_{1}(N)=G(N)$
J = 2 | $F_{2}(N)=\sum_{M=1}^{N}{N\choose M}G(M)F_{1}(N-M)$
J = J | $F_{J}(N)=\sum_{M=1}^{N}{N\choose M}G(M)F_{J-1}(N-M)$

Figure 7(b) presents the number of configurations generated for graphs with N
= 5, 10, and 15 when J varies from 1 to N. Figure 7(c) shows the number of
graphs as a function of N and J. To find the total number of graphs with one
junction layer (J = 1 to N) , we need to sum up all the graphs generated for
different number of junctions; see Fig. 7(d) for an illustration of this
growth.

(a)  (b)  (c)  (d)

Figure 7:  (a) Number of single-split graphs vs. number of non-root nodes, (b)
The number of multi-split graphs for graphs with $N=$ 5, 10, and 15 nodes when
the number of junctions varies from 1 to $N$, (c) 3D surface illustrating the
number of graphs as a function of the number of vertices and junctions, and
(d) the total number of multi-split graphs with only one junction layer as a
function of $N$.

The addition of junction nodes and layers leads to a rapid growth in the
number of graphs. It is desirable to impose constraints in a way that
meaningfully limits this growth. One strategy is to recognize that these
graphs represent physical designs, and that spatial system information can be
leveraged to reduce problem complexity. Here we employ a spatial clustering
method to group the nodes in sub-domains and enumerate the sub-branches only
in the nearby neighborhoods; this procedure is detailed in Section 4.2. This
approach, however, cannot access all possible designs (as can enumeration).
Therefore, we examine a second strategy where the junction locations are
enumerated. Exploring the trade-offs between performance improvement and
reducing problem complexity is explored in this article as in some cases
enumeration may be worth the added computational cost.

###  4.2 Generation of Multi-split Spatial Graphs

Algorithms 1– 4 describe the multi-split graph generation framework of this
work. The generation procedure passively choose junction nodes during run-time
instead of employing all the nodes in enumerations from the beginning. It
performs recursion for refinement. The enumerations can cover (i) tank-
junctions, (ii) junction-CPHXs, and (iii) both, based on the specific
requirements of the optimization problem; this study focuses on class (ii).

Figure 8: Clustering and the creation of a sample configuration for a set of
nodes. Figure shows (a) the spatial locations of source (T) and 17 nodes, (b)
first-level clusters, (c) centroid of each cluster, (d) chosen junctions and a
routing from T to junctions, (e) second-level clusters, (f) the centroid and
the corresponding junctions of the second-level, (g) a routing from first-
level to second-level junctions, (h) a sample enumeration in the internal
cluster with more than 1 node. Figure (i) presents the generated sample graph.
Algorithm 1 Selecting Number of Clusters to form Super-Nodes

1:$N\leftarrow$ number of nodes

2:for $K$ in range(1, $N$) do

3: cluster nodes into $K$ clusters via K-means

4: if the clustering is stable, stop

5:return $K$ and the corresponding nodes of each cluster

Algorithm 2 Generation of Hierarchical Tree of Super-Nodes

1:$T\leftarrow$ tree

2:super-node[0]$\leftarrow$ $\\{P[0],D[0]\\}$ here $P[0]$ is Tanks and $D[0]$
is the spatial data of all nodes

3:$L[0]$ includes super-node[0]

4:$numL\leftarrow$ number of levels in the tree

5:$m\leftarrow$ number of super-nodes in $T$

6:for $l$ in range(1, $numL$) do

7: $Ids\leftarrow$ index of all super-nodes in $L[l-1]$

8: for $i$ in $Ids$ do

9: cluster super-node $i$ into $K[i]$ clusters using Algorithm 1

10: for $k$ in range(1,K[i]) do

11: $m=m+1$

12: $P_{k}\leftarrow$ the closest node to centroid of $D(k)$ ($P_{k}$ is
junction.)

13: $P[i]\leftarrow$ $[P_{1},P_{2},...,P_{n},P_{k}]$

14: remove $P_{k}$ from $D(k)$

15: $T[m]\leftarrow$ $\\{P[k],D[k]\\}$

16: $L[l]\leftarrow$ add $T[m]$

17:return Tree T with numL levels

Algorithm 3 Graph Enumeration

1:$N\leftarrow$ number of nodes

2:$eG[0]\leftarrow[[0]]$ $\triangleright$ each sub-list represents the
adjacency list of a graph

3:$eG[1]\leftarrow[[(0,1)]]$

4:for $n$ in [2, $N$] do

5: $Parents[n]\leftarrow[1,\ldots,n-1]$

6: for $g$ in $eG[n-1]$ do

7: for each $Node$ in $g$ do

8: if $Node$ belongs to $Parents[n]$ then

9: add edge ($Node$, $n$) to $g$ and add the adjacency list of the new graph
to $eG[n]$

10:return all enumerated graphs $eG[N]$

Algorithm 4 Graph Generation for a Selected Tree Level

1:$l\leftarrow$ Tree level selected for graph generation

2:$Ids\leftarrow$ index of all super-nodes in $L[l]$ for tree generated by
Algorithm 2

3:for $i$ in $Ids$ do

4: $Q\leftarrow$ parents of the super-node ($P[i]$)

5: $g1[1:j]$ $\leftarrow$ enumerate super-node$[i]$ using Algorithm 3 with its
junction as its root

6: $g2$ $\leftarrow$ the circular graph of $Q$ with $(P[-1],P[0])$ edge
removed

7: $g[i,1:j]$ $\leftarrow$ merge $g2$ with every $g1[1:j]$

8:combine a sub-graph selected from each super-nodes of level l to generate
all the graphs

Algorithm 1 initially creates trees using nodes spatial data. The nodes data
are clustered to form super-nodes in a recursive manner. Here, a level
represents the levels of splitting in an architecture tree. A junction node is
defined for each super-node by choosing the node with the smallest Euclidean
distance from the cluster’s centroid, see Figure 8 and Algorithm 2. Note that
a junction represents a CPHX node where the coolant flow splits into branches.
Next, we enumerate all sub-graphs of a selected tree-level. This means, for a
chosen super-node in that tree level, its nodes are enumerated to form all
possible sub-graphs; Algorithm 4 describes how the enumeration algorithm
works. Afterwards, we connect the root node (Tank) and the corresponding
junction nodes of the super-node to its sub-graphs, see Algorithm 3. This
process is performed for all super-nodes of the selected level. Finally, an
architecture graph is created by choosing and merging one sub-graph from the
pool of sub-graphs generated for each super-node in that level, refer to
Algorithm 4. Note that various sub-clusters can be defined within a cluster to
form multiple layers of junction nodes.

In addition to the above mentioned graph generation algorithm, we also
generated graphs using a second strategy where the locations of junctions in
the trees are enumerated. The main difference between the first and second
graph generation strategies is in topological locations of junction nodes in
the architectural graphs.

Prior research performed by some of the authors on holistic design for 3D
spatial packaging and routing of interconnected systems [28, 29, 30, 31, 32]
will be employed in future work to capture the spatial aspects of the multi-
split configurations via optimal placement of junction nodes and cooling
circuit components such as the CPHXs, pumps, valves, and the tank with
simultaneous optimization of 3D lengths of the branches (or pipe segments)
while satisfying volume and multi-physics constraints. This can lead to a more
complete definition of the real-world multi-split configuration fluid-based
thermal system problem.

We generate all configurations within a class, evaluate the performance of
each configuration, and choose the optimal one. Figure 9 shows all 13
configurations made by the single-split algorithm for 3 nodes. The figure
depicts the base graphs and the extended physics graphs generated for
simulations. Figure 10 shows two examples configurations and their physics
models for a multi-split system with 19 nodes (root, 6 junction CPHXs, and 12
CPHXs nodes).

Figure 9: (a) Base and (b) complete physics graphs for single-split systems
with 3 CPHXs Figure 10: Two sample configurations for a multi-split system
with 18 CPHXs. Figures (a) and (b) show the base graphs, and Figs. (c) and (d)
show the complete physics graphs for the base graphs (a) and (b),
respectively.

##  5 Optimal Flow Control Problem (OLOC)

The objective of the optimal control problem is to maximize the thermal
endurance while satisfying limits on temperatures and mass flow rates. When
any of the node’s temperatures (defined as states) reaches the upper bound,
the OLOC terminates, and the final time is recorded as the thermal endurance.
Figure 11 shows the OLOC structure. Here, we seek to determine a control
trajectory $u(t)$ for each independent flow that maximizes the objective
function while satisfying constraints. Within this figure, Eq. ($c_{1}$)
presents the system states ($\bm{\xi}$), encompassing the vector of
temperature nodes ($\bm{T}$), and the flow rate of independent branches
($\bm{\dot{m}}_{\mathrm{indp}}$). The dynamics that shows how these states
evolve over time are shown in Eq. ($c_{8}$). Equation ($c_{2}$) defines the
control signals as the rate of change of the valve flow rates in independent
branches. The flow rate in dependent branches is determined by an algebraic
equation that ensures input flow rate of each branch is equal to the output
flow rate of that branch. This equation can be represented as a matrix
multiplication, as shown in Eq. ($c_{3}$), where matrix $M$ maps the flow rate
of independent branches ($\bm{\dot{m}}_{\mathrm{indp}}$) to the flow rate of
dependent branches ($\bm{\dot{m}}_{\mathrm{dp}}$). A practical illustration of
this calculation is provided in Figure 5. In Eq. ($c_{4}$), the total flow
rate of all branches is presented, comprising both independent and dependent
flow rates. It is worth mentioning that the order of combining independent and
dependent flow rates depends on the graph structure. For simplicity, in this
case, we first incorporate the independent flow rates and then the dependent
flow rates. The initial conditions for the temperature nodes are given by Eq.
($c_{5}$). In this equation, $\bm{T}_{w}$ represents the wall temperature of
the CPHX states, $\bm{T}_{f}$ represents the fluid temperature of the CPHX
states, and $\bm{T}_{l}$ represents the temperature of the Tank and LLHX
states. Eq. ($c_{6}$) presents inequality path constraints that ensure the
operating temperature of each component remains within an upper bound
throughout the entire time horizon. The first two terms in Eq. ($c_{7}$)
represent the inequality path constraint that guarantees the flow rates of
both independent and dependent branches remain within a specified bound. This
equation, when combined with Eq. ($c_{3}$), ensures that the input flow rate
of each branch is equal to its output flow rate. It also ensures that the
maximum flow rate cannot exceed the pump flow rate ($\dot{m}_{p}$). The last
term in Eq. ($c_{7}$) represents the limit on the derivative of the flow rate,
capturing the physical limitations of the valves [14]. Equations $c_{8}$ and
$c_{9}$ show how states defined in Eq. ($c_{1}$) evolve over time. In Equation
($c_{10}$), the objective value is represented, aiming to maximize the thermal
endurance. Similar to the approach in Ref. [14], a penalty term is
incorporated to facilitate solution smoothness and enhance convergence. The
parameter $\lambda$ is selected such that the total penalty cost remains below
$1\%$ of $t_{\mathrm{end}}$. Table 2 shows the parameters used in the physical
simulations for the studies in this article.

Figure 11: OLOC components Table 2: Parameters used in the physics modeling of
the thermal systems [14].

Parameter | Value
---|---
LLHX wall mass | 1.2 kg
CPHX wall mass | 1.15 kg
Tank fluid mass | 2.01 kg
Thermal sink temperature $T^{t}$ |  $15^{o}$ C
Tank/LLHX initial temperature , $T_{l,0}$ |  $15^{o}$ C
CPHX initall wall temperature, $T_{w,0}$ |  $20^{o}$ C
CPHX inital fluid temperatue , $T_{f,0}$ |  $20^{o}$ C
Thermal sink mass flow rate, $\dot{m}_{t}$ | 0.2 kg/s
Pupmp mass flow rate $\dot{m}_{p}$ | 0.4 kg/s
Valve rate limit $\ddot{m}_{f,\mathrm{max}}$ |  0.05 kg/$\mathrm{s}^{2}$
Penalty parameter, $\lambda$ | $0.01/(N_{f}{\ddot{\bm{m}}}_{\mathrm{f,max}}^{2})$

There are two main approaches for solving OLOC problems: 1) Indirect (optimize
then discretize) and 2) Direct (discretize then optimize). In the indirect
method, a differential algebraic equation is derived using optimality
conditions (the calculus of variations or the Pontryagin minimum principle).
These equations should then be discretized and solved numerically [33]. In
contrast, the direct method first discretizes the problem so it can be
transformed into a nonlinear program (NLP), which can then be solved by a
nonlinear programming solver such as SNOPT [34] or IPOPT [35]. The indirect
method provides more information about the structure of the problem, but
solving constrained problem can be challenging with this method.
Alternatively, direct methods use NLP solvers to solve complex problems
successfully; some well-established OLOC software tools based on the direct
method are available, for example GPOPS [36] and Dymos [37]. This paper uses
Dymos to solve the OLOC problems, which is an open source program developed in
Python. The computation cost of solving the OLOC problem depends on the size
of the graph, but on average, it takes approximately two minutes to solve the
problem for each configuration. Additionally, it is important to highlight
that in this article, the evaluation of the nonlinear optimal control problems
for each of the architectures has been parallelized. This parallelization
approach significantly reduces the computational cost associated with solving
these problems. All reported computational costs were obtained using a
workstation with an AMD EPYC 7502 32-Core Processor @ 2.5 GHz, 64 GB DDR4-3200
RAM, LINUX Ubuntu 20.04.1, and Python 3.8.10.

##  6 Case Studies

We present three case studies to illustrate how this work can help engineers
design optimal thermal management systems. In the first case study (Sect.
6.1), the goal is to obtain the optimal structure for thermal management
systems having 3 and 4 CPHXs. Here we use enumeration and compare both single-
split and multi-split cases. For 3 and 4 CPHXs systems, the heat loads are
$[12,4,1]$ kW and $[12,4,1,1]$ kW, respectively. In Sect. 6.2, we have shown
the results for the multiple-split case under two disturbances:
$[5,5,5,5,5,5]$ kW and $[5,7,6,4,5,5]$ kW. This example indicates optimal
configuration can change when disturbances are different. OLOC signals are
also compared for some configurations. While the main focus of this article
revolves around the results obtained from small graphs for detailed analysis
and discussions, it is important to note that the code is applicable to graphs
of any size. This versatility is demonstrated in Sec. 6.3, where a graph
consisting of 18 nodes, with a Tank, 14 CPHXs, and 3 junctions CPHXs, is
studied.

###  6.1 Comparing single-split and multi-split cases with 3 and 4 CPHX-nodes

(a) Population ranking (b) Flow rate of the cases shown in Fig.12(a) (c) Wall
temperature of the cases shown in Fig.12(a)

Figure 12: Comparison of multi-split and single-split cases with 3 nodes
having $dist=[12,4,1]$ kW.

Here we aim to find the best architectures among single-split and multi-split
configurations with 3 and 4 nodes. The result for the first scenario with 3
nodes is shown in Fig. 12. In this illustration, only the fluid nodes of CPHX
and the tank are depicted, while the disturbance of each CPHX is displayed in
its respective node. Here, Multi represents multi-split cases and Single
represents single-split cases. The horizontal axis in Fig. 12(a) shows the
percentile score and the vertical axis shows thermal endurance. A percentile
score represents the relative position of a value within a dataset by
indicating the percentage of values that are lower than it. Thus, the best
case is at the top right, and the worst case is at the bottom left. As we see,
the multi-split architecture yield the best result. It should be mentioned
that the results obtained depend on the heat load. For example, here, we have
an extreme load (12) that is much larger than other loads. Here, the multi-
split configurations usually has a better result since the node with the
maximum disturbance is connected to the tank and receives the maximum
available flow-rate. However, in many other cases, the flow rate that this
node receives is a fraction of the pump flow rate.

In Fig. 12, three cases are denoted as $M$ (Multi-split), $S_{1}$ (single-
split-1), and $S_{2}$ (single-split-2). The flow rate and wall temperature of
these three cases are studies in Fig. 12(b), and Fig. 12(c). Among the three
cases, $S_{2}$ exhibits the highest control authority due to its maximum
parallel flows (3). This suggests that the best results can be expected from
this configuration. However, this assumption may not necessarily hold true. A
significant difference between $M$ and $S_{2}$ lies in how the node with a 12
kW heat load is handled. In the case of $M$, this node receives the maximum
flow rate (pump flow rate) as it is directly connected to the pump.
Conversely, in $S_{2}$, the flow received by this node is a fraction of the
pump flow rate. As depicted in Fig. 12(b), the flow rates of these nodes
differ across each graph. Consequently, this disparity impacts the objective
function value and results in varying temperatures. Additionally, Fig. 12(c)
demonstrates that the wall temperature of nodes 12 and 4 reached the upper
bounds. Notably, the case labeled as "$M$" achieves this upper bound at a
later stage compared to the other cases, indicating a better objective
function value (specifically, 69.9 ∘ C).

As mentioned earlier, one of the advantages of using a multi-split graph in
this case is that the node with the maximum heat load is directly connected to
the pump and receives the maximum flow rate. Therefore, one might expect
similar results if all nodes are connected in series directly to the pump.
However, as depicted in Fig. 12(a), when these nodes are arranged in series
with the pump and nodes with higher loads are positioned closer to the pump
($S_{3}$), the achieved result, although satisfactory, is inferior to the
multi-split case. This disparity arises because the dynamics involved in this
problem, such as convection, advection, and bidirectional advection, are also
dependent on the graph’s structure. Consequently, even though the node 12 in
both $M$ and $S_{3}$ receive the pump flow rate, their distinct dynamics lead
to different objective function values. It should also be noted that even with
a fixed structure, altering only the load locations yields different outcomes.
For instance, as illustrated in Fig. 12(a), cases $S_{3}$ and $S_{4}$ possess
the same structure but exhibit substantially different objective values
($t_{f}$). Furthermore, simply changing the positions of the 12 and 4 heat
loads results in significantly different objective values for the multiple
split cases ($M$ and $M_{2}$). In these cases, we generally observe that when
the structure remains fixed, the objective value tends to be better when nodes
with higher heat loads are positioned closer to the pump. This can be
attributed to the fact that nodes with high heat loads require cooler fluid to
dissipate the heat. When these nodes are closer to the pump, the fluid
reaching them is relatively cooler. However, if these nodes are located far
away from the pump, the fluid reaching them is already hot as it has absorbed
heat from other nodes. Consequently, the objective value (in this case,
thermal endurance) decreases. Such studies provide valuable insights for
engineering purposes, enabling us to extract knowledge from optimization data.
In our future work, we intend to expand upon this idea to extract
interpretable knowledge that is understandable to humans.

(a) Population ranking (b) Flow rate of the cases shown in Fig.13(a) (c) Wall
temperature of the cases shown in Fig.13(a)

Figure 13: Comparison of multi-split and single-split cases with 4 nodes,
having $dist=[12,4,1,1]$ kW.

Figure 13 presents a similar study, but this time for graphs with 4 CPHXs. In
the case where we have two nodes with the same heat load ($1$), denoted as
$b_{1}$ in the first branch and $b_{2}$ in the second branch of case $M$, we
observe that the multi-split case ($M$) achieves a higher objective function
value compared to the single-split cases. The optimal solution is found in the
multi-split case ($M$), where the node with the maximum load is connected to
the pump and then divides into two branches. Among the single-split cases, the
best solution is achieved when all nodes are arranged in series and the nodes
with higher heat loads are positioned closer to the tank ($S_{1}$). On the
other hand, if the order of heat loads is reversed ($S_{3}$), the objective
value worsens. This is because in this scenario, the node with the highest
heat load receives fluid that is already hot since it has absorbed heat from
other nodes, resulting in a decrease in thermal endurance.

The flow rate and wall temperature characteristics of three graphs, namely
$M$, $S_{1}$, and $S_{2}$, are investigated in Figure 13(b) and Figure 13(c).
One notable difference between $M$ and $S_{1}$ is that, in $S_{1}$, the flow
rate in all nodes is the same as the pump flow rate, whereas in $M$, it is
not. As depicted in Fig. 13, the flow rate in node $1_{b2}$ of $M$ is nearly
zero initially, increases, and then decreases again. An interesting
observation is the change in the wall temperature of node $1_{b2}$. As shown
in Fig. 13(c), the wall temperature of this node increases as the flow rate in
that branch increases. This phenomenon can be explained by the fact that node
$1_{b2}$ has a lower heat load compared to other nodes, and the coolant fluid
entering this node is already warmer than its wall. Consequently, the flow
rate in this branch increases to allow the coolant fluid to dissipate some of
its heat to the CPHX of this node. Furthermore, as the flow rate in $1_{b2}$
surpasses zero, the coolant flow gets cooler and the rate of change in the
wall temperature of node $12$ decreases, resulting in a delayed approach to
the upper temperature bound (45 degrees Celsius) and thus increasing thermal
endurance. As a result, an optimal coordination of coolant fluid flow rate is
achieved, facilitating optimal heat transfer between different nodes and the
coolant, ultimately leading to the best objective value. A similar situation
was noticed in the case of 3 nodes (see Fig.12) around the 40.0 second mark.
At this point, the flow rate of node 1 was increased to help dissipate heat
from the fluid flow, resulting in cooler fluid. This, in turn, allowed the
other wall nodes to reach the upper bound at a later time, ultimately
improving thermal endurance. Gaining such intuitions through human experience
alone can be challenging. However, optimization studies like these assist
engineers in discovering optimal strategies to solve specific problems. By
leveraging these studies, engineers can avoid the need for trial and error,
thereby significantly reducing the time and resources required to achieve the
desired outcome.

###  6.2 Multi-split configurations with 6 CPHX-nodes and investigation of the
inner-loop results

In this section we define the locations of the CPHXs, and then, based on
spatial location, the junction nodes are produced. For this system, the
location of the CPHXs are defined as:
$[[2,0,0],[2,1,0],[3,1,0],[12,12,0],[15,10,0],[13,13,0]]$. Two sets of
disturbances are considered: case1 = [5,5,5,5,5,5] kW, case2 = [5,7,6,4,5,5]
kW. In this structure, the locations of junctions and their heat-load are
fixed but all other nodes will vary. Therefore, we have nine different
configurations in total, shown in Fig. 14. The result under these two
disturbance sets are shown in Fig. 15. As we see, the optimal configuration
depends on the disturbance values. For example, for the first scenario, the
best result is for the configuration 0, however, for the second scenario, the
best result is for configuration 6.

Figure 14: Different architectures produced by the code when the location of
CPHXs are defined as:
$[[2,0,0],[2,1,0],[3,1,0],[12,12,0],[15,10,0],[13,13,0]]$. Figure 15: Results
obtained for the multi-split configurations shown in Fig. 14 under two
different disturbances. The unit of disturbances is expressed in kW.

To understand what happens to each control and state signal when solving the
OLOC problem, consider Figs. 16 and 17. In the visualization, the fluid nodes
are represented by circles, and the corresponding wall temperatures are
indicated in the plot legend using the notation $w$. For instance, the wall
node linked to the fluid node $0$ is denoted as $w-0$. All nodes are
constrained to a maximum temperature of 45 degrees; if temperature constraints
are not active, it is possible to increase thermal endurance. Optimal thermal
endurance often occurs when all wall nodes reach the upper-bound at the same
time; if one of the nodes reaches the upper-bound sooner than the others, the
thermal endurance could be increased (the capacity of the system has not been
fully utilized).

Figure 16 compares the optimal control trajectories of the first scenario for
three cases: case 0 (best), case 4 (worst), and case 1 (in between). In the
best configuration, all wall temperature nodes ($w$) reach the upper-bound at
the same time. for configuration 1, only nodes 0, 1, and 3 reach the upper-
bound, and for configuration 4 only nodes 0 and 2 reached to the upper bound.
This is reasonable because, while all nodes have the same heat loads, the
nodes that are in series get the same flow rate. As a result, the last node
will have a higher temperature, as the fluid there has already absorbed the
heat from the previous nodes. The results obtained for the flow-rates and all
OLOC signals are shown in this figure.

Figure 16: Investigation of the OLOC signals for 3 cases presented in Fig. 14
and with heat loads = [5,5,5,5,5,5] kW. The multi-split configuration: 0 has
the maximum, 4 has the minimum, and 1 has a value in the mid-range of thermal
endurance.

Figure 17 compares the optimal control trajectories for the second scenario
where the disturbances are not the same for all nodes. Here, results are shown
for case 6 (best), case 4 (worst), and case 3 (in between). In this scenario,
the optimal solution is obtained from configuration 6 where junction 2 with
the maximum heat-load is in series with nodes 0, 1, see Fig. 15(b). In this
scenario, configuration 0 is the next optimal solution among the 9 cases.
Figure 17 also shows the control and state signals. Again for the optimal
solution we observe that all the wall nodes reach the upper-bound temperature
at the same times so this configuration has the maximum thermal endurance. By
comparing configurations 6 and 3 we notice that the only difference between
these architectures is the order of the nodes 0 and 1. In configuration 3, in
the first branch where all nodes are in series, node 0 is the last node and
has a larger disturbances than node 1. When the fluid reaches this node, it is
already warmer than the fluid in node 1; therefore, the temperature for that
node reaches the upper-bound temperature sooner than node 1. On the other
hand, in optimal configuration (6), node 0 is closer to the tank and receives
cooler fluid than does node 1.

In our future work, we plan to expand this study by considering a wider range
of graphs and conduct an in-depth analysis. Additionally, we aim to utilize
machine learning techniques to extract human interpretable data from the
optimization results. This approach will enable us to uncover valuable
insights and understand the underlying patterns and relationships in a more
systematic and interpretable manner. By leveraging machine learning methods,
we can enhance our understanding of the optimization process and gain
actionable knowledge that can be applied to further improve the design and
performance of such systems.

Figure 17: Investigation of the OLOC signals for 3 cases presented in Fig. 14
and with heat loads = [5,7,6,4,5,5] kW. The multi-split configuration: 6 has
the maximum, 4 has the minimum, and 3 has a value in the mid-range of thermal
endurance.

###  6.3 A multi-split configuration with 17 CPHX-nodes

In the previous section, smaller-size graphs were studied to facilitate an
easier discussion of the results. However, the generated code theoretically
has the capability of automatically generating and solving graphs of any size.
In this particular case, we considered a graph with 17 heat-exchanger nodes:
14 CPHXs, and 3 junction CPHXs which were added using the spatial locations of
the CPHXs. Figure 18 displays the generated graph. Some CPHXs are directly
connected to the tank, while others are connected to the created junctions.
All $C_{p}$ Fluid nodes have the same heat load of 4 kW, and the heat loads of
junctions $J_{2}$, $J_{3}$, and $J_{4}$ are 3 kW, 4 kW, and 5 kW,
respectively.

Figure 18 visualises the wall temperature (solid red line), fluid temperature
(dashed-dot green line), and flow rate (dotted-blue line) for each node.
Additionally, the dashed gray line represents the range of these variables,
with temperature ranging from 15 to 45 degrees and flow rate ranging from 0 to
0.4 kg/s. Observing the graph, we can see that the wall temperature of all
nodes at the end of each branch (nodes: 0, 1, 2, 3, 4, 6, 7, 8, 9, 10, 11, 12,
13) reaches the upper bound simultaneously. Similarly, the wall temperature of
$J_{4}$, which has the highest heat load among all three junctions, also
reaches the upper bound. Moreover, the flow rate of $J_{4}$ is higher than
that of $J_{2}$ and $J_{3}$. Additionally, since all nodes have the same heat
load, nodes connected to the same junction (or to the tank) exhibit almost
identical flow rate signals. As an example, the flow rate values of all nodes
within the following three categories are identical: $\\{6,7,8,0,1,2\\}$,
$\\{9,10,11\\}$, and $\\{12,13\\}$. The objective value of this problem
reaches 48.68, indicating the thermal endurance of this system under the given
heat load.

Figure 18: A graph with 18 nodes: 1 tank, 14 CPHXs, and 3 junction CPHXs.
Individual plots, adjacent to each node, display the wall temperature (solid
red line), fluid temperature (dashed-dot green line), and flow rate (dotted
blue line). The plots show temperature (left y-axis) in ${}^{\circ}\mathrm{C}$
and flow rate in $\mathrm{Kg/s}$ (right y-axis) vs time in $\mathrm{S}$
(x-axis). The dashed gray lines indicates the temperature values of 15
${}^{\circ}\mathrm{C}$ and 45 ${}^{\circ}\mathrm{C}$ , and flow rate of 0 and
0.4 $\mathrm{Kg/s}$.

##  7 Conclusion

This article presents the optimal flow control of fluid-based thermal
management systems with multi-split configurations. Graph-based modeling is
used to generate different configurations and to automatically construct their
dynamic equations. In addition, the spatial information of the heat-exchangers
is used to define junctions and determine the level of complexity of the
systems architecture. The presented generative algorithm can be used for
various applications in the domain of configuration design. Next, an open-loop
optimization problem is posed to solve the problem, and enumeration is
employed to determine the best structure.

The results include 3 parts. Part 1 compares the results of single and multi-
split systems composed of 3 or 4 CPHXs and discusses the results in detail.
Part 2 presents and compares the results for multi-split systems composed of 6
CPHXs and gives a comprehensive analysis on the inner-loop optimization.
Finally, part 3 studies a larger scale multi-split system made of 17 CPHXs; an
analysis of the results obtain is presented. The results show that multi-split
configurations result in a better configuration design in some cases. For many
optimal cases, all the wall nodes positioned at the end of CPHX system
branches reach the upper-bound temperature simultaneously. Additionally, we
observed that the optimal configuration depends on the disturbance values.
Finally, the results show the possibility of modeling and analysing large
system made of various CPHXs.

In the next step, we plan to use a population-based algorithm and compare its
results with an enumeration-based framework. Furthermore, our goal is to
utilize machine learning techniques to extract a knowledge base interpretable
by humans from the optimization data. This knowledge base can then be utilized
to guide designers in the process of designing similar systems. Some other
future work items include using more sophisticated hydraulic models,
performing 3D spatial optimization of the pipe network with simultaneous
energy loss minimization, and application of the proposed design framework to
some larger industry-relevant applications.

## References

  * [1] Rootzén, J., Wiertzema, H., Brolin, M., and Fahnestock, J., 2020, “Electrify everything! Challenges and opportunities associated with increased electrification of industrial processes,” .
  * [2] El-Refaie, A. M., 2016, “Growing role of electrical machines and drives in electrification,” 2016 XXII International Conference on Electrical Machines (ICEM), pp. 364–370.
  * [3] Nadel, S. M., 2019, “Electrification in the Transportation, Buildings, and Industrial Sectors: a Review of Opportunities, Barriers, and Policies,” Current Sustainable/Renewable Energy Reports.
  * [4] Mathew, J. and Krishnan, S., 2022, “A review on transient thermal management of electronic devices,” Journal of Electronic Packaging, 144(1).
  * [5] Bayat, S., Nejat Pishkenari, H., and Salarieh, H., 2021, “Observation of Stage Position in a 2-Axis Nano-positioner Using Hybrid Kalman Filter,” Scientia Iranica, 28(5), pp. 2628–2638.
  * [6] Smoyer, J. L. and Norris, P. M., 2019, “Brief historical perspective in thermal management and the shift toward management at the nanoscale,” Heat Transfer Engineering, 40(3-4), pp. 269–282.
  * [7] Wong, H.-S. P., Akarvardar, K., Antoniadis, D., Bokor, J., Hu, C., King-Liu, T.-J., Mitra, S., Plummer, J. D., and Salahuddin, S., 2020, “A density metric for semiconductor technology [point of view],” Proceedings of the IEEE, 108(4), pp. 478–482.
  * [8] Bayat, S., Pishkenari, H. N., and Salarieh, H., 2019, “Observer design for a nano-positioning system using neural, fuzzy and ANFIS networks,” Mechatronics, 59, pp. 10–24.
  * [9] Ye, P., Ernst, T., and Khare, M. V., 2019, “The last silicon transistor: Nanosheet devices could be the final evolutionary step for Moore’s Law,” IEEE spectrum, 56(8), pp. 30–35.
  * [10] He, Z., Yan, Y., and Zhang, Z., 2020, “Thermal management and temperature uniformity enhancement of electronic devices by micro heat sinks: A review,” Energy, p. 119223.
  * [11] Almubarak, A. A., 2017, “The Effects of Heat on Electronic Components,” International Journal of Engineering Research and Applications, 07, pp. 52–57.
  * [12] Mathew, J. and Krishnan, S., 2021, “A Review On Transient Thermal Management of Electronic Devices,” Journal of Electronic Packaging.
  * [13] Feng, C.-P., Chen, L.-B., Tian, G.-L., Wan, S.-S., Bai, L., Bao, R., Liu, Z., Yang, M., and Yang, W., 2019, “Multifunctional Thermal Management Materials with Excellent Heat Dissipation and Generation Capability for Future Electronics.” ACS applied materials & interfaces, 11 20, pp. 18739–18745.
  * [14] Peddada, S. R., 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,” Journal of Mechanical Design, 141(8).
  * [15] Panjeshahi, M., Ataei, A., Gharaie, M., and Parand, R., 2009, “Optimum design of cooling water systems for energy and water conservation,” Chemical Engineering Research and Design, 87(2), pp. 200–209.
  * [16] Muller, C. J. and Craig, I. K., 2016, “Energy reduction for a dual circuit cooling water system using advanced regulatory control,” Applied energy, 171, pp. 287–295.
  * [17] Ling, L., Zhang, Q., Yu, Y., and Shuguang, L., 2016, “Experimental investigation on the thermal performance of water cooled multi-split heat pipe system (MSHPS) for space cooling in modular data centers,” Applied Thermal Engineering, 107, pp. 591–601.
  * [18] Bennett, G., Elwell, C. A., Lowe, R., and Oreszczyn, T., 2016, “The Importance of Heating System Transient Response in Domestic Energy Labelling,” Buildings, 6, p. 29.
  * [19] Allison, J. T. and Herber, D. R., 2014, “Multidisciplinary Design Optimization of Dynamic Engineering Systems,” AIAA Journal, 52(4), pp. 691–710.
  * [20] Allison, J. T., Guo, T., and Han, Z., 2014, “Co-design of an Active Suspension Using Simultaneous Dynamic Optimization,” Journal of Mechanical Design, 136(8).
  * [21] Garcia-Sanz, M., 2019, “Control Co-Design: An Engineering Game Changer,” Advanced Control for Applications: Engineering and Industrial Systems, 1(1), p. e18.
  * [22] Bayat, S. and Allison, J. T., 2023, “Control Co-Design With Varying Available Information Applied to Vehicle Suspensions,” International Design Engineering Technical Conferences and Computers and Information in Engineering Conference, American Society of Mechanical Engineers.
  * [23] Herber, D. R., Guo, T., and Allison, J. T., 2017, “Enumeration of Architectures with Perfect Matchings,” Journal of Mechanical Design, 139(5).
  * [24] Herber, D. R. and Allison, J. T., 2019, “A Problem Class with Combined Architecture, Plant, and Control Design Applied to Vehicle Suspensions,” Journal of Mechanical Design, 141(10).
  * [25] Allison, J. T., Cardin, M.-A., McComb, C., Ren, M. Y., Selva, D., Tucker, C., Witherell, P., and Zhao, Y. F., 2022, “Artificial Intelligence and Engineering Design,” Journal of Mechanical Design, 144(2).
  * [26] Guo, T., Herber, D., and Allison, J. T., 2019, “Circuit Synthesis Using Generative Adversarial Networks (GANs),” AIAA Scitech 2019 Forum, p. 2350.
  * [27] Koeln, J. P., Williams, M. A., Pangborn, H. C., and Alleyne, A. G., 2016, “Experimental validation of graph-based modeling for thermal fluid power flow systems,” Dynamic Systems and Control Conference, Vol. 50701, American Society of Mechanical Engineers, p. V002T21A008.
  * [28] Peddada, S. R. T., Rodriguez, S. B., James, K. A., and Allison, J. T., 2020, “Automated Layout Generation Methods for 2D Spatial Packing,” International Design Engineering Technical Conferences and Computers and Information in Engineering Conference, Vol. Volume 11B: 46th Design Automation Conference (DAC), 10.1115/DETC2020-22627.
  * [29] Peddada, S. R. T., James, K. A., and Allison, J. T., 2020, “A Novel Two-Stage Design Framework for Two-Dimensional Spatial Packing of Interconnected Components,” https://doi.org/10.1115/1.4048817Journal of Mechanical Design, 143(3), 031706.
  * [30] Peddada, S., Zeidner, L., Ilies, H. T., James, K., and Allison, J. T., 2022, “Toward Holistic Design of Spatial Packaging of Interconnected Systems with Physical Interactions (SPI2),” J. Mech. Des, pp. 1–26, https://doi.org/10.1115/1.4055055.
  * [31] Peddada, S. R. T., James, K. A., and Allison, J. T., 2020, “A Novel Two-Stage Design Framework for 2D Spatial Packing of Interconnected Components,” International Design Engineering Technical Conferences and Computers and Information in Engineering Conference, Vol. Volume 11B: 46th Design Automation Conference (DAC), 10.1115/DETC2020-22695, V11BT11A032.
  * [32] Peddada, S. R. T., Dunfield, N. M., Zeidner, L. E., James, K. A., and Allison, J. T., 2021, “Systematic Enumeration and Identification of Unique Spatial Topologies of 3D Systems Using Spatial Graph Representations,” International Design Engineering Technical Conferences and Computers and Information in Engineering Conference, Vol. Volume 3A: 47th Design Automation Conference (DAC), 10.1115/DETC2021-66900, V03AT03A042.
  * [33] Bayat, S. and Allison, J. T., 2023, “SS-MPC: A user-friendly software based on single shooting optimization to solve Model Predictive Control problems,” Software Impacts, 17, p. 100566.
  * [34] Gill, P. E., Murray, W., and Saunders, M. A., 2005, “SNOPT: An SQP algorithm for large-scale constrained optimization,” SIAM review, 47(1), pp. 99–131.
  * [35] Biegler, L. T. and Zavala, V. M., 2009, “Large-scale nonlinear programming using IPOPT: An integrating framework for enterprise-wide dynamic optimization,” Computers & Chemical Engineering, 33(3), pp. 575–582.
  * [36] Patterson, M. A. and Rao, A. V., 2014, “GPOPS-II: A MATLAB software for solving multiple-phase optimal control problems using hp-adaptive Gaussian quadrature collocation methods and sparse nonlinear programming,” ACM Transactions on Mathematical Software (TOMS), 41(1), pp. 1–37.
  * [37] Falck, R., Gray, J. S., Ponnapalli, K., and Wright, T., 2021, “dymos: A Python package for optimal control of multidisciplinary systems,” Journal of Open Source Software, 6(59), p. 2809.