Energy conservation and indoor environment concerns have motivated extensive research on various aspects of control of Heating, Ventilating and Air-Conditioning (HVAC) and building systems. The study on optimal operation as well as modeling of HVAC and building systems is one of the fastest growing fields that contribute to saving energy and improving indoor environment.  This thesis is devoted to the development of a comprehensive modeling and optimization methodology for global multiple-stage optimal operation of HVAC and building systems. Two different dynamic models of a multizone variable air volume (VAV) system have been developed using (i) bottom-up and (ii) top-down approaches. The models take account of the dynamic interactions between building shell, VAV system components and control systems. The models describe the dynamics of fan, air distribution system, zone(s), cooling coil and primary plant (chiller) as one multivariable nonlinear system in a way that is useful for control analysis. Using the bottom-up approach a large-scale VAV system model has been developed. This model considers the interactions between flow field and thermal field via distributed capacity and variable air density considerations. An alternate model which is computationally more efficient was developed using the top-down approach. Model reduction techniques were applied to develop a reduced-order state space model of the VAV system. Results show that predictions from the reduced order model are within 5% of those from the large scale model.  Optimal control schemes are developed for the efficient operation of VAV systems. In the control scheme proposed it is necessary to compute optimal setpoint profiles for local controllers. The optimal control profiles so computed can be used as tracking signals for local controllers for moving the system states from one setpoint to another. In order to determine optimal setpoint profiles an optimization methodology for formulating and solving the multiple stage optimal operation problems has been developed. The methodology is based on the maximum principle of Pontryagin and perturbation method in order to deal with the multiple time-scale of the HVAC processes and building operating schedules. A solution methodology and the corresponding computer models have been developed for solving the multiple stage optimal operation problems.  The applications of the VAV model and the multistage optimization methodology have been demonstrated by considering several practical examples. The examples include (i) a comparison of optimal strategies for constant and variable air volume systems with and without time-of-day price structure for electrical energy, (ii) a two-zone VAV heating system and (iii) a five-zone VAV cooling system. Results showing the 24-hour optimal setpoint profiles, energy cost savings and the output responses such as zone temperatures and humidity ratios are given for different building operation schedules. These applications show that the developed models and optimization methodology can be used to determine energy efficient operating strategies for VAV systems without violating the thermal comfort in buildings.