A sequential programming method of first order for constrained optimization in infinite dimensional spaces is presented. It is referred to as the SP method. It relies on Lipschitz stability of the minimizers with respect to perturbations in an essential way, and is motivated by optimal control problems with partial differential equations as constraints. Convergence and convergence rate are analyzed based on solutions of first order optimality conditions and on Lagrange multiplier theory. The first order SP method does not rely on an approximation of the Hessian of the Lagrange functional and consequently it avoids instabilities due to possible indefiniteness far from a local minimum. It is especially well suited for bilinear control problems and can be extended to certain classes of nonsmooth and convex problems. It can efficiently be implemented by the use of saddle point solvers, with complexity which is between that of gradient and SQP methods. For the important class of bilinear control problems the method is stable when using damped updates. We also develop a globalization strategy, and a second order variant. The proposed methods are numerically tested for control in the coefficient problems or, equivalently, for bilinear optimal control problems.