The standard harmonic balance approach is applicable if the system under consideration can be split into a linear part $G(s)$ and a static nonlinear function $u = f(e) = f(r - y)$.
In this scenario it is always assumed that $f(e)$ is a function that takes a scalar argument.
However, is harmonic balance still applicable if $f:\mathbb{R}^n\rightarrow\mathbb{R}$ and $G(s)$ is a transfer function matrix with $n$ outputs?
For example, let $G(s)$ have two outputs, $y_1$ and $y_2$ and $f$ be a 2D function.
Can harmonic balance be somehow used to analyse such a system?