Economically Efficient Combined Plant and Controller Design Using Batch Bayesian Optimization: Mathematical Framework and Airborne Wind Energy Case Study

Combined Plant and Contr oller Design Using Batch Bayesian Optimization: A Case Study in Airborne W ind Ener gy Systems Ali Baheri Chris V ermillion Uni versity of North Carolina at Charlotte akhayatb@uncc.edu cvermill@uncc.edu Abstract W e present a nov el data-dri ven nested optimization framew ork that addresses the problem of coupling between plant and controller optimization. This optimization strat- egy is tailored towards instances where a closed-form ex- pression for the system dynamic response is unobtainable and simulations or experiments are necessary . Specifically , Bayesian Optimization, which is a data-dri ven technique for finding the optimum of an unknown and expensi ve-to- ev aluate objecti v e function, is employed to solve a nested optimization problem. The underlying objectiv e function is modeled by a Gaussian Process (GP); then, Bayesian Optimization utilizes the predicti v e uncertainty information from the GP to determine the best subsequent control or plant parameters. The proposed frame work dif fers from the majority of co-design literature where there exists a closed- form model of the system dynamics. Furthermore, we uti- lize the idea of Batch Bayesian Optimization at the plant optimization level to generate a set of plant designs at each iteration of the overall optimization process, recognizing that there will exist economies of scale in running multiple experiments in each iteration of the plant design process. W e validate the proposed frame work for Altaeros' Buoyant Airborne T urbine (BA T). W e choose the horizontal stabi- lizer area, longitudinal center of mass relativ e to center of buo yancy (plant parameters), and the pitch angle set-point (controller parameter) as our decision variables. Our results demonstrate that these plant and control parameters con- ver ge to their respecti ve optimal values within only a few iterations. Figure 1: Altaeros Buoyant Airborne T urbine (B A T), Image Credit: [1] 1 Introduction Airborne W ind Ener gy (A WE) systems are a ne w paradigm for wind turbines in which the structural elements of con ventional wind turbines are replaced with tethers and a lifting body (a kite, rigid wing, or aerostat) to harvest wind power from significantly increased altitudes (typically up to 600m). At those altitudes, winds are stronger and more con- sistent than ground-lev el winds. The v ast energy resource from high-altitude winds has attracted the attention of numerous research and commer- cial ventures over the past two decades ([1, 2, 3, 4, 5]). T o- date, many organizations in the A WE community ha ve fo- cused on optimizing operating altitude ([6, 7]) and cross- wind motion to maximize power output, while a limited number of studies have focused on combined plant and con- trol system designs, which hav e been sho wn in [8] to be coupled. Combined plant and controller problems consist of those scenarios in which the optimal controller depends on the physical system design (i.e., the plant) and vice versa. Com- bined plant and controller optimization (often termed co- design ) has been employed for a wide variety of systems, in- cluding automoti ve suspension systems (see [9], and [10]), elev ator systems [11], and the A WE application at hand (with initial studies reported in [8], and [12]). Broadly speaking, combined plant and controller strategies fall into four different cate gories: • Sequential : A sequential strategy completes the plant and controller optimization problem in successive or- der . For instance, authors in [13] use control proxy functions to augment the objective function of the plant, thereby separating the plant and and controller optimization problem into two sub-problems. • Iterative : The iterative approach fully optimizes the plant design for a giv en controller , then optimizes the controller design for a fixed plant, then repeats the cy- cle (see [14], [15]). • Nested : A nested optimization approach contains two loops: an inner loop that completes a full optimization of the controller and an outer loop that completes an it- eration of the plant optimization (see [8], and [9], [11], and [16]). • Simultaneous : In a simultaneous optimization strat- egy , both the plant and controller optimization prob- lems are carried out at the same time [17], [18]. An efficient decomposition-based variant of simultaneous optimization is proposed in [19]. Among these techniques for solving the co-design prob- lem, the nested co-design approach is unique in its abil- ity to le ver age critical dif ferences in control parameters (which can be modified during e xperiments) and plant pa- rameters (which, when experimental work is required, can generally only be modified between experiments). This sort of approach can therefore be extremely beneficial in complex systems where experiments will ultimately be re- quired. Howe ver , existing literature on nested co-design, which is summarized in T able 1, makes numerous simpli- fying assumptions (linearity , full state measurement, etc.). The tools used for the plant adjustment are typically local in nature and often unsuitable for complex systems. Fur- thermore, most existing literature on nested co-design uses continuous-time optimal control design techniques (such as LQR) for the controller design ([20]), without taking advan- tage of online adaptation capabilities. Recently , authors in [12] proposed a nested co-design approach in which the control parameter is adjusted dur- ing a simulation/experiment and plant parameters are opti- mized between simulations/e xperiments. The benefit here is that when the design optimization process inv olv es lengthy simulations or experiments, the ability to adjust controller parameters during the simulations/experiments can signifi- cantly reduce the time and cost of the optimization process. The framework in [12] used G-optimal design of Experi- ments (DoE) to select a batch of candidate plant parameters that cover the design space while maximizing a statistical information metric. Furthermore, extremum seeking (ES) control was utilized to adjust the control parameter(s) in real-time ov er the course of experiments/simulations. The authors recently extended their w ork in [22] to consider on- line adaptation mechanisms that are based on a global statis- tical information metric. While resolving some challenges, these approaches suffer from tw o main drawbacks: 1. Populating the design space with candidate points in order to merely gain the most information about the design space may lead to the ev aluation of candidate designs that have very poor associated performance. This will lead to significant effort expanded in charac- terizing re gions of the design space that are unlik ely to yield optimal design parameters. 2. ES only achieves local, rather than global optimality . One of the attracti ve problems in both the control and machine learning communities is optimizing system de- signs for real-world applications using scar ce data . This problem has been studied in the context of sequential decision-making problems aiming to learn the beha vior of an objectiv e function (called e xploration) while simultane- ously trying to maximize or minimize the objective (called exploitation). As an efficient and systematic approach for balancing exploration and e xploration in a partially observ- able en vironment, Bayesian Optimization has been applied to v arious real-world problems [23, 21, 6, 24, 25, 26, 27, 28]. In general, Bayesian Optimization aims to find the global optimum of an unknown, e xpensi ve-to-e v aluate, and black-box function within only a fe w e v aluations of the ob- jectiv e function at hand. The work presented here is an extension of the recent DSCC conference paper , [21], where we introduced a nov el machine learning variant of nested co-design frame work. In this frame w ork, Bayesian Optimization is used at both the plant and control design le vels of a nested co-design frame- work, in order to efficiently explore the design space. This framew ork is applicable to a v ast array of problems where there is no analytical expression for the performance metric in terms of decision variables. Furthermore, the proposed framew ork can easily support experiments in which control parameters are adapted during experiments. In this work, we extend our results from [21] in the follo wing ways: • W e employ Batch Bayesian Optimization at the plant optimization level to generate a set of plant designs at 2 T able 1: K e y existing literature on nested co-design A uthor Basic idea Optimization technique Restriction on system dynamic and objective function Fath y , et.al (2003) [20] Outer loop seeks to optimize the ov erall performance by varying the plant design. The inner loop generates the optimal controller for each plant parameter selected by the outer loop. Plant: Interior point method Controller: LQR design - The structure of system dynamics are known. - System of interest is L TI. Desse, et.al (2017) [12] The co-design problem is solved using optimal design of experiment for the plant optimization and extremum seeking for the control system optimization. Plant: G-optimal DoE Controller: extremum seeking - Extremum seeking may get stuck in a local optima. Baheri, et.al (2017) [21] A machine learning variant of nested co-design. Uses Bayesian Optimization as a global optimization tool for the black-box objectiv e function at both lev els. Plant and Controller: Bayesian Optimization - No restriction on the system or objectiv e function, b ut the framew ork requires simulation and/or experiments. each iteration of ov erall optimization process, recog- nizing that there will exist economies of scale in run- ning multiple experiments at once. • W e conduct a detailed economic assessment to quan- tify the cost associated with batch and non-batch Bayesian Optimization algorithms. T o v alidate our approach, we focus on the Altaeros Buoyant Airborne T urbine (B A T) (See Fig. 1). Among pa- rameters critical to system performance, we focus on hori- zontal stabilizer area and longitudinal center of mass rela- tiv e to center of buo yancy as plant parameters. W e take the trim pitch angle as the controller parameter . The rest of this paper is structured as follows: First, we present the fused plant and control optimization methodology that is used in the nested optimization frame work. Next, we summarize the plant and controller design optimization of the B A T . Fi- nally , we provide detailed results, demonstrating the effec- tiv eness of the proposed approach. 2 FUSED PLANT AND CONTR OLLER METHODOLOGY 2.1 PR OBLEM FORMULA TION The ultimate goal of this study is to solv e the following optimization problem: minimize p p , p c J ( p p , p c ) = Z T final 0 g  x ( t ); p p , p c  dt (1) subject to: ˙ x = f ( x , u , d ; p p , p c ) (2) p p ∈ P, p c ∈ C (3) where (1) describes the integral cost functional, (2) repre- sents a general dynamic model that governs the system dy- namic, and (3) presents hard constraints on plant and control parameters, which are denoted by p p and p c , respectiv ely . d is the state disturbance vector . This work will focus on co-design processes that are carried out in controlled en- vironments, where the en vironmental perturbation (mani- fested by d ) is consistent between simulations/experiments but nonetheless important. The nested optimization frame work consists of tw o main loops: • Outer loop : Plant parameters ( p p ) are adjusted in a direction that is chosen by one iteration of our chosen optimization tool (i.e., Bayesian Optimization). • Inner loop : For the candidate set of plant parameters generated by the outer loop, the inner loop completes a full optimization of the control parameter vector . Fig. 2 presents the proposed nested co-design frame work. A description of each variable within the general process is provided in T able 2. The process continues until con ver - gence to optimal parameters. Algorithm 1 also summarizes the process proposed in this work. 2.2 B A YESIAN OPTIMIZA TION The final goal is to optimize the plant and control pa- rameters, (i.e.,  p p , p c  ), according to some performance index, J . Bayesian Optimization is our chosen tool to solv e the aforementioned nested optimization problem. Specifi- cally , Bayesian Optimization aims to adjust p p in the outer loop, while it seeks the optimum of p c in the inner loop. T o illustrate the process of Bayesian Optimization, we take p as the general optimization variable in the rest of this section. Bayesian Optimization traditionally seeks to maximize the 3 T able 2: Description of k ey v ariables in volved in the combined plant and controller optimization V ariable Description p p Plant parameter(s) p c Control parameter(s) p ∗ c ( p p ) Optimal control parameter(s) for a candidate plant design B j Elements (candidate plant designs) in a batch at iteration j J ( p ∗ c ( p p ) , p p ) Integral performance v alue while operating at the optimal control parameter(s) and candidate plant design Bayesian Optimization Modify plant via Batch Bayesian Optimization j = j +1 i = i +1 Simulation/ Experiment p ⇤ c ( p p,j ) i = i +1 p c,i J ( p p,j , p c,i ) Continuous Controller Optimization (Section I I  D) (Section II  C) B j Figure 2: Machine learning v ariant of nested plant and con- troller co-design using Batch Bayesian Optimization value of some performance index. T o keep our exposition aligned with this traditional implementation, we define: R , − J ( p ) (4) as the rew ard to be maximized. Bayesian Optimization seeks to globally minimize an unknown, black-box, and expensi ve-to-ev aluate objectiv e function. Due to the computational costs associated with ev aluating this function and the possibility of the optimiza- tion getting stuck in a local optimum, it is crucial to select the location of each ne w ev aluation deliberately . T o ad- dress this, Bayesian Optimization is employed to find the global optimum of the objective function within few eval- uations on the real system. Broadly speaking, there are two main phases in Bayesian Optimization. First, the un- derlying function is modeled as a Gaussian Process (GP). Second, we choose an acquisition function that guides the optimization by determining the next point to e valuate. 2.2.1 Learning phase using Gaussian Processes (GPs) In this section, we introduce the basic properties of GPs. GPs are used to model the objective function of the main Algorithm 1 1: procedure F U S E D P L A N T A ND C O N T RO L L E R W I T H B O 2: while plant parameters not con verged do 3: Run one iteration of Bayesian Optimization in outer loop 4: for plant candidate parameters do 5: while control parameter not con verged do 6: Run full Bayesian Optimization in inner loop 7: Update control parameter optimization. As an attractiv e choice for non-parametric regression in machine learning, the ultimate goal of GPs is to find an ap- proximation of a nonlinear map, from an input vector to the function value. In general, GP models are able to char- acterize complex phenomena since they are able to handle the nonlinear ef fects and interaction between covariates in a systematic fashion. In general, a GP is fully specified by its mean function, µ ( p ) , and cov ariance function, k ( p , p 0 ) : R ( p ) ∼ G P  µ ( p ) , k ( p , p 0 )  , (5) The GP frame work is used to predict the inte gral cost func- tion , R ( p ) , at an arbitrary optimization variable(s) (i.e., plant or controller parameters), P =  p 1 , · · · , p t  , based on a set of t past observations, D 1: t =  p 1: t , R 1: t ( p )  . The function value, represented by R ( p ∗ ) in Eq. (6), for an un- observed input p ∗ , and the observed function values, is as- sumed to follow a multi variate Gaussian distrib ution [29]:  y t R ( p ∗ )  ∼ N  0 ,  K t + σ 2  I t k t k T t k ( p ∗ , p ∗ )   (6) where y t =  R ( p 1 ) , · · · , R ( p t )  is the vector of observed function values. The vector k t ( p ) = [ k ( p , p 1 ) , · · · , k ( p , p t )] encodes the relationship between the new input, p , and the past data points in P . The co- variance matrix has entries [ K t ] ( i,j ) = k ( p i , p j ) for i, j ∈ 4  1 , · · · , t  . The identity matrix is represented by I t and σ  represents the noise variance [29]. The entries in the covariance matrix are parameterized based on a kernel function. A common choice of kernel function is the squared exponential (SE) function, whose ev aluation between two inputs, p i and p j , is expressed as: k ( p i , p j ) = σ 2 0 exp  − 1 2 ( p i − p j ) T Λ − 2 ( p i − p j )  , (7) where γ =  σ 0 , Λ  are the hyper -parameters of the kernel function. W e select the hyper-parameters as the ones that maximize the marginal log-lik elihood of D [29]: γ ∗ = arg max log θ p ( y t | P , γ ) , (8) where log p ( y t | P , γ ) =  − 1 2 y T t K − 1 y t − 1 2 log | K | − t 2 log 2 π  . (9) Once the hyper-parameters are identified, the predictive mean and v ariance at p ∗ , conditioned on these past obser- vations, are e xpressed as: µ t ( p ∗ | D ) = k t ( p )  K t + I t σ 2   − 1 y T t , (10) σ 2 t ( p ∗ | D ) = k ( p , p ) − k t ( p )  K t + I t σ 2   − 1 k T t ( p ) , (11) In summary , the learning phase using GPs in volves two main steps: training and prediction. The training step con- sists of finding proper mean and cov ariance functions, as well as optimized hyper-parameters, in light of the data (Eq. 8). The prediction phase characterizes the objec- tiv e function value at an unobserv ed input in a probabilistic framew ork (Eqs. 10-11) [29]. These two equations implic- itly serve as surrogates for our unknown function and are used in the next phase to calculate the acquisition function (which is used to determine which point to ev aluate next). 2.2.2 Optimization phase As mentioned earlier, we need to choose an acquisition function, which guides the optimization by determining the next point to e valuate. Specifically , the acquisition function uses the predictiv e mean and variance (Eqs. 10-11) to bal- ance exploration of high-variance regions with exploitation of high-mean regions. The ultimate goal of Bayesian Optimization is to deter- mine the next optimization variable(s) based on past obser- vations. Among se veral choices of acquisition functions, we use an acquisition function belonging to the improv ement- based family [30]. More precisely , the next operating point is selected as the one that maximizes the expected improv e- ment: p t +1 = arg max p EI ( p t +1 | D 1: t ) (12) where: EI ( p t +1 | D 1: t ) = E  max  0 , R t +1 ( p ) − R ( p ) max  | D 1: t  . (13) Here, max  0 , R t +1 ( p ) − R ( p ) max )  represents the im- pr ovement to ward the best value of the objectiv e function so far , R ( p ) max . The improvement is positi ve when the predic- tion is higher than the best value of the objective function so far . Otherwise, it is set to zero. The inability of the acquisi- tion function to assume negati ve v alues reflects the fact that if the design worsens from one iteration to the next, then it is possible to simply rev ert to the pre vious best design. For - tunately , there exists a closed form expression for expected improv ement, given by [31]: EI ( p t +1 | D 1: t ) = (  µ t ( p ) − R ( p ) max  Φ( Z ) + σ t ( p ) φ ( Z ) , σ t ( p ) > 0 0 , σ t ( p ) = 0 (14) where Z = µ t ( p ) − R ( p ) max σ t ( p ) , (15) and Φ( . ) and φ ( . ) denote the cumulative and probability density function for the normal distribution, respecti vely . W e summarize the generic version of Bayesian Opti- mization in Algorithm 2. The algorithm is initialized by two pre viously-ev aluated optimization variables and corre- sponding costs (line 2). Then, at each step, a GP model is trained (line 4) to compute the predictive mean and v ari- ance (line 5). These statistical quantities are used to con- struct the acquisition function (line 6). Next, the point that maximizes the acquisition function is selected as the next candidate (line 7). At the next optimization instance, this point is added to the historical data (line 8), and the process repeats. 2.3 Bayesian Optimization-Based Co-Design: Batch Plant Design The algorithm introduced in our prior conference publi- cation [21] utilizes Bayesian Optimization to identify just one plant design candidate at each iteration of the outer loop. In this section, we examine the use of Batch Bayesian Optimization at the plant optimization level to generate a set of plant designs at each iteration of the o verall optimization 5 Algorithm 2 Bayesian Optimization (BO) 1: procedure G E N E R I C B AY E S I A N O P T I M I Z A T I O N 2: D ← Initialize :  p 1:2 , R ( p 1:2 )  3: for each iteration do 4: T rain a GP model from D 5: Compute mean and variance of GP  Eqs. (10- 11)  6: Compute acquisition function  Eq. (14)  7: Find p ∗ that optimizes acquisition function 8: Append  p ∗ , R ( p )  to D process, recognizing that there will exist economies of scale in running multiple experiments at once. T o mathematically introduce the idea of Batch Bayesian Optimization, we slightly modify the notation of acquisition function presented so far . Specifically , we take α ( p , I t,k ) as the acquisition function, where I represents the a vailable data set, D , plus the GP structure when n data points are av ailable. Consequently , subscripts t and k represent the k th element of t th batch, respectfully . The Local Penalization (LP) algorithm originally pre- sented in [32] is a heuristic approach for Batch Bayesian Optimization that works by iteratively penalizing the cur- rent peak in the acquisition function to find the next peak. Consequently , each successi ve element of the batch is cho- sen to maximize the modified acquisition function, which has been modified based on each pre vious element of the batch. According to the LP algorithm, e very element in t th batch is giv en by: p t,k = arg max p ∈ P n g  α ( p , I t, 0 )  Q k − 1 j =1 ϕ ( p ; p t,j ) o (16) where g is a differentiable transformation of α that keeps the acquisition function positiv e. ϕ ( p ; p t,j ) is the core com- ponent of LP algorithm and called the local penalizer cen- tered at p j (in the t th batch). If selected properly , the LP estimates the distance from p j to the true optimum of the cost function. If we believ e that the distance from p j is far from the true optimum, then a lar ge penalizer can discard a large portion of the feasible domain that should not be con- sidered in selecting one of the batch elements. On the other hand, if we belie ve that p j is close to the true optimum, a small penalizer keeps collecting elements in a close neigh- borhood of true optimum. Howev er , the main challenge lies in selecting an appropriate local penalizer . 2.3.1 Selecting a local penalizer As mentioned earlier , the local penalizer characterizes a be- lief about the distance from the batch locations to true opti- mum. Let R M = max p ∈ P R ( p ) . Consider the ball: B r j ( p j ) =  p ∈ P :   p j − p   ≤ r j  (17) where r j = R M − R ( p j ) L (18) and L is a valid Lipschitz constant of R . If R ( p ) ∼ G P  µ ( p ) , k ( p , p 0 )  , then we can choose a local penalizer , ϕ ( p ; p j ) , as the probability that p , an y point in P that is a potential candidate for batch elements, does not belong to B r j ( p j ) : ϕ ( p ; p j ) = 1 − p  p ∈ B r j ( p j )  (19) Howe ver , we need an analytical expression for ϕ ( p ; p j ) to compute each batch element. Proposition 1 provides an an- alytical form for the local penalizer [32]: Proposition 1. If R ( p ) ∼ G P  µ ( p ) , k ( p , p 0 )  , then ϕ ( p ; p j ) as defined in (18), is a local penalizer at p j such that: ϕ ( p , p j ) = 1 2 erfc( − z ) (20) wher e z = 1 p 2 σ 2 n ( p )  L   p − p j   − R M + µ n ( p j )  . (21) and erfc is the complementary err or function . Proposition 1 implies that if µ n ( p j ) is close to R M , then ϕ ( p , p j ) will hav e a more localized and small impact on α ( p ) . On the other hand, if µ n ( p j ) is far from R M , then ϕ ( p ; p j ) will hav e a big impact on α ( p ) . 2.3.2 Selecting the parameters L and R M An appropriate local penalizer relies on v alid choices for R M and L . Howe ver , the value of R M and L are unknown in general. T o approximate R M , one can take ˆ R M = max i  R i  (22) Regarding the parameter L we take ˆ L = max p ∈ P k ∇ µ ( p ) k (23) as a valid Lipschitz constant [32]. Finally , algorithm 3 summarizes the procedure presented for Batch Bayesian Optimization. 2.4 Bayesian Optimization-Based Co-Design: Continuous Contr oller Optimization Bayesian Optimization is an iteration-based optimization algorithm. Thus, in order to use Bayesian Optimization as a tool for adaptation of controller parameters in a continuous- time setting, it is crucial to make a clear connection be- tween the concept of discr ete iterations (used in Bayesian Optimization), and continuous time over the course of a 6 Algorithm 3 1: procedure B AT C H B AY E S I A N O P T I M I Z A T I O N W I T H L P 2: Inputs: batch size: n b 3: for t = 1 until con vergence do 4: Fit a GP to D t 5: Build the acquisition α ( p , I t, 0 ) using the cur- rent GP 6: ˜ α t, 0 ( p ) ← g  α ( p , I t, 0 )  7: ˆ L ← max p ∈ P k ∇ µ ( p ) k 8: for j = 1 to n b do 9: p t,j ← arg max ˜ α t,j − 1 ( p ) 10: ˜ α t,j ( p ) ← ˜ α t, 0 ( p ) Q k − 1 j =1 ϕ ( p ; p t,j ) 11: B n b t ←  p t, 1 , . . . , p t,n b  12: R t, 1 , . . . , R t,n b ← ev aluation of R at B n b t 13: D t +1 ← D t ∪  p t,j , R t,j  n b j =1 14: Fit GP to D n simulation/experiment. Each iteration in the Algorithm 1 corresponds to one window of time within the o verall sys- tem simulation. W e divide each of these windows into three different phases (See Fig. 3). Phase one represents the set- tling period. T o av oid unfairly introducing system transients from one iteration to another in our cost function calcula- tion, we allow the system to settle during the first period, then only compute the cost function v alue based on the sec- ond phase, which we refer to as the performance period. Fi- nally , the calculation of the subsequent decision v ariable(s) occurs in the third phase, using Bayesian Optimization. 2.5 Con vergence detection In any optimization problem, the stopping criteria is de- termined based on either a fixed iteration budget or certain con vergence criteria is met. Detecting con ver gence, rather than relying on a fixed iteration budget, is critical partic- ularly where experiments come into play , since ev ery ex- periment requires time and money . In this work we set the following stopping criteria for con vergence: | R ( p i ) − R ( p i − j ) | <  j = 1 , . . . , n (24) W e set n equal to 2 in this work. 3 PLANT AND CONTROLLER DESIGN OPTIMIZA TION of AL T AER OS B A T 3.1 Plant Model of Altaeros B A T The Altaeros B A T ([1]) features a horizontal axis turbine that is suspended within a buo yant shell. Unlike many A WE Iteration 1 Iteration 2 g ( x ; p p , p c ) t T Settle T Per f . T Calc . Figure 3: Clarification of iteration vs. time in the proposed framew ork concepts, the BA T is designed to remain substantially sta- tionary and passiv ely align with the wind. By accomplish- ing this, the B A T can achie ve secondary objectiv es, such as telecommunications, which Altaeros has publicly indicated interest in. T o that end, it is of great interest to design the combined plant and control system to achieve the steadiest possible flight under atmospheric disturbances, which is the focus of our case study . The dynamic model of Altaeros B A T system was orig- inally introduced in [33]. W e briefly introduce some of the important features of this model in this section. Fig. 4 shows the v ariables used in the dynamic model. The model, which was deriv ed using an Euler-Lagrange framew ork, describes the position and orientation of the B A T through six generalized coordinates: Θ , φ , ψ , L t , θ 0 , and φ 0 . Three of these generalized coordinates represent uncontrolled angles. Specifically , Θ is the azimuth angle (angle of the tether projection on the horizontal plane), Φ is the zenith angle (angle of tether with respect to vertical), and Ψ is the twist angle (about the tether axis). T able 3 represents the full state variables of A WE system. The three tethers are modeled as a single tether with length L t with a bridle joint at the top that splits in three attachment points on the BA T . The bridle joint is modeled through two angles, φ 0 and θ 0 , which are referred to as in- duced roll and induced pitch, respecti vely . The single tether approximation removes algebraic constraint equations, al- lowing the system to be described by ordinary dif ferential equations (ODEs) rather than differential algebraic equa- tions (D AEs). The center of mass location is modeled as a function of Φ (zenith angle), Θ (azimuth angle), and L t (av erage tether length). The induced angles are related to the tether lengths ( l 1 , l 2 , and l 3 ) through the following ex- 7  ⇥ l 2 l 1 l 3 x g z g y g y z ,r  ,p ✓ ,q x  0 ✓ 0 Figure 4: Ground-fixed and body-fix ed coordinates plus the key v ariables used in deriving Euler -Lagrangian dynamics. pressions: φ 0 = tan − 1  l 3 − l 2 l lat sep  (25) θ 0 = tan − 1  l 1 − 0 . 5( l 2 + l 3 ) l long sep  (26) where l long sep and l lat sep are longitudinal and lateral tether at- tachment separation distances, respectiv ely . l 1 , l 2 and l 3 are the distances between the bridal joint and the three tether attachment points. The control inputs are the tether release speeds, ¯ u i , which are giv en by , ¯ u i = d dt l i (27) Ultimately , the governing system equations gi ven below are deri ved using an Euler-Lagrange formulation and are expressed by: D ( Q ) ¨ Q + C ( Q, ˙ Q ) ˙ Q + g ( Q ) = τ  Q, ˙ Q, V wind , ψ wind  , (28) X = f ( Q, ˙ Q ) , (29) Ω = g ( Q, ˙ Q ) , (30) where: Q = [Φ Θ Ψ L t θ 0 φ 0 ] (31) X = [ x y z u v w ] (32) Ω = [ φ θ ψ p q r ] (33) Here, V wind is the wind speed, and ψ wind and τ represent the wind direction and vector of generalized aerodynamic forces and moments, respectively . Aerodynamic forces and moments are functions of α (angle of attack) and β (side slip angle), which describe the orientation of the apparent T able 3: Full state v ariables of an A WE system State variable Notation Zenith angle Φ Azimuth angle Θ T wist angle Ψ Zenith angle rate ˙ Φ Azimuth angle rate ˙ Θ T wist angle rate ˙ Ψ Unstreched tether length l Induced roll φ 0 Induced pitch θ 0 wind vector with respect to the body-fixed coordinates of the B A T . X and Ω represent the translational and rotational dynamics, respectively . Since we treat the horizontal stabi- lizer area, A H , as a design parameter to be optimized, the aerodynamic coefficients are modeled as explicit functions of the stabilizer areas as follows: C D,L,S ( α, β ) = C F D,L,S ( α, β ) + C H D,L,S ( α, β ) A H A ref + C V D,L,S ( α, β ) A V A ref (34) C M x ,M y ,M z ( α, β ) = C F M x ,M y ,M z ( α, β ) + C H M x ,M y ,M z ( α, β ) A H l H A ref l ref + C V M x ,M y ,M z ( α, β ) A V l V A ref l ref (35) Here, C D , C L , and C S represent the drag, lift, and side force coefficients, whereas C M x , C M y , and C M z represent the roll, pitch, and yaw moment coef ficients. Closed-Loop Contr oller The controller is designed to track three different set points, namely altitude ( z sp ), pitch ( θ sp ) and roll ( φ sp ). Due to the symmetrical configuration, φ sp is always set to zero. The block diagram of the controller is sho wn in Fig. 5. Bayesian Optimization is used to update θ sp , which is com- puted internally on the flight computer and is not a user- specified value (hence, θ sp represents a control parameter). T racking is achieved through the use of three lead filtered PD controllers, which indi vidually control altitude, pitch angle, and roll angle. The outputs of these three lead filters are denoted, respectiv ely , by ¯ v z , ¯ v θ , and ¯ v φ . These outputs represent the av erage tether release speed (used to regulate altitude), forward/aft tether speed difference (used to regu- late pitch), and port/starboard tether speed difference (used 8 ✓ sp z sp Bayesian Optimization Position, velocity feedback  sp =0 Altitude/attitude control ¯ u fwd ¯ u aft , fwd ¯ u aft , p ort J ⇣ p c ( p p ) , p p ⌘ Figure 5: Block diagram of closed loop flight controller for the B A T . z sp denotes a constant altitude set-point. W e choose p c = θ sp and p p = [ x cm − x cb A H ] T in our case study results. to regulate roll), respecti vely . These controller outputs are related to the tracking errors in altitude ( z e ), pitch ( θ e ) and roll ( φ e ) through: ¯ v z ( s ) = k z d s + k z p τ z s + 1 z e ( s ) , (36) ¯ v θ ( s ) = k θ d s + k θ p τ θ s + 1 θ e ( s ) , (37) ¯ v φ ( s ) = k φ d s + k φ p τ φ s + 1 φ e ( s ) . (38) The tether release speeds ¯ u center , ¯ u stbd , and ¯ u port serve as control inputs to three motors and are related to ¯ v z , ¯ v θ , and ¯ v φ through:   ¯ u center ¯ u stbd ¯ u port   =   1 − 1 0 1 1 1 1 1 − 1     ¯ v z ¯ v θ ¯ v φ   . (39) For our case study results, we focus on the follo wing plant parameters to be optimized: p p =  x cm − x cb , A H  (40) where x cm − x cb describes the longitudinal center of mass position relativ e to the center of buo yancy and A H repre- sents the horizontal stabilizer area. The controller parameter used in this work is gi ven by: p c = θ sp (41) where θ sp presents the trim pitch angle (which is pro- grammed on the flight computer and therefore constitutes a control parameter rather than a user input). 3.2 PERFORMANCE INDEX The performance index, to be minimized, tak es into ac- count two main system properties: 1. Gr ound footprint : It is of interest to minimize the land usage requirements for multiple systems. For this pur - pose, the horizontal projected area of land that BA T cov ers is used as a criterion for quantifying the ground footprint. This area is represented by A = π l 2 sin 2 Φ . Thus, as Φ decreases, the projected area decreases. 2. Quality of flight : A low v alue of heading angle typi- cally corresponds to few oscillations in the system and desirable direct-downwind operation. Furthermore, since we are focused on steady , level flight, we de- sire to have the B A T as stationary as possible. T o characterize the de gree to which we accomplish this goal, we penalize heading and roll angle tracking error ( ψ − ψ flow and φ − φ sp , respectively) in our perfor- mance metric. Ultimately , the performance index is denoted by: J  p c ( p p ) , p p  = Z T f T i  k 1 Φ 2 + k 2 ( ψ − ψ flow ) 2 + k 3 ( φ − φ sp ) 2 )  dt, (42) 3.3 SIMULA TION SETUP T o excite the system with a wind en vironment that is consistent across simulations, we implement a frequency approximation of vortex shedding of flow o ver a cylinder . The model perturbation is based on spectral analysis of flow ov er a cylinder in [34]. Each component of the velocity is approximated as a sinusoidal perturbation about a mean v e- locity found in [34]. Each of the velocity components is giv en by: v x = v base x + v x 0 sin( ω dist t ) (43) v y = v y 0 sin( ω dist t ) (44) v z = v z 0 sin( ω dist t ) (45) where v base = 0 . 606 m s , v x 0 = 0 . 0866 m s , v y 0 = 0 . 065 m s , v z 0 = 0 . 0087 m s , and ω dist = 1 H Z . This mechanism for perturbing the system is attractiv e because it (i) excites key lateral dynamics and (ii) can be implemented in later lab-scale experimental co-design us- ing the team's water channel setup described in [35] and pictured in Fig. 11. 4 RESUL TS W e ev aluated the proposed algorithm on the B A T nu- merical model. W e assessed the effecti veness of Batch Bayesian Optimization algorithm for 3 and 4 elements 9 2 4 6 8 Number of Iterations 1 1.5 2 2.5 3 Horizontal Stabilizer Area [m 2 ] 10 -3 2 4 6 8 Number of Iterations 0 2 4 6 Relative Center of Mass [m] 10 -3 2 4 6 8 Number of Iterations 0.05 0.1 0.15 0.2 0.25 0.3 Optimal Trim Pitch Angle [rad] 2 4 6 8 Number of Iterations 180 200 220 240 260 280 300 Integral Cost Function 1 2 3 4 5 6 Number of Iterations 1 1.5 2 2.5 3 Horizontal Stabilizer Area [m 2 ] 10 -3 1 2 3 4 5 6 Number of Iterations 3 4 5 6 Relative Center of Mass [m] 10 -3 1 2 3 4 5 6 Number of Iterations 0 0.1 0.2 0.3 Optimal Trim Pitch Angle [rad] 1 2 3 4 5 6 Number of Iterations 200 250 300 350 Integral Cost Function First element Second element Third element 1 2 3 4 5 Number of Iterations 1.5 2 2.5 3 Horizontal Stabilizer Area [m 2 ] 10 -3 1 2 3 4 5 Number of Iterations 1 2 3 4 5 6 Relative Center of Mass [m] 10 -3 1 2 3 4 5 Number of Iterations 0.21 0.22 0.23 0.24 Optimal Trim Pitch Angle [rad] 1 2 3 4 5 Number of Iterations 180 200 220 240 260 Integral Cost Function First element Second element Third element Forth element Figure 6: Conv ergence of plant parameters, control parameter , and integral cost function for 1, 3, and 4 batch sizes (from upper left to lower right) (plant designs) in each batch. W e also compare results from Batch Bayesian Optimization against results for a batch size of 1 (i.e., generic Bayesian Optimization). Fig. 6 shows the conv ergence of plant parameters (hor- izontal stabilizer area and relativ e center of mass), control parameter (trim pitch angle), and integral of cost function for different batch sizes. Each line structure in this figure represents the ev olution of a different batch element o ver the course of optimization. One can immediately see that fewer total iterations are required to con verge with larger batch sizes. Furthermore, as can be seen from this figure, the plant parameters con- ver ge after only a small number of iterations (each iteration corresponds to one round of system performance ev alua- tion). With the same simulation setup, comparing these re- sults with those reported in [12] rev eals that Bayesian Op- timization leads to a faster conv ergence than the optimal DoE proposed in that work. It should be emphasized that plant design changes are much more expensi ve compared to adjusting control parameters in instances when simula- tions are replaced with experiments (which is a long-term goal of the present work). Therefore, reducing the number of required plant reconfigurations in the design optimization process is a very important goal. Fig. 6 also illustrates the optimal control parameter in the inner loop for each plant design generated by the outer loop. Furthermore, Fig. 7 represents the sample ev olution of control parameter (trim pitch angle) and actual pitch an- gle in the inner loop ov er a simulation. Figs. 8, 9, and 10 illustrate the e volution of the individual instantaneous cost function components (i.e., roll tracking error , heading tracking error , and zenith angle) before and after the iterative optimization process to demonstrate the effecti veness of the proposed approach. 10 Figure 7: Sample ev olution of control parameter (i.e., trim pitch angle) and actual pitch angle in inner loop over the time 5 Economies of Scale So far , we ha ve explored the idea of Batch Bayesian Op- timization at the plant optimization lev el to generate a set of plant designs at each iteration of the o verall optimization process. In this section, we will assess, for the A WE system, whether this will result in economies of scale. T o assess the economies of scale that are introduced through a Batch Bayesian Optimization process, we focus our attention on a lab-scale e xperimental platform (depicted in Fig. 11) for closed-loop flight characterization of A WE systems, which is detailed in [35]. W ith this system, 3D printed models (depicted in Fig. 12) of A WE system lift- ing bodies are tethered and flo wn in the UNC-Charlotte wa- ter channel. Micro DC motors are used to regulate tether lengths, high-speed cameras are used for image capture, and a high-performance target computer is used for real-time image processing and control. Because this experimental platform represents the ultimate “end game” for the Batch Bayesian Optimization approach, we use it as the basis for the economies of scale analysis presented herein. In order to run an experimental batch (and non-batch) Bayesian Optimization in water channel, the following steps need to happen: • Prior to running any experiments: 3D print the main body . The cost associated with this step is giv en by: C 3 D =  m × c eng + c recharge  N T print (46) where c eng and c recharge represent the cost required for engineer(s) to conduct experiments and the f acility charge, respecti vely . T print denotes the time required to print the main body . m is the number of employees required to print the model (and fins). • After each batch is chosen, the following needs to hap- pen once per batch: (i) 3D print N sets of fins. This cost is gi ven by: C 3 D Fins =  m × c eng + c recharge  N T 3 D Fins (47) where T 3 D Fins represents the time required to print the fins. (ii) Schedule time in the water channel (this time may include lead time.) C lead = c lostTime T lead (48) where c lostTime and T lead represent the opportunity cost of lost time and the lead time, respectiv ely . (iii) Initialize the water channel: C Wchannel =  m 0 × c eng + c W recharge  T setup + N T exp + N T reconfig  (49) where T setup and T reconfig denote the time required to set up the experiment in the water channel and the time required to reconfigure the model, respectively . Fur - thermore, c W recharge and T exp represent the cost re- quired to run the water and time required to conduct an e xperiment in the inner loop, respecti vely . m 0 is the number of employees required to conduct experiments. • Between each experiment within a batch, the follo wing needs to happen N times per batch: (i) Sw ap out fins; (ii) Adjust center of mass; (iii) Ini- tialize image processing for a gi ven configuration; (i v) Run the experiment. T able 4 summarizes the design parameters for economic assessment. The total cost can be ultimately e xpressed as: C total = N conv ergence ×  C 3 D + C 3 D Fins + C lead + C Wchannel  (50) where N conv ergence denotes the number of iterations re- quired for the conv ergence based on the number elements per batch. As can be seen from C total , some terms depend on N , while some terms do not. This indicates that there exist economies of scale in running Batch Bayesian Opti- mization. T able 5 shows the cost associated with dif ferent batch sizes according to our economic assessment. One can conclude from this table that as the number of elements per batch increases, the associated cost decreases. 11 5 10 15 Time [sec] -0.2 -0.1 0 0.1 0.2 Roll Tracking Error [rad] Optimal Sub-optimal Figure 8: (Zoomed) roll tracking error before and after the optimization 5 10 15 Time [sec] -0.2 -0.1 0 0.1 0.2 Heading Tracking Error [rad] Optimal Sub-optimal Figure 9: (Zoomed) heading tracking error before and after the optimization 0 20 40 60 80 100 Time [sec] 0 0.2 0.4 0.6 0.8 1 Zenith Angle [rad] Optimal Sub-optimal Figure 10: Zenith angle before and after the optimization T able 4: Design parameters for economic assessment Parameter Description V alue c eng cost required to hire employee(s) to conduct experiments $30 hours c recharge cost required for equipment recharge $240 day c W recharge cost required for water channel rechar ge $2400 day c lostTime opportunity cost of lost time $1200 day T print time required to print the main body 12 hours T 3DFins time required to print the fins 4 hours T lead lead time 5 days T setup time required to setup experiment 30 min T reconfig time required to reconfigure the model 5 min T exp time required to conduct one experiment in inner loop 1200 sec N number of elements per batch 1, 3, 4 m number of employee(s) required to print the model and fins 1 m 0 number of employee(s) required to conduct experiments 2 T able 5: Cost in $ associated with running experiments with different batch sizes # of elements per batch Cost ( $ ) 1 54293 3 49200 4 44533 6 Conclusion This work presented a nested co-design frame work where Batch Bayesian Optimization replaced local tech- niques at both the plant and controller optimization lev els. Our results using an A WE system demonstrate that both the plant and controller parameters conv erge to their respectiv e optimal values within only a few iterations. Furthermore, the results confirm that economies of scale e xist with Batch Bayesian Optimization approach. References [1] Altaeros Energies, inc., http://www .altaerosenergies.com. [2] Ampyx Power , http:www .ampyxpower .com. [3] Makani Power , inc., http://makanipower .com. [4] SkySails gmbH and co., http://www .skysails.info. [5] Chris V ermillion, Tre y Grunnagle, Ronny Lim, and Ilya K olmanovsky . Model-based plant design and hi- erarchical control of a prototype lighter-than-air wind energy system, with experimental flight test results. IEEE T ransactions on Contr ol Systems T echnology , 22(2):531–542, 2014. [6] Ali Baheri and Chris V ermillion. Altitude optimiza- tion of airborne wind energy systems: A Bayesian Optimization approach. In American Contr ol Confer- ence , Seattle, US, 2017. 12 Figure 11: W ater channel experimen- tal setup at UNC Charlotte Figure 12: 1 100 -scale B A T model [7] Shamir Bin-Karim, Alireza Baf andeh, Ali Baheri, and Christopher V ermillion. Spatiotemporal optimization through gaussian process-based model predictiv e con- trol: A case study in airborne wind energy . IEEE T ransactions on Contr ol Systems T echnology , 2017. [8] Parvin NikpoorParizi, Nihar Deodhar, and Christo- pher V ermillion. Combined plant and controller per - formance analysis and optimization for an energy- harvesting tethered wing. In American Contr ol Con- fer ence (ACC), 2016 , pages 4089–4094. American Automatic Control Council (AA CC), 2016. [9] Hosam K. Fath y , Panos Y . Papalambros, and A. Galip Ulsoy . Integrated plant, observer , and controller opti- mization with application to combined passive/acti ve automotiv e suspensions. In ASME 2003 International Mechanical Engineering Congress and Exposition , pages 225–232. American Society of Mechanical En- gineers, 2003. [10] James T Allison, Tinghao Guo, and Zhi Han. Co- design of an acti ve suspension using simultaneous dy- namic optimization. Journal of Mechanical Design , 136(8):081003, 2014. [11] Hosam K. Fathy , Scott A. Bortoff, G. Scott Copeland, Panos Y . Papalambros, and A. Galip Ulsoy . Nested optimization of an elev ator and its gain-scheduled lqg controller . In ASME 2002 International Me- chanical Engineering Congr ess and Exposition , pages 119–126. American Society of Mechanical Engineers, 2002. [12] Joe Deese, Nihar Deodhar , and Chris V ermillion. Nested plant / controller co-design using g-optimal design and extremum seeking: Theoretical frame- work and application to an airborne wind energy sys- tem. In International F ederation of Automatic Con- tr ol , T oulouse, FR, 2017. [13] Diane L. Peters, PY Papalambros, and A G Ulsoy . Control proxy functions for sequential design and con- trol optimization. Journal of Mechanical Design , 133(9):091007, 2011. [14] Kamal Y oucef-T oumi. Modeling, design, and con- trol integration: a necessary step in mechatronics. IEEE/ASME T ransactions on Mechatr onics , 1(1):29– 38, 1996. [15] Julie A. Reyer and Panos Y . Papalambros. Optimal design and control of an electric dc motor . In Pr oceed- ings of the 1999 ASME Design Engineering T echnical Confer ences , pages 85–96. Citeseer, 1999. [16] Hosam K. Fathy , Panos Y . Papalambros, A. Galip Ul- soy , and Davor Hro vat. Nested plant/controller opti- mization with application to combined passive/acti ve automotiv e suspensions. In American Contr ol Confer - ence, 2003. Pr oceedings of the 2003 , volume 4, pages 3375–3380. IEEE, 2003. [17] Timothy W ard Athan and Panos Y Papalambros. A note on weighted criteria methods for compromise so- lutions in multi-objectiv e optimization. Engineering Optimization , 27(2):155–176, 1996. [18] Indraneel Das and John E. Dennis. A closer look at drawbacks of minimizing weighted sums of objectiv es for pareto set generation in multicriteria optimization problems. Structur al optimization , 14(1):63–69, 1997. [19] James T . Allison and Sam Nazari. Combined plant and controller design using decomposition-based de- sign optimization and the minimum principle. In ASME 2010 International Design Engineering T ech- nical Conferences and Computers and Information in Engineering Confer ence , pages 765–774. American Society of Mechanical Engineers, 2010. [20] Hosam K. Fathy , Julie A. Reyer , Panos Y . Papalam- bros, and A G Ulsov . On the coupling between the 13 plant and controller optimization problems. In Amer- ican Contr ol Conference , 2001. Proceedings of the 2001 , volume 3, pages 1864–1869. IEEE, 2001. [21] Ali Baheri, Joe Deese, and Chris V ermillion. Com- bined plant and controller design using Bayesian op- timization: A case study in airborne wind energy sys- tems. In ASME 2017 Dynamic Systems and Contr ol Confer ence , T ysons Corner, US, 2017. [22] Joe Deese and Chris V ermillion. Nested plant/controller co-design using G-optimal de- sign and continuous time adaptation la ws: Theoretical framew ork and application to an airborne wind energy system. Submitted to ASME Journal on Dynamic Systems, Measur ement, and Contr ol . [23] Ali Baheri, Shamir Bin-Karim, Alireza Bafandeh, and Christopher V ermillion. Real-time control us- ing bayesian optimization: A case study in airborne wind ener gy systems. Contr ol Engineering Practice , 69:131–140, 2017. [24] Ali Baheri and Chris V ermillion. Context-dependent bayesian optimization in real-time optimal control: A case study in airborne wind energy systems. In NIPS W orkshop on Bayesian Optimization, Long Beach, CA , 2017. [25] Ali Baheri, Praveen Ramaprabhu, and Chris V ermil- lion. Iterativ e in-situ 3D layout optimization of a re- configurable ocean. In ASME 2017 Dynamic Systems and Contr ol Confer ence , T ysons Corner, US, 2017. [26] Ali Baheri, Pra veen Ramaprabhu, and Christopher V ermillion. Iterativ e 3d layout optimization and para- metric trade study for a reconfigurable ocean current turbine array using bayesian optimization. Renewable Ener gy , 127:1052–1063, 2018. [27] Hany Abdelrahman, Felix Berkenkamp, Jan Poland, and Andreas Krause. Bayesian optimization for max- imum power point tracking in photov oltaic power plants. In Pr oc. Eur opean Contr ol Confer ence , Aal- borg, DK, 2016. [28] Roman Garnett, Michael A. Osborne, and Stephen J. Roberts. Bayesian optimization for sensor set selec- tion. In Proceedings of the 9th ACM/IEEE Interna- tional Confer ence on Information Processing in Sen- sor Networks , Stockholm, SE, 2010. [29] Carl Edward Rasmussen and Christopher K. I. W illiams. Gaussian processes for machine learning . Number ISBN 0-262-18253-X. The MIT Press, 2006. [30] Eric Brochu, Vlad M. Cora, and Nando De Fre- itas. A tutorial on bayesian optimization of expensi ve cost functions, with application to activ e user mod- eling and hierarchical reinforcement learning. arXiv pr eprint arXiv:1012.2599 , 2010. [31] A. Zilinskas J. Mockus, V . T iesis. The application of Bayesian methods for seeking the extremum. T owar ds Global Optimization , 2:117–129, 1978. [32] Javier Gonzalez, Zhenwen Dai, Philipp Hennig, and Neil Lawrence. Batch bayesian optimization via lo- cal penalization. In Arthur Gretton and Christian C. Robert, editors, Pr oceedings of the 19th International Confer ence on Artificial Intelligence and Statistics , volume 51 of Pr oceedings of Machine Learning Re- sear ch , pages 648–657, Cadiz, Spain, 09–11 May 2016. PMLR. [33] Chris V ermillion, Ben Glass, and Sam Greenwood. Evaluation of a water channel-based platform for char- acterizing aerostat flight dynamics: A case study on a lighter-than-air wind energy system. In Lighter -Than- Air Systems Confer ence, AIAA (P art of A viation 2014) , 2014. [34] Mustafa Sario ˘ glu and T ahir Y avuz. V ortex shedding from circular and rectangular cylinders placed hori- zontally in a turbulent flo w . T urkish Journal of Engi- neering and En vir onmental Sciences , 24(4):217–228, 2000. [35] Nihar Deodhar , Alireza Bafandeh, Joe Deese, Brian Smith, T im Muyimbwa, Christopher V ermillion, and Peter Tkacik. Laboratory-scale flight characterization of a multitethered aerostat for wind energy generation. AIAA Journal , pages 1–11, 2017. 14