This paper concerns optimization and analysis of telescope-coronagraph systems for direct imaging of exoplanets. In this paper, the coronagraph system, with arbitrary pupil geometry, is theoretically considered, and the governing equation for the pupil design is derived. The method of moments is applied to solve the generalized energy-concentration eigenvalue problem to obtain the optimal pupil apodization and complete sets of orthonormal basis functions for arbitrary pupil geometries. The method yields eigenvalues indicating the fraction of starlight energy encircled in the area of the focal-plane mask (FPM), where starlight can be occulted and/or nulled. In other words, a higher eigenvalue implies less leakage/spillover of light outside of the FPM region and into the planet-discovery zone. Thus, a higher eigenvalue supports better starlight suppression for a given type of coronagraph. This methodology is useful for semi-quantitatively ranking different modes of perturbation with respect to energy spillage in the focal plane independent of coronagraph design details. A model-order-reduction-based sensitivity analysis is conducted to investigate the coupling between different pupil modes induced by aberrations. A pupil mode recovery scheme is presented to offer a theoretically rigorous and computationally efficient approach to reconstruct the optimal pupil mode under an arbitrary phase perturbation. The reconstruction coefficients and recovery-effectiveness factors are derived theoretically and demonstrated numerically. Several numerical examples, including the LUVOIR A and B pupils, are provided to validate and demonstrate the applicability of the proposed methods. The reported methodology enables model-order reduction based on degree of focal-plane energy concentration and reconstruction of optimal pupil apodization vis-á-vis phase aberrations using a precomputed basis set. These features should enhance computational efficiency for coronagraph design and sensitivity analysis.