I'm looking forward to your comments to the following solution.
We will solve the problem for the constraints $\mathbf{1}^T\mathbf{x}=\mathbf{1}^T\mathbf{y}_i=\mathbf{1}^T\mathbf{x}=1$ and for a general matrix $\mathbf{A}$ (not necessarily non-negative or symmetric).
First of all, since $f$ is linear with respect to $\mathbf{x}$, an optimal solution is attained at an extreme point, thus we can replace the constraint $\mathbf{x}\in [0,1]^n$ by $\mathbf{x}\in\left\{0,1\right\}^n$. (Of course we have the same statement for $\mathbf{y}$ and $\mathbf{z}$ but we do not need that).
Since $\mathbf{x}\in\left\{0,1\right\}^n$ and $\mathbf{1}^T\mathbf{x}=1$, $\mathbf{x}$ has the form: one component of $\mathbf{x}$ is $1$ and all the others are $0$. Thus, $\mathbf{x}$ can only take $n$ different values. The idea is that we iterate over these $n$ values of $\mathbf{x}$ and minimize the function with respect to $(\mathbf{y},\mathbf{z})$ (which is efficient because we have an analytical solution), then we compare all the $n$ obtained optimal values and pick the best one.
Now, suppose that $\mathbf{x} = \mathbf{x}_0$ where $x_{i_0} = 1$ and $x_i=0$ for any $i\neq i_0$. We have
\begin{align*}f(\mathbf{x}_0,\mathbf{y},\mathbf{z}) &=\mathbf{z}^T\mathbf{A}\mathbf{y}_{i_0} + b_{i_0} + \sum_{1\le i\le n}\mathbf{c}_i^T\mathbf{y}_i +\mathbf{d}^T\mathbf{z} \\
&= \left( \mathbf{z}^T\mathbf{A}\mathbf{y}_{i_0} + \mathbf{c}_{i_0}^T\mathbf{y}_{i_0} +\mathbf{d}^T\mathbf{z} \right) + \sum_{1\le i\le n,i\neq i_0}\mathbf{c}_i^T\mathbf{y}_i + b_{i_0}.
\end{align*}
Denote $\mathbf{y}_i = (y^i_1,y^i_2,\ldots,y^i_p)$. Clearly, the minimum value of $\mathbf{c}_i^T\mathbf{y}_i$ is the smallest component of $\mathbf{c}_i$, and if $j_0(i)$ is the position of this component, then the minimum is attained when $y^i_{j_0(i)} = 1$ and $y^i_j = 0 \forall j\neq j_0(i).$
It remains to minimize $g(\mathbf{y}_{i_0},\mathbf{z}) = \mathbf{z}^T\mathbf{A}\mathbf{y}_{i_0} + \mathbf{c}_{i_0}^T\mathbf{y}_{i_0} +\mathbf{d}^T\mathbf{z}$.
Since $\mathbf{1}^T\mathbf{y}_{i_0} = \mathbf{1}^T\mathbf{z}=1$ we have
$$g(\mathbf{y}_{i_0},\mathbf{z}) = \mathbf{z}^T\mathbf{A}\mathbf{y}_{i_0} + \mathbf{z}^T\mathbf{1}\cdot\mathbf{c}_{i_0}^T\mathbf{y}_{i_0} +\mathbf{z}^T\mathbf{d}\cdot \mathbf{1}^T\mathbf{y}_{i_0} = \mathbf{z}^T \left( \mathbf{A} + \mathbf{1}\cdot\mathbf{c}_{i_0}^T + \mathbf{d}\cdot \mathbf{1}^T \right) \mathbf{y}_{i_0}.$$
According to this lemma: Find $\min x^TAy$ subject to $1^Tx=1^Ty=1,x\ge 0,y\ge 0$, the minimum value of $g$ is the smallest component of the matrix $\mathbf{A} + \mathbf{1}\cdot\mathbf{c}_{i_0}^T + \mathbf{d}\cdot \mathbf{1}^T$, and if $(k_0,j_0)$ is the position of this component, then the minimum is attained when $z_{k_0} =1$ and $y^{i_0}_{j_0}=1$.
Finally, we have the following algorithm to solve the problem:
Initialize the best-so-far minimum: $f_0=+\infty$.
Pre-compute the position $j_0(i)$ of the smallest components of the vectors $\mathbf{c}_i$. Compute the sum of these components: $S_c = \sum_{1\le i\le n}c^i_{j_0(i)}$.
For $i=1,2,\ldots,n$:
+) Find the smallest component $m_{k_0,j_0}$ of the matrix $\mathbf{M}=\mathbf{A} + \mathbf{1}\cdot\mathbf{c}_{i}^T + \mathbf{d}\cdot \mathbf{1}^T$.
+) Compute the current minimum $f=m_{k_0,j_0} + (S_c - c^i_{j_0(i)}) + b_{i}$
+) Compare with the best-so-far minimum: if $f \le f_0$ then $f_0 = f$ and set the corresponding optimal solution.
