Linearized models of power systems are often desirable to formulate tractable control and optimization problems that still reflect real-world physics adequately under various operating conditions. In this paper, we propose an approach that can make use of known data concerning the distribution of demand, and/or intermittent supply, to minimize expected model inconsistency with respect to the original non-linear model. The optimal linearization is obtained by approximating a generalized moment problem with a hierarchy of sparse semi-definite relaxations. The output is a linearization point that minimizes the expected absolute constraint violation with respect to the uncertain supply and demand. Numerical results for different power systems networks demonstrate the accuracy and scalability of our linearization method.