arXiv CSopen access

Scalable Physics-Informed Neural Differential Equations and Data-Driven Algorithms for HVAC Systems

Scalable Physics-Informed Neural Differential Equations and Data-Driven Algorithms for HVAC Systems Hanfeng Zhai1 *, Hongtao Qiao2†, Hassan Mansour2 , Christopher Laughman2

arXiv:2604.18438v1 [cs.LG] 20 Apr 2026

1 Department of Mechanical Engineering,

Stanford University, Stanford, CA 2 Mitsubishi Electric Research Laboratories, Cambridge, MA

April 21, 2026 Abstract We present a scalable, data-driven simulation framework for large-scale heating, ventilation, and air conditioning (HVAC) systems that couples physics-informed neural ordinary differential equations (PINODEs) with differential-algebraic equation (DAE) solvers. At the component level, we learn heat-exchanger dynamics using an implicit PINODE formulation that predicts conserved quantities (refrigerant mass Mr and internal energy Ehx ) as outputs, enabling physics-informed training via automatic differentiation of mass/energy balances. Stable long-horizon prediction is achieved through gradient-stabilized latent evolution with gated architectures and layer normalization. At the system level, we integrate learned components with DAE solvers (IDA and DASSL) that explicitly enforce junction constraints (pressure equilibrium and mass-flow consistency), and we use Bayesian optimization to tune solver parameters for accuracy–efficiency trade-offs. To reduce residual system-level bias, we introduce a lightweight corrector network trained on short trajectory segments. Across dual-compressor and scaled network studies, the proposed approach attains multi-fold speedups over highfidelity simulation while keeping errors low (MAPE below a few percent) and scales to systems with up to 32 compressor–condenser pairs.

Keywords: HVAC systems, Neural ODE, Gated recurrent unit, Vapor compression system, Differential algebraic equation

* Work done during internship at MERL; E-mail: [email protected] † Corresponding author; E-mail: [email protected]

1

Highlights • Implicit physics-informed neural ODE formulation that treats conserved quantities as outputs, enabling physics-informed training via automatic differentiation. • Gradient-stabilized latent evolution with gated architectures ensures stable training and long-horizon predictions for stiff thermo-fluid systems. • Lightweight corrector network compensates for system-level errors, improving mass and energy prediction accuracy. • DAE solvers (IDA and DASSL) with adaptive high-resolution stepping explicitly enforce algebraic constraints for pressure equilibrium and mass flow conservation. • Scalable to systems with 32 compressor–condenser pairs, representing one of the largest data-driven HVAC surrogates reported. • Achieves 4–9× speedup with MAPE < 2.5% compared to high-fidelity simulators (best: 2.04% MAPE, 58.91 s). • Bayesian optimization automatically tunes solver parameters across multiple solver types for optimal accuracy–efficiency trade-offs.

2

Contents 1

Introduction

7

2

Methodologies 2.1 Physics-Informed Neural ODE . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.1.1 Implicit Formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.1.2 Model Framework . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.1.3 Physics-Informed Loss . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.1.4 Training Algorithms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 Corrector Network for Cyclic Systems . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2.1 Network Architecture and Formulation . . . . . . . . . . . . . . . . . . . . . . 2.2.2 Training Regime Selection . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2.3 Deployment to HVAC Systems . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3 System Solver . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3.1 General System Solver Framework . . . . . . . . . . . . . . . . . . . . . . . . . 2.3.2 Algebraic Solver . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3.3 DAE Solver and Formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3.4 Integration with PINODE Framework . . . . . . . . . . . . . . . . . . . . . . . 2.3.5 Bayesian Optimization for Parameter Tuning . . . . . . . . . . . . . . . . . . . 2.3.6 Adaptive High-Resolution Stepping . . . . . . . . . . . . . . . . . . . . . . . .

9 9 9 11 13 14 16 16 16 17 19 20 21 23 24 25 26

3

Results & Discussions 3.1 Training PINODE Models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Corrector Network for System Simulations . . . . . . . . . . . . . . . . . . . . . . . . 3.3 Generalizability in Dual-Compressor Systems . . . . . . . . . . . . . . . . . . . . . . 3.3.1 Bayesian Optimization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3.2 Simulation Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

27 27 28 32 32 36

4

Scalability of HVAC Systems 41 4.1 Scaling Up HVAC . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 4.2 Stability and Solver Parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41

5

Conclusions

44

3

Nomenclature State Variables Mr Ehx p T h y yd ya ẏ

Refrigerant mass within heat exchanger control volume (kg) Internal energy within heat exchanger (J) Pressure (Pa) Temperature (K) Specific enthalpy (J/kg) System state vector Differential component of state vector Algebraic component of state vector Time derivative of state vector

Rate Variables Ṁr Ėhx ṁ ṁr,in ṁr,out Q̇ a

Rate of change of refrigerant mass (kg/s) Rate of change of internal energy (J/s) Mass flow rate (kg/s) Inlet refrigerant mass flow rate (kg/s) Outlet refrigerant mass flow rate (kg/s) Heat transfer rate to air stream (W)

DAE Solver Variables F G J hn ∆t hmax hmin ∆tIDA out ∆tDASSL out,min DASSL Nmax ∆ttarget tn k αi,k ATOL RTOL ϵ∆t ϵsoln τ npre npost

DAE residual function Nonlinear system function Jacobian matrix Time step size at step n Time step size Maximum time step size Minimum time step size IDA solver output interval DASSL minimum output interval DASSL maximum number of internal steps Target output interval for adaptive stepping Time at step n BDF order BDF coefficients (subscript i for coefficient index, subscript k for order) Absolute tolerance Relative tolerance Time step tolerance Solution tolerance for root-finding Control jump detection threshold Number of steps before control jump (high-res) Number of steps after control jump (high-res)

4

Neural Network Variables h( t ) h0 Xenc Senc Sdec x(t) s( t ) ŷ W θ dlatent dx ds Tenc Tdec Φenc ψ ζ fθ ϕcorr zin ϕcorr raw ϕcorr GP ϕcorr smooth ϕcorr m̂pred m̂bench m̂corr din

Latent state vector in neural ODE Initial latent state Encoder input sequence Encoder state sequence Decoder state sequence Time-dependent feature vector Time-dependent state vector Predicted output vector Neural network weight matrices Neural network parameters Latent dimension Feature dimension State dimension Encoder sequence length Decoder sequence length Encoder function (GRU) Decoder function (GRU) Latent projection network Neural ODE dynamics function Corrector network function Corrector network input vector Corrector network output (correction term) Raw correction from network GP-smoothed correction Final smoothed correction Normalized predicted mass/energy values Normalized benchmark mass/energy values Normalized corrected mass/energy values Corrector network input dimension

Loss Functions Ltotal Ldata Lphys Lcons Lcorr λphys λcons N

Total training loss Data fidelity loss Physics-informed loss Conservation loss Corrector network loss Physics loss weight Conservation loss weight Batch size

5

System Parameters nc nv ncond nevap np n Nsteps pliq psuct G(nc , nv )

Number of compressors (or compressor–condenser pairs) Number of valves (or valve–evaporator pairs) Number of condensers Number of evaporators Number of pressure variables (junction pressures) Total system dimension Number of simulation time steps Liquid manifold pressure node Suction manifold pressure node System topology graph with nc compressors and nv evaporators

Mathematical Symbols σ ⊙ R ∥·∥ ∂ ∇ MSE MAPE arg min D M X

Sigmoid activation function Element-wise multiplication Set of real numbers Euclidean norm Partial derivative Gradient operator Mean squared error Mean absolute percentage error Argument of the minimum Dataset Model Parameter space

Abbreviations PINODE DAE ODE DASSL BDF GRU RK4 HVAC GP EI

Physics-Informed Neural ODE Differential-Algebraic Equation Ordinary Differential Equation Differential-Algebraic System Solver Backward Differentiation Formula Gated Recurrent Unit Fourth-order Runge-Kutta method Heating, Ventilation, and Air Conditioning Gaussian Process Expected Improvement

6

1

Introduction

Heating, ventilation, and air conditioning (HVAC) systems account for a substantial portion of global energy consumption, with residential and commercial buildings responsible for approximately 40% of total energy usage [1]. As the demand for energy-efficient building operations grows, the need for accurate, real-time modeling and control of HVAC systems has become increasingly critical [2]. Large-scale HVAC networks, comprising multiple compressors, condensers, evaporators, and valves, exhibit complex thermo-fluid dynamics governed by nonlinear differential-algebraic equations (DAEs) that couple mass, energy, and momentum conservation laws [3]. Traditional approaches to HVAC modeling and simulation face a fundamental trade-off: high-fidelity physics-based models implemented in commercial software (e.g., Modelica [4], EnergyPlus [5]) achieve high accuracy but require computationally expensive numerical solvers that preclude real-time applications, while simplified reduced-order models sacrifice physical consistency for computational efficiency. Current physics-based simulators, despite their accuracy, suffer from significant computational bottlenecks that limit their practical deployment. High-fidelity simulations of large-scale HVAC systems are computationally intensive [6], make it challenging for real-time control, optimization, or co-simulation with building energy management systems [7]. The computational burden stems from the need to solve stiff, nonlinear DAEs with adaptive time-stepping, handle complex thermodynamic property calculations, and resolve fast transients in multi-component networks [3, 8, 9]. This research gap between accuracy and computational efficiency has motivated the development of data-driven surrogates that can approximate the behavior of physics-based simulators at a fraction of the computational cost [10, 11, 12, 13, 14]. To bridge this gap, we propose a hybrid modeling framework that combines the expressiveness of neural networks with the physical consistency of conservation laws, enabling both accurate predictions and efficient computation. Our approach integrates physics-informed neural ordinary differential equations (PINODEs) [15] for component-level dynamics with specialized DAE solvers [16, 17] for system-level constraint satisfaction. This dual-level strategy allows us to leverage the flexibility of machine learning for capturing complex, data-driven patterns while explicitly enforcing fundamental physical principles, such as mass and energy conservation, that are essential for thermodynamic validity [18, 19]. Existing machine learning surrogates for HVAC systems face several critical limitations that compromise their reliability and applicability. Conventional neural network approaches, such as feedforward networks or recurrent architectures (LSTMs, GRUs), could produce predictions that violate conservation laws, leading to unphysical results that accumulate errors over long prediction horizons [20]. Similarly, standard neural ODE formulations [21] struggle with numerical stability during training and inference, particularly when dealing with stiff dynamics or long time horizons [22, 23]. On the system-level side, most existing solvers for HVAC networks either simplify the algebraic constraints to make the system easier to solve (losing accuracy) or employ generic ODE integrators that cannot properly handle the differential-algebraic structure, which may lead to constraint violations and pressure imbalances [10]. Our framework addresses these limitations through gradient-stabilized latent evolution in the neural ODE component and explicit algebraic constraint enforcement in the DAE solver, ensuring both stability and physical consistency. To solve these challenges, we draw upon recent advances in continuous-time neural models and differential-algebraic equation solvers. Neural ODEs [21] and physics-informed neural networks [20, 24, 25, 26] have demonstrated how to incorporate physical laws, though

7

existing work often combines these ideas superficially without deeply embedding conservation laws into the latent dynamics or addressing stability issues in stiff systems. For the system-level solver, DASSL (Differential-Algebraic System Solver) [17] provides a robust foundation for handling index-1 DAEs, but its application to large-scale HVAC networks requires careful parameter tuning and adaptive time-stepping strategies. To ensure fair and systematic comparison across different solver types (algebraic, DAE-IDA, and DAE-DASSL), we employ Bayesian optimization [27, 28, 12, 3] to systematically explore the parameter space and characterize each solver’s performance, enabling objective evaluation of accuracy and computational efficiency. The key novelty of our proposed methodology consists of four core innovations: • Implicit PINODE formulation with conservation-aware states: Conserved quantities (refrigerant mass Mr and internal energy Ehx ) are treated as model outputs rather than prescribed inputs, enabling direct use of conservation structure in learning and downstream system simulation. • Stable long-horizon latent dynamics: A gradient-stabilized evolution scheme with gated architectures and layer normalization [29] improves stability for stiff thermo-fluid dynamics and supports accurate long-horizon prediction. • DAE-based system integration with explicit constraint enforcement: A specialized systemlevel DAE formulation explicitly enforces junction constraints (e.g., pressure equilibrium and mass flow consistency), improving thermodynamic consistency across coupled components. • Automated solver tuning and scalability assessment: Bayesian optimization is used to tune solver parameters for fair accuracy–efficiency comparison, and the framework is demonstrated on large-scale HVAC configurations. These innovations enable scalable and thermodynamically consistent modeling of large HVAC configurations, while achieving 4–9× computational speedup with MAPE < 2.5% relative to highfidelity simulators. The remainder of this paper is organized as follows. Section 2.1 presents the physics-informed neural ODE framework for heat exchanger modeling, including the encoder-decoder architecture, latent dynamics formulation, and physics-informed training strategy. Section 2.3 describes the DAE solver for large-scale HVAC systems, detailing the differential–algebraic formulation, DASSL integration, and Bayesian optimization for parameter tuning. Section 3 presents experimental results on real HVAC datasets, demonstrating accuracy, efficiency, and scalability. Section 4 discusses the scalability analysis and computational performance. Finally, Section 5 summarizes the contributions and discusses future research directions.

8

2

Methodologies

2.1

Physics-Informed Neural ODE

We develop a physics-informed neural ordinary differential equation (PINODE) framework for modeling heat exchanger dynamics in HVAC systems. The approach combines the expressiveness of neural networks with the physical consistency of conservation laws, enabling accurate longhorizon predictions while maintaining thermodynamic validity. 2.1.1

Implicit Formulation

A key design decision in our framework is the treatment of conserved quantities. Specifically, the refrigerant mass Mr and internal energy Ehx are treated as model outputs rather than known inputs. This implicit formulation contrasts with the explicit approach in [10], where mass and energy are provided as inputs to the heat exchanger model, following a discrete-time dynamical system form with fixed time steps. The heat exchanger model maps an 8-dimensional input vector xin ∈ R8 to a 9-dimensional output vector ŷ ∈ R9 : ŷ = MPINODE ( xin , sin ), (1) where the input vector xin contains the boundary conditions: xin = [ Ta,in , ϕa,in , ṁ a , Pamb , ṁr,in , hr,in , hr,out , Pr,out ]⊤ ,

(2)

comprising air-side properties (inlet temperature Ta,in , humidity ratio ϕa,in , mass flow rate ṁ a , ambient pressure Pamb ) and refrigerant-side properties (inlet mass flow rate ṁr,in , inlet enthalpy hr,in , outlet enthalpy hr,out , outlet pressure Pr,out ). The output vector ŷ includes both observable quantities and the conserved state variables:  ⊤ ŷ = p1 , p N , h1 , h N , Ta,out , Q̇ a , Mr , Ehx , Q̇lat ,

(3)

where p1 and p N are refrigerant pressures at the first and last nodes, h1 and h N are corresponding enthalpies, Ta,out is the air outlet temperature, Q̇ a is the heat transfer rate to air, Mr and Ehx are the conserved quantities (refrigerant mass and internal energy), and Q̇lat is the latent heat transfer rate. The key distinction from the explicit formulation in [10] lies in the treatment of Mr and Ehx . In explicit the explicit approach, these quantities are provided as inputs: xin = [ xin , Mr , Ehx ]⊤ ∈ R10 , and the model predicts only the observable outputs. In contrast, our implicit formulation treats Mr and Ehx as outputs, enabling direct computation of their time derivatives through automatic differentiation: d d pred pred (4) Ṁr = M̂r , Ėhx = Êhx , dt dt where the derivatives are computed with respect to time through the neural network’s computational graph. This capability enables the incorporation of gradient information into physicsinformed loss terms that enforce conservation laws during training, as the model must predict rates of change that are consistent with mass and energy balance: Ṁrtrue = ṁr,in − ṁr,out ,

true Ėhx = ṁr,in hr,in − ṁr,out hr,out − Q̇ a .

9

(5)

Di!erential Algebraic Equation Di!erential: F(t, y, ẏ) = ẏ → f (t, y), Algebraic: g(p, y, t) = 0, Di!erential variables: ! " (n +n ) → (1) (2) Mr , Mr , . . . , Mr cond evap y= . (n +n ) (1) (2) Ehx , Ehx , . . . , Ehxcond evap

Figure 1: Heat exchanger model inputs and outputs for the implicit PINODE formulation. The top panel shows the 8 input variables (air and refrigerant boundary conditions), while the bottom panel shows the 9 output variables, including the conserved quantities Mr and Ehx . The figure displays a representative 100-step segment (steps 100–210) of a simulation; in practice, the model can simulate thousands of time steps.

1

States advanced by solver: Mr , Ehx Air out

Air in Heat exchanger (control volume)

Refrigerant in

Refrigerant out

Figure 2: Simplified heat exchanger schematic (control volume view) illustrating the coupling between air-side and refrigerant-side flows. The conserved quantities Mr and Ehx are the differential states advanced by the system-level solver.

10

This physics-informed training approach is essential for learning dynamics that respect therpred pred modynamic constraints, ensuring that the predicted rates Ṁr and Ėhx align with the true true . conservation-based rates Ṁrtrue and Ėhx Figure 1 illustrates the input–output structure of the implicit PINODE formulation for heat exchanger modeling. The model receives 8 input variables that characterize the boundary conditions: air-side properties (inlet temperature, humidity ratio, mass flow rate, ambient pressure) and refrigerant-side properties (inlet and outlet mass flow rates, enthalpies, and outlet pressure). The PINODE model processes these inputs through its encoder–neural ODE–decoder architecture to predict 9 output variables, including the conserved quantities Mr and Ehx that are treated as differential states in the system-level DAE formulation. The figure shows a representative 100-step segment (steps 100–210) of a simulation trajectory; in practice, the model can simulate thousands of time steps to capture long-horizon system dynamics. The implicit formulation is fundamentally more physically meaningful: in real-world HVAC simulations, the mass and energy within each heat exchanger are unknown state variables that emerge from solving the system dynamics, not known inputs that can be directly measured or specified. Although this introduces training challenges due to potential error accumulation, we address these through careful architecture design, including gradient-stabilized latent evolution and regularization strategies detailed in the following sections. 2.1.2

Model Framework

The PINODE architecture consists of three main components: an encoder, a neural ODE, and a decoder, as illustrated in Figure 3. The encoder processes historical input sequences to extract a compact latent representation, which is then evolved forward in time through a continuoustime neural ODE that captures the underlying dynamics. The decoder maps the evolved latent trajectory back to the observable output space, enabling predictions of heat exchanger states and thermodynamic properties. Given encoder sequences Xenc ∈ RTenc ×dx and decoder sequences Sdec ∈ RTdec ×ds , where d x and ds denote feature and state dimensions respectively, the model operates as follows. Encoder: A GRU encoder Φenc processes the input sequence to extract a latent representation that captures the temporal dependencies in the input history. The encoder takes as input the feature sequence Xenc and state sequence Senc , processing them through recurrent layers to produce a summary representation: zT = Φenc ( Xenc , Senc ), (6) where Senc represents the encoder state sequence and zT is the final hidden state after processing all Tenc time steps. The encoder output is then projected into a lower-dimensional latent space via a feedforward network ζ to initialize the latent trajectory: ζ 0 = ζ (zT ) ∈ Rdlatent ,

(7)

where dlatent is the latent dimension (typically dlatent = 8), and ζ 0 serves as the initial condition for the neural ODE integration. Neural ODE: The latent trajectory ζ (t) evolves continuously in time according to a neural ODE, which learns the underlying dynamics of the heat exchanger system in a low-dimensional latent space. The evolution is governed by: dζ = f θ (ζ (t), x(t), s(t)), dt 11

(8)

Algorithm 1 PINODE Forward Pass for Heat Exchanger Prediction Require: Encoder sequences Xenc , Senc ; decoder sequences Xdec , Sdec ; time step ∆t Ensure: Predicted outputs ŷ for each decoder time step 1: Step 1: Encoding - Extract latent representation from input history 2: zT ← Φenc ( Xenc , Senc ) ▷ GRU encoder processes input sequence 3: ζ 0 ← ζ (zT ) ▷ Project to latent space Rdlatent 4: Step 2: Latent Evolution - Integrate neural ODE over decoder horizon 5: for each decoder time step t = 1, 2, . . . , Tdec do 6: Extract feature and state vectors: xt ← Xdec [t], st ← Sdec [t] 7: Compute latent derivative: dζ ▷ Neural ODE with gated architecture dt ← f θ ( ζ t−1 , xt , st ) 8: Integrate latent trajectory: ζ t ← RK4( f θ , ζ t−1 , xt , st , ∆t) ▷ 4th-order Runge-Kutta 9: end for 10: Step 3: Decoding - Map evolved latent states to predictions 11: for each decoder time step t = 1, 2, . . . , Tdec do 12: ŷt ← ψ(ζ t ) ▷ GRU decoder maps latent to observable space 13: end for 14: Step 4: Physics Constraints - Compute derivatives for physics loss pred pred 15: Ṁr , Ėhx ← dtd [ M̂r , Êhx ] ▷ Automatic differentiation 16: return Predicted outputs ŷ = [ ŷ1 , . . . , ŷTdec ] and derivatives where f θ is a neural network parameterized by θ that defines the vector field in the latent space, and x(t) and s(t) are the time-dependent feature and state vectors at the decoder input that provide external forcing and boundary conditions. The neural ODE uses a gated architecture inspired by GRU cells to ensure stable gradient flow and prevent numerical instabilities during long-horizon integration: r t = σ ( Wr [ x t , s t ] + Ur ζ t ) , zt = σ (Wz [ x t , st ] + Uz ζ t ), (9) ζ̃ t = tanh(Wh [ xt , st ] + Uh (rt ⊙ ζ t )), dζ = (1 − zt ) ⊙ (ζ̃ t − ζ t ), dt where σ denotes the sigmoid function, ⊙ is element-wise multiplication, and Wr , Wz , Wh , Ur , Uz , Uh are learnable weight matrices. The reset gate rt controls how much of the previous latent state is retained, while the update gate zt balances between the previous state and the candidate update ζ̃ t . This gated formulation mitigates vanishing and exploding gradients by constraining the rate of change in the latent space, enabling stable training and accurate long-horizon predictions. The latent trajectory is integrated forward in time using a fourth-order Runge-Kutta (RK4) scheme, which provides high accuracy while maintaining computational efficiency: ζ t+∆t = RK4( f θ , ζ t , xt , st , ∆t),

(10)

where ∆t is the time step size. The RK4 method evaluates the vector field f θ at four intermediate points within each time step, combining these evaluations to produce a fourth-order accurate approximation of the latent state at the next time point. Our work is inspired by [15]. Decoder: A GRU decoder ψ maps the evolved latent trajectory ζ (t) back to the observable output space, transforming the low-dimensional latent representation into predictions of heat exchanger 12

states and thermodynamic properties: ŷ = ψ(ζ (t)),

(11)

where ŷ contains the predicted outputs including pressures, enthalpies, temperatures, heat transfer rates, and the conserved quantities (mass and energy). To improve training stability, we employ layer normalization [29] at multiple stages: after the encoder, in the latent space, and before the final output layer. Dropout regularization is applied to prevent overfitting during training.

x

φ

h

ζ=

Z T

H(ζ) dt

ψ

y

0

GRU Encoder: φ(x)

MLP Decoder: h(φ)

Latent Dynamics: ζ̇ = H(ζ)

GRU Decoder: ψ(ζ)

Figure 3: Schematic of the physics-informed neural ODE (PINODE) architecture for heat exchanger modeling. The model encodes input sequences into a latent space, evolves the latent dynamics using a neural ODE, and decodes back to observable states. Physics constraints are enforced through a physics-informed loss term that penalizes violations of mass and energy conservation. 2.1.3

Physics-Informed Loss

The training objective combines data fidelity with physics constraints to ensure the learned dynamics respect thermodynamic conservation laws. For a heat exchanger, the key conserved quantities are the refrigerant mass Mr and internal energy Ehx within the control volume. The total loss function is defined as:

Ltotal = Ldata + λphys Lphys + λcons Lcons ,

1 N ∥ŷi − yi ∥2 , N i∑ =1 i 1 N h pred pred true 2 true 2 ∥ Ṁ − Ṁ ∥ + ∥ Ė − Ė ∥ Lphys = , r r hx hx N i∑ =1

Ldata =

Lcons =

(12)

 1 N  true 2 ∥ M̂r − Mrtrue ∥2 + ∥ Êhx − Ehx ∥ , ∑ N i =1

where λphys and λcons are weighting hyperparameters (typically set to 0.5), N is the batch size, and ŷi and yi are the predicted and ground-truth state vectors, respectively.

13

The physics loss Lphys enforces consistency between the model’s predicted rates of change and the true rates computed from conservation laws. For a heat exchanger control volume, the mass and energy rates are: Ṁrtrue = ṁr,in − ṁr,out , (13) true Ėhx = ṁr,in hr,in − ṁr,out hr,out + Q̇ a , where ṁr,in and ṁr,out are the inlet and outlet mass flow rates, hr,in and hr,out are the corresponding specific enthalpies, and Q̇ a is the heat transfer rate to the air stream. pred

pred

The predicted rates Ṁr and Ėhx are computed via automatic differentiation of the model output with respect to time, leveraging the continuous-time formulation of the neural ODE: dŷ ∂ψ dζ = , dt ∂ζ dt

(14)

where the derivatives are evaluated using the chain rule through the decoder and neural ODE. This automatic differentiation capability is a key advantage of the continuous-time formulation, as it enables direct computation of time derivatives without finite difference approximations, ensuring smooth and accurate gradient information for the physics-informed loss terms. Since Mr and Ehx are treated as algebraic variables in the implicit formulation, the conservation loss Lcons provides additional supervision to ensure these quantities match the ground truth. 2.1.4

Training Algorithms

The training procedure alternates between validation and training phases within each epoch, as detailed in Algorithm 2. During the validation phase, the model is evaluated on held-out data to monitor generalization performance and adjust the learning rate. The training phase performs gradient-based optimization with physics-informed loss terms. Key implementation details: • Optimizer: Adam optimizer with an initial learning rate of 10−3 , reduced by a factor of 0.5 when validation loss plateaus (patience of 25 epochs). • Gradient clipping: Gradients are clipped to a maximum norm of 1.0 to prevent instability during training. • Physics weight scheduling: The physics loss weight λphys is set to 0.5 throughout training, balancing data fidelity with physical consistency. • Batch processing: Sequences are processed in batches, with encoder sequences of length Tenc and decoder sequences of length Tdec . The algorithm uses an implicit problem formulation, where algebraic constraints are enforced through the conservation loss term, which penalizes deviations in the conserved quantities Mr and Ehx . The PINODE framework provides a principled approach to learning HVAC dynamics that respects physical constraints while maintaining the flexibility of neural network models. By embedding conservation laws directly into the training objective, the model learns representations that generalize better to unseen operating conditions and maintain thermodynamic consistency over long prediction horizons. The stabilized latent dynamics, achieved through gated architectures 14

Algorithm 2 Physics-Informed Neural ODE Training for Heat Exchanger Require: Training data Dtrain , validation data Dval , model M (PINODE) Ensure: Trained model M∗ with optimized parameters 1: Initialize M with random weights W, optimizer (Adam), learning rate scheduler 2: Set λphys = 0.5, λcons = 0.5 3: for epoch e = 1, 2, . . . , E do 4: Validation: Set model to evaluation mode 5: for each batch ( Xenc , Senc , Xdec , Sdec , ∆t, Y ) ∈ Dval do 6: Ŝ ← M( Xenc , Senc , Xdec , Sdec , ∆t) 7: Ldata ← MSE(Ŝ, Y ) true from physics; extract Ṁpred , Ėpred from model 8: Compute Ṁrtrue , Ėhx r hx pred

pred

true ) Lphys ← MSE( Ṁr , Ṁrtrue ) + MSE( Ėhx , Ėhx 10: Lcons ← MSE(Ŝ[ Mr ], Y [ Mr ]) + MSE(Ŝ[ Ehx ], Y [ Ehx ]) 11: Lval ← Ldata + λphys · Lphys + λcons · Lcons 12: end for ∗ , save M∗ ← M 13: Update learning rate; if Lval < Lval 14: Training: Set model to training mode 15: for each batch ( Xenc , Senc , Xdec , Sdec , ∆t, Y ) ∈ Dtrain do 16: Ŝ ← M( Xenc , Senc , Xdec , Sdec , ∆t) 17: Compute Ldata , Lphys , Lcons as in validation 18: Ltotal ← Ldata + λphys · Lphys + λcons · Lcons total 19: ∇W ← ∂L∂W ; clip if needed; update W ← W − η · ∇W 20: end for 21: end for 22: return Best model M∗ with lowest validation loss

9:

15

and layer normalization, enable stable gradient flow during backpropagation through the neural ODE, addressing a key challenge in training continuous-time neural models. In the following sections, we demonstrate how this framework integrates with DAE solvers to handle algebraic constraints in large-scale HVAC systems. Compressor and expansion valve In addition to the PINODE heat exchanger models, the compressor and expansion valve are modeled as static (memoryless) mappings implemented with lightweight multilayer perceptrons (MLPs). These models take local inlet/outlet thermodynamic conditions and actuation signals (e.g., compressor speed or valve opening) as inputs and output the corresponding mass flow rates and outlet enthalpies needed by the system solver. We do not elaborate these static models in detail here; the training setup and feature definitions follow [10].

2.2

Corrector Network for Cyclic Systems

While physics-informed neural ODEs (PINODEs) provide accurate predictions for individual heat exchangers, system-level simulations involving multiple interconnected components can accumulate errors due to algebraic coupling constraints and numerical integration inaccuracies. To address this, we introduce a corrector neural network that learns to compensate for systematic biases in mass and energy predictions at the system level. The corrector network is trained on a short initial segment of the simulation trajectory and then applied throughout the entire simulation to improve prediction accuracy. 2.2.1

Network Architecture and Formulation

The corrector network ϕcorr : Rdin → R4 is a fully connected neural network that takes as input the concatenated normalized state and output vectors from both indoor heat exchangers (for dualcompressor simulation):   x̂HEX,1 ŷHEX,1  din  zin =  (15)  x̂HEX,2  ∈ R , ŷHEX,2

where x̂HEX,i and ŷHEX,i are the normalized state and output vectors for the i-th heat exchanger, i h (1) (2) (1) (2) ⊤ for the four respectively. The network outputs a correction term ϕcorr = ϕE , ϕE , ϕM , ϕM target variables: internal energy E and refrigerant mass M for both compressors. The network architecture consists of three fully connected layers with sigmoid activations: ϕcorr (zin ) = ConstrainTanh (W3 σ(W2 σ(W1 zin + b1 ) + b2 ) + b3 ) ,

(16)

where σ denotes the sigmoid activation function, Wi and bi are the weight matrices and bias vectors, and ConstrainTanh applies a scaled hyperbolic tangent to constrain the output within a specified range. The hidden layer dimension matches the outdoor heat exchanger model’s hidden size to maintain consistency with the overall architecture. 2.2.2

Training Regime Selection

The corrector network is trained on a short initial segment of the simulation trajectory, typically spanning 850 time steps starting from step 950. This training regime is selected to capture the 16

system’s transient behavior and initial coupling dynamics while keeping the computational cost minimal. During training, the dual-compressor system simulation is run for this initial segment, collecting pairs of predicted and benchmark (Dymola) values for the mass and energy variables. The training objective minimizes the mean squared error between the corrected predictions and the benchmark values:

Lcorr =

o 2 1 T n (t) (t) (t) m̂ + ϕ ( z ) − m̂bench , corr in ∑ pred T t =1

(17)

h i (t) (t) (t) (t) (t) ⊤ where m̂pred = Ê1 , Ê2 , M̂1 , M̂2 are the normalized predicted mass and energy values at (t)

time step t, m̂bench are the corresponding normalized benchmark values from Dymola, and T is the number of training time steps. The loss function penalizes deviations between the corrected predictions and the reference benchmark, encouraging the network to learn systematic correction patterns that compensate for accumulated errors.

The training procedure, detailed in Algorithm 3, begins by running the dual-compressor system simulation for the selected training segment (typically 850 steps starting from step 950) using the physics-informed algebraic solver. This generates a sequence of system states, from which we extract the normalized state and output vectors from both condensers to form the input vectors (t) zin . Simultaneously, we collect the corresponding benchmark mass and energy values from the high-fidelity Dymola simulation, which serve as the ground truth targets. The network is trained using the Adam optimizer with a learning rate of 10−3 for 500,000 epochs, ensuring convergence to a stable correction mapping. The training process iteratively adjusts the network weights to minimize the correction loss, learning to predict corrections that, when added to the raw predictions, bring them closer to the benchmark values. This supervised learning approach allows the corrector network to capture systematic biases that arise from component coupling and numerical integration errors, which are difficult to eliminate through solver parameter tuning alone. 2.2.3

Deployment to HVAC Systems

After training, the corrector network is deployed throughout the entire simulation to improve prediction accuracy. Figure 4 shows the schematic of applying the corrector neural network to the dual-compressor system simulation. raw ( t ). To At each simulation time step, the corrector network generates a raw correction term ϕcorr ensure smooth temporal evolution and reduce high-frequency noise, the corrections are postprocessed using a two-stage smoothing procedure:

1. Gaussian Process smoothing: A Gaussian process regressor with a radial basis function (RBF) kernel and white noise kernel is fitted to the raw correction time series for each output dimension:   ( t − t ′ )2 ′ k (t, t ) = C exp − + σn2 δ(t − t′ ), (18) 2ℓ2 where C = 1.0, ℓ = 2000 s, and σn2 = 0.3 are the kernel hyperparameters. The GP provides a GP ( t ) with uncertainty quantification. smoothed correction ϕcorr

17

Algorithm 3 Corrector Network Training Require: Dual-compressor system with PINODE models, benchmark Dymola data, training segment [tstart , tend ] ∗ Ensure: Trained corrector network ϕcorr 1: Initialize corrector network ϕcorr with random weights 2: Initialize Adam optimizer with learning rate η = 10−3 3: Data Collection Phase 4: Run system simulation for segment [tstart , tend ] using algebraic solver 5: for each time step t ∈ [tstart , tend ] do (t) (t) (t) (t) 6: Extract normalized state/output vectors: x̂HEX,1 , ŷHEX,1 , x̂HEX,2 , ŷHEX,2 7: 8:

(t)

(t)

(t)

(t)

(t)

Form input: zin ← [ x̂HEX,1 ; ŷHEX,1 ; x̂HEX,2 ; ŷHEX,2 ] (t)

(t)

(t)

(t)

(t)

Extract predicted values: m̂pred ← [ Ê1 , Ê2 , M̂1 , M̂2 ] (t)

Extract benchmark values: m̂bench ← from Dymola simulation 10: end for 11: Training Phase 12: for epoch e = 1, 2, . . . , Emax where Emax = 500, 000 do 13: for each time step t ∈ [tstart , tend ] do (t) (t) 14: Compute correction: ϕcorr ← ϕcorr (zin ) 9:

(t)

(t)

(t)

Compute corrected prediction: m̂corr ← m̂pred + ϕcorr

15:

(t)

(t)

(t)

Compute loss: Lcorr ← ∥m̂corr − m̂bench ∥2 17: end for (t) 18: Compute average loss: Lcorr ← T1 ∑t Lcorr 19: Compute gradients: ∇θ ← ∂L∂θcorr 20: Update weights: θ ← θ − η · ∇θ using Adam optimizer 21: end for ∗ 22: return Trained network ϕcorr with optimized parameters θ∗ 16:

Training NN: L = MSE(Cpred + ϕcorr , Csys )

State

xHEX

Normalize

Output

yHEX

Normalize

Extract: C ≡ [Mr , Ehx ]

⊕

Xnorm ∥Ynorm

Corrector NN ϕcorr Post-processing: GP smoothing + EMA

Cpred

+

Ccorr

Update yHEX

corr yHEX

Figure 4: Schematic of the corrector network deployment in the dual-compressor HVAC system simulation. The network takes normalized state and output vectors from both condensers as input and outputs corrections for the mass and energy predictions, which are then smoothed and applied to improve system-level accuracy.

18

2. Exponential moving average: The GP-smoothed corrections are further smoothed using an exponential moving average (EMA) with α = 0.95: smooth smooth GP ϕcorr (t) = αϕcorr (t − 1) + (1 − α)ϕcorr ( t ).

(19)

The final corrected predictions are computed as: (t)

(t)

smooth m̂corr = m̂pred + ϕcorr ( t ),

(20)

which are then denormalized and integrated back into the system state. The correction is only applied if the corrected values remain within the valid normalized range [−1, 1]; otherwise, the uncorrected predictions are used to maintain numerical stability. This approach enables the corrector network to learn and compensate for systematic biases that arise from the coupling between components, numerical integration errors, and model approximations, significantly improving the overall system-level prediction accuracy while maintaining computational efficiency.

2.3

System Solver

While the PINODE framework captures the dynamics of individual heat exchangers, large-scale HVAC systems with multiple interconnected components introduce algebraic constraints that must be satisfied simultaneously. These constraints arise from pressure equilibrium conditions at component junctions, mass and energy conservation across the entire system, and the coupling between differential states (mass charge Mr and internal energy Ehx ) and algebraic variables (pressures, temperatures, mass flow rates). Standard ODE solvers cannot directly handle these algebraic constraints, necessitating a differential-algebraic equation (DAE) formulation.

Figure 5: Schematic of the dual-compressor HVAC system topology. The system consists of two compressors, two condensers (indoor heat exchangers), one evaporator (outdoor heat exchanger), and one expansion valve. The components are interconnected through junctions where pressure equilibrium must be maintained, creating algebraic constraints that require a DAE formulation rather than a standard ODE solver. Figure 5 summarizes the dual-compressor cycle used in our system-level studies. Refrigerant flows through two parallel compressor–condenser branches (indoor heat exchangers), merges at a high-pressure liquid region, expands through a single valve, passes through the outdoor evaporator, and returns through a low-pressure suction path that feeds both compressors. Each compressor sets mass flow and work input into its condenser; the shared liquid and suction regions 19

couple the branches so that junction pressures and split mass flows must be solved together with the differential states Mr and Ehx in each heat exchanger. Those shared manifolds and tees are where mass-flow balance and pressure equilibrium appear as algebraic relations, while Mr and Ehx evolve according to the PINODE and conservation laws described above—precisely the mixed differential–algebraic structure a DAE solver is meant to treat. The same topology is driven in simulation by time-varying actuator commands: compressor speed (or capacity) for each compressor and opening of the expansion valve. Figure 6 plots these actuation signals for the system in Figure 5 over the simulation horizon, so the curves should be read as the boundary inputs applied to the components labeled in the schematic (two compressors, one valve, with dynamics coupled through the interconnecting lines and junctions).

Figure 6: Time-varying actuation signals for the dual-compressor HVAC system of Figure 5. Compressor speeds and valve opening are shown versus time. Rapid changes in these inputs (jumps with |∆u| > 5.0) trigger adaptive high-resolution time stepping in the DAE integration, as detailed in Section 2.3. Together, Figures 5 and 6 specify which components and couplings the solver resolves and which excitations drive the transient: the schematic fixes the network structure and algebraic unknowns, while the actuation traces determine how aggressively that structure is perturbed over time. 2.3.1

General System Solver Framework

The system solver framework provides a unified approach for advancing the HVAC system state in time, regardless of whether an algebraic solver or a DAE solver is employed. The key insight is that both approaches follow the same fundamental sequence: solve algebraic constraints (pressures), evaluate component models (PINODEs), compute state derivatives, and advance the solution. The primary difference lies in how the time advancement is performed: algebraic solvers use explicit time stepping with fixed or adaptive intervals, while DAE solvers use implicit integration 20

with adaptive error control. The general framework operates as follows. At each time step, the system state y (containing mass and energy for each heat exchanger) is known, and the solver must determine the state at the next time point. The algebraic pressure system g(p, y, t) = 0 must be solved first, as the component models (compressors, valves, heat exchangers) require pressure information to compute mass flow rates and other thermodynamic properties. Once pressures are determined, each component model can be evaluated: compressors and valves compute mass flow rates based on pressure differences, while heat exchangers use PINODE models to predict outlet conditions and heat transfer rates. These predictions are then used to compute the time derivatives of mass and energy for each heat exchanger, which are integrated to advance the system state. The algorithm highlights the common structure shared by all solver types: the algebraic pressure system must be solved before component evaluation, component models (especially PINODEs) provide the necessary predictions, and state derivatives are computed from conservation laws. The key distinction is in Step 4: algebraic solvers use explicit time stepping with a fixed or adaptively chosen time step ∆t, while DAE solvers solve an implicit system that enforces the DAE residual F(t, y, ẏ) = 0 using adaptive BDF methods. This unified framework ensures that regardless of the solver choice, the system maintains thermodynamic consistency through proper handling of algebraic constraints and conservation laws. 2.3.2

Algebraic Solver

The default algebraic solver uses the Powell hybrid method (via scipy.optimize.root with method=’hybr’) to solve for junction pressures that satisfy mass flow conservation throughout the system. At each time step, the solver finds the junction pressures p = [ p1 , p2 , . . . , pn p ]⊤ by solving the nonlinear system of residual equations: r(p) = 0,

(21)

where each residual ri enforces mass flow balance at junction i: ri =

ṁin,i − ṁout,i , mscale

(22)

with mscale = 0.01 kg/s used for numerical scaling. The mass flow rates ṁin,i and ṁout,i are computed by evaluating each component’s mass flow model (compressor, valve, or heat exchanger) using the current junction pressures and upstream conditions. The Powell hybrid method combines Powell’s method with a quasi-Newton approach, making it robust for systems with n p unknowns where n p is the number of junction pressures. For large systems with more than 10 pressures, the implementation automatically switches to a bounded least-squares solver (scipy.optimize.least_squares with method=’trf’) with pressure bounds [2 × 105 , 6 × 106 ] Pa to ensure physical feasibility and improve convergence. The algebraic solver employs adaptive time stepping based on a Runge–Kutta–Fehlberg (RK45) method to determine the optimal time step size ∆t at each iteration. The adaptive stepping mechanism leverages the PINODE model’s latent dynamics to predict the time step that maintains numerical accuracy while maximizing computational efficiency. Specifically, the RK45 method computes both 4th-order (z4 ) and 5th-order (z5 ) solutions to the latent ODE: dh = f θ (h(t), x(t), s(t)), dt 21

(23)

Algorithm 4 General System Solver Framework for HVAC Systems Require: Initial state y0 , time horizon [t0 , tend ], actuation signals u(t), solver type (Algebraic: Powell hybrid, DAE-IDA, or DAE-DASSL) Ensure: Solution trajectory {yn , tn }nN=0 1: Initialize: y ← y0 , t ← t0 , n ← 0: 2: Determine time step ∆t or hn based on solver type and adaptive control 3: while t < tend do 4: Step 1: Solve Algebraic Pressure System 5: Solve g(p, y, t) = 0 for junction pressures p 6: Use Powell hybrid method (algebraic solver) or Levenberg–Marquardt (DAE solvers) with tolerance ϵsoln 7: Step 2: Evaluate Component Models 8: for each component i in system do 9: Extract inlet conditions from upstream components and pressures p 10: if component i is heat exchanger then 11: ŷi , Q̇ a,i ← PINODEi ( xin,i , sin,i ) ▷ Predict outputs and heat transfer 12: else if component i is compressor or valve then 13: ṁi , hout,i ← ComponentModeli ( pin , pout , hin ) ▷ Compute mass flow and enthalpy 14: end if 15: end for 16: Step 3: Compute State Derivatives 17: for each heat exchanger i do 18: Ṁr,i ← ṁin,i − ṁout,i ▷ Mass conservation 19: Ėhx,i ← ṁin,i hin,i − ṁout,i hout,i − Q̇ a,i ▷ Energy conservation 20: end for 21: Assemble derivative vector: ẏ ← [ Ṁr,1 , Ėhx,1 , . . . , Ṁr,N , Ėhx,N ]⊤ 22: Step 4: Advance Solution 23: if solver type is algebraic (Powell hybrid) then 24: yn+1 ← yn + ẏ · ∆t ▷ Explicit Euler step 25: tn+1 ← tn + ∆t 26: else if solver type is IDA or DASSL (DAE) then 27: Solve implicit system: F(tn+1 , yn+1 , ẏn+1 ) = 0 ▷ DAE residual 28: Use BDF method with adaptive step size hn and order k 29: tn+1 ← tn + hn (adaptively determined) 30: end if 31: n ← n + 1, y ← yn+1 , t ← tn+1 32: end while 33: return Solution trajectory {(yn , tn )}nN=0

22

where h(t) is the latent state, x(t) and s(t) are the feature and state vectors, and f θ is the neural ODE function. The error estimate ϵ = ∥z4 − z5 ∥ is used to adaptively adjust the time step: ∆tnew = 2.5



0.6 · ϵ∆t · ∆t0 ϵ

1/4

,

(24)

where ϵ∆t is the time step tolerance (typically 1.1 × 10−3 or 0.3 for stability) and ∆t0 is the current time step. This adaptive mechanism ensures that the time step is automatically reduced when the system dynamics are rapidly changing (e.g., during control input discontinuities) and increased when the dynamics are smooth, balancing accuracy and computational efficiency throughout the simulation. 2.3.3

DAE Solver and Formulation

For a large-scale HVAC system with nc compressors, nv valves, ncond condensers, and nevap evaporators, we formulate the system as an index-1 DAE. The differential variables are the mass charge Mr and internal energy Ehx for each heat exchanger, while the algebraic variables are the junction pressures p j , which are solved separately using a root-finding algorithm. The state vector y ∈ R2(ncond +nevap ) contains only the differential variables: y=

"

(1)

(2)

(1)

(2)

(n

Mr , Mr , . . . , Mr cond (n

Ehx , Ehx , . . . , Ehxcond

# +nevap ) ⊤

+nevap )

,

(25)

where the first ncond + nevap elements correspond to refrigerant mass for each heat exchanger, and the remaining elements correspond to internal energy. The DAE residual function F(t, y, ẏ) = 0 is formulated as: F(t, y, ẏ) = ẏ − f(t, y), (26) where f(t, y) is the right-hand side function that computes the time derivatives of mass and energy for each heat exchanger: Ṁr,i = ṁin,i − ṁout,i ,

Ėhx,i = ṁin,i hin,i − ṁout,i hout,i − Q̇ a,i ,

(27) (28)

where ṁin,i and ṁout,i are the inlet and outlet mass flow rates for heat exchanger i, h denotes specific enthalpy, and Q̇ a,i is the heat transfer rate predicted by the PINODE model for that component. The mass flow rates and enthalpies depend on the current system state and the junction pressures, which are determined by solving the algebraic constraint system: g(p, y, t) = 0,

(29)

where p = [ p1 , p2 , . . . , pn p ]⊤ are the junction pressures. These constraints encode pressure equilibrium at component junctions, mass flow conservation at splitters and mergers, and componentspecific relationships (e.g., compressor and valve characteristics). The pressure system is solved using a Levenberg–Marquardt algorithm (via scipy.optimize.least_squares) at each time step to ensure fast convergence.

23

We implement two DAE solvers to solve the implicit system F(t, y, ẏ) = 0: the IDA solver and the DASSL solver. Both solvers use variable-order, variable-step backward differentiation formula (BDF) methods, but differ in their implementation details and numerical characteristics. At each time step tn , both solvers approximate the derivative using a k-th order BDF formula: ẏn ≈

1 k αi,k yn−i , hn i∑ =0

(30)

where hn = tn − tn−1 is the step size, and αi,k are BDF coefficients for orders k ∈ {1, 2, 3, 4, 5}. Substituting into the DAE residual yields a nonlinear system: ! 1 k G(yn ) = F tn , yn , αi,k yn−i = 0, (31) hn i∑ =0 which both solvers solve via modified Newton iteration. Both solvers adaptively select the step size hn and BDF order k based on local truncation error estimates, and use weighted root-meansquare error norms with absolute and relative tolerances (both set to ATOL = RTOL = ϵsoln , typically 10−6 ) to control accuracy. IDA Solver We implement the IDA (Implicit Differential-Algebraic) solver using the SUNDAE library (specifically sksundae), which provides a Python interface to the SUNDIALS IDA solver. The key distinguishing feature of IDA is its use of standard modified Newton iteration without explicitly forming the step-dependent Jacobian matrix. Additional parameters include maximum and minimum step sizes (hmax , hmin ) and the output interval (∆tIDA out ), which controls the temporal resolution of saved solution points. If the SUNDAE IDA solver is unavailable, the implementation falls back to CVODE (also from SUNDIALS) using the Adams–Moulton method, though this is less suitable for DAE systems. DASSL Solver We employ the DASSL (Differential-Algebraic System Solver) algorithm [17] via the PyDAS Python wrapper, which provides access to the original DASSL Fortran implementation. The main difference from IDA is that DASSL explicitly forms and uses a step-dependent Jacobian matrix in its Newton iteration: Jn =

∂F α0,k ∂F + . ∂y hn ∂ẏ

(32)

This explicit Jacobian formulation can provide better convergence properties for certain stiff systems. The implementation uses adaptive steps with intermediate target times, saving solution points only when a minimum output interval (∆tDASSL out,min ) has elapsed, allowing DASSL to choose its own internal time steps while maintaining control over output density. 2.3.4

Integration with PINODE Framework

The DAE solver integrates seamlessly with the PINODE models for each heat exchanger component through a tightly coupled iterative procedure. At each time step, the following sequence is executed: 1. Initialize state: The DAE solver starts with the current system state yn (mass and energy for each heat exchanger) and determines the time step hn based on adaptive error control. 24

2. Solve pressure system: The junction pressures p are solved using a Levenberg–Marquardt algorithm (via scipy.optimize.least_squares) to satisfy the algebraic constraints: g(p, yn , tn ) = 0.

(33)

This ensures pressure equilibrium at junctions and mass flow conservation throughout the network. 3. Predict heat exchanger outputs: For each heat exchanger, the PINODE model receives the current inputs (inlet conditions from upstream components, boundary pressures from the solved pressure system) and predicts the outlet states (pressures, enthalpies, temperatures) and heat transfer rate Q̇ a . 4. Compute mass and energy derivatives: The predicted Q̇ a and flow properties (mass flow rates ṁin , ṁout , enthalpies hin , hout ) are used to evaluate the right-hand side of the mass and energy conservation equations: Ṁr,i = ṁin,i − ṁout,i ,

Ėhx,i = ṁin,i hin,i − ṁout,i hout,i − Q̇ a,i .

(34) (35)

5. Advance DAE solution: The DAE solver (IDA or DASSL) advances to the next time step using the computed residual function F(t, y, ẏ) = ẏ − f(t, y), where f(t, y) contains the mass and energy derivatives computed in the previous step. This tight coupling ensures that the neural network predictions respect the system-wide algebraic constraints, maintaining thermodynamic consistency across all components. The iterative nature of this process (pressure solving → PINODE prediction → derivative evaluation → DAE advancement) ensures that the solution satisfies both the differential equations and the algebraic constraints simultaneously. 2.3.5

Bayesian Optimization for Parameter Tuning

For large-scale systems, the DAE solver performance depends critically on solver-specific parameters that control time stepping, tolerances, and output resolution. We employ Bayesian optimization [30] to automatically tune these parameters for each solver type. For the algebraic solver (Powell hybrid method), we optimize two parameters: θalg = [ϵ∆t , ϵsoln ], where ϵ∆t is the time step tolerance and ϵsoln is the solution tolerance for the root-finding algorithm. For the DAE-IDA solver, we optimize five parameters: θIDA = [ϵ∆t , ϵsoln , hmax , hmin , ∆tIDA out ], where hmax and hmin control the adaptive step size bounds, and ∆tIDA controls the output interval. out For the DAE-DASSL solver, we optimize six parameters: θDASSL = DASSL DASSL DASSL [ϵ∆t , ϵsoln , hmax , hmin , ∆tout,min , Nmax ], where ∆tout,min controls the minimum spacing between DASSL limits the maximum number of internal DASSL steps. saved solution points and Nmax The optimization objective minimizes a weighted combination of the mean absolute percentage error (MAPE) over all system outputs and the simulation time: θ∗ = arg min [wMAPE · MAPEall (θ) + wtime · tsimulation (θ)] , θ∈X

25

(36)

where wMAPE and wtime are weighting factors (in our implementation, wMAPE = wtime = 0.5), and X is the parameter space defined by log-uniform distributions over specified ranges. Using a Gaussian process surrogate model with expected improvement (EI) acquisition function, the algorithm efficiently explores the parameter space with typically 100 function evaluations, converging to near-optimal solver settings that balance accuracy and computational efficiency. Algorithm 5 Bayesian optimization of parameter optimization n

1: Initialize: Sample n0 points { θi }i=0 1 uniformly, evaluate yi = f (θi ) 2: for n = n0 + 1 to N do

Fit GP on Dn−1 to get µn−1 , σn−1 4: θn = arg maxθ EI(θ) 5: Evaluate yn = f (θn ) 6: Update Dn = Dn−1 ∪ {(θn , yn )} 7: end for 8: return θ∗ = arg min(θ,y)∈D N y

3:

2.3.6

▷ Find next point ▷ Run DAE simulation

Adaptive High-Resolution Stepping

To handle rapid changes in control inputs (e.g., compressor speed or valve opening), we implement adaptive high-resolution time stepping around control discontinuities for both IDA and DASSL solvers. The algorithm pre-computes a high-resolution mask by detecting control jumps where |∆u| > τ (typically τ = 5.0 units) in the actuation signals. For each detected jump, a high-resolution window is created spanning npre = 5 steps before and npost = 50 steps after the jump. The target output spacing is adaptively determined based on the current time index: ( 2.5 s if in high-resolution window, ∆ttarget = 7.5 s otherwise.

(37)

However, the implementation differs between IDA and DASSL solvers due to their different interfaces and adaptive stepping mechanisms. IDA Solver For the IDA solver (implemented via the SUNDAE library, which provides a Python interface to SUNDIALS IDA), the target spacing ∆ttarget directly controls the output interval ∆tIDA out . The solver pre-computes evaluation points teval = {t0 , t0 + ∆ttarget , t0 + 2∆ttarget , . . . , t1 } and passes them to the IDA solver via ida_solver.solve(teval , y0 , y0′ ). The IDA solver then performs adaptive BDF stepping internally between these output points, automatically adjusting its internal time step size hn based on local truncation error estimates while ensuring outputs are provided at the specified evaluation points. DASSL Solver For the DASSL solver (implemented via the PyDAS Python wrapper for the DASSL Fortran library), the adaptive stepping mechanism is more incremental. The solver uses ∆ttarget to determine an intermediate step increment ∆tincrement = 0.5 · ∆ttarget and a minimum output spacing ∆tDASSL out,min (typically 0.1 · ∆ttarget or from the pydas_min_output_dt parameter). During simulation, the solver calls dae_solver.step(tnext ) repeatedly with tnext = min(tcurrent + 26

∆tincrement , t1 ), allowing DASSL to adaptively choose its internal step size. Outputs are saved only when the time since the last saved output exceeds ∆tDASSL out,min , ensuring appropriate temporal resolution while allowing DASSL to take smaller internal steps when needed for stability. Both approaches ensure accurate resolution of transient dynamics during control changes while maintaining computational efficiency during steady-state operation. The implementation also includes output caching to avoid redundant computations when transitioning from high-resolution to low-resolution regions.

3

Results & Discussions

3.1

Training PINODE Models

The physics-informed neural ODE (PINODE) models for individual heat exchangers form the foundation of our system-level simulation framework. We train separate PINODE models for indoor (condenser) and outdoor (evaporator) heat exchangers using the training procedure described in Section 2.1.

Figure 7: Parity plot comparing predicted versus true outputs for the training and testing sets of the condenser (indoor heat exchanger) PINODE model. The model demonstrates excellent agreement with the reference data across both datasets. The PINODE models demonstrate excellent performance on both training and testing datasets, as shown in Figure 7. The parity plots reveal strong linear correlation between predicted and true values, with data points closely aligned along the diagonal, indicating high prediction accuracy. The loss history in Figure 8 indicates stable convergence with minimal overfitting, as evidenced by the parallel trajectories of training and testing losses. The physics-informed loss component Lphys , which penalizes violations of conservation laws through gradient information, and the conservation loss component Lcons , which directly constrains mass and energy predictions, operate on similar data scales, enabling effective multiobjective optimization. For the validation set, the overall loss is typically smaller because the

27

Figure 8: Loss history for the training and testing sets of the condenser (indoor heat exchanger) PINODE model. The stable convergence with minimal gap between training and testing losses indicates good generalization and minimal overfitting. mean squared error contributions from individual components (data, physics, and conservation terms) are balanced, demonstrating that the model successfully learns to satisfy both data fidelity and physical constraints simultaneously. This validates the model architecture and training procedure, ensuring that the well-trained PINODE models serve as reliable building blocks for the system-level simulations discussed in the following sections.

3.2

Corrector Network for System Simulations

While individual PINODE models provide accurate predictions for isolated heat exchangers, system-level simulations involving multiple interconnected components can accumulate errors due to algebraic coupling constraints, numerical integration inaccuracies, and model approximations. To address this, we employ a corrector network (described in Section 2.2) that learns to compensate for systematic discrepancies between predicted and true mass-energy (ME) terms at the system level. We begin with manually inspected, physically justifiable parameters (e.g., ϵ∆t = 10−3 ) for the system solver. The corrector network is trained on a short initial segment of the simulation trajectory (850 time steps) and then deployed throughout the entire simulation to improve prediction accuracy. This lightweight neural network takes normalized state and output predictions from both condensers as input and outputs correction terms for the four ME components (E1 , E2 , M1 , M2 ). The training loss of the corrector network is shown in Figure 9. The network converges rapidly, with training completing in approximately 30 seconds, demonstrating the computational efficiency of this approach. The monotonic decrease in loss indicates stable learning, and the final low loss value suggests that the network successfully captures the systematic biases present in the system-level predictions. The effectiveness of the corrector network is immediately apparent when comparing system-level predictions. Figure 10 shows a significant improvement in the parity plot after applying the corrector network, with predictions much closer to the reference Dymola simulation results. The 28

Figure 9: Training loss history of the corrector neural network. The network converges rapidly, demonstrating efficient learning of systematic biases in mass-energy predictions. baseline predictions (without corrector) show noticeable scatter and deviation from the line of perfect prediction, particularly for extreme values, while the corrected predictions exhibit much tighter clustering around the ideal one-to-one correspondence, indicating substantially improved accuracy. The time-stepping behavior of the dual-compressor cycle is illustrated in Figure 11, highlighting the adaptive time step selection throughout the simulation. The observed step size variation is closely tied to the system dynamics: the solver automatically reduces the time step during periods of rapid transients (e.g., control input changes and strong component interactions), while increasing it during quasi-steady regimes where the solution evolves more smoothly. This behavior reflects effective local error control, allowing the solver to allocate computational effort where it is most needed without uniformly refining the entire trajectory. Importantly, the presence of the corrector network does not introduce instability or excessive stiffness; instead, it preserves stable integration while enabling accurate tracking of system dynamics with a balanced computational cost. Overall, the adaptive stepping pattern provides indirect evidence that the combined PINODE–corrector framework produces numerically well-behaved trajectories suitable for efficient time integration. To further quantify the improvement, we examine the normalized mass-energy (ME) terms directly. Figure 12 shows the evolution of the four ME components (E1 , E2 , M1 , M2 ) before and after applying the corrector network. In the baseline case (Figure 12(a)), the predicted trajectories exhibit systematic deviations from the reference solution, most notably during transient periods where coupling effects and integration errors are more pronounced. After applying the corrector network (Figure 12(b)), these discrepancies are significantly reduced, with the predicted trajectories more closely aligned with the reference across all components. The improvement is particularly evident in the reduction of transient overshoots and bias, indicating that the corrector effectively compensates for accumulated system-level errors arising from component coupling and numerical integration. While Figure 12 focuses on component-level mass-energy variables, we further examine how these improvements propagate to the system-level energy evolution.

29

Before applying corrector NN

After applying corrector NN

(a)

(b)

Figure 10: Comparison of parity plots (a) before and (b) after applying the corrector neural network. The corrector network significantly improves the agreement between PINODE predictions and the reference Dymola simulation results, with predictions clustering much closer to the diagonal.

Figure 11: Adaptive time step evolution during the dual-compressor cycle simulation with the corrector network. The solver adjusts the time step size based on local error estimates, using smaller steps during transients and larger steps during steady operation.

30

Before applying corrector NN

(a)

After applying corrector NN

(b)

Figure 12: Comparison of normalized mass-energy (ME) terms for the dual-compressor system: (a) before and (b) after applying the corrector network. The corrector network reduces systematic deviations between predicted and reference ME values, with improvements observed across both energy terms (E1 , E2 ) and mass terms (M1 , M2 ).

Before applying corrector NN

(a)

After applying corrector NN

(b)

Figure 13: Refrigerant energy evolution in the dual-compressor HVAC system: (a) without and (b) with the corrector network. Improvements in local mass-energy variables translate to more consistent system-level energy predictions, with reduced discrepancies relative to the reference solution, particularly during transient regimes.

31

Figure 13 shows the refrigerant energy evolution of the overall system. In the baseline case (Figure 13(a)), noticeable deviations from the reference persist, particularly during transient regimes. After applying the corrector network (Figure 13(b)), the predicted energy trajectory more closely tracks the reference across the entire time horizon, with reduced transient discrepancies. This behavior is consistent with the correction of local mass-energy states leading to improved global energy consistency.

3.3

Generalizability in Dual-Compressor Systems

3.3.1

Bayesian Optimization

To systematically identify optimal solver parameters, we apply Bayesian optimization to characterize the design space for solver tolerances across three different solver types: the algebraic (Powell hybrid) solver, the DAE-IDA solver, and the DAE-DASSL solver. The objective function balances simulation accuracy (measured by overall MAPE) and computational efficiency (simulation time), allowing us to find Pareto-optimal configurations that trade off these competing objectives. For the algebraic solver, we optimize two key parameters: the time step tolerance (ϵ∆t ) and the solution tolerance (ϵsoln ). Figure 14 shows the evolution of the objective function during the optimization process, demonstrating convergence to improved solutions over 100 iterations. The objective function decreases monotonically, indicating that the Bayesian optimization successfully explores the design space and identifies progressively better parameter combinations. The corresponding design space contours in Figure 15 reveal the sensitivity of both MAPE and simulation time to different tolerance combinations, identifying regions of optimal performance. The contour plots show that tighter tolerances generally improve accuracy but increase computational cost, creating a clear trade-off that the optimization must navigate.

Figure 14: Evolution of the objective function during Bayesian optimization for the algebraic solver. The objective function decreases over 100 iterations, demonstrating convergence to improved parameter configurations that balance accuracy and computational efficiency.

32

Figure 15: Design space characterization in Bayesian optimization for the algebraic solver. The contour plots show the sensitivity of (left) overall MAPE and (right) simulation time to different combinations of ϵ∆t and ϵsoln . The optimal configuration (marked with a yellow star) achieves a balance between accuracy and efficiency. For the two DAE solvers, we conduct a high-dimensional design space search since these solvers inherently require more parameters than the algebraic solver. The optimization results are shown in Figures 16 and 17, revealing distinct optimal regions compared to the algebraic solver. Unlike the algebraic solver, which optimizes over a 2-dimensional parameter space (ϵ∆t , ϵsoln ), the DAE-IDA solver explores a 5-dimensional design space (ϵ∆t , ϵsoln , hmax , hmin , ∆tIDA out ), as described in Section 2.3. The objective function evolution (Figure 16) demonstrates superior convergence behavior compared to the algebraic solver, with the Bayesian optimization algorithm efficiently navigating the higher-dimensional space to identify optimal configurations that balance accuracy and computational efficiency. The 2D design space projections (shown for the two primary tolerances) reveal different sensitivity patterns, indicating that the DAE-IDA solver has distinct optimal operating characteristics that require careful multi-dimensional tuning. Similarly, the DAE-DASSL solver optimization involves six parameters: ϵ∆t , ϵsoln , hmax , hmin , DASSL ∆tDASSL out,min , and Nmax . The optimization results are presented in Figures 18 and 19. The design space exploration reveals different optimal parameter regions for DASSL compared to both the algebraic solver and the IDA solver, reflecting the unique numerical properties and adaptive timestepping strategies of each solver. To compare the overall performance across all three solvers, we examine the Pareto fronts in Figure 20. The Pareto front represents the set of non-dominated solutions, where no solution can improve one objective without worsening the other. From the Pareto front plots, we observe that the DAE-IDA solver generally outperforms the other solvers, achieving better trade-offs between accuracy and computational efficiency across the design space. The DAE-IDA solver’s Pareto front lies closer to the origin (lower MAPE and lower simulation time), indicating superior overall performance. The algebraic solver shows competitive performance but with a less favorable Pareto front, while the DAE-DASSL solver exhibits a wider spread of solutions, suggesting more variability in its performance characteristics. 33

Figure 16: Evolution of the objective function during Bayesian optimization for the DAE-IDA solver. The optimization explores a 5-dimensional parameter space to find optimal configurations balancing accuracy and computational efficiency.

Figure 17: Design space characterization in Bayesian optimization for the DAE-IDA solver. The contour plots show the sensitivity of (left) overall MAPE and (right) simulation time to different combinations of ϵ∆t and ϵsoln , with other parameters fixed at their optimal values.

34

Figure 18: Evolution of the objective function during Bayesian optimization for the DAE-DASSL solver. The optimization explores a 6-dimensional parameter space to identify optimal configurations.

Figure 19: Design space characterization in Bayesian optimization for the DAE-DASSL solver. The contour plots show the sensitivity of (left) overall MAPE and (right) simulation time to different combinations of ϵ∆t and ϵsoln , with other parameters fixed at their optimal values.

35

Figure 20: Comparison of Pareto fronts for the three solvers on a log-log scale. The DAE-IDA solver achieves the best overall trade-off between accuracy (MAPE) and computational efficiency (simulation time), with its Pareto front lying closest to the origin. 3.3.2

Simulation Results

Using the optimal parameters identified through Bayesian optimization, we evaluate the performance of each solver on the dual-compressor system. Figure 21 compares the overall prediction accuracy for all three solvers, demonstrating that both DAE solvers achieve excellent agreement with the reference Dymola simulation, with the DAE-IDA solver showing particularly strong performance. All three solvers produce predictions that cluster closely around the diagonal, indicating high accuracy, but the DAE-IDA solver shows the tightest clustering and best overall agreement with the reference data. The quantitative performance metrics are summarized in Tables 1 and 2. Table 1 shows the best overall objective (weighted combination of MAPE and simulation time), where the DAE-IDA solver achieves the lowest MAPE (2.04%) and fastest simulation time (58.91 s) simultaneously, representing the best balance between accuracy and efficiency. The algebraic solver achieves a slightly higher MAPE (2.20%) with a longer simulation time (73.33 s), while the DAE-DASSL solver shows comparable accuracy (2.19% MAPE) but slightly better efficiency (62.30 s) than the algebraic solver. Table 2 focuses on the fastest simulation times, where the DAE-IDA solver again excels with a simulation time of 15.43 s, though with a slightly higher MAPE of 3.56%. This demonstrates the DAE-IDA solver’s ability to achieve both high accuracy and high computational efficiency, depending on the optimization objective. Table 1: Best Overall Objective (Weighted Combination of MAPE and Simulation Time)

MAPEall [%] tsimulation [s]

Algebraic

DAE-IDA

DAE-DASSL

2.2045% 73.33

2.0366% 58.91

2.1888% 62.30

36

Algebraic solver

(a) DAE-IDA solver

(b)

DAE-DASSL solver

(c)

Figure 21: Comparison of parity plots for (a) algebraic solver, (b) DAE-IDA solver, and (c) DAEDASSL solver using optimal parameters identified through Bayesian optimization. All three solvers demonstrate excellent agreement with the reference Dymola simulation, with the DAEIDA solver showing the tightest clustering around the diagonal.

Table 2: Best Simulation Time (Minimum Time with Corresponding MAPE)

MAPEall [%] tsimulation [s]

Algebraic

DAE-IDA

DAE-DASSL

2.2486% 29.14

3.5551% 15.43

7.3143% 35.93

37

The computational speedups achieved by our framework are substantial when compared to highfidelity simulators. The DAE-IDA solver achieves a simulation time of 15.43 s, representing approximately a 9.5× speedup compared to the Dymola-based simulator reported in [10], which requires 146.8 s for similar HVAC system simulations. The algebraic solver achieves a 5.0× speedup (29.14 s), while the DAE-DASSL solver achieves a 4.1× speedup (35.93 s). These 4–9× computational speedups, combined with MAPE errors below 2.5%, demonstrate that our hybrid framework successfully bridges the accuracy–efficiency trade-off that has limited the practical deployment of physics-based HVAC simulators. The optimal time step sizes selected by each solver reflect their different numerical strategies and adaptive time-stepping mechanisms. Figure 22 compares the optimal ∆t evolution for the three solvers, revealing distinct temporal resolution patterns. The DAE-IDA solver achieves smaller time steps on average, enabling higher temporal resolution and better capture of transient dynamics, which contributes to its superior accuracy. The algebraic solver operates with larger time steps for computational efficiency, while the DAE-DASSL solver shows intermediate behavior, balancing between the two extremes. The time step evolution patterns reflect each solver’s internal error control mechanisms and their ability to adapt to local solution characteristics, as discussed in Section 2.3. For the algebraic solver, the RK45 adaptive stepping mechanism (see Section 2.3) automatically adjusts the time step based on local truncation error estimates, reducing ∆t when dynamics are rapidly changing and increasing it during smooth transients. The DAE solvers (IDA and DASSL) employ variable-order, variable-step BDF methods with adaptive error control, using weighted root-mean-square error norms to adaptively select both step size hn and BDF order k, enabling fine-grained temporal resolution control that is evident in the time step evolution patterns. The observed behavior aligns with the theoretical framework for adaptive time stepping, where local error estimates drive automatic step size selection to maintain accuracy while maximizing computational efficiency. Finally, Figure 23 shows the refrigerant energy predictions using the optimal solver parameters for all three solvers. All three solvers demonstrate good agreement with the reference Dymola simulation, validating the effectiveness of the Bayesian optimization approach in identifying highperforming parameter configurations. The DAE-IDA solver consistently shows the best overall performance, achieving high accuracy while maintaining computational efficiency. The energy predictions track the reference trajectory closely across all time periods, including rapid transients and steady-state operation, demonstrating the robustness of the optimized solver configurations. Based on the dual-compressor simulations, the DAE-IDA solver emerges as a strong candidate for general-purpose HVAC dynamical simulations for several reasons. First, it consistently achieves favorable computational efficiency, often matching or exceeding the fastest configurations identified for the other solvers. Second, it maintains low prediction error across a wide range of operating conditions, with MAPE values comparable to or slightly better than those of the algebraic and DAE-DASSL solvers. Third, and most importantly, the DAE-IDA solver demonstrates superior performance on the Pareto front, indicating a more robust trade-off between accuracy and computational cost across the explored parameter space. This suggests that, beyond individual optimal points, the IDA solver provides more reliable performance under varying tolerance and stepping configurations. However, it is important to note that for smaller-scale systems, such as the dual-compressor case considered here, the performance gap between solvers is relatively modest. In particular, the 38

Figure 22: Comparison of optimal time step sizes (∆t) for the three solvers obtained through Bayesian optimization. The DAE-IDA solver achieves the smallest time steps, enabling higher temporal resolution, while the algebraic solver operates with larger time steps for computational efficiency. algebraic solver remains competitive in terms of both runtime and accuracy when appropriately tuned, and in some cases can achieve comparable best-case performance. Therefore, while DAEIDA offers better overall robustness and scalability, the choice of solver for small systems may still depend on specific optimization objectives, implementation simplicity, and parameter tuning considerations.

39

Algebraic

(a) DAE-IDA

(b)

DAE-DASSL

(c)

Figure 23: Refrigerant energy predictions using optimal solver parameters for (a) algebraic solver, (b) DAE-IDA solver, and (c) DAE-DASSL solver. All three solvers demonstrate good agreement with the reference Dymola simulation, with the DAE-IDA solver showing the best overall performance in terms of accuracy and computational efficiency.

40

4

Scalability of HVAC Systems

4.1

Scaling Up HVAC

We scale up the HVAC systems by simultaneously increasing the number of compressor– condenser pairs (nc ) and valve–evaporator pairs (nv ), demonstrating the scalability of our algorithms. This represents the first demonstration of such system scaling in the literature. Compared with the previous work by Ma et al. [10], we show that with our new system solver and effective data-driven components, large-scale systems can be solved efficiently, which was indeed a significant challenge for traditional physics-based solvers. The large-scale system topology follows a parallel–merge architecture that generalizes the dualcompressor configuration: 1. Parallel compression stage: Each of the nc compressors discharges into its own dedicated condenser, operating in parallel. 2. Liquid manifold: All condenser outlets merge into a single liquid manifold node, creating a common high-pressure liquid reservoir. 3. Parallel expansion stage: The liquid manifold feeds nv parallel valve–evaporator pairs, each operating independently. 4. Suction manifold: All evaporator outlets merge into a single suction manifold node, creating a common low-pressure vapor reservoir. 5. Compressor inlets: The suction manifold feeds all nc compressor inlets, completing the cycle. This topology ensures that each compressor–condenser pair operates independently while sharing common liquid and suction manifolds, enabling efficient parallel operation and load distribution across multiple compressors. The system solver automatically adapts to the increased number of junction pressures, where n p scales with nc and nv according to the system topology. For systems with n p > 10, the algebraic solver automatically switches from Powell hybrid method to bounded least-squares optimization (Levenberg–Marquardt) to ensure robust convergence. For DAE solvers, the state vector dimension scales linearly with the number of heat exchangers: y ∈ R2(ncond +nevap ) , where ncond = nc and nevap = nv for the parallel–merge architecture, and each heat exchanger contributes two differential variables (mass charge Mr and internal energy Ehx ).

4.2

Stability and Solver Parameters

To evaluate the stability and performance of our solvers on large-scale systems, we conduct simulations for systems ranging from nc = 2 to nc = 32 compressor–condenser pairs, with an equal number of valve–evaporator pairs (nv = nc ). Each simulation runs for Nsteps = 500 time steps to assess both short-term stability and computational efficiency. The solvers are configured with optimized parameters identified through Bayesian optimization (Section 2.3), ensuring fair comparison across different system sizes. The optimization parameters include tolerances (ϵ∆t , ϵsoln ), step DASSL DASSL for DASSL). size bounds (hmax , hmin ), and output intervals (∆tIDA out for IDA, ∆tout,min and Nmax Figure 25 compares the time step evolution for a system with nc = 2 compressors and nv = 2 evaporators across all three solver types. The DAE-IDA solver achieves the finest temporal 41

Algorithm 6 Large-scale HVAC system construction Require: Number of compressor–condenser pairs nc , number of valve–evaporator pairs nv Ensure: Configured system with nc compressors and nv evaporators 1: Initialize: Create nc compressors C = {C1 , . . . , Cnc }, nc condensers Hc = { Hc,1 , . . . , Hc,nc }, nv valves V = {V1 , . . . , Vnv }, nv evaporators He = { He,1 , . . . , He,nv } 2: Connect parallel compression: 3: for k = 1 to nc do 4: Connect Ck .out → Hc,k .in ▷ Compressor discharge to condenser inlet 5: end for 6: Create liquid manifold: 7: Set pliq = Hc,1 .out ▷ Liquid manifold pressure node 8: for k = 2 to nc do 9: Connect Hc,k .out → pliq ▷ Merge condenser outlets 10: end for 11: Connect parallel expansion: 12: for k = 1 to nv do 13: Connect pliq → Vk .in ▷ Liquid manifold to valve inlet 14: Connect Vk .out → He,k .in ▷ Valve outlet to evaporator inlet 15: end for 16: Create suction manifold: 17: Set psuct = He,1 .out ▷ Suction manifold pressure node 18: for k = 2 to nv do 19: Connect He,k .out → psuct ▷ Merge evaporator outlets 20: end for 21: Connect compressor inlets: 22: for k = 1 to nc do 23: Connect psuct → Ck .in ▷ Suction manifold to compressor inlet 24: end for 25: return Configured system with topology G(nc , nv )

42

…

… …

… …

…

Figure 24: Scaling topology for large-scale HVAC systems. The system scales from a dualcompressor configuration (left) to configurations with nc compressors and nv evaporators (right). The parallel–merge architecture ensures independent operation of compressor–condenser pairs while maintaining thermodynamic consistency through shared liquid and suction manifolds. resolution with smaller, more adaptive time steps (∆t ≈ 0.1–1.0 s), enabling better capture of transient dynamics. The DAE-DASSL solver shows intermediate behavior with time steps in the range ∆t ≈ 0.5–2.0 s, while the algebraic solver operates with larger, more uniform time steps (∆t ≈ 1.0–5.0 s) for computational efficiency. All three solvers successfully complete the Nsteps = 500 simulation, demonstrating stability for this system size (nc = nv = 2). Figure 26 shows the computational scaling behavior as the number of components increases from nc = 2 to nc = 32. Across all solvers, the simulation time tsimulation exhibits clear superlinear growth with system size. Over the range considered, the increase in runtime is broadly consistent with polynomial scaling, with an effective exponent between quadratic and cubic, although the limited number of data points precludes a precise characterization. This trend is expected given the increasing complexity of the coupled algebraic pressure system with n p = O(nc ) junction pressures, as well as the growth of the state space y ∈ R2(nc +nv ) for the DAE solvers. Among the methods, the DAE-IDA solver demonstrates comparatively favorable scaling behavior, maintaining moderate simulation times even for larger systems (e.g., tsimulation ≲ 200 s at nc = 32). The algebraic solver performs competitively for small to medium system sizes (nc ≤ 8), but its runtime increases more rapidly at larger scales, indicating reduced scalability. The DAE-DASSL solver exhibits intermediate behavior, with consistently higher computational cost than DAE-IDA but more stable scaling than the algebraic formulation across the tested range. The time step evolution patterns observed in Figure 25 directly reflect the adaptive time-stepping mechanisms described in Section 2.3: the algebraic solver’s RK45 method adaptively adjusts ∆t based on local error estimates, while the DAE solvers’ BDF methods use adaptive error control to

43

Figure 25: Time step (∆t) evolution comparison for a large-scale system with nc = 2 compressors and nv = 2 evaporators. The DAE-IDA solver achieves the finest temporal resolution with smaller, adaptive time steps, while the algebraic solver uses larger, more uniform steps. All three solvers successfully complete the Nsteps = 500 simulation, demonstrating stability for medium-scale systems. select both step size and order, resulting in the distinct temporal resolution behaviors observed across system sizes. The results demonstrate that our framework successfully scales to large HVAC systems, with the DAE-IDA solver providing the best balance between accuracy and computational efficiency across all system sizes tested. This scalability is crucial for practical applications where HVAC systems may contain dozens of compressors and heat exchangers, enabling system-level optimization and control design for real-world installations.

5

Conclusions

This work presents a comprehensive framework for modeling large-scale HVAC systems using physics-informed neural ODEs (PINODEs) integrated with differential-algebraic equation (DAE) solvers. The framework combines data-driven component models with physics-based constraints, enabling accurate and efficient simulation of complex multi-compressor systems. Our key contributions include: (1) physics-informed neural ODE models for individual heat exchangers that capture transient dynamics while respecting thermodynamic principles, (2) a corrector neural network that compensates for system-level errors by learning from short training segments, (3) integration of DAE solvers (IDA and DASSL) to handle algebraic constraints arising from pressure equilibrium and mass flow conservation, and (4) Bayesian optimization for automatic parameter tuning across multiple solver types. The results demonstrate that the DAE-IDA solver achieves superior performance across all metrics, achieving an overall MAPE of 2.04% with a simulation time of 58.91 seconds for the dualcompressor system, outperforming both the algebraic solver and DAE-DASSL solver. The correc44

Figure 26: Computational scaling comparison for HVAC systems with increasing numbers of compressor–condenser pairs (nc ranging from 2 to 32). The simulation time tsimulation exhibits clear superlinear growth with system size. Over the tested range, the increase in runtime is broadly consistent with polynomial scaling, with an effective exponent between quadratic and cubic, although the limited number of data points precludes a precise characterization. This trend reflects the increasing complexity of the coupled algebraic pressure system with n p = O(nc ) junction pressures, as well as the growth of the state space y ∈ R2(nc +nv ) for DAE solvers. Among the methods, the DAE-IDA solver demonstrates comparatively favorable scaling behavior, maintaining moderate simulation times even for larger systems (e.g., nc = 32). Failed simulations are marked with red crosses.

45

tor network significantly improves prediction accuracy, reducing systematic biases in mass and energy predictions and bringing the system-level results into close agreement with reference Dymola simulations. Bayesian optimization successfully identifies optimal solver parameters, revealing distinct optimal regions for each solver type and enabling automated tuning that balances accuracy and computational efficiency. The framework demonstrates excellent scalability, successfully simulating systems with up to 32 compressor–condenser pairs while maintaining computational efficiency. The DAE-IDA solver exhibits superior scaling behavior, maintaining reasonable simulation times even for large systems, while the algebraic solver shows steeper computational scaling due to the increasing complexity of the pressure system. This scalability is crucial for practical applications, as real-world HVAC installations often contain dozens of compressors and heat exchangers, requiring efficient system-level simulation for optimization and control design. The combination of physics-informed neural networks, system-level correction, and robust DAE solvers provides a powerful approach for modeling complex HVAC systems that was previously challenging with traditional physics-based methods. This work opens new possibilities for system-level optimization, predictive control, and design exploration of large-scale HVAC installations, with potential applications in building energy management, fault detection, and real-time control systems.

Declaration of Competing Interests None.

CRediT authorship contribution statement Hanfeng Zhai: Conceptualization, Data curation, Formal analysis, Funding acquisition, Investigation, Methodology, Software, Validation, Visualization, Writing - original draft, Writing - review & editing. Hassan Mansour: Conceptualization, Data curation, Formal analysis, Funding acquisition, Investigation, Methodology, Software, Validation, Visualization, Writing - original draft, Writing review & editing. Hongtao Qiao: Conceptualization, Data curation, Formal analysis, Funding acquisition, Investigation, Methodology, Software, Validation, Visualization, Writing - original draft, Writing - review & editing. Christopher Laughman: Conceptualization, Data curation, Formal analysis, Funding acquisition, Investigation, Methodology, Software, Validation, Visualization, Writing - original draft, Writing review & editing.

Data availability The data that has been used is confidential.

46

Appendix This appendix provides additional results for the outdoor heat exchanger (evaporator) PINODE model, complementing the indoor heat exchanger results presented in Section 3.1. The outdoor heat exchanger model follows the same training procedure and architecture as the indoor model, ensuring consistency across both components used in the system-level simulations.

Figure 27: Parity plot comparing predicted versus true outputs for the training and testing sets of the evaporator (outdoor heat exchanger) PINODE model. The model demonstrates excellent agreement with the reference data across both datasets, with data points closely aligned along the diagonal, indicating high prediction accuracy for the outdoor heat exchanger component. Figure 27 shows the parity plot for the outdoor heat exchanger model, demonstrating excellent performance on both training and testing datasets. The strong linear correlation and tight clustering around the diagonal indicate that the PINODE model successfully captures the dynamics of the evaporator component, providing accurate predictions for use in system-level simulations. The loss history in Figure 28 confirms stable convergence with minimal overfitting, as evidenced by the parallel trajectories of training and testing losses throughout the training process. This validates the model architecture and training procedure for the outdoor heat exchanger, ensuring that both indoor and outdoor components are well-trained and ready for integration into the system-level simulation framework.

References [1] U.S. Energy Information Administration. Annual energy outlook 2021. Technical report, U.S. Energy Information Administration, 2021. [2] Abdul Afram and Farrokh Janabi-Sharifi. Review of modeling methods for hvac systems. Applied Thermal Engineering, 67(1–2):507–519, June 2014.

47

Figure 28: Loss history for the training and testing sets of the evaporator (outdoor heat exchanger) PINODE model. The stable convergence with minimal gap between training and testing losses indicates good generalization and minimal overfitting, validating the model architecture and training procedure for the outdoor heat exchanger component. [3] Ankush Chakrabarty, Emilio Maddalena, Hongtao Qiao, and Christopher Laughman. Scalable bayesian optimization for model calibration: Case study on coupled building and hvac dynamics. Energy and Buildings, 253:111460, December 2021. [4] Modelica Association. Modelica–a unified object-oriented language for physical systems modeling: Language specification version 3.0. Technical report, Modelica Association, 2005. [5] Drury B Crawley, Linda K Lawrie, Frederick C Winkelmann, W F Buhl, Y Joe Huang, Curtis O Pedersen, Richard K Strand, Richard J Liesen, Daniel E Fisher, Michael J Witte, et al. Energyplus: creating a new-generation building energy simulation program. Energy and buildings, 33(4):319–331, 2001. [6] Hongtao Qiao, Vikrant Aute, and Reinhard Radermacher. Transient modeling of a flash tank vapor injection heat pump system–part i: Model development. International journal of refrigeration, 49:169–182, 2015. [7] Yixing Chen, Tianzhen Hong, and Mary Ann Piette. A review of data-driven building energy consumption prediction studies. Renewable and Sustainable Energy Reviews, 81:1192–1205, 2018. [8] Per Sahlin, Lars Eriksson, Pavel Grozman, Hans Johnsson, Alexander Shapovalov, and Mika Vuolle. Whole-building simulation with symbolic dae equations and general purpose solvers. Building and Environment, 39(8):949–958, August 2004. [9] Ranjan Anantharaman. Approximation of Large Stiff Acausal Models. PhD thesis, Massachusetts Institute of Technology, Cambridge, MA, 2023. [10] Jiacheng Ma, Yiyun Dong, Hongtao Qiao, and Christopher R. Laughman. A physicsconstrained deep learning framework for dynamic modeling of vapor compression systems. Applied Thermal Engineering, 254:123734, October 2024.

48

[11] Ankush Chakrabarty, Emilio Maddalena, Hongtao Qiao, and Christopher Laughman. Datadriven calibration of joint building and hvac dynamic models using scalable bayesian optimization. In Proceedings of Building Simulation 2021: 17th Conference of IBPSA, BS 2021. KU Leuven, September 2021. [12] Dafang Zhao, Daichi Watari, Yuki Ozawa, Ittetsu Taniguchi, Toshihiro Suzuki, Yoshiyuki Shimoda, and Takao Onoye. Data-driven online energy management framework for hvac systems: An experimental study. Applied Energy, 352:121921, December 2023. [13] Vincent Taboga, Clement Gehring, Mathieu Le Cam, Hanane Dagdougui, and Pierre-Luc Bacon. Neural differential equations for temperature control in buildings under demand response programs. Applied Energy, 368:123433, August 2024. [14] Zhelun Chen, Jin Wen, Anthony Kearsley, and Amanda Pertzborn. Smoothing techniques in dynamic building system simulation. In 2021 International Conference on Instrumentation, Control, and Automation (ICA), page 156–161. IEEE, August 2021. [15] Aleksei Sholokhov, Yuying Liu, Hassan Mansour, and Saleh Nabi. Physics-informed neural ode (pinode): embedding physics into models using collocation points. Scientific Reports, 13(1), June 2023. [16] Ida – sundials: An open-source software library of time integrators and nonlinear solvers. https://computing.llnl.gov/projects/sundials/ida, 2025. Accessed: 2025-12-10. [17] Linda R Petzold. A description of dassl: A differential/algebraic system solver. Technical Report SAND82-8637, Sandia National Labs., 1982. [18] Evan Fleming, Shaoyi Wen, Li Shi, and Alexandre K. da Silva. Thermodynamic model of a thermal storage air conditioning system with dynamic behavior. Applied Energy, 112:160–169, December 2013. [19] A. F. I. Mamadou, M. O. K. Idrissou, S. A. O. Sanya, J. V. C. Vargas, and A. R. Alexandra. Intelligent optimal control of endoreversible single-effect hvac-ar system using machine learning. International Journal of Air-Conditioning and Refrigeration, 33(1), October 2025. [20] Maziar Raissi, Paris Perdikaris, and George Em Karniadakis. Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. Journal of Computational Physics, 378:686–707, 2019. [21] Ricky TQ Chen, Yulia Rubanova, Jesse Bettencourt, and David K Duvenaud. Neural ordinary differential equations. Advances in neural information processing systems, 31, 2018. [22] Edward De Brouwer, Jaak Simm, Adam Arany, and Yves Moreau. Gru-ode-bayes: Continuous modeling of sporadically-observed time series. Advances in neural information processing systems, 32, 2019. [23] Tobias Thummerer and Lars Mikelsons. Eigen-informed neuralodes: Dealing with stability and convergence issues of neuralodes. arXiv preprint arXiv:2302.10892, 2023. [24] Xuhui Meng, Zhen Li, Dongkun Zhang, and George Em Karniadakis. Ppinn: Parareal physics-informed neural network for time-dependent pdes. Computer Methods in Applied Mechanics and Engineering, 370:113250, 2020. [25] Hanfeng Zhai, Quan Zhou, and Guohui Hu. Predicting micro-bubble dynamics with semiphysics-informed deep learning. AIP Advances, 12(3), 2022. 49

[26] Lu Lu, Xuhui Meng, Zhiping Mao, and George Em Karniadakis. Deepxde: A deep learning library for solving differential equations. SIAM review, 63(1):208–228, 2021. [27] Peter I Frazier. Bayesian optimization. In Recent advances in optimization and modeling of contemporary problems, pages 255–278. Informs, 2018. [28] Hanfeng Zhai and Jingjie Yeo. Computational design of antimicrobial active surfaces via automated bayesian optimization. ACS Biomaterials Science & Engineering, 9(1):269–279, 2022. [29] Jimmy Lei Ba, Jamie Ryan Kiros, and Geoffrey E Hinton. Layer normalization. arXiv preprint arXiv:1607.06450, 2016. [30] Jonas Mockus, Vytautas Tiesis, and Antanas Zilinskas. On bayesian methods for seeking the extremum. Optimization techniques IFIP technical conference, pages 400–404, 1975.

50