Data Center Chiller Plant Optimization via Mixed-Integer Nonlinear Differentiable Predictive Control

Data Cen ter Chiller Plan t Optimization via Mixed-In teger Nonlinear Differen tiable Predictiv e Con trol J´ an Boldo c k ´ y ∗ , ∗∗ Cary F aulkner ∗∗ Elad Mic hael ∗∗ Martin Gulan ∗ Aaron T uor ∗∗ J´ an Drgo ˇ na ∗∗∗ , ∗∗ ∗ Slovak University of T e chnolo gy, Br atislava, Sl ovakia (e-mail: { jan.b oldo cky, martin.gulan } @stub a.sk) ∗∗ Pacific Northwest National L ab or atory, R ichland, W A 99352 USA (e-mail: { c ary.faulkner, elad.michael, aar on.tuor } @pnnl.gov) ∗∗∗ Johns Hopkins University, MD 21218, USA, (e-mail: jdr gona1@jh.e du) Abstract: W e presen t a computationally tractable framew ork for real-time predictive control of multi-c hiller plan ts that inv olve both discrete and contin uous control decisions coupled through nonlinear dynamics, resulting in a mixed-in teger optimal con trol problem. T o address this challenge, w e extend Differentiable Predictiv e Con trol (DPC)—a self-sup ervised, mo del- based learning methodology for approximately solving parametric optimal con trol problems—to accommo date mixed-in teger control p olicies. W e b enchmark the prop osed framework against a state-of-the-art Model Predictive Con trol (MPC) solv er and a fast heuristic Rule-Based Con troller (RBC). Simulation results demonstrate that our approac h ac hieves significan t energy sa vings ov er the RBC while maintaining orders-of-magnitude faster computation times than MPC, offering a scalable and practical alternativ e to con ven tional combinatorial mixed-in teger con trol formulations. Keywor ds: control and managemen t of energy systems, mixed-integer con trol, differentiable programming, energy systems optimization 1. INTR ODUCTION The rapid expansion of artificial intelligence and digital infrastructure has led to a sharp increase in global data cen ter electricit y demand. In the United States alone, data centers curren tly consume ab out 150 TWh of elec- tricit y annually and are pro jected to grow to betw een 200 TWh and 500 TWh by 2030 (Aljb our et al., 2024). With this expansion comes an accompan ying demand for high-efficiency cooling systems in data centers. Chiller plan ts are commonly deploy ed to p erform this task. A c hiller plant t ypically includes heat pumps (often referred to as c hillers), pumps, and cooling tow ers, with the c hillers b eing the most energy-in tensiv e comp onents within the plan t. The con trols of chiller plan ts in v olve b oth contin u- ous decision v ariables (e.g., chilled water supply temp era- ture, differen tial pressure) and discrete decision v ariables (e.g., n umber of c hillers activ e, pump mode). Optimizing these coupled decisions can significantly reduce electricity consumption, ensure reliable op eration, and extend the ⋆ This research was supp orted b y the U.S. Department of Energy Building T echnologies Office (BTO), USA within the Office of Energy Efficiency and Renew able Energy (EERE), USA under Contract No. DE-AC05-76RL01830. J´ an Drgoˇ na was supported by the Ralph O’Connor Sustainable Energy Institute at Johns Hopkins Universit y . J.B. and M.G. also ac knowledge the contribution of the Slov ak Research and Developmen t Agency under the pro ject APVV-22- 0436, and the Scientific Grant Agency of the Slov ak Republic under the grant 1/0401/26. lifespan of costly mec hanical assets. Ho wev er, this requires the solution to a high-dimensional mixed-in teger nonlinear optimal control problem, which is computationally chal- lenging to solve in real time. The most widely used industrial strategy is Rule-Based Con trol (RBC), v alued for its in terpretability and low implemen tation cost. Several works hav e focused on re- fining RBC through heuristic tuning and performance optimization. F or example, Ma and W ang (2011) devel- op ed a sup ervisory RBC framework that stages chillers based on the part-load ratio (PLR) to improv e efficiency; Huang et al. (2016) prop osed adaptive PLR thresholds to balance energy use and load following; Liu et al. (2017) form ulated rule-based optimization of chiller sequencing using plant-lev el p erformance metrics; and F aulkner et al. (2025) recently performed simulations optimizing con trol parameters of RBC strategies to balance chiller, pump, and to wer loads. While RBC metho ds offer fast and in ter- pretable op eration, they rely on fixed heuristics and t ypi- cally fail to accoun t for future disturbances, often resulting in sub optimal energy efficiency or constraint violations. T o ov ercome these limitations, researchers hav e inv esti- gated predictiv e control strategies for chiller plants (Sala- Cardoso et al., 2020; T erzi et al., 2020; Chan et al., 2022; F an and Zhou, 2023; He et al., 2023), that leverage dy- namic mo dels and load forecasts to proactively schedule c hillers and optimize contin uous setp oin ts. F or example, Mo del Predictive Con trol (MPC) has been sho wn to im- pro ve chiller staging under v arying ambien t and load con- ditions (Sala-Cardoso et al., 2020; He et al., 2023), while data-driv en or learning-based predictive control v ariants adapt to plant-specific nonlinearities (T erzi et al., 2020; Chan et al., 2022; F an and Zhou, 2023). How ever, most existing metho ds handle either contin uous (Pan and Li, 2024) or discrete con trol actions in isolation, and few address the mixed-in teger nature of chiller plant op eration. A principled wa y to unify discrete and contin uous decisions under constraints is Mixed-In teger Mo del Predictive Con- trol (MI-MPC) (McAllister and Ra wlings, 2022; Dua et al., 2002; Ric hards and How, 2005; Kirches, 2011; T ak ap oui et al., 2020). MI-MPC provides a systematic framework for incorporating nonlinear dynamics and hard constraints while pro ducing high-quality solutions. How ev er, its com- putational complexit y grows com binatorially with the horizon length and the n umber of integer v ariables, mak- ing real-time operation in tractable, particularly for large- scale systems suc h as multi-c hiller plants. Approximation strategies such as mo ve blo c king, surrogate linearizations, or w arm-starts help alleviate but do not eliminate these scaling limitations, esp ecially when tight sampling perio ds and long prediction horizons are required. A representa- tiv e example in the con text of cen tral chiller plants is the decomposition-based MPC framework of Deng et al. (2015), which combines dynamic programming and mixed- in teger linear programming to obtain tractable sub optimal solutions; how ever, the approach still requires rep eated online optimization, limiting the computational scalability . These challenges hav e motiv ated a growing b o dy of re- searc h on learning-based alternatives that av oid online com binatorial search. In imitation-based or approximate MPC approac hes, machine learning (ML) policies are trained to mimic optimal or near-optimal mixed-integer solutions offline and then deplo y ed as explicit con trollers, offering fast real-time ev aluation while maintaining reason- able constraint satisfaction(Karg and Lucia, 2018; Cauligi et al., 2020; Domahidi et al., 2014). In parallel, warm- start tec hniques utilize ML on previously solv ed MIP in- stances to initialize subsequent optimizations with go od primal and dual guesses, reducing branch-and-bound ex- ploration and improving real-time feasibility in receding- horizon settings (Masti and Bemp orad, 2019; Marcucci and T edrake, 2021; Reiter et al., 2024). Beyond imitation learning, the broader field of learning to optimize (L2O) has also explored data-driv en acceleration of MIP solvers through learned branching, no de selection, cutting, and heuristic policies that reduce solv e times b y orders of mag- nitude (Bertsimas and Stellato, 2022; T ang et al., 2025; He et al., 2014; Zarpellon et al., 2021). These dev elopments motiv ate a self-sup ervised Differentiable Predictiv e Con- trol framework, which merges the interpretabilit y of MPC with the scalability of L2O. T o address these challenges, w e build on our prior work on Differen tiable Predictive Con trol (DPC) (Drgo ˇ na et al., 2022; Drgo ˇ na et al., 2024), which form ulates parametric MPC problems as differen tiable programs, enabling offline self-sup ervised learning of explicit neural control poli- cies through gradien t-based learning. Building up on this foundation, our recen t work on Mixed-In teger DPC (MI- DPC) (Boldo c k ´ y et al., 2025) extended this idea to handle mixed-in teger optimal control problems (MI-OCP) with discrete decision v ariables. In this w ork, w e demonstrate the application of MI-DPC to a c hallenging m ulti-chiller plan t optimization problem c haracterized b y nonlinear dy- namics, mixed-in teger decisions, sub ject to input and state constrain ts. The prop osed MI-DPC approac h provides a scalable self-sup ervised learning algorithm for synthesizing real-time capable explicit neural con trol p olicies that yield high-qualit y appro ximate solutions to the underlying MI- OCP , without online combinatorial search. The contributions of this work include: (i) a nonlinear con trol-oriented multiple chiller system model; (ii) exten- sion of the MI-DPC framework to mixed-in teger nonlinear optimal control problems arising in chiller plant op era- tions; (iii) extension of the MI-DPC framework to han- dle binary decision v ariables and the introduction of the binary-v ariance regularization tec hnique, necessary to pre- v ent high-frequency switc hing b ehavior; (iv) a systematic p erformance comparison of the prop osed MI-DPC against an RBC and an MPC solv ed with a state-of-the-art mixed- in teger solver; and (v) the release of op en-source co de 1 enabling replication of the numerical exp erimen ts. 2. PR OBLEM FORMULA TION W e consider a simplified m ulti-chiller plan t, illustrated in Figure 1, with c hillers op erating in parallel. In the follo wing, we use a sup erscript to denote an asso ciation of a v ariable with chiller lo op i , i.e. ˙ m (1) represen ts the mass flo w through the first chiller. W e further denote quantities sharing the same units b y a common symbol with a subscript, such as T e , T s , T r for the ev aporator, supply and return temp eratures, resp ectively . In the follo wing, we use M to denote the n umber of chillers in the system. Pump 1 Pump M m ( M ) m (1) Chiller M ( T e ( M ) ) Chiller 1 ( T e (1) ) T s ( M ) T s (1) T r Thermal Load (Q load ) ( i ) ( i ) Fig. 1. Sc hematic of a multiple-c hiller plant arranged in a parallel configuration. 2.1 System Mo del The con tin uous-time thermal dynamics of the m ulti-c hiller system are 1 Code is av ailable at github.com/pnnl/MI-DPC. Q ( i ) ( t ) = η r ˙ m ( i ) ( t ) c p δ ( i ) ( t )  T r ( t ) − T ( i ) s ( t )  , (1a) ˙ T r ( t ) = 1 C r ˜ Q load ( t ) − M X i =1 Q ( i ) ( t ) ! , (1b) ˙ T ( i ) s ( t ) = − c p η s C ( i )  δ ( i ) ( t ) ˙ m ( i ) ( t )  T ( i ) s ( t ) − T ( i ) e ( t )  , (1c) where Q ( i ) ( t ) ∈ R ≥ 0 is the delivered co oling from chiller i at time t , δ ( i ) ( t ) ∈ { 0 , 1 } denotes chiller’s on/off state, c p ∈ R > 0 is the sp ecific heat capacity of the co olant, and C r , C ( i ) ∈ R > 0 are the thermal capacitances of the return and chiller lo ops. Moreo ver, η r , η s ∈ (0 , 1) denote the heat-exc hanger efficiency co efficien ts, and ˜ Q load ( t ) ∈ R ≥ 0 is the effective thermal load, estimated using a low-pass, discrete-time finite impulse resp onse filter ˜ Q load ( t ) = L X l =0 h l Q load ( t − l ∆ t ) , (2) where Q load ∈ R ≥ 0 is the d irect server load—an exogenous v ariable in this context, ∆ t denotes the length of a discrete- time step, L ∈ N ≥ 1 is the filter order, and ( h l ) L l =0 ∈ R > 0 are the filter co efficien ts, with P l h l = 1. Note that Q load is treated as a deterministic known v ariable in our study . W e consider the FIR filter to capture higher-order heat-transfer dynamics, thereby enhancing the mo del’s expressiv e capacity . Naturally , the filter co efficien ts can b e a sub ject of system identification along with other system parameters. The p ow er consumption of chillers P ( i ) chiller ( t ) ∈ R ≥ 0 and pumps P ( i ) pump ( t ) ∈ R ≥ 0 are mo deled as P ( i ) chiller ( t ) = Q ( i ) ( t ) COP ( i ) ( t ) + ρ ( i ) δ ( i ) ( t ) , (3a) COP ( i ) ( t ) = a ( i ) 0 + a ( i ) 1 Q ( i ) ( t ) Q ( i ) max + a ( i ) 2 Q ( i ) ( t ) Q ( i ) max ! 2 , (3b) P ( i ) pump ( t ) = γ ( i )  ˙ m ( i ) ( t ) δ ( i ) ( t )  3 , (3c) for the system constants a ( i ) 0:2 , γ ( i ) , ρ ( i ) , Q ( i ) max ∈ R . The constan t ρ represen ts the base pow er required for chiller op eration. In this formulation, the chiller’s co efficient of p erformance COP ( i ) ( t ) ∈ R ≥ 0 is a function of the part load ratio (PLR), the ratio of curren t co oling to maxim um co oling capacity . Ensuring that each chiller is op erating near its highest efficiency load ratio is the core of the underlying optimal control problem. In the follo wing, we discretize the contin uous time dynamics given b y (1b) and (1c) using a fourth-order Runge-Kutta method with a time-step length ∆ t and discrete time index k . 2.2 Optimal Contr ol Pr oblem Let u ( i ) k := { δ ( i ) ( t ) , T ( i ) e ( t ) , ˙ m ( i ) ( t ) } denote the set of con- trol v ariables and x ( i ) k := { T ( i ) s ( t ) , T r ( t ) } denote the set of state v ariables at discrete time step k . Then the optimal con trol problem (OCP) can b e formulated as follows: min { u ( i ) k } M i =1 , N − 1 k =0 N − 1 X k =0 M X i =1 P ( i ) chiller ,k + P ( i ) pump ,k + ∥ ∆ δ ( i ) k ∥ 2 R + λ W ( Q k − Q load ,k ) 2 (4a) s.t. x ( i ) k +1 = f  x ( i ) k , u ( i ) k , ˜ Q load ,k , ∆ t  , (4b) Q k = X i Q ( i ) k , ∀ k ∈ { 0 , . . . , N − 1 } , (4c) X i δ ( i ) k ≥ 1 , ∀ k ∈ { 0 , . . . , N − 1 } , (4d) u ( i ) min ≤ u ( i ) k ≤ u ( i ) max , ∀ k ∈ { 0 , . . . , N − 1 } , (4e) x ( i ) min ≤ x ( i ) k ≤ x ( i ) max , ∀ k ∈ { 0 , . . . , N } , (4f ) Q ( i ) k ≤ Q ( i ) max , ∀ k ∈ { 0 , . . . , N − 1 } , (4g) (3a) − (3c) , ∀ k ∈ { 0 , . . . , N − 1 } , (4h) δ ( i ) k ∈ { 0 , 1 } , ∀ k ∈ { 0 , . . . , N − 1 } , (4i) x ( i ) (0) = x ( i ) ( t ) . (4j) Here, the ob jective function (4a) p enalizes the chiller and pump p ow er consumption, switching of chiller on/off sta- tus ∆ δ ( i ) k := δ ( i ) k +1 − δ ( i ) k ∀ k ∈ { 0 , . . . , N − 2 } —weigh ted b y a di- agonal matrix R of size M , and the squared error b etw een total delivered co oling Q ∈ R ≥ 0 and the thermal load de- mand, w eighted b y coefficient λ W ∈ R > 0 . W eighting matri- ces R and λ W directly gov ern the trade-off b et ween integer switc hing frequency and load tracking error minimization. Larger R suppresses frequent switching b ehavuour, while larger λ W v alues prioritizes minimizing the load tracking error. The dynamics (4b) is solved using RK4 integration of the thermal equations (1b)–(1c). Constraint (4d) en- sures that at least one chiller remains activ e at all times, serving as a pro cess safet y constraint. Finally , constraints (4e)–(4h) enforce op erational and physical systems limits. Note that the OCP (4) is non-trivial, primarily due to system dynamics, where the binary decision v ariable is coupled with tw o bilinear terms in (1b) and (1c). Addi- tionally , the OCP includes a set of p olynomial constraints in (4h), making it a mixed-integer nonlinear program. 2.3 Rule-Base d Contr ol Due to the complexit y of the given OCP , a rule-based con- trol (RBC) is the most common control strategy used for impro ving energy efficiency of the chiller plan ts. Assuming the PLR–COP c haracteristic is known, a typical approach for multi-c hiller systems is to stage the c hillers based on system’s PLR, defined as PLR k = P M i =1 Q ( i ) k [ Q (1) max , . . . , Q ( M ) max ][ δ (1) k , . . . , δ ( M ) k ] ⊤ . (5) In cases where the chiller plant consists of units with similar PLR–COP characteristics, a simple sequen tial threshold-based staging strategy can b e employ ed: s k +1 =    min ( M , s k +1 ) , if PLR k > t, max (1 , s k − 1 ) , if PLR k < t, s k , otherwise, (6) where s k := P M i =1 δ ( i ) k ∈ Z > 0 denotes the n umber of activ e chillers at time k , and t, t ∈ (0 , 1) are the stag- ing thresholds—selected to a void high-frequency switc hing while maintaining the effective op erating range of the system. In addition, the con tinuous control v ariables are t ypically fixed at constant v alues that preven t pro cess constrain t violations. 3. METHODOLOGY In this section, we describ e the metho dology designed to reform ulate the defined OCP (4) as a differen tiable program suitable for the DPC framew ork. W e preserv e the differen tiable con trol ob jectiv e (4a) and enforce con- strain ts (4b)–(4h) using either pro jections, where the con- strain ts are enforced within a differen tiable computational graph by pro jecting the v ariables onto their feasible sets, or relaxations using p enalty terms, where the constraints are relaxed and incorp orated within a loss function. The tec hnique used to ensure the integralit y of the discrete- v alued decision v ariable is discussed in Section 3.2. 3.1 Pr oblem R eformulation T o cast (4) as a differentiable program, we relax the constrain ts (4e)–(4f) using the p enalt y metho d. Moreov er, the constraint (4d) is enforced by fixing δ (2) ( k ) = 1, whereas (4b), (4g), and (4h) can b e incorporated within the differentiable computational graph and thus do not require any relaxation. Considering the ab ov e, the OCP can b e rec ast as min θ 1 ,θ 2 ,θ 3 E ξ k ∼ P ξ  N − 1 X k =0 M X i =1 P ( i ) chiller ,k + P ( i ) pump ,k + ∥ ∆ δ ( i ) k ∥ 2 R + λ W ( Q k − Q load ,k ) 2 + λ x q ( x k ) + λ u p ( u k ) + ∥ ˜ δ ( i ) (1 − ˜ δ ( i ) ) ∥ 2 Λ  (7a) s.t. x ( i ) k +1 = f  x ( i ) k , u ( i ) k , ˜ Q load ,k , ∆ t  , (7b) ˙ m ( i ) k = π θ 1 ( ξ k ) , (7c) T ( i ) e , k = π θ 2 ( ξ k ) , (7d) ˜ δ ( i ) k = π θ 3 ( ξ k ) , (7e) δ ( i ) k ∈ { 0 , 1 } , ∀ k ∈ { 0 , . . . , N − 1 } , (7f ) δ (2) k =1 , ∀ k ∈ { 0 , . . . , N − 1 } , (7g) ξ 0 := [ x ( i ) ( t ) , ˜ Q load ( t ) , Q load ( t ) , . . . , Q load ( t + N ∆ t )] ⊤ . (7h) Here, we optimize the parameters θ ( · ) of the control p oli- cies π θ 1 , π θ 2 , and π θ 3 that map the v ector of control param- eters ξ , defined in (7h), to contin uous-v alued control v ari- ables ˙ m , T e , and a slack v ariable ˜ δ ∈ R M − 1 —represen ting the relaxed v alue of the binary control v ariable δ . Notice that by fixing one of the binary v ariables, the dimension of ˜ δ is reduced to M − 1. T o compute control v ariables, we use separate neural mo dules that are represented by fully connected deep neural netw orks. Moreov er, we sample the con trol paramete rs from a known probability distribution P ξ to emulate scenarios that are lik ely to o ccur during the system’s op eration, which defines a self-sup ervised nature of the DPC metho dology . F unctions q ( · ) and p ( · ) denote state and input p enalt y functions w eighted by λ x and λ u , resp ectiv ely . Equation (7g) states that one chiller in the system m ust alwa ys remain in op eration. 3.2 Differ entiable R ounding T o facilitate gradient-based optimization, each comp o- nen t of the problem (7) m ust b e differentiable with well- defined gradients. This requirement is satisfied for all terms in (7) except for constraint (7f), whic h enforces the integralit y of the binary v ariable δ , thereby introduc- ing non-differentiabilit y . W e enforce the binary integralit y through a threshold-based discretization scheme with a fixed threshold of 0.5, given as δ ( i ) k = ( 1 , if ˜ δ ( i ) k > 0 . 5 , 0 , otherwise. (8) This function b ehav es analogously to a Heaviside (unit step) function exhibiting a discontin uity at 0 . 5. T o en- able gradient-based optimization, we leverage a Straight- through estimation heuristic in tro duced b y Bengio et al. (2013). W e approximate the gradient of (8) using its dif- feren tiable surrogate, given by a scaled Sigmoid function σ ( x ) = 1  1 + e − µ ( x − 0 . 5)  , (9) where µ ∈ R > 0 is a slop e parameter determining the sharpness of the transition, and the function is centered at 0 . 5. Accordingly , the deriv ative is approximated as ∂ δ ( i ) k ∂ ˜ δ ( i ) k ≈ ∂ σ  ˜ δ ( i ) k  ∂ ˜ δ ( i ) k . (10) 3.3 Binary-V arianc e R e gularization T o promote binary polarity of the relaxed v ariable ˜ δ and to complemen t the ∆ δ p enalty in suppressing high-frequency switc hing, we in tro duce a binary-v ariance regularization (BVR) defined as ∥ ˜ δ ( i ) (1 − ˜ δ ( i ) ) ∥ 2 Λ , where Λ ∈ R M − 1 × M − 1 denotes the weigh ting matrix. As the weigh t increases, the regularization p enalizes v alues of ˜ δ near the discretization threshold (0 . 5), as well as those outside the [0 , 1] range. Based on extensive numerical exp eriments, this regular- ization has prov en crucial to the prop osed framew ork, as its omission leads to ˜ δ v alues clustering around the discretization threshold, thereby increasing the likelihoo d of high-frequency switc hing behavior. The influence of this regularization on the relaxed binary v ariable for different p enalt y magnitudes is illustrated in Figure 3. 3.4 Mixe d-Inte ger Policy Optimization The computational graph, illustrated in Figure 2, com- prises four comp onents: a contin uous neural p olicy , a dif- feren tiable rounding lay er, a differentiable ordinary differ- en tial equation (ODE) solver (numerical integrator), and a mo del predictive loss function. Given syn thetically sam- pled con trol parameters that include the initial conditions and an N -long vector of exogenous v ariables (serv er load), a closed-lo op rollout of length N is p erformed through the integrated system dynamics, constituting the single sho oting forward pass. Consequently , the loss function is ev aluated based on the resulting closed-lo op tra jectories. control parameters i n i t i a l c o n d i t i o n s r e f e r e n c e s c o n s t r a i n t s continuous neural policy e n e r g y m i n i m i z a t i o n s w i t c h i n g p e n a l t y chiller plant predictive control loss function r e f e r e n c e t r a c k i n g c o n s t r a i n t s p e n a l t i e s dif ferentiable rounding ... ... chiller dynamics ODE solver f Fig. 2. Conceptual diagram of the mixed-in teger differentiable predictiv e control for nonlinear c hiller plan t optimization. Green dashed arrows represent the forward pass, while red dashed arrows represent the backw ard pass. 0 5 10 15 20 Time step — k 0 . 0 0 . 5 1 . 0 ˜ δ ( k ) Λ = 0 Λ = 1 Λ = 10 Λ = 100 Fig. 3. Effect of binary-v ariance regularization on relaxed in teger v alues, showing switching b eha vior under a step change in co oling load from 150kW to 500kW o ccurring at time k = 20 (setup from Sec. 4, N = 20). The dotted black line indicates the rounding threshold. The backw ard pass is used to enable the p olicy opti- mization by propagating gradients through the compu- tational graph. Sp ecifically , gradients of the loss func- tion with resp ect to p olicy parameters are obtained by bac k-propagating through the integrated system dynam- ics using either the backpropagation through time algo- rithm (Puskorius and F eldk amp, 1994) or the adjoint sen- sitivit y metho d (Griewank and W alther, 2008). The dif- feren tiability of the discretization lay er and the numerical in tegrator ensures that the gradient signals are transmit- ted contin uously through closed-lo op tra jectories. These gradien ts are then used to adjust the p olicy parameters in the direction of decreasing loss. Consequently , the policy parameters are iteratively up dated. The high-dimensional and non-conv ex nature of the underlying parameter space necessitates the use of sto c hastic gradient descent (SGD), for whic h global optimality cannot b e guaranteed. 4. NUMERICAL EXPERIMENTS In this section, we present and discuss the results of a nu- merical study in which w e test the three con trol strategies. The control strategies are tested for systems comprising differen t n um b ers of c hillers and v arious prediction horizon lengths. Eac h system consists of iden tical c hiller units with parameters listed in T ab. 1. 4.1 Exp erimental Setup Throughout the study , the system dynamics w ere sim- ulated using a constant time step length of ∆ t = 180 s and a finite impulse resp onse filter with co efficien ts h l = [0 . 45 , 0 . 2 , 0 . 15 , 0 . 1 , 0 . 05 , 0 . 05]. The load signal for eac h exp erimen t w as algorithmically generated to emulate daily v ariations in server load. Sp ecifically , the nighttime load amplitude was sampled from a uniform distribution U (100kW , 350kW ), while the daytime amplitude was sam- pled from U (300kW , 0 . 75 P i Q ( i ) max ) with additive white noise and a transition p eriod of 4 h. The system states w ere uniformly initialized within their respective op erating ranges T r ∼ U  T min r , T max r  and T s ∼ U  T min s , T max s  . Each neural control p olicy mo dule employ ed an identical archi- tecture, consisting of three hidden lay ers of size 200 with ReLU activ ations and [0 , 1] normalization on the input la yer. The control loss w eights were manually finetuned to Λ = 200, λ x , λ u = 10, λ W = 0 . 001, R = 20. It is imp ortant to emphasize that the magnitude of Λ critically influences the b ehavior of the prop osed framework. Sp ecifically , lo w er v alues of Λ led to instances of high-frequency switch- ing, whereas excessively large v alues resulted in training instabilit y . F urthermore, the hyperparameters were also man ually adjusted a learning rate of 0 . 006, a batc h size of 10 000, with training and developmen t dataset sizes of 30 000 and 10 000, resp ectively . The learning rate w as selected to ensure stable optimization, while the relatively large batch size w as c hosen to provide reliable gradient es- timates. T o further mitigate the risk of gradient explosion during optimization, the gradient norms w ere clipp ed with a maximum allo wed norm of 100, suc h that all gradients { g i } w ere uniformly rescaled to satisfy ∥ g ∥ 2 ≤ 100. The sigmoid slope co efficient was k ept at the default v alue µ = 1, as with scaling v alues, an excessively conserv ative binary b ehavior was exhibited, thereby limiting actuation flexibilit y and reducing p otential energy sa vings. In infer- ence the contin uous con trol inputs of MI-DPC policies T able 1. Chiller system parameters. Symbol Description V alue Unit c p Specific heat of water 4 . 184 kJ / kg ◦ C C Thermal capacitance 14 644 kJ / ◦ C C r Thermal capacitance 29 288 kJ / ◦ C ρ Base chiller p ow er 10 kW Q max Maximum rated cooling 500 kW a 0 Efficiency co efficien t 1 − a 1 Efficiency co efficien t 19 . 33 − a 2 Efficiency co efficien t − 18 . 33 − η s Heat exchanger efficiency 0 . 7 − η r Heat exchanger efficiency 0 . 75 − γ Pump p o wer coefficient 9.62e-4 kWs 3 / kg 3 ˙ m min Minimum mass flow rate 5 kg / s ˙ m max Maximum mass flow rate 20 kg / s T min s , T min e , T min r Minimal temp eratures 8 ◦ C T max s , T max e Maximal temp eratures 12 ◦ C T max r Maximal return temp. 40 ◦ C w ere pro jected onto their feasible sets via element wise clipping. The n umerical experiments were conducted using the Neu- roMANCER v1.5.6 library (Drgona et al., 2023) run- ning with PyT orch v2.8.0 and CUDA 12.8 with MI-DPC. The MI-MPC formulation was implemented in Pyomo v6.9.4 (Bynum et al., 2021) and solved using Gurobi v12.0.3 (Gurobi Optimization, LLC, 2024). Similarly to (P an and Li, 2024), the nonconv exit y arising from bilinear coupling betw een contin uous and binary decision v ariables within the heat-transfer equations, the problem was refor- m ulated using McCormic k env elope relaxations—resulting in a mixed-integer quadratic programming appro ximation solv able b y Gurobi. In addition, higher-order nonlineari- ties, such as cubic term in (3c), were similarly reform ulated using slack v ariables and McCormic k relaxation. Never- theless, the MI-MPC form ulation remains computationally in tractable for real-time applications, as discussed in the follo wing Section. The closed-lo op simulation of MI-DPC, implemented in a receding-horizon manner, is shown in Figure 6. The observ ed increase in chiller p ow er consumption b et ween 80h and 90h is attributed to the PLR-COP characteristic. During this interv al, the chiller’s op erating p oint shifts to ward a less efficien t PLR range, resulting in a lo w er COP and, consequently , increased p ow er consumption despite comparable thermal load v alues. In addition, the con troller sho ws a tendency to activ ate an additional chiller within this time in terv al, as reflected by the rising v alue of the relaxed integer v ariable ˜ δ . How ever, ˜ δ remains b elow the rounding threshold, as activ ating an additional unit w ould lik ely increase the load-trac king error, thereby outw eighing the p otential improv ement of energy efficiency within the ob jective function. 4.2 Computational Sc alability T able 2 summarizes the simulation results obtained with the RBC and MI-DPC policies obtained from a sim ulation length of 7 days. Note that we do not include MI-MPC results due to computational in tractability; instead, w e rep ort the mean inference time, with the maximum solver time limited to 180 s—a practically tractable upper b ound for real-time implementation, equiv alen t to the control sampling p erio d. This highlights the notorious scalability limitations of hard mixed-in teger problems. In contrast, the MI-DPC approac h alleviates computational scalabilit y issues in both offline and online stages. As shown in Fig. 4, the mean inference time remains nearly constant ev en as the prediction horizon increases, since the same net work architecture is used across all scenarios. The n umber of trainable parameters increases only with the input dimension, which grows with the prediction horizon to accommo date the longer preview of the load v ariable. When it comes to the training time (TT), scalability with resp ect to the prediction horizon length exhibits an appro ximately linear trend. As evidenced b y the simulation results, the principal adv an tage of MI-DPC lies in its computational efficiency during b oth the controller synthesis phase and inference. The metho d ac hieves inference times several orders of magnitude faster than the implicit MI-MPC formulation, 100 150 200 TT [s] M = 2 M = 3 M = 4 M = 5 5 10 15 20 30 40 50 60 70 80 90 100 Prediction horizon length — N 1 . 8 1 . 9 2 . 0 MIT [s] × 10 − 4 Fig. 4. Computational scalabilit y of MI-DPC across differ- en t num b er of chillers ( M ) and horizon lengths ( N ). The top panel illustrates the total training time (TT), while the b ottom panel rep orts the mean inference time (MIT). thereb y p otentially enabling solving the mixed-integer nonlinear optimal con trol problems on computationally constrained edge devices in real time—a capability not attainable with conv en tional approaches. 4.3 Performanc e Comp arison In comparison to the RBC, the prop osed MI-DPC frame- w ork ac hieved up to 11 . 17 % energy sa vings. These sa vings resulted primarily from more adaptiv e and frequent c hiller staging—alb eit at the exp ense of a bigger relative load trac king error, as well as from reduced water pump p o wer consumption. With MI-DPC we are able to sim ultane- T able 2. Summary of sev en-day simulation re- sults for the RBC, MI-DPC, and MI-MPC im- plemen tations for different num b ers of chillers ( M ) and prediction horizon lengths ( N ). Method Metric M = 2 M = 3 N = 5 N = 10 N = 15 N = 5 N = 10 N = 15 RBC EC [MWh] 13.52 - - 17.53 - - EC Chillers [MWh] 13.23 - - 17.14 - - EC Pumps [MWh] 0.28 - - 0.39 - - COP [-] 4.08 - - 3.88 - - Num. of switches [-] 9 - - 22 - - Mean RCE [%] 5.38 - - 5.90 - - MIDPC EC [MWh] 12.39 12.16 12.24 15.96 15.64 15.57 Savings [%] 8.38 10.01 9.45 8.95 10.79 11.17 EC Chillers [MWh] 12.34 12.13 12.20 15.87 15.55 15.52 EC Pumps [MWh] 0.04 0.04 0.04 0.09 0.08 0.05 COP [-] 4.44 4.54 4.51 4.26 4.37 4.38 Num. of switches [-] 15 14 12 30 30 28 Mean RCE [%] 7.21 8.49 8.83 7.39 7.72 9.31 MIT [s] 1.9e-04 1.9e-04 1.9e-04 1.9e-04 1.9e-04 1.9e-04 TT [s] 117.36 115.53 121.47 108.46 113.10 117.40 NTP [-] 248205 251205 254205 249408 252408 255408 MIMPC Inference Time 157.44 181.33 181.70 182.77 181.74 180.98 Used acronyms: EC—Energy Consumption, EC-COP—Effective Chiller Coefficient of P erformance, MIT—Mean Inference Time, NTP—Number of T rainable P arameters, TT—T raining Time, RCE— Relative Control Error. 1 Measured using a single thread of Intel i9-14900KF. 2 T rained using an NVIDIA GeF orce R TX 5090. 3 Solved with the gurobi solver. ously sc hedule chiller activ ation and adjust the con tinuous con trol v ariables in an adaptive manner. The RBC was implemen ted with staging thresholds of t = 0 . 15, t = 0 . 6 and constant mass flows and ev ap oration temp eratures of ˙ m = 10 kg / s and T e = 10 ◦ C, resp ectiv ely . In T able 2, the Effectiv e Chiller COP (EC-COP) is com- puted as EC-COP = P N s k =0 P M i =1 Q ( i ) k P N s k =0 P M i =1 P ( i ) chiller,k , where N s is the length of sim ulation. The Relativ e Con trol Error (R CE) is computed as R CE k = | Q load, k − P M i =1 Q i k | Q load ,k · 100% . The significance of the predictiv e capability of MI-DPC b ecame eviden t through extensiv e numerical experiments. In addition to energy savings, the MI-DPC framework out- p erformed the RBC in terms of pro cess constrain ts. Sp ecif- ically , when considering cooling ramp-rate constraints, i.e., ˙ Q min ≤ ˙ Q ≤ ˙ Q max , whic h hav e been observed in real c hiller plan t op erations and reported in (Shan et al., 2021); the RBC may fail to activ ate chillers in time when en- coun tering a steep increase in co oling demand. As illus- trated in Figure 5, the RBC stages chillers solely based on the staging threshold. How ev er, due to the co oling ramp limitation, the return temp erature can exceed the pro cess safety b ounds. In contrast, the MI-DPC frame- w ork anticipates future load changes and stages chillers proactiv ely based on the known future co oling demand. In the exp eriment, during the time interv al from 3h to 7h, the results demonstrate that MI-DPC effectively adapts to abrupt c hanges in the co oling load profile, even when suc h changes w ere not represented in the training dataset. 0 1000 2000 Q, Q load [kW] Q load RBC MI-DPC 0 50 T r [ ° C] T min r , T max r bounds 0 5 10 15 20 25 Time [h] 1 2 3 s [-] Fig. 5. Closed-lo op simulation results of a chiller plant ( M = 3, Q max = 1 MW) with MI-DPC and RBC p olicies, highligh ting the imp ortance of predictive action for stable op eration, when considering co oling ramp-rate constrain ts. Moreo ver, extensive simulation exp eriments indicate that extending the prediction horizon b eyond N = 15 do es not provide additional energy savings for the considered system. This b ehavior can b e attributed primarily to the relativ ely fast thermal dynamics of the modeled c hiller plan t, which allows for effective con trol within relativ ely shorter horizon lengths. Extending the prediction horizon length may b e adv antageous for systems with slow er dy- namics, such as those including thermal storage tanks or larger thermal masses. F rom a computational standp oint, the MI-DPC framework maintains scalability even as the prediction horizon length increases, as illustrated in Fig- ure 4. 4.4 Pr actic al Chal lenges and Limitations In this study , we assume a fully white-b ox, deterministic system mo del for p olicy synthesis and testing, pro viding a con trolled setting to ev aluate and b enc hmark the MI-DPC framew ork. In practical applications, mo del parameters m ust b e identified from real plant data and mo del fidelity v alidated to ensure robustness and generalizability . This can b e addressed through grey-b ox parameter estimation, where the system parameters are estimated using, for instance, a nonlinear least-squares, or physics-informed mo deling strategies (Rack auc k as et al., 2021; Ko c h et al., 2025). How ever, practical MPC deploymen ts are further c hallenged by mo del–plant mismatch, unmo deled dynam- ics, sensor and actuator noise, and actuator delays, which ma y degrade closed-lo op p erformance and compromise constrain t satisfaction if not explicitly accounted for. Cur- ren t MPC frameworks address these issues through dis- turbance modeling, state estimation, uncertain t y propaga- tion, and constraint tigh tening strategies (Qin and Badg- w ell, 2003; Bemp orad and Morari, 2007). It is further recognized that the problem is solved approx- imately , without a guaranteed global optimum. Moreov er, while in tegrality and input constraints satisfaction are en- forced through pro jection-based metho ds, strict feasibility for state constraints is not guaranteed. These asp ects will b e addressed in future w ork. 5. CONCLUSIONS W e presen ted a mixed-integer nonlinear Differentiable Pre- dictiv e Control (MI-DPC) framework for real-time op- timization of multi-c hiller plants, unifying discrete stag- ing and contin uous setp oin t control within a single dif- feren tiable program. The approach leverages (i) differen- tiable rounding with straight-through estimator gradients to handle binary decisions, (ii) binary-v ariance regulariza- tion to suppress chattering, (iii) a nonlinear chiller plant dynamics solver with a differen tiable ODE solver, and (iv) a mo del-predictiv e loss ev aluated ov er rollout tra jectories, yielding an explicit neural p olicy with constant-time in- ference. In simulations on multi-c hiller systems, MI-DPC ac hieved up to 11% energy savings relativ e to an indus- trial RBC baseline while preserving pro cess constraints, and deliv ered orders-of-magnitude faster inference than MI-MPC, indicating strong potential for a practical and scalable deplo yment on large-scale HV AC op erations. 0 500 Q [kW] Q load Q 0 50 100 P ( i ) chiller [kW] i = 1 i = 2 20 40 T r , T ( i ) s [ ° C] T r i = 1 i = 2 0 . 0 0 . 5 1 . 0 ˜ δ , δ ( i ) [-] i = 1 i = 2 ˜ δ 0 25 50 75 100 125 150 175 Time [h] 0 10 20 δ ( i ) ˙ m ( i ) [kg/s] i = 1 i = 2 0 25 50 75 100 125 150 175 Time [h] 0 2 4 COP [-] Fig. 6. Closed-lo op seven-da y simulation results of a tw o-chiller system obtained with MI-DPC for a prediction horizon of N = 15. Building on the findings of this study , several extensions are envisioned. The c hiller plant mo del is to be extended to include co oling tow ers, allowing for further energy savings p oten tial through optimization of condenser-water lo op setp oin ts (Huang et al., 2017; W ang et al., 2019). The in tegration of water-side economizers can also be explored, as these hav e b een shown to reduce energy through free- co oling (F aulkner et al., 2025) and introduce additional discrete con trol v ariables for economizer activ ation. F uture w ork will also incorp orate probabilistic co oling-load fore- casting to more accurately capture real-w orld uncertaint y . Finally , the framew ork may b e adopted for experimental v alidation in real-world chiller plants. DECLARA TION OF GENERA TIVE AI AND AI-ASSISTED TECHNOLOGIES IN THE WRITING PR OCESS No generative AI and AI-assisted technologies were used in the writing pro cess. REFERENCES Aljb our, J., Wilson, T., and P atel, P . (2024). Po w ering in telligence: Analyzing artificial intelligence and data cen ter energy consumption. EPRI White Pap er no. 3002028905 . Bemp orad, A. and Morari, M. (2007). Robust mo del pre- dictiv e control: A surv ey . In R obustness in identific ation and c ontr ol , 207–226. Springer. Bengio, Y., L´ eonard, N., and Courville, A. (2013). Es- timating or propagating gradien ts through sto chastic neurons for conditional computation. arXiv pr eprint arXiv:1308.3432 . Bertsimas, D. and Stellato, B. (2022). Online mixed- in teger optimization in milliseconds. INFORMS Journal on Computing , 34(4), 2229–2248. Boldo c k ´ y, J., Jav an, S.D., Gulan, M., M¨ onnigmann, M., and Drgo ˇ na, J. (2025). Learning to solve parametric mixed-in teger optimal control problems via differentiable predictive control. arXiv pr eprint arXiv:2506.19646 . Byn um, M.L., Hack ebeil, G.A., Hart, W.E., Laird, C.D., Nic holson, B.L., Siirola, J.D., W atson, J.P ., and W o o druff, D.L. (2021). Pyomo–optimization mo deling in python , v olume 67. Springer Science & Business Media, third edition. Cauligi, A., Culb ertson, P ., Stellato, B., Bertsimas, D., Sc hw ager, M., and Pa v one, M. (2020). Learning mixed- in teger conv ex optimization strategies for rob ot plan- ning and control. In 59th IEEE Confer enc e on De cision and Contr ol , 1698–1705. Chan, K., W ong, V.T., Y ow, A.K., Y uen, P ., and Chao, C.Y. (2022). Dev elopment and p erformance ev aluation of a chiller plant predictive op erational control strategy b y artificial intelligence. Ener gy and Buildings , 262. Deng, K., Sun, Y., Li, S., Lu, Y., Brouw er, J., Mehta, P .G., Zhou, M., and Chakrab orty , A. (2015). Mo del predictiv e control of cen tral c hiller plant with thermal energy storage via dynamic programming and mixed- in teger linear programming. IEEE T r ansactions on A utomation Scienc e and Engine ering , 12(2), 565–579. doi:10.1109/T ASE.2014.2352280. Domahidi, A., Ullmann, F., Morari, M., and Jones, C.N. (2014). Learning decision rules for energy efficient building con trol. Journal of Pr o c ess Contr ol , 24(6), 763– 772. Energy Efficient Buildings Sp ecial Issue. Drgo ˇ na, J., Ki ˇ s, K., T uor, A., V rabie, D., and Klauˇ co, M. (2022). Differentiable predictive control: Deep learning alternativ e to explicit mo del predictiv e con trol for un- kno wn nonlinear systems. Journal of Pr o c ess Contr ol , 116, 80–92. Drgona, J., T uor, A., Ko c h, J., Shapiro, M., Jacob, B., and V rabie, D. (2023). NeuroMANCER: Neural Mo dules with Adaptive Nonlinear Constraints and Ef- ficien t Regularizations. URL https://github.com/ pnnl/neuromancer . Drgo ˇ na, J., T uor, A., and V rabie, D. (2024). Learning constrained parametric differentiable predictive control p olicies with guarantees. IEEE T r ansactions on Sys- tems, Man, and Cyb ernetics: Systems , 54(6), 3596–3607. Dua, V., Bozinis, N.A., and Pistikopoulos, E.N. (2002). A multiparametric programming approac h for mixed- in teger quadratic engineering problems. Computers & Chemic al Engine ering , 26(4), 715–733. F an, C. and Zhou, X. (2023). Mo del-based predictive con trol optimization of c hiller plan ts with w ater-side economizer system. Ener gy and Buildings , 278, 112633. F aulkner, C.A., Ho, J., Shi, C., F an, C., Ildiri, N., and Zuo, W. (2025). Developmen t of near-optimal adv anced con trol sequences for chiller plants with w ater-side econ- omizers in US climates (ashrae rp-1661). Scienc e and T e chnolo gy for the Built Envir onment , 31(1), 18–35. Griew ank, A. and W alther, A. (2008). Evaluating Deriva- tives . So ciety for Industrial and Applied Mathematics. Gurobi Optimization, LLC (2024). Gurobi Optimizer Reference Man ual. URL https://www.gurobi.com . He, H., Daume I I I, H., and Eisner, J.M. (2014). Learning to search in branc h and b ound algorithms. A dvanc es in neur al information pr o c essing systems , 27. He, K., F u, Q., Lu, Y., W ang, Y., Luo, J., W u, H., and Chen, J. (2023). Predictive control optimization of chiller plants based on deep reinforcement learning. Journal of Building Engine ering , 76, 107158. Huang, S., Zuo, W., and Sohn, M.D. (2016). Amelioration of the co oling load based chiller sequencing control. Applie d Ener gy , 168, 204–215. Huang, S., Zuo, W., and Sohn, M.D. (2017). Improv ed co oling tow er control of legacy chiller plan ts by opti- mizing the condenser w ater set point. Building and Envir onment , 111, 33–46. Karg, B. and Lucia, S. (2018). Deep learning-based em b edded mixed-integer mo del predictive con trol. In 2018 Eur op e an Contr ol Confer enc e , 2075–2080. Kirc hes, C. (2011). F ast Numeric al Metho ds for Mixe d- Inte ger Nonline ar Mo del-Pr e dictive Contr ol . Vieweg. Ko c h, J., Shapiro, M., Sharma, H., V rabie, D., and Drgona, J. (2025). Learning neural differen tial algebraic equa- tions via op erator splitting. In Confer enc e on De cision and Contr ol (CDC) . Liu, Z., T an, H., Luo, D., Y u, G., Li, J., and Li, Z. (2017). Optimal chiller sequencing control in an office building considering the v ariation of chiller maximum co oling capacit y . Ener gy and Buildings , 140, 430–442. Ma, Z. and W ang, S. (2011). Sup ervisory and optimal con trol of central c hiller plants using simplified adaptiv e mo dels and genetic algorithm. Applie d Ener gy , 88(1), 198–211. Marcucci, T. and T edrake, R. (2021). W arm start of mixed-in teger programs for model predictive control of h ybrid systems. IEEE T r ansactions on Automatic Contr ol , 66(6), 2433–2448. Masti, D. and Bemp orad, A. (2019). Learning bi- nary w arm starts for m ultiparametric mixed-in teger quadratic programming. In 18th Eur op e an Contr ol Con- fer enc e , 1494–1499. McAllister, R.D. and Ra wlings, J.B. (2022). Adv ances in mixed-integer mo del predictive control. In 2022 A meric an Contr ol Confer enc e , 364–369. P an, C. and Li, Y. (2024). Nonlinear mo del predictive con trol of chiller plant demand resp onse with koopman bilinear mo del and krylov-subspace mo del reduction. Contr ol Engine ering Pr actic e , 147, 105936. Pusk orius, G. and F eldk amp, L. (1994). T runcated back- propagation through time and Kalman filter training for neuro control. In 1994 International Confer enc e on Neur al Networks , volume 4, 2488–2493. Qin, S.J. and Badgw ell, T.A. (2003). A survey of industrial mo del predictive control techn ology . Contr ol engine er- ing pr actic e , 11(7), 733–764. Rac k auck as, C., Ma, Y., Martensen, J., W arner, C., Zub ov, K., Sup ek ar, R., Skinner, D., Ramadhan, A., and Edel- man, A. (2021). Universal differen tial equations for scien tific mac hine learning. URL abs/2001.04385 . Reiter, R., Quirynen, R., Diehl, M., and Di Cairano, S. (2024). Equiv ariant deep learning of mixed-integer optimal con trol solutions for vehicle decision making and motion planning. IEEE T r ansactions on Contr ol Systems T e chnolo gy . Ric hards, A. and How, J. (2005). Mixed-in teger program- ming for con trol. In 2005 Americ an Contr ol Confer enc e , 2676–2683. Sala-Cardoso, E., Delgado-Prieto, M., Kamp ouropoulos, K., and Romeral, L. (2020). Predictive chiller op eration: A data-driven loading and scheduling approac h. Ener gy and Buildings , 208, 109639. Shan, K., W ang, S., and Zhuang, C. (2021). Con trolling a large constant speed centrifugal c hiller to provide grid frequency regulation: A v alidation based on onsite tests. Applie d Ener gy , 300, 117359. T ak ap oui, R., Mo ehle, N., Boyd, S., and Bemp orad, A. (2020). A simple effective heuristic for embedded mixed- in teger quadratic programming. International Journal of Contr ol , 93(1), 2–12. T ang, B., Khalil, E.B., and Drgo ˇ na, J. (2025). Learning to optimize for mixed-in teger non-linear programming. URL . T erzi, E., Bonetti, T., Saccani, D., F arina, M., F agiano, L., and Scattolini, R. (2020). Learning-based predictive con trol of the co oling system of a large business cen tre. Contr ol Engine ering Pr actic e , 97, 104348. W ang, L., Lee, E.W.M., Y uen, R.K., and F eng, W. (2019). Co oling load forecasting-based predictive optimisation for c hiller plants. Ener gy and Buildings , 198, 261–274. Zarp ellon, G., Jo, J., Lo di, A., and Bengio, Y. (2021). P arameterizing branch-and-bound searc h trees to learn branc hing p olicies. In 35th AAAI Confer enc e on Artifi- cial Intel ligenc e .