Mathematica cannot handle this, as there is no algorithm that will deterministically achieve the global optimum. What you can do is do local optimization for several choices of initial points, and from that hopefully get the global optimum.
The algorithm normally used in this case is the see-saw. The idea is that your probability of success is a bilinear function of the states and observables. For a fixed choice of observables, the probability of success is a linear function of the states, which is trivial to optimize, and for a fixed choice of states the probability of success is a linear function of the observables, which is again trivial to optimize. You choose then randomly some initial observables, get the optimal states for them, then for these states obtain the optimal observables, then for these states get the optimal observables, and so on, until you converge.
More precisely, the probability of success is given by
$$p_\text{succ} := \sum_{x,y,b} c_{xyb} \operatorname{tr} \rho_x E^b_y,$$ where $\rho_x$ will be the 8 different states you use to encode the bits, and $\{E^0_y,E^1_y\}$ the 7 different POVMs that should get the functions you want. The coefficients $c_{xyb}$ are then determined by the functions.
Define then
$$ A^b_y := \sum_{x} c_{xyb} \rho_x $$
and
$$ B_x:= \sum_{y,b} c_{xyb} E^b_y,$$
so that
$$p_\text{succ} = \sum_{x} \operatorname{tr} \rho_x B_x = \sum_{yb} \operatorname{tr} A^b_y E^b_y. $$
Then a given choice of observables determines $B_x$, and you get the optimal states for them by maximizing $\sum_{x} \operatorname{tr} \rho_x B_x $, and a given choice of states determines $A^b_y$, and you get the optimal observables for them by maximizing $\sum_{yb} \operatorname{tr} A^b_y E^b_y$.