## Abstract

Increasing performance demands and constraints are necessitating the design of highly complex, integrated systems across multiple sectors, including transportation and energy. However, conventional design approaches for such systems are largely siloed and focused on steady-state operation. To accommodate tightening operating envelopes, new design paradigms are needed that explicitly consider system-component interactions and their implications on transient performance at the system design stage. In this work, we present a model fidelity-based decomposition (MFBD) hierarchical control co-design (HCCD) algorithm designed to optimize system performance characteristics, with an emphasis on robustness to transient disturbances during real-time operation. Our framework integrates system level control co-design (CCD) with high-fidelity component design optimization in a computationally efficient manner for classes of highly coupled systems in which the coupling between subproblems cannot be fully captured using existing analytical relationships. Our algorithm permits scalable decomposition of computationally intensive component models and addresses coupling issues between subproblems in part by introducing an intermediate optimization procedure to solve for reduced-order model parameters that maximize the accuracy of the lumped-parameter control model required in the CCD algorithm. We demonstrate the merits of the MFBD HCCD algorithm, in comparison to an all-at-once (AAO) CCD approach, through a case study on aircraft dynamic thermal management. Our results show that our decomposition-based solution matches the AAO optimal cost to within 2.5% with a 54% reduction in computation time.

## 1 Introduction and Background

Performance demands on a wide range of engineered systems are increasing rapidly. For example, increased electrification in the aircraft industry is placing a premium on rapid transient cooling [1,2]. Unfortunately, design of thermal systems, for integration with larger systems, is largely siloed and focused on steady-state operation. More specifically, thermal systems are often designed using offline static optimization routines that size the system and its components to meet relevant constraints under worst-case loading scenarios. A control systems engineer is then tasked with designing the control architecture and/or control policy for the already-designed plant to yield transient (dynamic) performance which meets a set of desired specifications. These types of design methods generally result in over-designed systems with excessive mass or volume that then undermines other performance objectives (e.g., efficiency) in certain types of systems. While this example is specific to challenges in dynamic thermal management, similar problems are occurring across many engineering sectors. Simply put, as advances in technology create the need for more complex and integrated systems with tightening operating envelopes, new design paradigms are needed to meet increased performance demands.

Fortunately, combined plant (system) and control design, also referred to as control co-design (CCD), provides a pathway towards reimagining the way in which we design such systems. CCD is a process wherein a system’s control elements and/or control policy, and thus its dynamic performance, is considered at the plant design stage [311]. In other words, CCD algorithms consider both plant and control variables together throughout the design process. The merits of CCD are rooted in the fact that there are fundamental limitations to what can be achieved with the addition of any control policy for a defined system design [12]. Leveraging CCD to exercise both plant and control degrees-of-freedom earlier in the design process thereby gives design engineers the ability to identify system designs with poor transient performance characteristics prior to a system build.

A comprehensive review of past and present CCD literature can be found in Ref. [13]. A common class of CCD algorithms in the recent literature involves applying a technique called direct transcription, or DT, wherein a system’s dynamic state trajectory and open-loop control signal are discretized and treated as decision variables in an all-at-once (AAO) CCD optimization algorithm. The DT CCD problem can then be solved as a nonlinear program where the system’s discretized dynamic state equations are applied as equality constraints. However, algorithms optimizing open-loop control policies have distinct limitations for other classes of systems. More specifically, open-loop control signals are ill-equipped to handle uncertain or unpredictable exogenous disturbances. For example, next generation aircraft thermal management systems (TMSs) will be required to dissipate time-varying heat load profiles that are either unknown or unpredictable. Therefore, CCD algorithms for these types of thermal systems must include the design of feedback control policies, rather than open-loop control signals, to help ensure the system is designed robustly with respect to the unknown disturbance profiles.

In Ref. [14], we demonstrated the performance benefits of a CCD algorithm that included the design of a feedback control policy with robustness guarantees for TMSs. Results showed that our CCD algorithm could generate system designs with improved transient robustness and little compromise to system design elements, such as physical mass, when compared to more conventional thermal system design methods. However, certain elements of our system model, while physics-based, were derived with a lumped-parameter modeling approach to allow for feedback control design. A general requirement for model-based feedback control synthesis, and thus CCD algorithms optimizing feedback control policies, is a significantly reduced-order model of the system dynamics. However, reduced-order modeling can sometimes present challenges in optimizing physically meaningful high-fidelity design features in a CCD algorithm. Using CCD models with the appropriate fidelity is especially important in the design of thermal-fluid systems, where small changes in specific physical component design variables can have large effects on the highly nonlinear, semi-empirical correlations that capture the heat transfer characteristics of the overall system. Researchers have explored using a technique called balanced co-design to address the issue of design-appropriate models in CCD [13]. In balanced co-design, equality constraints are used to relate high-fidelity physical design variables to intermediate lumped-parameter quantities such as spring constants and damping coefficients used in reduced-order dynamic models [8]. Unfortunately, balanced co-design is not easily amenable to certain classes of systems, including thermal-fluid systems, due to the computationally intensive iterative models required to characterize realizable component performance. Moreover, in some cases, a direct analytical relationship between high-fidelity physical variables and related coupled quantities that are needed to derive a control-oriented model does not exist. This is especially true in systems such as TMSs where heat exchangers with semi-empirical heat transfer correlations influence system dynamics as shown in Fig. 1.

Figure 1(b) depicts a notional TMS where a shell and tube heat exchanger (HX) is used to transfer heat from the system’s hot working fluid to a secondary cold fluid. Figure 1(a) illustrates the complex geometric structure of an actual shell and tube HX. Accurately defining the design and performance of this type of HX requires a high-fidelity design model that relates intricate design features such as tube pitch and shell diameter to semi-empirical heat transfer correlations that ultimately dictate the outlet temperatures of each fluid. This high-fidelity model itself requires a computationally intensive iterative solution where fluid properties are evaluated at the average of respective inlet and outlet temperatures at each successive iteration until the outlet temperatures converge. Figure 1(c), conversely, depicts a typical reduced-order modeling strategy for a HX in a system level context where the component dynamics are lumped into a small number of states that can be measured. This lower level of model fidelity is essential for a dynamic model suitable for feedback control design where model complexity hinders the ability to solve for an optimal control policy. Unfortunately, there is often not an exact analytical relationship between various coupled quantities in the high-fidelity physical design model and the reduced-order model required for dynamic analysis due to the nature of heat transfer correlations and the simplifications required in the reduced-order model. In other words, reduced-order thermal system models used in feedback control analysis, while physics-based, are generally tuned iteratively against experimental data for an already-designed system.

Fig. 1
Fig. 1
Close modal

The computationally-intensive iterative model required to accurately relate physical design variables to performance characteristics of thermal-fluid components, combined with the complicated coupling structure between this model and the reduced-order model required to solve the system level CCD problem, motivates the need for a hierarchical, or decomposition-based, CCD optimization algorithm for classes of thermal-fluid systems. We propose decomposing high-fidelity component models requiring iterative solutions into their own subproblems. By leveraging a hierarchical decomposition-based structure, design decisions with differing fidelity requirements can be coordinated across a hierarchy in a computationally efficient, and scalable, manner.

Hierarchical optimization methods have been researched across various disciplines in both the design optimization and CCD communities with applications ranging from structural engine design to electrical circuit design and multidisciplinary mechatronics design [1524]. The CCD community has also investigated integrating architectural design optimization with co-design for vehicle suspension systems [25]; a similar framework was applied to the architectural design and open-loop control of a single-split architecture for fluid-based thermal management in Ref. [26].

Analytical target cascading (ATC), and more specifically the generalized augmented Lagrangian coordination (ALC), is a popular decomposition-based method for hierarchical design optimization in the literature [13]. In ATC algorithms, a system is decomposed into several subsystems at multiple levels. A single element at the highest level, called the system level subproblem, is used to coordinate the whole system. Target values for higher level elements are cascaded to lower level subsystems and copies of responses, produced by the lower level subproblems, are passed to the upper levels. Variations between targets and responses are then minimized in each decomposed subproblem. ATC has been studied extensively in the context of design optimization of electric vehicles (PHEVs) [27,28] where functional mappings can relate coupling variables across the different levels; recently, authors integrated ATC concepts with a direct transcription CCD problem for a PHEV [29].

ATC and other existing decomposition-based CCD methods are readily applicable for certain types of coupled systems where the various subproblems have directly shared variables or sparse coupling patterns related by easily defined functions or mappings. However, ATC is difficult to apply for CCD of TMSs due to the modeling requirements explained previously and the absence of an exact analytical relationship between certain quantities in the highly coupled system level control-oriented model and decomposed high-fidelity component model(s). Moreover, even if a more defined relationship were to exist, the high level of coupling between the HX and the overall TMS would require several additional decision variables in each ATC subproblem, thereby negating the computational advantage gained via the ATC decision variable coordination method. This is especially true given the highly nonlinear nature of thermal management design problems where multistart optimization techniques are generally required to avoid converging to suboptimal local solutions. Absent from the literature is a hierarchical CCD algorithm that allows computationally efficient, scalable decomposition of subproblems that rely on semi-empirical modeling elements and are strongly coupled to the system level.

In this work, we present a model fidelity-based decomposition (MFBD) hierarchical control co-design (HCCD) algorithm designed to optimize system design variables and feedback control policies to achieve desired steady-state and transient performance objectives. Our algorithm permits scalable decomposition of computationally intensive component models and addresses coupling issues in part by introducing an intermediate optimization step to solve for reduced-order model parameters which maximize the accuracy of the lumped-parameter control-oriented model needed in the system level CCD algorithm. Key contributions and elements of the algorithm are summarized as follows:

1. The MFBD HCCD algorithm integrates CCD with high-fidelity component design optimization by decomposing computationally intensive component models into standalone design subproblems.

2. The top (system) level CCD algorithm optimizes a feedback control policy to ensure controllability and robust transient response for systems subjected to exogenous disturbances with unknown time-varying profiles.

3. Standalone component decomposition accepts desired system level setpoints and allows parallel optimization of high-fidelity physical component designs.

4. An intermediate coupling parameter optimization step maximizes the accuracy of lumped parameters needed in the system level reduced-order feedback control model.

5. A constraint mapping procedure in the intermediate coupling step ensures that system level design updates are made without affecting component level feasibility.

6. A sensitivity-based metric ensures convergence for variables that dominate changes in objective functions while avoiding unnecessary iteration due to changes in design variables that have a weaker effect on the objective functions.

We demonstrate the efficacy of the proposed framework with a case study on dynamic thermal management. We compare our results with an analogous solution to an all-at-once system level nested CCD problem where all decision variables are optimized at a system level.

The rest of the paper is organized as follows. In Sec. 2, we describe the proposed MFBD HCCD algorithm. In Sec. 3, we present the system and component level modeling equations needed to apply the MFBD HCCD algorithm to a notional aircraft thermal management system. In Secs. 4 and 5, we formalize two CCD problems for a case study on dynamic thermal management: one using the proposed MFBD HCCD algorithm and a second using a comparable all-at-once nested CCD algorithm. The results are discussed in Sec. 6 and conclusions from the work are summarized in Sec. 7.

## 2 Model Fidelity-Based Decomposition Hierarchical CCD Algorithm

The top level of the MFBD HCCD framework is the system level, responsible for optimizing all system design and control objectives. The system level is itself a nested optimization with dedicated outer and inner loops. In other words, the outer loop of the system level optimizes a set of system (plant) design variables dj as shown in Eq. (1), subject to system level plant and control constraints along with component feasibility constraints Ci generated from the optimal solutions to the decomposed subproblem(s). At each iteration of the system level outer loop, the solution to the inner loop optimal feedback control problem (Eq. (2)), ϕ(uj*), is included in the objective function. The inner loop uses a linearization of the reduced-order dynamic system model to optimize a set of feedback control design variables uj subject to various controller feasibility constraints; the inner loop solution quantifies the candidate system design’s potential for robust response to unknown time-varying disturbances. We note that for the feedback control architecture considered in this work, uj is not a control input signal or trajectory; rather, the variables in the set uj include a static gain matrix K that is used to compute an optimal feedback control signal during real-time system operation based on system measurements such as lumped fluid temperatures in a TMS.

The system level outputs a set S = {s1, s2, …, sm} of required setpoint vectors si∈[1,m] to be embedded as requirements in the decomposed component level subproblems. The parallel component level optimization, described analytically in Eq. (3), is then used to optimize the design of those components whose model descriptions require computationally expensive iterative solutions. The ith component subproblem takes required setpoints from the system level and optimizes a set of design variables di local to that component, subject to its component-specific constraints. For example, a decomposed HX optimization might take fluid inlet temperatures and mass flowrates, along with a required heat transfer rate, as system level setpoints and optimize a high-fidelity physical HX design geometry to match those requirements. By embedding the system level setpoints in the component level constraint vectors, we can ensure each optimal component design will adhere to necessary constraints needed at the system level. A set D* = {d1*, d2*, …, dm*} of optimal component design vectors di* from the decomposed component level optimization is passed to the subsequent system level iteration.

An intermediate step links the system level with the decomposed component level. At this intermediate step, lumped coupling parameters required to fully define the reduced-order dynamic model used at the system level, xt, are optimized such that the characteristics of the reduced-order model match those of the high-fidelity component models optimized at the last iteration, thereby improving the likelihood of system level feasibility. As an example, in the context of a TMS, these lumped coupling parameters would include heat transfer coefficients that are used in the reduced-order model; the coupling parameter optimization would then be used to minimize outlet temperature differences between the decomposed high-fidelity HX model and those predicted by the reduced-order model.

The intermediate step is also used to create feasibility maps, defined as Ci, for each decomposed component. This involves creating a binary (Boolean) set, denoted by $B$, that characterizes how changes to each setpoint vector si affect component feasibility. For this process, a discretized region Pi(si) around each n-dimensional vector si is mapped to a set of zeroes and ones, where a zero represents a combination of setpoints for which a given component design remains feasible and a one represents a combination of setpoints for which that component design becomes infeasible. For example, consider a HX with four setpoints passed down from the system level, defined as $s∈R4$: hot/cold fluid inlet temperatures and hot/cold fluid flowrates. A feasibility map of zeros and ones would be created based on evaluating a set of HX design feasibility constraints for all combinations of a user-defined range of these setpoints. The resulting binary map Ci would then be applied as a single equality constraint in the next iteration of the system level. Albeit dependent on the choice of Pi, this technique improves the likelihood that the optimal setpoints determined by the system level at iteration k do not cause the component subproblems to become infeasible at iteration k + 1. We also note that while problem-dependent, the combinations of setpoints that comprise Pi should be chosen to cover the allowable design space for a given setpoint. For example, in the context of a TMS, the range of HX mass flowrate setpoints in Pi should cover the full range of mass flowrates that can be achieved by the fluid pump in the system. A specific example of the feasibility mapping procedure will be detailed in Sec. 4.

Figure 2 depicts the operational structure of the algorithm. The MFBD HCCD algorithm starts at the lower component level at iteration index k = 0. The set S(0) of setpoint vectors si is initialized by the user for the first iteration. The decomposed component subproblems are solved in parallel, and the set of optimal component design variables $D*(k)$ is sent to the intermediate coupling step. The intermediate step outputs a vector of optimized reduced-order modeling parameters $xt*(k)$ to link the system level control model with the high-fidelity component models and a set of i = 1, 2, …, m constraint maps $Ci*(k)$ such that system level iteration k chooses plant and control variables that preserve component feasibility. The optimal component variables, lumped coupling parameters, and feasibility maps are then sent to the kth system level iteration. The system level nested CCD algorithm then optimizes a larger set of design variables $dj*(k)$ and control variables $uj*(k)$ that includes a static feedback control policy.

Fig. 2
Fig. 2
Close modal

The iteration index is updated (k = k + 1) and new setpoints S(k) are sent to the component subproblems. The MFBD HCCD algorithm iterates between the system and component levels until convergence is achieved according to Eq. (6), where J represents the Jacobian of a given objective function f. Taking the dot product of each objective function’s Jacobian and the vector of changes in its design variables between successive iterations of the algorithm provides a sensitivity-based convergence criterion. The metric captures the dominant decision variables for a given objective and ensures these variables converge to optimal values before the algorithm is terminated. Consider the term ∂fj/∂di in the Jacobian of a given objective function, fj. This term quantifies the sensitivity of objective function fj with respect to decision variable di. If this sensitivity is low, di does not affect the objective to a high degree and the variable ’s contribution to Eq. (6) will be weighted less relative to changes in other decision variables. In this case, the algorithm will be allowed to converge without unnecessary iteration for less important decision variables. Conversely, if the objective function’s sensitivity to variable di is high, the partial derivative will have a higher value. The criteria in Eq. (6) will then ensure continued iteration until these more influential decision variables converge across multiple iterations.

System Level Control Co-Design
$dj*(k)=argmindjfj(dj,xt*,D*,ϕ(uj*))subjecttogj(dj,xt*,D*)≤0,∃uj*∈(0,∞),Ci(dj,xt*,D*)=0∀i$
(1)
$uj*(k)=argminujϕ(uj,dj,D*,xt*)subjecttogu(uj,dj,D*,xt*)≤0$
(2)
$S={s1(dj*),s2(dj*),…,sm(dj*)}$
Component Decomposition
$di*(k)=argmindifi(di,si)subjecttogi(di,si)≤0$
(3)
$D*={d1*,d2*,…,dm*}$
Coupling and Feasibility Mapping
$xt*(k)=argminxt‖q(dj*,xt)−q(D*)‖2subjecttogj(xt,dj*,D*)≤0,gu(xt,dj*,uj*,D*)≤0$
(4)
$Ci*(k):=gi(di*(k),Pi(si))≤0→B,i∈[1,m]$
(5)
Convergence Criteria
$|Jz(fz,dz*(k))|⋅|(dz*(k)−dz*(k−1))|≤ε∀z$
(6)

## 3 Preliminaries: Modeling a Thermal Management System

In this section, we present a reduced-order (control-oriented) model for a notional aircraft TMS and a high-fidelity design model for a shell and tube HX in the TMS. The models presented here will be used in a case study on MFBD HCCD in Secs. 46 where we optimize the design and control of the TMS.

An aircraft TMS’s (see Fig. 3) general purpose is to use a working fluid to extract heat as needed and safely transfer it to secondary fluids. Every TMS executes four essential processes: heat addition, heat rejection, thermal storage, and thermal transport. The working fluid is stored in the thermal storage component. Within the heat addition (cold plate) component, the working fluid absorbs heat from a source such as power electronics. A pump transports the hotter working fluid around the thermal loop where the absorbed heat is rejected either to ambient or a separate system via a heat rejection component such as a heat exchanger.

Fig. 3
Fig. 3
Close modal

### 3.1 Low-Fidelity SL TMS Model.

As described in Sec. 1, the model used in the system level CCD algorithm must capture the dominant dynamics of the system in a computationally efficient manner such that feedback controller(s) can be designed within the nested CCD inner loop. This is done in part by treating the HX component as a single tube with a lumped fluid capacitance and a lumped tube wall capacitance. While these simplifications allow us to leverage the system level model for optimal feedback control synthesis, they do introduce complications for semi-empirical heat transfer correlations used within the model. These complications will be discussed explicitly in Sec. 3.3. Our system level TMS model was initially derived in Ref. [14] and is summarized here for completeness.

Fluid with mass mt and temperature Tt exits a storage tank at a mass flowrate denoted as $m˙c$. A heat load $Q˙a$ is impinged on the surface of a cold plate heat sink (heat addition component) which has lumped fluid temperature Ta and thermal fluid capacitance Ca. The surface (or wall) of the cold plate is at temperature Tw,a and has capacitance Cw,a. After flowing through the cold plate heat sink, fluid can exit the cycle with mass flow rate $m˙e$. In practice, this might represent fuel exiting the TMS to power the aircraft’s engine. The remaining working fluid enters a shell and tube HX (heat rejection component) with lumped fluid temperature Tr, wall temperature Tw,r, fluid capacitance Cr, and wall capacitance Cw,r. Secondary fluid used for heat transfer within the shell and tube HX has mass flowrate $m˙s$ and inlet temperature Ts. After passing through the HX, the fluid returns to the thermal storage tank. We note that this work will not consider fluid exiting the TMS; in other words, we assume $m˙e=0$.

Six dynamic states describe the first law dynamics of the TMS: mt, Tt, Tr, Ta, Tw,r, and Tw,a. Three quantities—$Q˙a$, Ts, and $m˙e$—are treated as measurable input disturbances. The mass flowrates $m˙c$ and $m˙s$ are the two control inputs. The equations governing the system’s dynamics, derived from first principles, are given in Eq. (7), and the heat transferred from the HX tubes to the colder secondary fluid is computed at each time-step according to Eq. (8). In Eqs. (7) and (8), As is a heat transfer surface area, α is a lumped heat transfer coefficient, and C and cp are fluid capacity and specific heat capacity, respectively. The subscript r refers to a quantity associated with the HX, subscript a refers to a quantity associated with the cold plate, and subscript s refers to a quantity associated with the secondary fluid in the HX component.

Pumps are used to transport both the primary and secondary fluids; the flowrates that can be produced by either pump are functions of that pump’s impeller diameter. Therefore, to incorporate pump design elements in the TMS, we use a static model that maps impeller diameter to relevant specifications such as overall pump mass, achievable mass flowrate, and pump efficiency. Given a baseline, or reference, pump design with an impeller diameter Dref and maximum mass flowrate $m˙¯ref$, we can compute a maximum achievable mass flowrate associated with any candidate pump ($m˙¯cand$) from general pump affinity laws [30] given in Eq. (9a), where Dcand is the impeller diameter required to generate $m˙¯cand$.
$dmtdt=−m˙e$
(7a)
$mtdTtdt=(m˙c−m˙e)(Tr−Tt)$
(7b)
$CrdTrdt=As,rαr(Tw,r−Tr)+(m˙c−m˙e)cp(Ta−Tr)$
(7c)
$CadTadt=As,aαa(Tw,a−Ta)+m˙ccp(Tt−Ta)$
(7d)
$Cw,rdTw,rdt=As,rαr(Tr−Tw,r)+Q˙s$
(7e)
$Cw,adTw,adt=As,aαa(Ta−Tw,a)+Q˙a$
(7f)
$Q˙s=m˙scp,s(Ts−Ts′)Ts′=Tw,r+(Ts−Tw,r)e−αsAs,s/m˙scp,s$
(8)
$m˙¯candm˙¯ref=(DcandDref)3$
(9a)
$mcand=ρπ4wimpDcand2$
(9b)
Assuming a disc-shaped impeller, the mass of the candidate pump can then be computed as shown in Eq. (9b) where wimp is the width of the impeller and ρ is the density of the material from which the impeller is constructed. We also model efficiency curves for each candidate pump. Figure 4 illustrates the effect of mass flowrate on the efficiency curve for two differently sized pumps.
Fig. 4
Fig. 4
Close modal

In Fig. 4, we see that the larger pump (Pump 1) can command a higher mass flowrate ($m˙¯1$) than the smaller pump (Pump 2). Pump 1 also has a higher efficiency, η1, throughout its operating range. However, Pump 1 requires a larger impeller diameter to achieve the higher flowrates and thus comes with additional mass as compared to Pump 2. Figure 4 also illustrates the effect of mass flowrate on a pump’s efficiency. The mass flowrate corresponding to the highest efficiency of a given pump is referred to as the “best efficiency point” or BEP. The shaded regions around each BEP are preferred operating zones. Operating a pump at flowrates far from its BEP can cause wear and fatigue over time, thereby shortening the life span of the pump. Therefore when selecting a pump, it is preferable to choose a design for which the BEP is near the expected nominal operating point for a given application. We note that the pump component designs could also be decomposed and solved as subproblems in our MFBD HCCD framework; however, for demonstration purposes, we only decompose the HX as a subproblem.

The system level model derived here, albeit physics-based, is intentionally designed with lumped-parameter techniques. As described in Sec. 1, the heat transfer characteristics of HXs are very nonlinear functions of the geometry of the component itself, including parameters such as the configuration of tubes within a shell for shell and tube type HXs. In the system level model, however, we intentionally treat the HX component as a single long tube with two states whose dynamics are described by Eqs. (7c) and (7e). By doing so, the system level model can more easily be used for designing feedback controllers in the inner loop of the system level CCD algorithm discussed later in Sec. 4. We can then use a higher-fidelity shell and tube HX model at the decomposed component level, derived in the next subsection, to optimize its detailed design characteristics.

### 3.2 High-Fidelity Shell and Tube HX Model.

To optimize a shell and tube HX design for the aircraft TMS, we use a high-fidelity HX model originally derived in Ref. [31] and summarized here for the reader’s convenience. The high-fidelity model captures the key relationships between the HX geometry and resulting heat transfer correlations. The model provides a physically realizable means of designing a component with the heat transfer characteristics desired at the system level.

Shell and tube HXs contain two different fluid paths: tubeside, or primary, fluid, and shellside, or secondary, fluid, as shown in Fig. 5. The tubeside (primary) fluid enters the tube bundle of the HX via an inlet plenum and is routed through a set of parallel tubes before exiting through an outlet plenum. The shellside (secondary) fluid enters the HX and flows internally relative to the device’s “shell” but outside of the tubes. A set of evenly spaced baffles routes the secondary fluid through the shell as it absorbs heat from the primary fluid within the tubes.

Fig. 5
Fig. 5
Close modal
Figure 5 illustrates several geometric characteristics of a shell and tube HX. The distance between any two tube centers is defined as the tube pitch PT while the distance between any two outer tube edges is referred to as the tube clearance C. Each tube in the tube bundle has an identical inner diameter (Di), outer diameter (Do), and wall thickness (wt). The shell, with diameter Ds, houses a number (NT) of tubes in the bundle. Finally, adjacent baffles are separated by baffle spacing B. This set of variables completely defines the geometry of the HX and allows us to compute heat transfer coefficients for both fluid paths as well as a steady-state rate of heat transfer for any particular geometry. The steady-state heat transfer rate in the HX is given by
$Q˙hx=(UA)hxΔTlmtd$
(10a)
$=(m˙cp)t(Tin−Tout)t$
(10b)
$=(m˙cp)s(Tout−Tin)s$
(10c)
where T represents temperature, subscripts in and out represent inlet and outlet, and subscripts s and t signify shellside fluid and tubeside fluid, respectively. The term ΔTlmtd is a log mean temperature difference [32] computed as
$ΔTlmtd=(Tin,t−Tin,s)−(Tout,t−Tout,s)ln(Tin,t−Tin,sTout,t−Tout,s)$
(11)
and the overall heat transfer coefficient (UA)hx, or overall HTC, is a highly nonlinear function of the intricate tube and shell geometries along with the inlet temperatures and mass flowrates. This overall HTC is computed from standard semi-empirical correlations that are functions of both the tube and shell geometries, along with the inlet temperature and flowrate of each fluid. We note that the overall HTC includes contributions from both tubeside and shellside convection coefficients as detailed in Ref. [31]. With known mass flowrates and inlet temperatures for both fluids, Eq. (10) contains three unknowns: outlet temperatures for the tubeside and shellside fluids (Tt,out, Ts,out) and the nominal (steady-state) rate of heat transfer $Q˙hx$. An iterative process is needed to solve the set of equations due to the fluid property lookups required. At the initial iteration, Eq. (10) is solved with fluid properties evaluated at the known inlet temperatures. The process is repeated iteratively with fluid properties for each fluid evaluated at the average of its inlet temperature and outlet temperature from the previous iteration until the HX outlet temperatures converge to within a defined tolerance.

### 3.3 System Level TMS and HX Component Model Coupling.

As discussed in the Introduction, thermal systems such as the aircraft TMS considered here are characterized by strong coupling among the system components. The lumped-parameter system level TMS model (shown in Fig. 3), and the detailed HX model (shown in Fig. 5) are coupled by both geometric quantities like HX surface area (and thus HX mass) and by model inputs and outputs such as the mass flowrates, inlet temperatures, and outlet temperatures of each fluid. In particular, while the overall HTC in the component level model is computed using an iterative solution that depends on the detailed geometry of the HX, the system level model assumes lumped heat transfer coefficients αr and αs, based on semi-empirical heat transfer correlations, to compute heat transfer. Unfortunately, a fully physics-based analytical relationship between the system level HTCs and the coupled HX outlet temperatures Tr and $Ts′$ does not exist. Therefore, traditional decomposition-based strategies such as ATC are not easily amenable to design optimization for this class of thermal systems. Instead, the intermediate coupling parameter step (Eq. (4)) in our proposed MFBD HCCD algorithm can be used to generate correction factors for the lumped parameter HTCs αs and αr that minimize the difference between coupled parameters, namely the outlet temperatures, of the two models. This method will be explicitly detailed in Sec. 4.

## 4 MFBD HCCD for an Aircraft TMS

In this section, we apply the MFBD HCCD algorithm described in Sec. 2 to the notional aircraft TMS described in Sec. 3. While the analysis here focuses on dynamic thermal management, the framework derived in Sec. 2 can be applied to other types of complex systems with highly-coupled structures and subproblems. The MFBD HCCD algorithm applied to the TMS is depicted visually in Fig. 6, and the corresponding optimization problems for each level are given in Eqs. (12)(16). For this case study, our overall goal is to minimize the total mass of the TMS while ensuring we can safely absorb $Q˙a=5$ kW of heat in the cold plate component without exceeding a nominal (steady-state) cold plate surface temperature Tw,a of 60 °C. We note that the algorithm can be used to optimize other objectives, such as minimizing the overall system footprint. Moreover, other performance objectives can be enforced as constraints by the user, such as constraining the total system footprint to be less than an allowable threshold. Here, however, we focus on minimization of system and component mass to easily compare results with an all-at-once nested CCD algorithm which will be described in Sec. 5.

Fig. 6
Fig. 6
Close modal
System Level Control Co-Design
$dsl*(k)=argmindslλ1⋅M(dsl,xt*,dhx*)+λ2⋅γ*(usl*,dsl,xt*,dhx*)subjecttodsl,min≤dsl≤dsl,maxTw,ae(dsl,dhx*,xt*))≤TmaxTse(dsl,dhx*,xt*)=TcoldS˙gene(dsl,dhx*,xt*)≥0m˙ce(dsl)≤r1⋅m˙¯p,1(dsl)m˙se(dsl)≤r2⋅m˙¯p,2(dsl)ηp,1e(dsl)≥s1⋅η1BEPηp,2e(dsl)≥s2⋅η2BEPChx(dsl,dhx*,xt*)=0∃usl*∈(0,∞)$
(12)
$usl*(k)=argminuslγ(usl,dsl,xt*,dhx*)$
subject to Eqs. (18a)(19b) (13)
$γ>0,ρ(XY)≤γ2$
(13)
$dsl=[mt,m˙c,m˙s,Dchan,Lchan,Nchan,Dp,1,Dp,2]T$
$usl={γ,X,Y}$
$shx={Tae(dsl*),m˙c(dsl*),m˙s(dsl*)}$
HX Decomposition
$dhx*(k)=argmindhxmhx(dhx,shx)subjecttodhx,min≤dhx≤dhx,maxQ˙hx(dhx,shx)≥Q˙desBhx≤Dshell≤5BhxTt,out(dhx,shx)≥Ts,out(dhx,shx)Ac,t(dhx)≤Ac,s(dhx)C(dhx),De(dhx)>0ΔPt(dhx,shx),ΔPs(dhx,shx)≤δPmax$
(14)
$dhx=[Pt,Do,tube,Ltube,Ntube,Dshell,Bhx]T$
$ghx:=Eq.(14)constraints$
Coupling and Feasibility Mapping
$xt*(k)=argminxt=[λs,λr]T‖q(xt,dsl*,dhx*)‖2$
(15)
subject to Eq. (12)(14) constraints
$αs*=λsαs,αr*=λrαr$
$q=[Tre(xt,dsl*)−Tt,o(dhx*),Ts′e(xt,dsl*)−Ts,o(dhx*)]T$
$Chx*(k):=ghx(dhx*,Phx(shx+))≤0→B$
(16)

To initialize the algorithm, we first declare initial setpoint mass flowrates for the primary fluid ($m˙c$) and secondary fluid ($m˙s$), as well as an initial value for the hot fluid inlet temperature for the shell and tube HX. The cold fluid inlet temperature is a constant-valued input to the algorithm. We note that these setpoints are functions of the optimal system level design in subsequent iterations. The following subsections describe the optimization algorithms at each level and their resulting integration into the MFBD HCCD framework introduced in Sec. 2.

### 4.1 Decomposed HX Component Optimization.

An extended version of the HX algorithm presented here was derived in Ref. [31]. Here, we present an amended version of that algorithm to focus on minimization of mass; we also update the heat exchanger decision variable vector and the constraint vector to integrate the component and system levels. More specifically, the decomposed component level accepts system level optimal setpoints shx and optimizes a high-fidelity shell and tube HX design to meet the desired specifications. The full statement of the optimization problem is given in Eq. (14).

The component level HX objective function, given in Eq. (14), is a function of the HX decision variables and the setpoints received from the system level. The decision variables dhx for the HX algorithm are geometric quantities that dictate the heat transfer characteristics of the component, including tube pitch, tube outer diameter, length of a single tube, number of tubes, inner shell diameter, and baffle spacing (see Table 1). We note that the original algorithm in Ref. [31] also treated the tubeside (primary) and shellside (secondary) mass flowrates as decision variables; here, these variables are optimized at the system level and passed down to the component level as setpoints. In other words, the system optimizes the flowrates for the system level objective(s), and the decomposed component level uses the setpoints to design a shell and tube HX to realize the desired system performance. At each iteration of the HX algorithm, Eq. (10) is solved using the HX decision variables and system level setpoints until the shellside and tubeside outlet temperatures (Tout,s and Tout,t) converge to within a tolerance of 10−6 °C.

Table 1

Decision variables for shell and tube HX optimization in MFBD HCCD algorithm

VariableDescription
PtTube pitch
Do,tubeTube outer diameter
LtubeLength of tubes in tube bundle
NtubeNumber of tubes in tube bundle
DshellShell inner diameter
VariableDescription
PtTube pitch
Do,tubeTube outer diameter
LtubeLength of tubes in tube bundle
NtubeNumber of tubes in tube bundle
DshellShell inner diameter

The component level HX optimization problem contains several high-fidelity constraints. First, each design variable is both lower bounded (dhx,min) and upper bounded (dhx,max). Second, we constrain the HX design to achieve a minimum rate of heat transfer $Q˙des$. This ensures that the HX design optimization results in a physical component capable of achieving the heat removal rate, and thus the temperature characteristics, desired by the system level TMS. Additionally, in accordance with typical design practices, the shell diameter is constrained to be larger than the baffle spacing but no larger than five times the baffle spacing. The shell diameter is further constrained to ensure it is large enough to accommodate all NT tubes. The tube clearance C and shell equivalent diameter De, both functions of several decision variables, are required to be positive values. We constrain the tubeside fluid outlet temperature to be greater than the shellside fluid outlet temperature to prevent heat transfer crossover during system operation. Finally, we add constraints to ensure that the pressure drops for the tubeside and shellside fluids are below an allowable tolerance.

### 4.2 Coupling Step and Feasibility Mapping.

As discussed in Sec. 3.3, the semi-empirical nature of HTCs, combined with the lumped-parameter modeling structure required for feedback control synthesis, can introduce error in the system level prediction of the HX dynamics. Therefore, it is not possible to explicitly pass HTCs between levels. Instead, before the system level CCD algorithm can be solved, we must update adjustment factors within the reduced-order model to ensure the lumped parameter HTCs accurately define the behavior of the optimized HX from the component decomposition. The problem statement for our coupling optimization procedure applied to a TMS is defined in Eq. (15).

The coupling procedure takes the last iteration’s optimal HX design $dhx*$ and optimal system design $dsl*$ and uses adjustment factors λs and λr to update the nominal values of the lumped system level HTCs αs and αr such that the differences in HX output temperatures (q) between the two models are minimized. The optimized adjustment factors are used to scale αs and αr, which are themselves functions of lumped HX quantities and the system level temperatures and flowrates in the subsequent iteration of the system level CCD algorithm. We note that a suboptimal initial condition $dsl*(0)$ is used in the coupling parameter optimization preceding the initial system level CCD iteration.

We also use the intermediate step to perform a mapping procedure that ensures HX feasibility is preserved at each iteration of the system level algorithm. This is critical in the context of our framework wherein the system level model does not have access to the high-fidelity component level model and constraint set. For our case study, we create a set $Phx(shx+)$ of setpoint combinations where each vector $shx+$ is a combination of possible values for the n = 3 system level setpoints Ta, $m˙c$, and $m˙s$. We discretize allowable ranges for each setpoint into vectors of m = 50 points. For the inlet temperature Ta, the temperature range considered spans the ambient temperature to the maximum allowable cold plate surface temperature. The setpoint range for each flowrate spans all achievable mass flowrates for the respective primary and secondary fluid pumps. Each of the mn combinations of setpoints is evaluated against the HX constraints in Eq. (14) using the optimal HX geometry that will be passed to the system level. We then map the setpoint combinations to a Boolean set, denoted by $B$, with mn elements to characterize the feasibility of each setpoint combination. A feasible combination of setpoints receives a zero value in the Boolean map and an infeasible combination receives a value of one. The resulting map, defined as Chx, is used with n-dimensional nearest neighbor interpolation to formulate an equality constraint at the system level to ensure optimal system level setpoints (Ta, $m˙c$, and $m˙s$) at iteration k do not violate the HX feasibility constraints for the optimized HX design from that same iteration.

### 4.3 System Level Nested CCD Algorithm.

The full statement of the system level CCD algorithm is given in Eqs. (12) and (13). Recall that the CCD algorithm is itself a nested optimization wherein an outer plant loop optimizes an overall system objective with a contribution from a nested inner control loop that synthesizes an optimal feedback control law to achieve transient performance specifications. In this case study, we use the system level algorithm to minimize a combination of the overall TMS mass M and a robust control metric γ* that characterizes how well a candidate TMS design can respond to time-varying exogenous heat load profiles. More specifically, we minimize the overall mass of the system such that we can maintain a cold plate surface temperature less than Tw,a = 60 °C while absorbing a heat load no less than $Q˙a=5$ kW in magnitude; the inner loop optimizes an H robust feedback control law to regulate (control) the nominal cold plate surface temperature during transient (dynamic) heat loading.

The outer loop decision variables (Table 2) include design variables for each component not optimized at the decomposed level, such as the initial (nominal) amount of fluid in the storage tank (mt) and geometric variables for the channels in the cold plate (Dchan, Lchan, Nchan). Additionally, nominal (steady-state) operational variables that dictate how the components operate together in an integrated system, such as mass flowrates $m˙c$ and $m˙s$, are treated as design variables. The candidate plant design must adhere to a number of constraints. First, each design variable is both upper and lower bounded. We also constrain the nominal cold plate surface temperature $Tw,ae$ to be less than our maximum acceptable 60 °C value during nominal, or steady-state, operation. The nominal value of the secondary fluid inlet temperature $Tse$ of the shell and tube HX is constrained to equal a known value Tcold. In practice, the secondary fluid in the HX may not be controllable; therefore, we treat it as a known disturbance. We constrain the pump impeller diameters (Dp,1 and Dp,2) to ensure the pumps meet constraints related to efficiency and operating range (Fig. 4). In other words, each pump must nominally operate within a user-defined range of its best efficiency point to avoid unnecessary wear on the component. The optimal steady-state mass flowrates $m˙c$ and $m˙s$ are constrained to be less than the maximum achievable pump flowrates by factors of safety r1, r2 < 1 to ensure our feedback controller can modulate each mass flowrate as needed to absorb higher levels of heat during transients. The HX feasibility map described in Sec. 4.2 enforces a component constraint on each system design. Finally, we constrain each system design such that a feasible solution to the inner loop control problem exists.

Table 2

Decision variables for system level CCD optimization in MFBD HCCD algorithm

VariableDescription
mtNominal amount of working fluid in tank
$m˙c$Primary fluid mass flowrate in system
$m˙s$Secondary fluid mass flowrate in HX
DchanDiameter of fluid channels in cold plate
LchanLength of fluid channels in cold plate
NchanNumber of fluid channels in cold plate
Dp,1Impeller diameter of primary fluid pump
Dp,2Impeller diameter of secondary fluid pump
XSolution to Matrix Riccati equation in Eq. (18a)
YSolution to Matrix Riccati equation in Eq. (18b)
γrobust control metric for transient system capabilities
VariableDescription
mtNominal amount of working fluid in tank
$m˙c$Primary fluid mass flowrate in system
$m˙s$Secondary fluid mass flowrate in HX
DchanDiameter of fluid channels in cold plate
LchanLength of fluid channels in cold plate
NchanNumber of fluid channels in cold plate
Dp,1Impeller diameter of primary fluid pump
Dp,2Impeller diameter of secondary fluid pump
XSolution to Matrix Riccati equation in Eq. (18a)
YSolution to Matrix Riccati equation in Eq. (18b)
γrobust control metric for transient system capabilities
The inner loop control decision variables are non-physical variables used to synthesize a state-feedback H controller to meet the control objectives discussed previously [33]. At each iteration of the outer plant loop, our candidate TMS dynamics in Eq. (7) are linearized to produce a state-space linear time-invariant model representation of the form
$x˙=Ax+B1w+B2uz=C1xy=C2x$
(17)
where x is a vector of the system’s dynamic states, u is a vector of control inputs, w is a vector of exogenous disturbances (uncontrollable quantities), y is a vector of physically measurable states for feedback control, and z is a vector of performance outputs to be regulated via feedback control. The matrices A, B1, B2, C1, and C2 characterize the dynamic performance of the system near its nominal operating point. We note that the vector of control inputs u in Eq. (17) is different from the control decision variables usl defined in the system level-nested CCD algorithm. Instead, usl represents the nonphysical decision variables in the CCD inner loop that are used, in part, to optimize a static gain matrix K that is then used to compute the physical control input u in real time.
The H controller optimizes a state-feedback control law to minimize adverse effects of unknown disturbance profiles w on regulated signals z. In other words, H optimal feedback control formulations quantify the best theoretical robust performance, or disturbance rejection capability, attainable for a given system design. In this case study, we would like to regulate the cold plate surface temperature in response to transient heat loading while assuming full-state feedback; in other words, we assume C1 = [0, 0, 0, 0, 0, 1] and C2 is an identity matrix. The optimal control problem is stated fully in Eq. (13). We seek the smallest positive value of decision variable γ for which the feedback control constraints are satisfied. In Eq. (13), the matrices X and Y are constrained to be symmetric matrices that are solutions to the following matrix algebraic Riccati equations:
$ATX+XA+C1TC1+X(γ−2B1B1T−B2B2T)X=0$
(18a)
$AY+YAT+B1B1T+Y(γ−2C1TC1−C2TC2)Y=0$
(18b)
Additional constraints are placed such that the eigenvalues λ of the following matrices each have negative real values:
$Reλi[A+(γ−2B1B1T−B2B2T)X]<0∀i$
(19a)
$Reλi[A+Y(γ−2C1TC1−C2TC2)]<0∀i$
(19b)
Finally, the spectral radius ρ(XY), or the eigenvalue with the largest absolute value, of the product XY, is constrained to be less than γ2. The control optimization uses a bisection method to find the smallest value γ for which this feedback controller is realizable. The optimal value γ*, which is embedded in the outer loop objective function, quantifies the best achievable robust performance for a given plant. In other words, a smaller γ* guarantees better robust performance. An optimal feedback control gain matrix, computed as $K=−B2TX$, is then guaranteed to achieve optimal transient performance of the TMS during operation. A full discussion and derivation of the optimal H control policy used here can be found in Chapter 9 of Ref. [33].

Equations (12) and (13) illustrate that the system level objective function is not only a function of the system level design and control variables dsl and usl but also of the optimal solutions from the previous iteration of the decomposed HX algorithm $dhx*(k)$ and the optimized coupling parameters $xt*(k)$. For example, the total mass of the optimal HX from Eq. (14) is included in the system level mass calculation (M) in Eq. (12). Once an optimal plant and controller are designed at the system level, the new optimal mass flowrates ($m˙c$, $m˙s$) and nominal HX inlet temperature (Ta) for that system design are passed down to the component level as new setpoints for the next iteration of the HX optimization algorithm.

While the system level algorithm in our MFBD HCCD framework does not have full access to an entire set of design variables for the HX, we will show that we can generate optimal system and control designs with results very close to those generated with a true all-at-once nested CCD algorithm. In the following section, we present an analogous all-at-once CCD problem formulation that will be used as a benchmark for our algorithm.

## 5 All-at-Once Nested CCD for an Aircraft TMS

To better characterize the efficacy of our decomposition-based framework, we formulate and solve an analogous AAO nested CCD algorithm wherein all design and control variables are included in a single system level optimization. The full statement of the AAO CCD problem is given in Eq. (20).

All-at-Once Nested Control Co-Design
$d¯*=argmind¯λ1⋅M(d¯,xt*)+λ2⋅γ*(u*,d¯,xt*)subjecttod¯min≤d¯≤d¯maxgsl(d¯,xt*)ghx(d¯)∃u¯*∈(0,∞)$
(20)
$u¯*=argminu¯γ(u¯,d¯,xt*)subjecttogu(d¯,xt*,u¯)$
$xt*=argminxt=[λs,λr]T‖[q1(xt,d¯),q2(xt,d¯)]T‖2d¯=[dsl,dhx]Tu¯={γ,X,Y}dsl=[mt,m˙c,m˙s,Dchan,Lchan,Nchan,Dp,1,Dp,2]dhx=[Pt,Do,tube,Ltube,Ntube,Dshell,Bhx]αs*=λs*αs,αr*=λr*αrq1=Tre(xt,dsl)−Tt,o(dhx)q2=Ts′e(xt,dsl)−Ts,o(dhx)gsl:=Eq.(12)constraintsghx:=Eq.(14)constraintsgu:=Eq.(13)constraints$

In addition to the system design variables, the AAO problem also includes the full set of heat exchanger design variables. At each iteration of the AAO algorithm, the HX model described in Sec. 3.2 is used to iteratively solve for the HX outlet temperatures and heat transfer rate resulting from a given set of design variables. An unconstrained parameter adjustment optimization is nested within the overall algorithm in Eq. (20) to generate lumped HTC adjustment factors that guarantee the linearized reduced-order control model required for the controller optimization is accurate with respect to the current iteration’s HX design. The AAO nested inner loop control optimization in Eq. (20) is the same as derived in Sec. 4.

While the AAO CCD algorithm does contain a single level with access to the full design space, it is extremely computationally intensive due to the iterative HX model now included at each iteration as well as the need for the nested coupling parameter optimization within each iteration. Because of this, the AAO algorithm is not easily scalable for larger and more complex systems. It is also important to note that this AAO CCD algorithm does not allow for an exact one-to-one comparison with the MFBD HCCD algorithm. This is primarily due to the additional nested coupling parameter optimization required in the AAO CCD algorithm to relate the high-fidelity HX design to an analogous reduced-order dynamic model within a given iteration. However, it still serves as a reasonable benchmark for evaluating the proposed decomposition-based strategy.

## 6 Design Optimization Results: MFBD HCCD vs. AAO Nested CCD

The results presented in this section directly compare our proposed MFBD HCCD framework with the analogous single-level AAO nested CCD algorithm. Each optimization problem is solved using the NOMAD nonlinear optimization algorithm which interfaces with matlab through the open-source OPTI toolbox [34]. A multistart approach was applied in each case to avoid convergence to suboptimal local minima. The multistart approach involves sampling the bounded design space at 2n points, where n is the number of decision variables, and then solving the given optimization problem starting from the 10 best feasible points. Each simulation was performed on a 32-core machine with a 3.0 GHz processor and 32 GB RAM. Values for various parameters and weightings used in the case study are shown in Table 3.

Table 3

Optimization parameters and weightings used in comparative CCD case study

Parameter (units)ValueParameter (units)Value
λ1 (–)1λ2 (–)0.01
r1 (–)0.90r2 (–)0.90
s1 (–)0.80s2 (–)0.80
Tcold (C)10$Q˙des$ (kW)5
Tmax (C)60δPmax (kPa)30
Parameter (units)ValueParameter (units)Value
λ1 (–)1λ2 (–)0.01
r1 (–)0.90r2 (–)0.90
s1 (–)0.80s2 (–)0.80
Tcold (C)10$Q˙des$ (kW)5
Tmax (C)60δPmax (kPa)30

We assume water to be both the primary fluid for the TMS loop and the secondary fluid in the HX component. The optimal design variables for both the MFBD HCCD approach and the AAO nested CCD approach are given in Table 4, along with the initial conditions and upper/lower bounds used in the case study. We note that the MFBD HCCD optimization problem took four iterations across each level to converge; more specifically, both levels converged to within a tolerance ɛ = 0.1 (see Eq. (6)) after the fourth iteration across the hierarchy.

Table 4

Comparison of optimal system designs and characteristics: MFBD HCCD versus AAO Nested CCD

DesignUnitsInitialLowerUpperMFBD HCCDAAO nested CCDPercent
variableconditionboundboundoptimal valueoptimal valuedifference
mtkg5.01.01001.001.000.00
$m˙c$kg/s0.100.0100.500.04860.041517.0
$m˙s$kg/s0.200.0101.00.1920.1872.54
Dchanm6.3 × 10−32.0 × 10−40.151.50 × 10−31.30 × 10−311.2
Lchanm0.100.0130.610.02650.038431.0
Nchan5020100022246.69
Dp,1m0.300.0760.460.1870.1775.62
Dp,2m0.300.0760.460.2980.2950.847
Ptm6.3 × 10−33.2 × 10−30.0154.16 × 10−34.23 × 10−31.71
Do,tubem3.2 × 10−32.4 × 10−30.0142.40 × 10−32.40 × 10−30.00
Ltubem0.910.0253.70.5040.5080.740
Ntube12011000100981.71
Dshellm0.0760.0500.300.05080.05080.00
Bhxm0.0250.0250.300.02540.02540.011
M* (kg)15.97715.5962.45
γ* (–)0.3800.3605.55
System level objective function value (–)15.98115.5992.45
Run time (h)3.7318.13554.1
DesignUnitsInitialLowerUpperMFBD HCCDAAO nested CCDPercent
variableconditionboundboundoptimal valueoptimal valuedifference
mtkg5.01.01001.001.000.00
$m˙c$kg/s0.100.0100.500.04860.041517.0
$m˙s$kg/s0.200.0101.00.1920.1872.54
Dchanm6.3 × 10−32.0 × 10−40.151.50 × 10−31.30 × 10−311.2
Lchanm0.100.0130.610.02650.038431.0
Nchan5020100022246.69
Dp,1m0.300.0760.460.1870.1775.62
Dp,2m0.300.0760.460.2980.2950.847
Ptm6.3 × 10−33.2 × 10−30.0154.16 × 10−34.23 × 10−31.71
Do,tubem3.2 × 10−32.4 × 10−30.0142.40 × 10−32.40 × 10−30.00
Ltubem0.910.0253.70.5040.5080.740
Ntube12011000100981.71
Dshellm0.0760.0500.300.05080.05080.00
Bhxm0.0250.0250.300.02540.02540.011
M* (kg)15.97715.5962.45
γ* (–)0.3800.3605.55
System level objective function value (–)15.98115.5992.45
Run time (h)3.7318.13554.1

The results in Table 4 show that the percent difference in the optimal system level objective function value between the two cases was less than 2.5%, while the MFBD HCCD algorithm produced a solution in approximately 54% less time. When comparing the two solutions, we see that the percent difference in optimal values between the two cases is less than 6% for ten of the 14 decision variables. The largest differences are in the variables pertaining to the cold plate component’s channels; this is likely due to the fact that the cold plate contributes a small amount of mass to the overall system when compared with the HX and the two pumps. Another difference between the two cases is in the primary fluid flowrate $m˙c$, where the optimal value was 0.0486 kg/s for the MFBD HCCD algorithm compared with 0.0415 kg/s for the AAO CCD case. The difference in these values is likely due to the fact that the system level does not have access to the higher-fidelity HX design space in the MFBD HCCD algorithm, so a flowrate that satisfies the latest decomposed HX design constraints must be chosen in accordance with the feasibility mapping procedure outlined in Eq. (16). Nevertheless, the proposed MFBD HCCD algorithm produces very similar optimal design results to the AAO nested CCD case in approximately half the computation time.

Designing for optimal transient performance in the presence of unknown time-varying disturbances is a critical aspect of the proposed framework. The inner control loop in our system level nested CCD algorithm generates an optimal robustness metric γ* along with a corresponding set of feedback control gains to optimally control the system under transient loading. Figures 7 and 8 compare the simulated transient performance of each optimal system. In Fig. 7(a), each system is subjected to a series of transient heat pulses up to approximately 15% above the nominal 5 kW value. Figure 7(b) shows the regulated cold plate surface temperature resulting from each optimal design and Fig. 7(c) shows the HX outlet temperatures increase as more heat is absorbed by each system’s primary fluid. Figure 8 shows each controller’s response in modulating the primary and secondary fluid mass flowrates. We see that the cold plate surface temperature is regulated to within one-half degree Celsius during the transient loading in both the MFBD HCCD and AAO CCD cases; this is done primarily by modulating the primary flowrate $m˙c$ higher during the transients.

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

While we do not formally prove convergence of the proposed MFBD HCCD algorithm, the results of this case study verify that convergence of the algorithm itself is achievable, and moreover that the solution approximates that of a comparable non-decomposed CCD algorithm. Future work will investigate specific conditions under which convergence can be guaranteed.

## 7 Conclusions

In this work, we presented a MFBD HCCD algorithm designed to achieve desired steady-state and transient performance characteristics via optimization of the system design and an associated feedback control policy. Our algorithm permits scalable decomposition of computationally intensive component models and addresses coupling issues in part by introducing an intermediate optimization step to solve for reduced-order model parameters that maximize the accuracy of the lumped-parameter control-oriented model needed for feedback control synthesis in the system level CCD algorithm. We demonstrated the efficacy of the proposed algorithm with a case study on dynamic thermal management. More specifically, we compared our results to the solution of a comparable all-at-once (AAO) system level nested CCD problem. The results showed that our proposed MFBD HCCD algorithm matched the AAO CCD solution to within 2.5% while locating a solution in roughly 54% less time. Future work will investigate formulating a rigorous proof of convergence for the proposed algorithm as well as conducting experimental validation of the design approach for integrated thermal systems.

## Acknowledgment

This work is supported by the U.S. Office of Naval Research Thermal Science and Engineering Program under contract number N00014-17-1-2333.

## Conflict of Interest

There are no conflicts of interest.

## References

1.
Doty
,
J.
,
Yerkes
,
K.
,
Byrd
,
L.
,
Murthy
,
J.
,
Alleyne
,
A.
,
Wolff
,
M.
,
Heister
,
S.
, and
Fisher
,
T. S.
,
2017
, “
Dynamic Thermal Management for Aerospace Technology: Review and Outlook
,”
J. Thermophys. Heat Transfer
,
31
(
1
), pp.
86
98
. 10.2514/1.T4701
2.
Bodie
,
M.
,
Russell
,
G.
,
McCarthy
,
K.
,
Lucas
,
E.
,
Zumberge
,
J.
, and
Wolff
,
M.
,
2010
, “
Thermal Analysis of an Integrated Aircraft Model
,”
48th AIAA Aerospace Sciences Meeting Including the New Horizons Forum and Aerospace Exposition
,
Orlando, FL
, AIAA, p.
288
.
3.
Fathy
,
H. K.
,
Papalambros
,
P. Y.
,
Ulsoy
,
A. G.
, and
Hrovat
,
D.
,
2003
, “
Nested Plant/Controller Optimization With Application to Combined Passive/Active Automotive Suspensions
,”
Proceedings of the 2003 American Control Conference
,
Denver, CO
,
June 4–6
, Vol.
4
,
IEEE
, pp.
3375
3380
.
4.
Allison
,
J. T.
, and
Nazari
,
S.
,
2010
, “
Combined Plant and Controller Design Using Decomposition-Based Design Optimization and the Minimum Principle
,”
ASME 2010 International Design Engineering Technical Conferences and Computers and Information in Engineering Conference
,
,
Aug. 15–18
,
American Society of Mechanical Engineers
, pp.
765
774
.
5.
Hung
,
Y.-H.
, and
Wu
,
C.-H.
,
2015
, “
A Combined Optimal Sizing and Energy Management Approach for Hybrid In-Wheel Motors of EVs
,”
Appl. Energy
,
139
(
1
), pp.
260
271
. 10.1016/j.apenergy.2014.11.028
6.
Youcef-Toumi
,
K.
,
1996
, “
Modeling, Design, and Control Integration: A Necessary Step in Mechatronics
,”
IEEE/ASME Trans. Mechatron.
,
1
(
1
), pp.
29
38
. 10.1109/3516.491407
7.
Peters
,
D. L.
,
Papalambros
,
P. Y.
, and
Ulsoy
,
A. G.
,
2013
, “
Sequential Co-Design of an Artifact and its Controller Via Control Proxy Functions
,”
Mechatronics
,
23
(
4
), pp.
409
418
. 10.1016/j.mechatronics.2013.03.003
8.
Allison
,
J. T.
,
Guo
,
T.
, and
Han
,
Z.
,
2014
, “
Co-Design of an Active Suspension Using Simultaneous Dynamic Optimization
,”
ASME J. Mech. Des.
,
136
(
8
), p.
081003
. 10.1115/1.4027335
9.
Jiang
,
Y.
,
Wang
,
Y.
,
Bortoff
,
S. A.
, and
Jiang
,
Z.-P.
,
2015
, “
Optimal Codesign of Nonlinear Control Systems Based on a Modified Policy Iteration Method
,”
IEEE Trans. Neural Netw. Learn. Syst.
,
26
(
2
), pp.
409
414
. 10.1109/TNNLS.2014.2382338
10.
Patil
,
R.
,
Filipi
,
Z.
, and
Fathy
,
H.
,
2010
, “
Computationally Efficient Combined Design and Control Optimization Using a Coupling Measure
,”
IFAC Proc. Vol.
,
43
(
18
), pp.
144
151
. 10.3182/20100913-3-US-2015.00126
11.
Abzug
,
M. J.
, and
Larrabee
,
E. E.
,
2002
,
Airplane Stability and Control: a History of the Technologies that Made Aviation Possible
(
2nd ed. No. 14 in Cambridge Aerospace Series
),
Cambridge University Press
,
Cambridge, UK
.
12.
Goodwin
,
G.
,
Graebe
,
S.
, and
,
M.
,
2001
,
Control System Design
,
Prentice Hall
,
.
13.
Allison
,
J. T.
, and
Herber
,
D. R.
,
2014
, “
Special Section on Multidisciplinary Design Optimization: Multidisciplinary Design Optimization of Dynamic Engineering Systems
,”
AIAA J.
,
52
(
4
), pp.
691
710
. 10.2514/1.J052182
14.
Nash
,
A. L.
, and
Jain
,
N.
,
2019
, “
Combined Plant and Control Co-Design for Robust Disturbance Rejection in Thermal-Fluid Systems
,”
IEEE Trans. Control Syst. Technol.
, pp.
1
8
15.
Renzi
,
C.
,
2016
, “
A Genetic Algorithm-Based Integrated Design Environment for the Preliminary Design and Optimization of Aeronautical Piston Engine Components
,”
,
86
(
9–12
), pp.
3365
3381
. 10.1007/s00170-016-8433-7
16.
Zhao
,
W.
,
Yang
,
Z.
, and
Wang
,
C.
,
2018
, “
Multidisciplinary Hybrid Hierarchical Collaborative Optimization of Electric Wheel Vehicle Chassis Integrated System Based on Driver’s Feel
,”
Struct. Multidiscip. Optim.
,
57
(
3
), pp.
1129
1147
. 10.1007/s00158-017-1801-6
17.
Liu
,
Y.
,
Yin
,
Y.
, and
Guo
,
Z.
,
2012
, “
Static and Dynamic Design Based on Hierarchical Optimization for Materials and Structure of Porous Metals
,”
Sci. China Technol. Sci.
,
55
(
10
), pp.
2808
2814
. 10.1007/s11431-012-4917-3
18.
Martins
,
R.
,
Lourenço
,
N.
, and
Horta
,
N.
,
2015
, “
,”
Expert Syst. Appl.
,
42
(
23
), pp.
9137
9151
. 10.1016/j.eswa.2015.08.020
19.
Zheng
,
C.
,
Hehenberger
,
P.
,
Le Duigou
,
J.
,
Bricogne
,
M.
, and
Eynard
,
B.
,
2017
, “
Multidisciplinary Design Methodology for Mechatronic Systems Based on Interface Model
,”
Res. Eng. Des.
,
28
(
3
), pp.
333
356
. 10.1007/s00163-016-0243-2
20.
Takata
,
S.
, and
Umeda
,
Y.
,
2007
,
,
Springer
,
London
.
21.
Jansen
,
S.
, and
Welp
,
E. G.
,
2005
, “
Model-based Design of Actuation Concepts: a Support for Domain Allocation in Mechatronics
,”
Proceedings of the 15th International Conference on Engineering Design
,
Melbourne, Australia
, The Design Society.
22.
,
Z.
,
Singule
,
V.
,
Vechet
,
S.
, and
Ondrusek
,
C.
,
2010
, “
Development of Energy Harvesting Sources for Remote Applications as Mechatronic Systems
,”
Proceedings of the 14th International Power Electronics and Motion Control Conference
,
Ohrid, Macedonia
,
Sept. 6–8
,
IEEE
, pp.
T10
T13
.
23.
Wang
,
L.
, and
Nee
,
A. Y.
, eds.
2009
,
Collaborative Design and Planning for Digital Manufacturing
,
Springer
,
London
.
24.
Liu
,
T.
,
Azarm
,
S.
, and
Chopra
,
N.
,
2017
, “
A Decentralized Approach for Multi-Subsystem Co-Design Optimization Using Direct Collocation Method
,”
Volume 2A: 43rd Design Automation Conference
,
Cleveland, OH
,
Aug. 6–9
,
American Society of Mechanical Engineers
, p. V02AT03A004.
25.
Herber
,
D. R.
, and
Allison
,
J. T.
,
2019
, “
A Problem Class With Combined Architecture, Plant, and Control Design Applied to Vehicle Suspensions
,”
ASME J. Mech. Des.
,
141
(
10
), p.
101401
. 10.1115/1.4043312
26.
,
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
,”
ASME J. Mech. Des.
,
141
(
8
), p.
083401
. 10.1115/1.4043203
27.
Alexander
,
M. J.
,
Allison
,
J. T.
,
Papalambros
,
P. Y.
, and
Gorsich
,
D. J.
,
2011
, “
Constraint Management of Reduced Representation Variables in Decomposition-Based Design Optimization
,”
ASME J. Mech. Des.
,
133
(
10
), p.
101014
. 10.1115/1.4004976
28.
Alexander
,
M. J.
,
Allison
,
J. T.
, and
Papalambros
,
P. Y.
,
2011
, “
Reduced Representations of Vector-Valued Coupling Variables in Decomposition-Based Design Optimization
,”
Struct. Multidiscip. Optim.
,
44
(
3
), pp.
379
391
. 10.1007/s00158-011-0636-9
29.
Behtash
,
M.
, and
Alexander-Ramos
,
M. J.
,
2018
, “
Decomposition-Based MDSDO for Co-Design of Large-Scale Dynamic Systems
,”
Volume 2A: 44th Design Automation Conference
,
,
Aug. 26–29
,
American Society of Mechanical Engineers
, p. V02AT03A003.
30.
Fox
,
R. W.
, and
McDonald
,
A. T.
,
1994
,
Introduction to Fluid Mechanics
,
John Wiley and Sons
,
New York
.
31.
Nash
,
A. L.
, and
Jain
,
N.
,
2019
, “
Dynamic Design Optimization for Thermal Management: A Case Study on Shell and Tube Heat Exchangers
,”
ASME 2019 Dynamic Systems and Control Conference
,
Park City, UT
,
Oct. 8–11
,
American Society of Mechanical Engineers Digital Collection
, p. V002T23A004.
32.
Incropera
,
F. P.
,
Lavine
,
A. S.
,
Bergman
,
T. L.
, and
DeWitt
,
D. P.
,
2007
,
Fundamentals of Heat and Mass Transfer
,
Wiley
,
Hoboken, NJ
.
33.
,
S.
, and
Postlethwaite
,
I.
,
2005
,
Multivariable Feedback Control Analysis and Design
, 2nd ed.,
Wiley
,
West Sussex, England
.
34.
Currie
,
J.
, and
Wilson
,
D. I.
,
2012
, “
OPTI: Lowering the Barrier Between Open Source Optimizers and the Industrial MATLAB User
,”
Foundations of Computer-Aided Process Operations
,
N.
Sahinidis
, and
J.
Pinto
, eds.