Quantifying effects of system-wide uncertainties (i.e., affecting structural, piezoelectric, and/or electrical components) in the analysis and design of piezoelectric vibration energy harvesters have recently been emphasized. The present investigation proposes first a general methodology to model these uncertainties within a finite element model of the harvester obtained from an existing finite element software. Needed from this software are the matrices relating to the structural properties (mass, stiffness), the piezoelectric capacitance matrix as well as the structural-piezoelectric coupling terms of the mean harvester. The thermal analogy linking piezoelectric and temperature effects is also extended to permit the use of finite element software that do not have piezoelectric elements but include thermal effects on structures. The approach is applied to a beam energy harvester. Both weak and strong coupling configurations are considered, and various scenarios of load resistance tuning are discussed, i.e., based on the mean model, for each harvester sample, or based on the entire set of harvesters. The uncertainty is shown to have significant effects in all cases even at a relatively low level, and these effects are dominated by the uncertainty on the structure versus the one on the piezoelectric component. The strongly coupled configuration is shown to be better as it is less sensitive to the uncertainty and its variability in power output can be significantly reduced by the adaptive optimization, and the harvested power can even be boosted if the target excitation frequency falls into the power saturation band of the system.