We propose a method for reducing the order of dynamic models of temperature and humidity in multizone buildings. Low-order models of building thermal dynamics are useful for model-based HVAC control techniques, especially to computationally intensive ones such as Model Predictive Control (MPC). Even a lumped parameter model for a multi-zone building, which is a set of non-linear coupled ordinary differential equations, can have large state-space dimension. Model reduction techniques are useful to simplify such models. Although there are a number of well-developed techniques for model reduction of linear systems, techniques available for non-linear systems are limited. The method we propose exploits the linear portion of the model to compute a transformation (by using balanced realization) and a specific sparsity pattern of the non-linear portion that building thermal models possess to obtain the reduced order model. Simulations show that the prediction of the zone temperatures and humidity ratios by the reduced model is quite close to that from the full-scale model, even when substantial reduction of model order is specified that reduces computation time by a factor of six or more. (C) 2011 Elsevier B.V. All rights reserved.