I present an approach for constructing self-consistent dynamical models of our Galaxy based on distribution functions (DF) in the action space. The parameters of the DF are adjusted to match observed rotation curve and distribution of stars in the Solar neighborhood, and the gravitational potential is computed from the DF iteratively. With a dozen or so parameters, we have a complete description of an axisymmetric approximation of the Galaxy in dynamical equilibrium, which may be used to construct N-body models, serve as the basis for perturbation theory for non-axisymmetric modes, and much more.