Multiple polynomial regression of third order in v≥2 real variables is considered, where the experimental region is a ball centered at zero or a symmetric cube centered at zero. Besides the full cubic regression the authors consider also symmetric submodels which are between the full quadratic and the full cubic models. The problem is to find (numerically) an optimal approximate design under a given convex and differentiable optimality criterion. The large dimension of the problem is considerably reduced by restricting to designs which are invariant w.r.t. the group of permutations and sign changes of coordinates. This is justified by the fact that many (convex) optimality criteria are invariant w.r.t. the induced matrix group. The authors develop the mathematical tools to apply algorithms from Gaffke & Mathar to compute (nearly) optimal designs. Numerical results are presented for mixtures of Kiefer’s ℝlsquo;Åp-criteria under a fixed regression model, and for mixtures over competing models. These mixture criteria aim at designs which are robust under change of criterion or change of regression model, respectively, and are thus of considerable practical interest.