A stochastic modeling and uncertainty analysis methodology for energy system synthesis/design is proposed in this paper and applied to the development of the fuel processing subsystem (FPS) of a proton exchange membrane fuel cell (PEMFC) system. The FPS consists of a steam methane reformer, both high and low temperature water-gas shift reactors, a CO preferential oxidation reactor, a steam generator, a combustor, and several heat exchangers. For each component of the system, detailed thermodynamic, geometric, chemical kinetic, and cost models are developed and integrated into an overall model for the subsystem. Conventionally, in energy system synthesis/design, such models are treated deterministically, using a specific set of non-probabilistic input variable values that produce a specific set of non-probabilistic output variable values. Even though these input values, which include the specific load profile (i.e. electrical, thermal, and/or aerodynamic) for which the system or subsystem is synthesized/designed, can have significant uncertainties that inevitably propagate through the system to the outputs, such deterministic approaches are unable to quantify these uncertainties and their effect on the final synthesis/design and operation/control. This deficiency can, of course, be overcome by treating the inputs and outputs probabilistically. The difficulty with doing this, particularly when large-scale dynamic optimization with a large number of degrees of freedom is being used to determine the optimal synthesis/design and operation/control of the system, is that the traditional probabilistic approaches (e.g., Monte Carlo Method) are so computationally intensive that combined with large-scale optimization it renders the problem computationally intractable. This difficulty can be overcome by the use of approximate approaches such as the response sensitivity analysis (RSA) method based on Taylor series expansion. In this study, RSA is employed and developed by the authors for use with dynamic energy system optimization. Load profile and cost models are treated as probabilistic input values and uncertainties in output results investigated. The results for the uncertainty analysis applied to the optimization of the FPS synthesis/design and operation/control are compared with those found using a Monte Carlo approach with good results. In this paper, the FPS synthesis/design and operation optimization is treated as a multi-objective optimization problem to minimize the capital cost and operating cost simultaneously, and uncertainty effects on the optimization are assessed by taking uncertainties into account in the objectives and constraints. Optimization results show that there is little effect on the objective (the operating cost and capital cost), while the constraints (e.g., that on the CO concentration) can be significantly affected during the synthesis/design and operation/control optimization.