Suppose $P(y,z)$ is a integer coefficient polynomial represented entirely in terms of the "variables" (though not linearly independent) of the form $y_i+z_j$ as $i$ and $j$ range over all positive integers with integer coefficients. The coefficients are not necessarily nonnegative, however we assume that there is an expression in terms of these types of monomials for $P(y,z)$ that has nonnegative integer coefficients. For example, we could have $$P(y,z)=(y_1+y_3+z_3+z_5)(y_2+y_4+z_6+z_7)-(y_1+z_5)(y_2+z_7)$$ and it should be easy to see that this does have some positive expression.
I'm trying to find any expression for a given $P(y,z)$ that has nonnegative integer coefficients, it can be any such representation if there is more than one (I just need one) and I have infinitely many polynomials $P(y,z)$ that I want to do this with.
The way I'm doing this is generating a set of $y_i+z_j$ monomials that can possibly occur in such an expression, then using integer programming (where we vectorize it by expanding all products and taking the coefficients of the resulting ordinary monomials). I have some ways of culling elements from this set, for example if one of the coefficients of the ordinary monomials in the expansion of a given $y_i+z_j$ monomial exceeds a coefficient of the final expression, we can remove that monomial from the set.
The most efficient way I found to reduce the size of the set of possible monomials is by looking at each ordinary monomial containing only $y$ variables, say $y_1^{a_1}\cdots y_n^{a_n}$, and, assuming without loss of generality that $a_1\neq 0$, look at monomials of the form $y_2^{a_2}\cdots y_n^{a_n}z_1^{b_1}\cdots z_n^{b_n}$ that occur, and adding factors of $y_1+z_i$ $b_i$ times, and doing this for each $y$ variable. In practice, it actually works for my application to assume that $b_i\in \{0,1\}$ for all $i$. Regardless, the number of basis elements ends up being huge when the degree exceeds $6$. Is there a better way?