We propose a level set-based methodology to design optimal shell structures accounting for uncertainties in the material properties and the applied loads. Our approach aims to minimize the expectation and the variance of the compliance function. The method yields an optimized material distribution within a fixed optimization domain, defined by an isogeometric linear elastic Kirchhoff-Love shell \cite{HubnerScherer2024}. The Level Set Method (LSM) is adopted to describe the material properties and the topology of the domain, iteratively adapting the shape based on the computed descent direction of the cost function, given by a distributed shape derivative in a probabilistic setting. The structures are optimized in two distinct random modeling scenarios. First, real random variables are used to model uncertainties in the orientation of the applied loads. Second, Matérn-type Gaussian random fields are generated to model continuous distributions of material imperfections within the domain. These probabilistic effects are captured in the computation of the shape derivative. The Gaussian random fields are sampled using the Whittle Stochastic Partial Differential equation in an isogeometric discretization \cite{Duswald2024}. We furthermore analyze the resulting shape under different penalizations of the variance. The implementation is based on the open source tIGAr library with a penalty based approach to handle non-matching patches based on the PENGoLINS library \cite{Zhao2022}. The inclusion of uncertainties in the optimization ensures a robust result, such that the developed methodology can serve as an optimization setting for industrial applications.