The numerical solution of the electronic structure problem in Kohn-Sham density functional theory may in certain cases yield fractional occupancy of the single-particle orbitals. In this paper, we propose a quadratically convergent approach for simultaneous optimization of orbitals and occupancies in systems with fractional occupation numbers (FONs). The starting guess for orbitals and FONs is obtained via the relaxed constraint algorithm. Numerical results are presented for benchmark cases. (C) 2003 American Institute of Physics.