In the analysis of constrained holonomic systems, the Lagange multiplier method yields a system of second-order ordinary differential equations of motion and algebraic constraint equations. Conventional holonomic or nonholonomic constraints are defined as geometric constraints in this paper. Previous works concentrate on the geometric constraints. However, if the total energy of a dynamic system can be computed from the initial energy plus the time integral of the energy input rate due to external or internal forces, then the total energy can be artificially treated as a constraint. The violation of the total energy constraint due to numerical errors can be used as information to control these errors. It is a necessary condition for accurate simulation that both geometric and energy constraints be satisfied. When geometric constraint control is combined with energy constraint control, numerical simulation of a constrained dynamic system becomes more accurate. A new convenient and effective method to implement energy constraint control in numerical simulation is developed based on the geometric interpretation of the relation between constraints in the phase space. Several combinations of energy constraint control with either Baumgarte's Constraint Violation Stabilization Method (CVSM) are also addressed.