Cannot multiply a quadratic expression by an affine expression

Does that mean I can’t define matrix product at all?

Correct. In the current version of JuMP, you cannot create expressions like X * Y * Z where X, Y and Z are arrays of JuMP variables.

One work-around is:

using JuMP, LinearAlgebra
n = 5
model = Model()
@variable(model, R[1:n, 1:n], Symmetric)
@variable(model, Q[1:n, 1:n], Symmetric)
@variable(model, RQ[1:n, 1:n], Symmetric)
@constraint(model, RQ .== R * Q)
@expression(model, P, Q' * QR)

This limitation is fixed in the work-in-progress nonlinear rewrite: https://github.com/jump-dev/JuMP.jl/pull/3106.