The problem of controlling the total power and power distribution in a large pressurized water reactor (PWR) core to follow a known time-varying load schedule has been formulated as a multistage optimization problem. The control problem is solved subject to hard constraints, which can be applied on total power, control variables and their rate of change, local power densities and their rate of change, and on more global power distribution measures such as axial and quadrant offsets. Based on a three-dimensional linearized nodal core model with some slightly nonlinear features, the optimal control problem is solved by quadratic programming. The method, called multistage mathematical programming, has been studied in simulations. A large PWR core, which was unstable with respect to both axial and azimuthal xenon oscillations, was represented by a simplified three-dimensional nonlinear nodal core simulator model. The three-dimensional oscillations were successfully damped at constant load, and an efficient anticipatory control was obtained for load cycling operation.