We present a framework to solve a finite-time optimal control problem for parabolic partial differential equations (PDEs) with diffusivity-interior actuators, which is motivated by the control of the current density profile in tokamak plasmas. The proposed approach is based on reduced order modeling (ROM) and successive optimal control computation. First we either simulate the parabolic PDE system or carry out experiments to generate data ensembles, from which we then extract the most energetic modes to obtain a reduced order model based on the proper orthogonal decomposition (POD) method and Galerkin projection. The obtained reduced order model corresponds to a bilinear control system. Based on quasi-linearization of the optimality conditions derived from Pontryagin’s maximum principle, and stated as a two-boundary-value problem, we propose an iterative scheme for suboptimal closed-loop control. We take advantage of linear synthesis methods in each iteration step to construct a...