The open-source PlanetProfile framework was developed to investigate the interior structure of icy moons based on self-consistency and comparative planetology. The software, originally written in Matlab, relates observed and measured properties, assumptions such as the type of materials present, and laboratory equation-of-state data through geophysical and thermodynamic models to evaluate radial profiles of mechanical, thermodynamic, and electrical properties, as self-consistently as possible. We have created a Python version of PlanetProfile. In the process, we have made optimization improvements and added parallelization and parameter-space search features to utilize fast operation for investigating unresolved questions in planetary geophysics, in which many model inputs are poorly constrained. The Python version links to other scientific software packages, including for evaluating equation-of-state data, magnetic induction calculations, and seismic calculations. Physical models in PlanetProfile have been reconfigured to improve self-consistency and generate the most realistic relationships between properties. Here, we describe the software design and algorithms in detail, summarize models for major moons across the outer solar system, and discuss new inferences about the interior structure of several bodies. The high values and narrow uncertainty ranges reported for the axial moments of inertia for Callisto, Titan, and Io are difficult to reconcile with self-consistent models, requiring highly porous rock layers equivalent to incomplete differentiation for Callisto and Titan, and a high rock melt fraction for Io. This effect is even more pronounced with the more realistic models in the Python version. Radial profiles for each model and comparison to prior work are provided as Zenodo archives.