Peak power consumption is a universal problem across energy control systems in electrical grids, buildings, and industrial automation where the uncoordinated operation of multiple controllers result in temporally correlated electricity demand surges (or peaks). While there exist several different approaches to balance power consumption by load shifting and load shedding, they operate on coarse grained time scales and do not help in de-correlating energy sinks. The Energy System Scheduling Problem is particularly hard due to its binary control variables. Its complexity grows exponentially with the scale of the system, making it impossible to handle systems with more than a few variables. We developed a scalable approach for fine-grained scheduling of energy control systems that novelly combines techniques from control theory and computer science. The original system with binary control variables are approximated by an averaged system whose inputs are the utilization values of the binary inputs within a given period. The error between the two systems can be bounded, which allows us to derive a safety constraint for the averaged system so that the original system's safety is guaranteed. To further reduce the complexity of the scheduling problem, we abstract the averaged system by a simple single-state single-input dynamical system whose control input is the upper-bound of the total demand of the system. This model abstraction is achieved by extending the concept of simulation relations between transition systems to allow for input constraints between the systems. We developed conditions to test for simulation relations as well as algorithms to compute such a model abstraction. As a consequence, we only need to solve a small linear program to compute an optimal bound of the total demand. The total demand is then broken down, by solving a linear program much smaller than the original program, to individual utilization values of the subsystems, whose actual schedule is then obtained by a low-level scheduling algorithm. Numerical simulations in Matlab show the effectiveness and scalability of our approach.