We study the spin S = 1 antiferromagnetic Heisenberg model on the square lattice with, in addition to the nearest-neighbor interaction, a three-site interaction of the form (S-i . S-j)(S-j . S-k) + H.c. This interaction appears naturally in a strong coupling expansion of the two-orbital, half-filled Hubbard model. For spin 1/2, this model reduces to a Heisenberg model with bilinear interactions up to third neighbors, with a second-neighbor interaction twice as large as the third-neighbor one, a very frustrated model with an infinite family of helical classical ground states in a large parameter range. Using a variety of analytical and numerical methods, we show that the spin-1 case is also very frustrated, and that its phase diagram is even richer, with possibly the succession of seven different phases as a function of the ratio of the three-site interaction to the bilinear one. The phases are either purely magnetic phases with collinear order or of mixed magnetic and quadrupolar character with helical order.