A new methodology has been developed for determining the optimal operation of water distribution systems for water quality purposes. The methodology is based upon describing the operation as a discrete time optimal control problem that can be used to determine the optimal operation schedules of the pumps in distribution systems. The solution methodology is based upon a mathematical programming approach resulting in a large scale nonlinear programming problem which can not be solved by using existing nonlinear codes. The solution of the optimization problem is obtained by interfacing a simulation code, EPANET, which solves the hydraulic and water quality constraints with a nonlinear optimization code, GRG2. Bound constraints on the state variables are incorporated into the objective function using the augmented Lagrangian and the bracket penalty methods. Three objective functions can be used in the model: (1) the minimization of the deviations of actual substance concentrations from the desired concentration values, (2) the minimization of the total pump duration times, and (3) the minimization of the total energy cost. The effectiveness of methodology has been tested using both hypothetical and existing water distribution systems.