In this paper, we develop an effective scheme to estimate stability regions by using an energy function that is a generalization of the Lyapunov functions. It is shown that the scheme can optimally estimate stability regions. A fairly comprehensive study for the structure of the constant energy surface lying inside the quasi-stability region is presented. A topological characterization, as well as a dynamical characterization for the point on the quasi-stability boundary and the point on the stability boundary with the minimum value of an energy function are derived. These characterizations are then used in the development of a computational scheme to estimate quasi-stability regions. By utilizing an energy function approach (or Lyapunov function approach), this scheme can significantly reduce conservativeness in estimating the stability region, because the estimated stability region characterized by the corresponding energy function is the largest one within that stability region.