In some multitubular packed bed catalytic reactors, a graded reactor activity profile in the axial direction is created by loading the reactor with mixtures of active and inert pellets. This graded bed technique is developed to balance the reaction heat generation with reactor heat removal capacities and potentially improve reactor performance. In this study, we propose a model-based mathematical optimization approach to systematically design the optimal catalytic activity profile, taking the partial oxidation of ortho-xylene as the example. To carry out the task, we propose a multizone optimization formulation where the constraints are complicated by differential-algebraic equations (DAEs) from the reactor model. The reactor is described by one-dimensional mass and energy balance equations adopted from the so-called a model. We discretize the DAEs with orthogonal collocation over finite elements, which gives rise to a translated nonlinear optimization problem that can be built in generic algebraic modeling system (GAMS) and solved with off-the-shelf nonlinear programming (NLP) solvers. Our optimization results show that for ortho-xylene oxidation, compared with uniform catalyst distribution, the addition of a second optimized activity zone increases production by 26%, and further improvements are available with more zones. The methodology presented is generic and can accommodate different chemistries and constraints, which may result in significantly different optimization results and grading patterns.