Achieving robust and scalable convergence for simulation of realistic power flow cases can be challenging. One specific issue relates to the disconnected solution space that is created by the use of piecewise-discontinuous models of power grid devices that perform control mechanisms. These models are generally resolved by outer iteration loops around power flow, which can result in solution oscillations, increased iteration count, divergence or even convergence to a solution in an unstable operational region. This paper introduces a continuously differentiable model for device control mechanisms that is incorporated within the power flow formulation. To ensure robust power flow convergence properties, recently introduced homotopy methods are extended to include these continuous models. The scalability and efficacy of the proposed formulation is demonstrated on several large-scale test cases that represent the US Eastern Interconnect network, the Synthetic USA, and the Nigerian grid.