Shape optimization models with one or more shapes are considered in this chapter. Of particular interest for applications are problems in which where a so-called shape functional is constrained by a partial differential equation (PDE) describing the underlying physics. A connection can made between a classical view of shape optimization and the differential-geometric structure of shape spaces. To handle problems where a shape functional depends on multiple shapes, a theoretical framework is presented, whereby the optimization variable can be represented as a vector of shapes belonging to a product shape space. The multi-shape gradient and multi-shape derivative are defined, which allows for a rigorous justification of a steepest descent method with Armijo backtracking. As long as the shapes as subsets of a hold-all domain do not intersect, solving a single deformation equation is enough to provide descent directions with respect to each shape. Additionally, a framework for handling uncertainties arising from inputs or parameters in the PDE is presented. To handle potentially high-dimensional stochastic spaces, a stochastic gradient method is proposed. A model problem is constructed, demonstrating how uncertainty can be introduced into the problem and the objective can be transformed by use of the expectation. Finally, numerical experiments in the deterministic and stochastic case are devised, which demonstrate the effectiveness of the presented algorithms. * This work has been partly supported by the state of Hamburg within the Landesforschungsförderung under project "Simulation-Based Design Optimization of Dynamic Systems Under Uncertainties" (SENSUS) with project number LFF-GK11, and by the German Academic Exchange Service (DAAD) within the program "Research Grants-Doctoral Programmes in Germany, 2017/18."