Methods based on the chain-of-states approach are very useful in bridging timescale gaps for the studies of rare events. However, the requirements of ensuring equal distances between replicas along a chain significantly affect the optimization of objective functions and lead to slow convergence and numerical instability during path optimization. In this work, we present a new algorithm that represents the equal-distance requirements as holonomic constraints. This scheme thus decouples path optimization from equal-spacing replicas under the framework of constrained optimization using Lagrange multipliers. Straightforward implementation of super-linear minimization schemes such as the Adopted Basis Newton Raphson (ABNR) method can thus be achieved. As a result, the quality of convergence and stability of path optimization is greatly improved compared to conventional methods such as the nudged elastic band (NEB) method. In addition to minimum energy paths (MEP), we also apply this method to find the minimum free energy paths (MFEP) in molecular systems. Moreover, a novel hyper plane constraint is also developed to compute the potential of mean force along an optimized path. This set of new methods is applied to study the isomerization of alanine dipeptide and the helix-to-hairpin transition of an amyloid beta peptide. To model the helix-to-hairpin structural transition, as high as 513 replicas were used, resulting in more than 60,000 atoms in the chain for path optimization. This example thus demonstrates the robustness and efficiency of our scheme of constrained path optimization. Finally, outlook of applying the methods developed in this work to simulate transitions in complicated molecular systems will be given.