We investigate numerical methods for solving large-scale saddle point systems which arise during the feedback control of flow problems. We focus on the instationary Stokes equations that describe instationary, incompressible flows for moderate viscosities. After a mixed finite element discretization we get a differential-algebraic system of differential index two [J. Weickert, Navier-Stokes Equations as a Differential-Algebraic System, Preprint SFB393/96-08, Department of Mathematics, Chemnitz University of Technology, Chemnitz, Germany, 1996]. To reduce this index, we follow the analytic ideas of [J.-P. Raymond, SIAM J. Control Optim., 45 (2006), pp. 790-828] coupled with the projection idea of [M. Heinkenschloss, D. C. Sorensen, and K. Sun, SIAM J. Sci. Comput., 30 (2008), pp. 1038-1063]. Avoiding this explicit projection leads to solving a series of large-scale saddle point systems. In this paper we construct iterative methods to solve such saddle point systems by deriving efficient preconditioners based on the approaches of Wathen and colleagues, e. g., [M. Stoll and A. Wathen, J. Comput. Phys., 232 (2013), pp. 498-515]. In addition, the main results can be extended to the nonsymmetric case of linearized Navier-Stokes equations. We conclude with numerical examples showcasing the performance of our preconditioned iterative saddle point solver.