This paper proposes and experimentally validates a hierarchical control framework for fluid flow systems performing thermal management in mobile energy platforms. A graph-based modeling approach derived from the conservation of mass and energy inherently captures coupling within and between physical domains. Hydrodynamic and thermodynamic graph-based models are experimentally validated on a thermal-fluid testbed. A scalable hierarchical control framework using the graph-based models with model predictive control (MPC) is proposed to manage the multidomain and multi-timescale dynamics of thermal management systems. The proposed hierarchical control framework is compared to decentralized and centralized benchmark controllers and found to maintain temperature bounds better while using less electrical energy for actuation.

## Introduction

Modern aircraft and large ground vehicles are complex machines consisting of interconnected subsystems that interact in multiple energy domains with a variety of dynamic timescales. A longstanding trend of these systems is increasing electrical power requirements, necessarily accompanied by increased thermal loading [1,2]. This additional burden on thermal management systems is compounded by the continual desire to reduce the size and weight of mobile platforms, and further exacerbated in aviation by a decreased ability to remove thermal energy due to the use of composite skin materials with high thermal resistance and reduction in ram air heat exchanger cross sections [3]. Many approaches have been explored to meet this challenge, with a focus on employing transient analyses to improve upon traditional steady-state analyses [4]. These approaches include integrated aircraft modeling [3], characterization of component performance under transient operation [5,6], design optimization [7], improved topologies for thermal management [8], enhanced mission planning [9], and advanced control approaches [10–13].

Due to the inherent complexity and multidomain nature of thermal management systems, modular modeling frameworks have been developed to support system and control design. A graph-based approach to modeling not only provides this modularity and scalability but can also readily facilitate model-based control design, as shown in Ref. [14] for building thermal systems, [15] for process systems, and [12] and [16] for mobile thermal energy management systems.

To validate both modeling and control approaches for energy systems, experimental testbeds have been developed across a range of application areas. Examples include the vapor compression refrigeration testbeds of Refs. [17] and [18], the hydraulic hybrid vehicle testbed of Ref. [19], the aircraft fuel thermal management system (FTMS) testbed of [12], and the shipboard chilled water distribution system testbed of Ref. [20].

This paper proposes a hierarchical framework to address the challenge of controlling the multiple timescales and energy domains present in fluid-based thermal management systems. Hierarchical control approaches have been widely investigated in the literature, including [21–25]. While much of this work focuses on deriving analytical guarantees of stability and/or robustness under appropriate assumptions, practical demonstration of modeling and hierarchical control frameworks in application to physical systems is not readily found in the literature. Therefore, the purpose of this paper is to provide such a demonstration. The graph-based model validation of Ref. [26] is first extended to a more complex testbed configuration. Then, this paper builds upon the initial simulation study of Ref. [13] by experimentally demonstrating the hierarchical control approach and comparing this to several baseline control approaches. In comparison to these benchmarks, the hierarchical controller is found to maintain temperature bounds better while using less energy for actuation.

The remainder of this paper is organized as follows. Section 2 presents the graph-based modeling framework. Model validation on an experimental testbed is shown in Sec. 3. Section 4 presents the proposed hierarchical control framework, and Sec. 5 discusses the baseline decentralized and centralized control approaches. Section 6 experimentally demonstrates and compares the three control approaches. Conclusions are provided in Sec. 7.

*Notation:* A vector *x* of elements *x _{i}* is denoted by $x=[xi]$, while $M=[mij]$ denotes a matrix

*M*of elements

*m*. Lower case superscripts are used throughout this paper in the naming of variables, while upper case superscripts indicate mathematical functions, such as a transpose.

_{ij}## Graph-Based Dynamic Modeling

### Generic Graph Formulation.

The graph-based models of this paper are derived by applying conservation equations to a component or system, inherently capturing the storage and transport of mass or energy. When modeling a system as a graph, capacitive elements that store energy, or mass, are represented as vertices and the paths for the transport of energy, or mass, between storage elements are represented as edges. A key feature of graph-based models that makes them attractive for model-based hierarchical control is that a graph-based model of a complete system can easily be partitioned into separate models for individual subsystems by clustering the graph into subgraphs based on an analysis of the edges and vertices. Edges that are cut because of this partitioning represent coupling terms between the subsystems, for which controllers of the hierarchy can account by exchanging information. This ease of partitioning reduces the labor required to generate the many models used throughout a model-based hierarchical control framework.

Let the oriented graph $G=(v,e)$ represent the structure of storage and exchange of a conserved quantity throughout a system $S$. Graph $G$ is said to be of order *N _{v}* with vertices $v=[vi],i\u2208[1,Nv]$, and of size

*N*with edges $e=[ej],j\u2208[1,Ne]$. As shown in the notional graph example of Fig. 1, each edge

_{e}*e*is incident to two vertices and indicates directionality from its

_{j}*tail*vertex $vjtail$ to its

*head*vertex $vjhead$. The set of edges directed into vertex

*v*is given by $eihead={ej|vjhead=vi}$, while the set of edges directed out of vertex

_{i}*v*is given by $eitail={ej|vjtail=vi}$ [27].

_{i}*v*of graph $G$ has an associated dynamic state in system $S$, denoted as

_{i}*x*, representing an amount of stored energy or mass. Similarly, each edge

_{i}*e*has an associated value

_{j}*y*describing the transfer rate of energy (i.e., power flow) or mass between adjacent vertices. The orientation of each edge indicates the direction of positive flow, from $vjtail$ to $vjhead$. Therefore, the dynamics of each vertex satisfy the conservation equation

_{j}where *C _{i}* is the storage capacitance of the

*i*th vertex.

*y*along each edge is a function of the states of the vertices to which it is incident and may also be a function of an input signal

_{j}*u*. The transfer rate along each edge is, therefore, given by

_{j}In addition to capturing the exchange of energy or mass within the graph, the modeling framework must account for exchange with entities external to the graph. Sources to graph $G$ are modeled by source edges $ein=[ejin],j\u2208[1,Nin]$ with associated power flows $yin=[yjin]$, which are treated as disturbances to the system that may come from neighboring systems or the environment. Therefore, edges belonging to *e*^{in} are not counted among the edges *e* of graph $G$, and transfer rates in *y*^{in} are not counted among the internal transfer rates *y* of system $S$.

Sinks of graph $G$ are modeled by sink vertices $vout=[vjout],j\u2208[1,Nout]$ with associated states $xout=[xjout]$. The sink vertices are counted among the vertices *v* of graph $G$ but the sink states *x*^{out} are not included in the state vector *x* of system $S$. Instead, the sink states *x*^{out} are treated as disturbances to the system associated with neighboring systems or the environment.

*M*can then be partitioned as

where the indexing of edges is assumed to be ordered such that $M\xaf$ is a structural mapping from power flows *y* to states *x*, and $M\xaf$ is a structural mapping from *y* to sink states *x*^{out}.

*D*, the dynamics of all states in system $S$ are given by

where $C=diag([Ci])$ is the diagonal matrix of capacitances.

*y*in $S$ is given by

### Domain-Specific Modeling.

An exposition of how conservation of mass and energy can be applied to derive lumped parameter graph-based models for a variety of fluid-thermal components is found in Ref. [16]. These components include a fluid reservoir, a flow split/junction, a pump, a pipe, a cold plate (CP) heat exchanger, and a liquid-to-liquid brazed-plate heat exchanger. Figure 2 shows the mass conservation and thermal energy conservation graphs for each component. For graphs in the hydraulic domain, vertices represent dynamic states of pressure, while edges represent the rate of mass transfer between vertices. For graphs in the thermal domain, vertices represent dynamic states of temperature, while edges represent thermal power flow between vertices due to convection in heat exchangers (HXs) or fluid transport. Dashed lines, indicating sources or sinks of each component, consist of variables determined by neighboring components or disturbances.

### Hydraulic Graph Modeling.

For notation purposes, a superscript *m* is used to denote capacitances, functions, and inputs associated with hydraulic graphs. The reader is referred to Ref. [16] for a detailed derivation of the model equations that follow.

For all hydraulic vertices except those of a reservoir, the hydraulic capacitance is given by $Cim=Vi\rho /E$, where *V* is the fluid volume in the component and both the fluid density *ρ* and the bulk modulus *E* are assumed to be constant. For reservoirs, $Cim=Ac,i/g$, where *A _{c}* is the reservoir cross-sectional area and

*g*is the gravitational constant.

*L*,

*D*, and

*A*are the fluid flow length, diameter, and cross-sectional area, respectively, $\Delta h$ is the height difference between the inlet and outlet flow,

_{c}*f*is the friction factor, and

*K*is the minor loss coefficient. For pumps, the mass flow rate is given by

_{L}*H*is determined using an empirical map as a linear function of pump speed

*ω*and pressure differential across the pump

where $\alpha 1,\u2009\alpha 2$, and $\alpha 3$ are constants.

### Thermal Graph Modeling.

For notation purposes, a superscript *t* is used to denote capacitances and functions associated with thermal graphs. The reader is referred to Ref. [16] for a detailed derivation of the model equations that follow.

For all vertices associated with a fluid temperature, the thermal capacitance is given by $Cit=\rho Vicp$, where the specific heat capacitance of the fluid *c _{p}* is assumed to be constant. For all vertices associated with heat exchanger wall temperatures, $Cit=Mw,icp,w$, where

*M*is the mass of the wall and $cp,w$ is the specific heat capacitance of the wall material.

_{w}*A*is the convective surface area and

_{s}*h*is the heat transfer coefficient, approximated in this paper as a function of mass flow rate and temperature as $h\u2248\beta 1+\beta 2m\u02d9T$, where $\beta 1$ and $\beta 2$ are constants.

### Multigraph System Representation.

Table 1 summarizes the quantities associated with each element of the graph-based framework, including in the generic sense of Sec. 2.1 and in the specific hydraulic and thermal domains of Secs. 2.3 and 2.4, respectively. The combined hydraulic and thermal system dynamics may be represented as two coupled graphs, as shown in Fig. 3. The hydraulic graph, denoted as $Gm$ in Fig. 3, is governed by mass conservation laws, while the thermal graph, denoted as $Gt$, is governed by energy conservation laws. It is assumed in this paper that the coupling between these two graphs is limited to the unidirectional influence of hydrodynamics on thermodynamics.

where $Z\u2208\mathbb{R}Net\xd7(Nem+Next)$.

To capture the dynamics of pumps, including rate limits and time delays between each pump command $uip,i\u2208[1,Np]$ and the actual pump effort $uim$ affecting the hydraulic graph, each $uim$ is paired with a single-input single-output (SISO) system $Sip$ as shown in Fig. 3. Each $Sip$ models the state of the *i*th pump $uim$ as a function of its commanded value $uip$.

### Hydraulic Graph Linearization.

In general, the graph-based models have a nonlinear form but satisfy the generic relationships of Eqs. (1) and (2) for each vertex and edge. For control design, in particular, it is often useful to use a linear representation of the system dynamics. A benefit of the graph-based approach is that a linear model of the full system can be generated by individual linearization of each edge relationship.

*i, j*. Linearizing this expression about an equilibrium operating condition using a first-order Taylor Series gives linear mass flow rate equations of the form

where the $\Delta $ denote deviations from the linearization point and *a _{j}*,

*b*are constant coefficients.

_{j}$M\xaf\u0303m$ represents the columns of $M\xafm$ corresponding to edges associated with pumps, and $b\u0303m$ is the vector of input coefficients for edges associated with pumps (i.e., edges *j* for which $c4,j\u22600$ in Eq. (19)).

### Thermal Graph Linearization.

*i, j*. Linearizing this expression about an equilibrium operating condition using a first-order Taylor Series gives linear power flow equations of the form

## Modeling Example and Validation

### Example Configuration Description.

Appendix A describes the experimental testbed used for validation in this paper. The testbed is pictured in Fig. 4 in an example system configuration used for demonstration throughout this paper. The corresponding schematic is shown in Fig. 5. This configuration is notionally representative of a simplified aircraft FTMS in which heat from actuators, generators, engine oil, and other transient loads is absorbed, stored in liquid fuel, and rejected through transfer to neighboring systems [28].

The sample configuration has eight pumps arranged in four sets of two in series. In this paper, the *a* and *b* pumps of each set receive the same commands. Therefore, for notational convenience, the two pumps in each set are referred collectively. For example, pumps 1a and 1b are collectively termed as “pump 1.”

The secondary loop (identified as the left half of the system in Fig. 5) absorbs thermal loads from the heaters mounted to cold plate 1 (CP 1), through which fluid is driven by pump 1. This loop has dedicated thermal storage available in reservoir 2, and the ability to exchange thermal energy with the main loop across brazed-plate heat exchanger 1 (HX 1) with fluid driven by pump 2.

The primary loop (identified as the right half of the system in Fig. 5) includes two parallel fluid flow paths out of reservoir 2. The path driven by pump 3 passes through heat exchanger 1, exchanging thermal energy with the secondary loop. The path driven by pump 4 passes through CP 2 and CP 3 absorbing thermal energy produced by their heaters. The two flow paths then junction and pass through HX 2, by which thermal energy is transferred out of the system to the thermal sink.

### Graph-Based Representation of Example Configuration.

The hydrodynamics of the example testbed configuration in Fig. 5 are represented by the system graph shown in Fig. 6, formed by the interconnection of the individual hydraulic component graphs from Fig. 2. This hydraulic graph consists of 32 vertices and 34 edges, which in turn set the number of pressure states and mass flow rates in the corresponding graph-based hydraulic model.

where each *Q* is the heat load to the corresponding cold plate heat exchanger, $m\u02d9ext$ is the mass flow rate of chilled fluid through the right side of HX 2, and *T _{c}* is the temperature of the fluid exiting the chiller and entering the right side of HX 2. The thermal graph consists of 39 vertices (one of which is a sink vertex), 41 edges, and 4 source edges. This results in a corresponding graph-based thermal model with 38 temperature states and 41 thermal power flows.

### Experimental Validation of Example Configuration.

Figure 8 shows the commands and disturbances applied to the experimental system and models for validation. The linearization point used for the linear models is the steady-state operating condition of the nonlinear models subject to commands and disturbances that fall approximately in the middle of the operating range. To demonstrate the repeatability of the system across multiple runs, five experimental trials were conducted with the same commanded sequence. The traces for the chiller outlet temperature of Fig. 8 show the envelope between the maximum and minimum values measured at each time among the five trials.

The heat loads plotted in Fig. 8 are translated into a reference current for each bank of resistive heaters using an empirical map between the applied electrical current and the achieved thermal load. Each reference current is tracked by proportional–integral (PI) control of the corresponding solid-state relay.

The chiller is set to track a temperature set point of 20 °C. Figure 8 shows that error from this set point of about $0.5$° C on average is present due to measurement and tracking error within the chiller unit's internal controller.

Figures 9(a) and 9(b) show a selection of hydraulic and thermal signals, respectively, that result from applying the inputs and disturbances of Fig. 8. All experimental traces plotted show the envelope between the maximum and minimum values measured at each time among five experimental trials. To make the width of these envelopes more clear, a closer view of several signals is provided in Fig. 10, which demonstrates that the testbed exhibits a high degree of repeatability.

Figure 9(a) demonstrates a close matching between the experimental data and the hydraulic graph-based models. While offset occurs at times between the models and data, this is generally small relative to the magnitude of the gains when commands change. Where differences exist between the two models, especially in the traces for the pumps 2b and 3b mass flow rate, the nonlinear model is more accurate. This is due to the error incurred by linearization of the terms under the square root in Eqs. (8) and (9).

Figure 9(b) similarly demonstrates a high degree of accuracy in the nonlinear thermal graph-based model. The discrepancies that occur can be attributed to unmodeled friction/losses and errors in heat transfer coefficients. However, in the interval from 750 to 850 s significant error occurs for the linear model in signals pertaining to cold plates of the primary loop. This is largely due to the combination of a low speed in pump 3 and a high speed in pump 4, which falls far from the linearization conditions. Linearization of the inherently bilinear fluid-thermal power flow equation $P=m\u02d9cpT$ decouples the relationship between mass flow rates and temperatures, and this can result in large error under some operating conditions. However, the linear thermal model still preserves the correct gains during this time interval, as is critical to the design of stabilizing model-based controllers. The accuracy of the linear model at most other times across the 1000 s mission is close to that of the nonlinear model. While these models could be improved at the cost of increased complexity, the accuracy of the proposed graph-based models is sufficient for their intended use in closed-loop control.

## Graph-Based Hierarchical Control

### Overview of Hierarchical Control Framework.

Figure 11 diagrams the hierarchical control framework developed for the example testbed configuration. This framework consists of four layers, as enclosed by the dashed box. The update period of each control layer is included in Fig. 11 below its name. Lower layers have faster update rates because they are responsible for governing dynamics of faster timescales, progressing from slow thermal dynamics, to faster thermal dynamics, to hydrodynamics, and finally to pump dynamics in the pump control layer.

The thermal system control layer at the top of the hierarchy is responsible for coordinating the overall long-term thermal behavior of the system, leveraging preview of expected thermal disturbances for the mission at hand. Using model predictive control (MPC) with an update interval of 80 s, this controller is designed to primarily manage slow thermal dynamics, such as reservoir fluid temperatures, cold plate wall temperatures, and heat exchanger wall temperatures.

While the thermal system control layer plans for the coordination of the full thermal system over a long-time horizon, its slow update rate prohibits governing faster system dynamics or quickly compensating for model error or differences between preview of expected future disturbances and the true disturbances that affect the system. However, simply increasing the update rate of the thermal system control layer while maintaining the same time horizon may not be computationally tractable for complex systems, as this would increase the number of prediction steps to be solved by the optimization program while decreasing the time between consecutive updates in which the program must be solved. This motivates the introduction of the faster updating thermal subsystem control layer.

The thermal subsystem control layer consists of two MPC controllers governing the thermal dynamics of the primary and secondary loops. The graph of the full thermal system in Fig. 7 is partitioned into separate graphs for each loop. Using the same general control formulation as for the thermal system controller, graph-based thermal controllers of a 2 s update period are designed for each of the thermal subsystems.

References to be tracked for select temperatures and thermal power flows are communicated from the thermal system control layer to the thermal subsystem control layer. This leverages the system-level coordination and long-time horizon in the thermal system controller, while the faster update rate of the thermal subsystem layer compensates for model and preview error to achieve thermal objectives. Because long-horizon planning is performed by the thermal system controller, the thermal subsystem controllers can be constructed with a relatively small number of steps in their prediction horizons.

The third layer in the hierarchy is the hydraulic subsystem control layer, which like the thermal subsystem control layer is partitioned into subsystems corresponding to each fluid loop. References for mass flow rates are sent from each controller in the thermal subsystem control layer to the corresponding controller in the hydraulic subsystem control layer. In this paper, we assume that only pressure measurements are available as hydraulic feedback, as pressure sensors may be preferred over mass flow rate sensors in implementation for their reduced cost, increased accuracy, and faster response time.

The controllers in the hydraulic subsystem control layer determine references for the pump states. In the pump control layer, a set of decoupled SISO controllers track the desired pump states by commanding the pump inputs, compensating for dynamics and time delays within each pump.

### Thermal System Control Layer Formulation.

*ith*temperature. However, a secondary objective to maintaining temperature constraints is keeping the mass flow rates required to achieve the references small when possible, reducing the pump effort required for the system. The MPC in the thermal control layer solves the following constrained nonlinear program:

In Eqs. (32)–(33), $Nht$ is the number of steps in the prediction horizon of the thermal system control layer MPC. The subscript “*preview*” is included in some variables to indicate that preview of expected values for these signals is assumed to be available over the time horizon of the controller as part of the system's planned mission. Equation (32) defines the cost function as minimizing the thermal slack variable $st=[sit],i\u2208[1,Nvt]$ with weighting $\lambda st$, the mass flow rates under control $m\u02d9$ with weighting $\lambda ut$, the tracking error of a subset of system temperatures *T*^{track} from references *T*^{ref} with weighting $\lambda vt$, the tracking error of a subset of thermal power flows *P*^{track} from references *P*^{ref} with weighting $\lambda et$, and the sum of the differences between mass flow rates in consecutive steps with weighting $\lambda dt$. The last of these objectives is used to ensure that control decisions remain smooth in time across the prediction horizon.

Equation (33*a*) sets the temperatures of the first step in the horizon equal to the current estimated values *T*_{est}. Equations (33*b*) and (33*c*) define *s ^{t}* as measuring the extent by which temperature constraints are violated and impose that $sit(k)\u22650\u2009\u2200i,k$. Equation (33

*d*) imposes a discrete form of the linear thermal graph-based model from Eq. (28). Equation (33

*e*) imposes the nonlinear thermal power flow equation from Eq. (16) to compute those power flows with tracking objectives

*P*

^{track}. While the linearized thermal power flow equation of Eq. (27) could be used instead, the algebraic nature of this equation means that large error can occur at any steps in the horizon when values of $m\u02d9t$ are far from their linearization points, and so when computational resources are sufficient to employ a nonlinear solver, the nonlinear power flow model is preferred.

Equations (33*f*) and (33*g*) constrain the mass flow rates of the system. From Eqs. (17), (33*f*) defines the relationship between the mass flow rates that affect the thermal system $m\u02d9t$, the mass flow rates controlled by the hydraulic system $m\u02d9$, and the mass flow rates that are external disturbances to the system $m\u02d9ext$.

Equation (33*g*) defines the envelope of allowable mass flow rates. The presence of flow splits and junctions in the fluid loops means that it is possible for fluid to flow in the reverse direction from the orientation indicated by the arrows in Fig. 5. As reverse flow would typically be undesirable in an aircraft FTMS, constraints must be included in the controllers to avoid this behavior. In addition, constraints are required to define the envelope of achievable mass flow rates in each loop. The approach used to determine these constraints is detailed in Appendix B. For the thermal system model, $m\u02d9=[m\u02d9prim\u02d9\u2009sec\u2009],\u2009H=[Hpri00H\u2009sec\u2009]$, and $z=[zpriz\u2009sec\u2009]$ from Eqs. (B1) and (B2) of Appendix B.

Finally, Eq. (34) dictates that the mass flow rates chosen for the first step in the prediction horizon $m\u02d9(0)$ be equal to those for the corresponding time chosen in the previous iteration of the controller $m\u02d9last(1)$. This is performed under the assumption that the controller is afforded one step in duration to calculate a solution, after which the controller output decisions corresponding to the *k* = 1 step are applied.

### Thermal Subsystem Control Layer Formulation.

To create a model for each subsystem controller, the system graph in Fig. 7 is partitioned into subsystem graphs corresponding to the primary and secondary loops. The vertex corresponding to the wall temperature of HX 1, which couples the thermal dynamics of the two loops, is represented as a sink vertex in each subsystem graph.

Select temperatures and thermal power flows determined in the thermal system control layer are passed as references to the thermal subsystem control layer. Temperature references are downsampled to the update rate of the latter using a first-order hold, while power flow references are downsampled using a zero-order hold because they are calculated with an algebraic equation and can change quickly due to the dependence on mass flow rates.

The MPC formulation for each of the thermal subsystem controllers is identical to that in Eqs. (32)–(34) subject to use of signals and model parameters corresponding to the thermal subsystem graph-based models diagrammed in Fig. 7. For the primary thermal subsystem $m\u02d9=m\u02d9pri,H=Hpri,$ and $z=zpri$, and for the secondary thermal subsystem $m\u02d9=m\u02d9\u2009sec\u2009,H=H\u2009sec\u2009,$ and $z=z\u2009sec\u2009$.

### Hydraulic Control Layer Formulation.

Mass flow rates determined in the thermal subsystem control layer are passed as references to the hydraulic subsystem control layer. These references are downsampled to the update rate of the hydraulic subsystem controllers using a zero-order hold. Each of the MPC controllers in the hydraulic control layer solves a constrained quadratic program using a discrete form of the linear hydraulic graph-based model from Eqs. (21) and (23). The formulation for the hydraulic controllers is found in Ref. [13].

### Pump Control Layer Formulation.

The pump control layer consists of MPC controllers for each pump. These track pump state references from the hydraulic subsystem control layer *u ^{m}* by commanding the pump inputs

*u*, compensating for pump dynamics and delays. These references are downsampled to the update rate of the pump controllers using a zero-order hold. Each of the MPC controllers in the pump control layer solves a constrained quadratic program using a discrete form of the pump model from Eq. (18). The formulation for the pump controllers is found in Ref. [13].

^{p}## Benchmark Controllers

### Decentralized Benchmark Control.

Figure 12 diagrams the decentralized benchmark controller for the example system configuration of this paper. Each set of pumps is paired with a PI loop that seeks to bring a temperature measurement from the system to track a constant reference. Table 2 lists the signal used as feedback for each PI loop.

The choice of the feedback signal for each loop reflects the primary purpose of the overall system to manage the temperature of the cold plate walls. PI 1 tracks a reference temperature for the CP 1 wall. PI 2 and PI 3 track reference temperatures for reservoir 1 and the HX 1 primary side outlet, respectively, governing the exchange of thermal energy between the secondary and primary loops by controlling fluid flow on either side of HX 1. Because pump 4 moves fluid through both CP 2 and CP 3, using the wall temperature of just one of these cold plates as the feedback signal to PI 4 could result in an inability to properly manage the temperature of the other cold plate. Therefore, the CP 3 outlet fluid temperature, which is affected by the wall temperature of both CP 2 and CP 3, is used as the feedback signal to control pump 4.

### Centralized Benchmark Control.

## Closed-Loop Experiments

While heat loads in the experimental testbed are generated by resistive heaters, in application these would be generated by high-power electrical equipment, such as electrical actuators and batteries, which may have tightly constrained operating temperatures. In this paper, temperature limits are embedded into thermal MPC formulations as soft constraints $T\xaf$ and $T\xaf$, parametrized in accordance with Table 3.

To evaluate the ability of the controllers to manage tight constraints, only 5 °C separate the minimum and maximum temperature constraints for CP 1. Managing the temperature of CP 1 not only involves proper control of the pumps in the secondary subsystem but also coordination with the primary subsystem to transfer thermal energy across heat exchanger 1, providing cooling to the secondary subsystem.

The heat load to each cold plate, the temperature of the chiller outlet fluid, and the mass flow rate of the chiller outlet fluid all serve as disturbances to the system. The chiller is given a constant temperature set point of 20 °C and mass flow rate of 0.35 kg/s. These are assumed to be known to the thermal MPC controllers as preview information across their prediction horizon.

Figure 14 shows the cold plate heat loads for a 4800 s mission, which serves as a case study to evaluate closed-loop performance. Each heat load nominally performs several large step changes on the order of 1 kW, as shown in the top subplot of Fig. 14. However, the applied loads also include additive disturbances consisting of higher frequency uniformly distributed noise stepped every 40 s of up to 200 W in each direction, as shown in the bottom subplot of Fig. 14. The heat loads are saturated to enforce a 1.7 kW maximum load to each cold plate. The nominal loading profile represents reconfiguration of the system as mission phases change over hundreds or thousands of seconds, which in application could coincide with changes in duty cycles of electrical equipment generating heat, charge and discharge of batteries, etc. As such, the nominal loads are assumed to be known to the thermal MPC controllers as mission preview information across their prediction horizon. However, the higher frequency noise added to the nominal heat loads is not included in this preview information and instead serves as an unknown thermal disturbance.

### Decentralized Benchmark Control.

Table 2 lists the specific reference temperatures and gains for each PI controller loop, corresponding to the input/output pairings shown in Fig. 12. These were designed manually by iterating over closed-loop simulations with the goal of bringing the system to achieve thermal objectives. The reference temperature for the CP 1 wall of 42.5 °C falls exactly in the middle of its constraints. The CP 3 fluid outlet temperature reference is 10 °C below the CP 2 wall temperature upper bound, and 5 °C below the CP 3 wall temperature upper bound.

Figure 15(a) shows the pump commands chosen by the decentralized PI controller in experimental application and Fig. 15(b) shows the resulting cold plate wall temperatures. While the temperatures of CP 1 and CP 3 are held closely within their constraints, the temperature of CP 2 is not, with multiple durations of violations more than 5 °C at peak. This emphasizes a drawback of employing a decentralized single-input single-output approach to control complex systems. One actuator may be responsible for managing many elements but the control approach does not incorporate a comprehensive awareness of the condition of all those elements. This is often the case in thermal management systems, where a single fluid line may be responsible for cooling multiple thermal loads [28].

### Centralized Benchmark Control.

The centralized controller is updated every 80 s. The number of steps in the prediction horizon is $Nht,sys=20$, therefore the time horizon is 1600 s. For this controller, *T*^{track} consists of only the temperature of the CP 1 wall, which tracks a constant reference of 42.5 C, and there are no thermal power flow tracking objectives. The following weightings are used in the objective function of Eq. (32): $\lambda st,sys=104,\lambda ut,sys=10\u22124,\u2009\lambda vt,sys=105,\u2009\lambda dt,sys=10\u22123$.

The contribution to the total objective function cost of each term of Eq. (32) is further scaled by the dimension and magnitude of the signals in each term. Therefore, the weightings above indicate that the highest priority is placed via $\lambda st,sys$ on minimizing violations of the temperature bounds. The weighting on minimizing mass flow rates $\lambda ut,sys$ is comparatively small, indicating that this should only be done when possible without violating the temperature constraints.

Figure 16(a) shows the pump commands chosen by the centralized controller in experimental application and Fig. 16(b) shows the resulting simulated cold plate wall temperatures. Large constraint violations occur in cold plate temperatures of the primary loop, which exceed constraints by up to 10 °C. This is due to model error in the slow updating controller. During periods of constraint violation, the controller's model predicts that the CP 2 and CP 3 wall temperatures will decrease to within constraints over the next update period, during which the mass flow rates applied are those chosen for the second step in the horizon of the previous iteration of the controller due to the one-step delay built into the control design to accommodate computation time. This limits the ability of the controller to leverage measurement feedback in compensating for model and preview errors.

### Hierarchical Control.

The thermal system control layer of the hierarchical framework in Fig. 11 is identical to the centralized benchmark controller of the previous section. The two controllers in the thermal subsystem control layer are updated every 2 s, with $Nht,sub=2$ steps in the prediction horizon. A long-time horizon is not critical for the thermal subsystem controllers because their primary objective is to track references from above in the hierarchy that have already been optimized over a long-time horizon.

For each thermal subsystem controller, *T*^{track} consists of references from the thermal system control layer for the temperature of all CP walls, HX walls, and reservoirs present in the subsystem. These represent the slowest dynamics, whose behaviors evolve over the timescale of the thermal system control layer. Because the HX 1 wall is treated as a sink temperature in both subsystem thermal graphs, the temperature for this vertex predicted by the thermal system control layer is included in the sink state preview $Tpreviewout$ to each subsystem controller. *P*^{track} for each thermal subsystem controller consists of references for the thermal power flow along edges that couple the subsystems to each other and to the chiller. In Fig. 7, these are labeled in the primary loop as edges 39 and 42, and in the secondary loop as edge 16. These power flow references ensure that the coordination among subsystems planned in the thermal system control layer is achieved by the thermal subsystem control layer with no requirement for direct communication between subsystem controllers.

The following weightings are used in the objective function of Eq. (32) for each thermal subsystem controller: $\lambda st,sub=107,\lambda ut,sub=10\u22126,\lambda vt,sub=105,\lambda et,sub=103,\lambda dt,sub,pri=108,\lambda dt,sub,\u2009sec\u2009=3\xd7109$. These weightings indicate that the highest priority is on tracking references from the layer above in the hierarchical framework.

The two controllers in the hydraulic subsystem control layer are updated every 1 s, with $Nhm,sub=3$ steps. For all pressures, $p\xafi=10\u2009kPa\u2009\u2200\u2009i$, and $\u2009pi\xaf=200\u2009kPa\u2009\u2200\u2009i$.

The four controllers in the pump control layer are updated every 0.25 s with $Nhp=10$ steps. For all pumps, $u\xafip=20%\u2200\u2009i,$ and $\u2009u\xafip=65%\u2009\u2200\u2009i$.

Figure 17(a) shows the pump commands chosen by the hierarchical controller in the experiment. Figure 17(b) shows the resulting cold plate wall temperatures, which exhibit significantly improved constraint management as compared to the centralized benchmark controller. This improvement is largely due to the inclusion of the thermal subsystem control layer, which leverages a faster update rate to compensate for thermal model error and reject much of the unknown noise added to the nominal heat load profile in tracking references from above in the hierarchy.

### Controller Comparison Summary.

Figure 18(a) shows the total temperature constraint violations under each controller, computed by integrating the magnitude of constraint violations over time for each cold plate. The integral of violations in CP 1 is small enough for all controllers to not be visible in this figure. Figure 18(b) shows the peak temperature violation across the mission in each cold plate under each controller.

From Figs. 18(a) and 18(b), it is clear that the hierarchical framework outperforms both of the baseline controllers, having 7.7% of the total violations of the centralized baseline controller, 21% of the total violations of the decentralized PI baseline controller, and greatly reduced peak violations in CP 2 and CP 3.

Figure 19 shows the total energy consumed by the pumps under each controller, calculated as a function of the electrical current measured for each pump during the experiments. Although the objective of minimizing mass flow rates is given a relatively small weight in the MPC formulations, the centralized and hierarchical approaches still require less total pump energy than the decentralized PI approach. Most importantly, Fig. 19 shows that differences in thermal performance between the three control approaches are not due to increasing the overall actuator effort of the system but instead due to a better allocation of actuation.

## Conclusions

This paper applies a graph-based modeling approach and hierarchical control framework for thermal management. Graph-based dynamic models are derived from conservation of mass and thermal energy, where vertices represent storage elements and edges capture the transport of mass and energy. Experimental validation with a testbed fluid-thermal system demonstrates the high accuracy of the modeling approach. The graph-based framework facilitates model-based control design and is especially well suited to hierarchical control designs, where the control structure should reflect the structure inherent in the graph framework. A scalable hierarchical framework is proposed to manage the multidomain and multi-timescale dynamics present in the system. The proposed hierarchical controller is experimentally demonstrated on a testbed system and compared to decentralized and centralized benchmark controllers, where it is found to perform significantly better in managing thermal objectives by compensating for model error and rejecting unknown disturbances.

In cases where a centralized controller can be solved using a complete system model at a fast update rate and long-time horizon, implementation of hierarchical control is likely not warranted. However, beyond a certain level of complexity executing such a centralized controller becomes intractable given the limited computational resources on board vehicle systems. Therefore, hierarchical control represents an enabling technology for achieving high performance in highly complex thermal management systems, where centralized control is not sufficiently scalable and decentralized control can result in poor performance or necessitate overconservative designs due to an inability to manage coupling between subsystems.

Future work will extend the model-based control design to better capture nonlinear regimes in the system dynamics by performing a local linearization at each controller update and/or representing the system as a set of switched linear models with modes specific to operating regions. Formal procedures will also be explored for choosing the number of layers of the hierarchy, the partitioning of dynamics within each layer, and the update rate of each layer. Future work will also include control of other energy domains present in vehicle energy systems, using the graph-based models for electrical and turbomachinery components presented in Ref. [32]. Lastly, analytical techniques for ensuring stability will be incorporated into the hierarchical control demonstration, such as the passivity approach for graph-based models in Ref. [33], which has been extended to switched graph-based models in Ref. [34].

## Funding Data

National Science Foundation Graduate Research Fellowship (Grant No. DGE-1144245).

National Science Foundation Engineering Research Center for Power Optimization of Electro-Thermal Systems (POETS) with cooperative agreement EEC-1449548.

Air Force Research Laboratory (Contract No. FA8650-14-C-2517).

### Appendix A: Experimental Testbed Overview

This experimental testbed was developed to emulate features of fluid-based thermal management systems while being rapidly reconfigurable to allow for numerous configurations. Table 4 and Fig. 20 contain specifications and images of the components and sensors currently included in the testbed. The working fluid is an equal parts mixture of propylene glycol and water. Components are connected via flexible tubing.

Centrifugal pumps are the primary fluid movers in the system. Speed is controlled via a PWM % duty cycle with less than 20% corresponding to a constant 1300 rpm, 65% and above corresponding to 4500 rpm, and a linear trend between. Peak power consumption of the pumps is 20 W with a peak efficiency of 35%.

Liquid-to-liquid brazed-plate HXs transfer thermal energy between fluid loops in either a parallel-flow or counter-flow configuration.

Each CP heat exchanger consists of several 47 $\Omega $ resistive heaters mounted to an aluminum cold plate that has copper tubing passing through. The heaters on each cold plate are wired to a solid-state relay actuating the heater power output. Up to four heaters can be mounted on each cold plate; however, in this paper just two are used, allowing a maximum heat load of 1.7 kW to be applied to each cold plate.

The reservoirs act as thermal storage elements. A liquid-level sensor inside each reservoir is used to calculate its liquid mass and therefore its thermal capacitance.

A 1.5 HP (1.12 kW) industrial chiller acts as a thermal energy sink (e.g., a vapor compression system). With variable temperature control from –10 °C to 70 °C, the chiller can emulate a wide range of sink conditions.

Infrared cameras were used to identify locations on the HX and CP walls that closely represent the average wall temperature, at which surface temperature sensors are affixed. The infrared image in Fig. 21 shows CP 1 and reservoir 1 of the example testbed configuration in Fig. 4. The cable for the CP 1 wall temperature sensor leads from the center of the plate across its left side.

Sensors and actuators are connected to a National Instruments CompactDAQ, exchanging sensor measurements and actuator commands with National Instruments LabVIEW software on a desktop computer at a rate of 10 Hz. From within LabVIEW, signals can be exchanged with matlab/simulink either by running the two programs simultaneously and communicating via the user datagram protocol or by embedding matlab code in LabVIEW using a matlab script node.

### Appendix B: Hydraulic Coupling Constraints

To determine the constraints on mass flow rate for closed-loop control, the nonlinear hydraulic model is simulated to steady-state at all combinations of pump speeds in the range of 20–65% PWM in increments of 0.25%. As a safety margin against complete flow reversal, any input combinations resulting in a mass flow rate less than 0.03 kg/s are excluded from the allowable hydraulic operating conditions. The resulting envelope of mass flow rates through pumps 3 and 4 is shown in Fig. 22(a).

Figure 22(b) shows the envelope of pump commands generating mass flow rates in the envelope of Fig. 22(a). This demonstrates that coupling between the two pump commands must be taken into account to avoid reverse flow.

where *H*_{pri} is a matrix of appropriate dimensions and *z*_{pri} is a vector. Vertices used to define this polyhedron are circled in Fig. 22(a).