We report a theoretical framework for the study of the optimal control of multisurface molecular systems via a set of nondegenerate excitation fields. The resulting control equations in the strong response regime are presented in terms of both the Liouville-space density matrix dynamics and the Hilbert-space wave function evolution. We further derive a pair of eigenequations for the optimal pump-pump fields in the pure-state control of three-surface molecular systems in the weak response regime. The globally optimal pair of pump-pump fields in this case are identified. Application to the control of a rovibronic level on the final excited surface reveals a symmetry relation within the optimal pair of pump-pump fields in the weak response regime. For numerical demonstrations, we consider the control of the I-2 molecular system involving the initial ground X, the intermediate B, and the final E surface. The target is chosen as an outgoing vibrational wave packet in the hound region of the final E electronic state. The optimal control fields in both the strong and weak response regimes are calculated and further parameterized to fit simple experimentally realizable laser pulses.