Motivated by several real-world applications in continuum mechanics, mathematical finance or life sciences, the talk highlights recent advances in the analysis and, in particular, the development of robust numerical solution schemes for mathematical programs with equilibrium or complementarity constraints (MPECs or MPCCs), generalized Nash equilibrium problems (GNEPs) and (other) quasi-variational inequalities (QVIs) involving partial differential operators, respectively. For MPECs and MPCCs the derivation of suitable stationarity conditions (such as appropriate variants of C-, M- and strong stationarity) is briefly addressed, and for GNEPs as well as for QVIs the existence of solutions and associated characterizations, which are suitable for numerical realization, are discussed. In particular, for GNEPs novel constraint qualifications are presented, and for QVIs a rather general approach to existence is introduced. For all problems addressed in the talk efficient numerical solution schemes relying on relaxation concepts, the Moreau-Yosida-regularization or semismooth Newton techniques are presented. Also, multilevel approaches relying on adaptive discretization are briefly mentioned, and the practical behavior of the presented solvers is highlighted.