Edit: Here is what I think is a counter-example.
Let $\phi$ be the indication function of the even natural numbers, let $\mathcal U$ be an ultrafilter supported on the even naturals, and let $\mathcal V$ be an ultrafilter defined on the even negative integers. Define $\mu\in M(\mathbb Z)$ by
$$\int f(x) \ d\mu(x) = \lim_{x\in\mathcal V} f(x).$$
Then $\mu\in I_\phi$ (as, in fact, $\int \phi(x+y) \ d\mu(y)=0$ for all $x$). If $\nu\in I_\phi$ then $\nu$ must assign the same measure to $2\mathbb N$ and $2\mathbb N+1$, say $\alpha\leq 1/2$. You also need to argue that $\nu$ must assign zero measure to any finite set (else it won't be $\phi$-invariant). So for any $y\in\mathbb Z$,
$$\int \phi(x+y) \ d\nu(x) = \nu(2\mathbb N-y) = \begin{cases}
\nu(2\mathbb N) &: y\in 2\mathbb Z, \\ \nu(2\mathbb N+1) &: y\in 2\mathbb Z+1, \end{cases}
= \alpha.$$
Thus
$$\int \int \phi(x+y) \ d\nu(x) \ d\mu(y) = \int \alpha \ d\mu(y) = \alpha,$$
as $\mu$ is a probability measure. By contrast, let $\nu$ be defined by
$$\int f(x) \ d\nu(x) = \lim_{y\in\mathcal U} f(y).$$
Then
$$\int \phi(x+y) \ d\nu(x) = \begin{cases} 1 &: y\in 2\mathbb Z, \\
0 &: y \in 2\mathbb Z+1, \end{cases}$$
and so
$$\int \int \phi(x+y) \ d\nu(x) \ d\mu(y) = \int \chi_{2\mathbb Z}(y) \ d\mu(y)
= 1.$$
So $F$ is not maximised on $I_\phi$.
In fact, by replacing $2\mathbb N$ by $k\mathbb N$, I think you get that $F$ has norm one, but $F(\nu)\leq 1/k$ for any $\nu\in I_\phi$.
But somehow, to my mind, what's wrong is that the $\mu\in I_\phi$ you choose is very poor. So here's a revised conjecture:
Let $\mu\in I_\phi$ maximise the integral $\int \phi(x) \ d\mu(x)$. Then $F$ attains its maximum on $I_\phi$.
Old post: (Explains my thinking).
I think of these questions using the Arens products, from abstract Banach algebra theory. So I work over the complex numbers; but this is not a problem.
Consider $A=\ell^1(\mathbb Z)$ with the convolution product, so $A$ is commutative. Then $A^*=\ell^\infty(\mathbb Z) = C(\beta\mathbb Z)$ is an $A$-module: $(a\cdot f)(b) = f(ba)$ for $a,b\in A,f\in A^*$. Then $A^{**}=M(\beta\mathbb Z)$ the space of finite Borel measures on the Stone-Cech compactification $\beta\mathbb Z$. Your space $M(\mathbb Z)$ is just the positive measures $\mu\in A^{**}$ with $\mu(1)=1$.
We try to extend the product of $A$ to $A^{**}$. Firstly we define a bilinear map $A^{**}\times A^*\rightarrow A^*$ by
$$(\mu\cdot f)(a) = \mu(a\cdot f) \qquad (\mu\in A^{**}, f\in A^*, a\in A).$$
But then we have two choices for the product on $A^{**}$:
$$(\mu \Box \lambda)(f) = \mu(\lambda\cdot f), \quad
(\mu\diamond\lambda)(f) = \lambda(\mu\cdot f)
\qquad (\mu,\lambda\in A^{**}, f\in A^*).$$
A little thought shows that $\mu\diamond\lambda = \lambda\Box\mu$.
So if $\phi\in A^*$ if positive then $\mu\in I_\phi$ if and only if $\mu\cdot\phi
= \mu(\phi) 1$. This follows, as writing $\delta_x\in A=\ell^1(\mathbb Z)$ for the point mass at $x\in\mathbb Z$, we have
$$(\phi\cdot\delta_x)(\delta_y) = \phi(\delta_{x+y}) \implies
(\mu\cdot\phi)(\delta_x) = \mu(\phi\cdot\delta_x)
= \int \phi(x+y) \ d\mu(y).$$
So the condition that $\mu\in I_\phi$ becomes that $(\mu\cdot\phi)(\delta_x)$ is constant in $x$, which is seen to be equivalent to $\mu\cdot\phi = \mu(\phi) 1$.
Similarly, your map $F$ is just $F(\nu) = (\mu\Box\nu)(\phi)$.
As you allude to, it's known that $\lambda\Box\mu \not= \mu\Box\lambda$ for arbitrary $\lambda,\mu$. However, we say that $f\in A^*$ is "weakly almost periodic" (WAP) if $(\lambda\Box\mu)(f) = (\mu\Box\lambda)(f)$ for all $\mu,\lambda\in A^{**}$. So if $\phi$ is WAP and $\mu\in I_\phi$ then for any $\nu\in M(\mathbb Z)$,
$$F(\nu) = (\mu\Box\nu)(\phi) = (\nu\Box\mu)(\phi) = \nu(\mu\cdot\phi)
= \nu(1) \mu(\phi) = \mu(\phi),$$
as $\nu$ is a probability measure. So actually $F$ is constant on $M(\mathbb Z)$ and so certainly attains its maximum at a point of $I_\phi$.
So, to be interesting, we need to ask the question for $\phi$ which are not WAP. An alternative characterisation of $\phi$ being in WAP is that the set of translates of $\phi$ in $\ell^\infty(\mathbb Z)$ forms a relatively weakly compact set. A nice characterisation of Grothendieck shows that this is equivalent to
$$\lim_n \lim_m \phi(x_n+y_m) = \lim_m \lim_n \phi(x_n+y_m)$$
whenever all the limits exist, for sequences $(x_n),(y_m)$ in $\mathbb Z$. If $\phi$ is the indicator function of $\mathbb N$, then it's not in WAP.
We may as well assume that $\|\phi\|_\infty=1$.
Another "easy" case is when we can find $\nu\in I_\phi$ with $\nu(\phi)=1$. Then $F(\nu) = \mu(\nu\cdot\phi) = \mu(1) \nu(\phi) = 1$; while for any $\lambda\in M(\mathbb Z)$, clearly $|F(\lambda)| = |\mu(\lambda\cdot\phi)| \leq 1$ as $\mu$ is a probability measure, and $\lambda\cdot\phi$ is bounded by $1$ (again, as $\lambda$ is a probability measure and $\phi$ is bounded by $1$). Notice that this case covers your example of when $\phi$ is the indicator function of $\mathbb N$.
So a test case is to find $\phi$ not in WAP and with $\nu(\phi)<\|\phi\|_\infty$ for all $\nu\in I_\phi$ (notice that $I_\phi$ is always non-empty, as $\mathbb Z$ is amenable). Do you have an example of such a $\phi$?
Actually, if $\phi$ is the indicator function of the even natural numbers, then that's an example. And that leads to my (hopeful) counter-example.
