I am currently facing the problem of interpolating between two unit quaternions $q\left(t=0\right)=q_0 \in SO\left(3\right)$ and $q\left(t=t_{\mathrm{target}}\right)q_1 \in SO\left(3\right)$, where $t$ indicates time. Unfortunately there is an additional constraint on the angular velocity that is $\omega\left(t=0\right) = \omega_0$ and $\omega\left(t=t_{\mathrm{target}}\right) = \omega_1$. All angular velocities are assumed to be vectors of dimension 3x1. I know that SLERP would be a possible approach if we weren't facing constraints on the initial and target velocity. Summing up, I seek to solve the following problem:
Find
$q\left(t\right) = f\left(q_0, q_1, \omega_0, \omega_1, t\right)$
such that:
$q\left(0\right) = f\left(q_0, q_1, \omega_0, \omega_1, 0\right) = q_0$
$q\left(t_{\mathrm{target}}\right) = f\left(q_0, q_1, \omega_0, \omega_1, t_{\mathrm{target}}\right) = q_1$
$\omega\left(t=0\right) = \omega_0$
$\omega\left(t=t_{\mathrm{target}}\right) = \omega_1$
$\vert\vert q \vert\vert = 1$
Are there existing techniques or procedures to tackle such kind of a problem?
Let's define the unit quaternion $\mathbf{q}$ as a function of time $t$, $0 \le t \le 1$, as $$\mathbf{q}(t) = \left ( \, x(t) sin(\varphi(t)) ,\;\, y(t) sin(\varphi(t)) ,\;\, z(t) sin(\varphi(t)) ,\;\, cos(\varphi(t)) \, \right ) \tag{1}\label{NA1}$$ where $$\hat{r}(t) = \left( x(t) , y(t) , z(t) \right) \tag{2}\label{NA2}$$ is the axis of rotation at time $t$. For $\mathbf{q}(t)$ to be an unit quaternion for all $t$, $$\left\lVert\hat{r}(t)\right\rVert = 1 \tag{3}\label{NA3}$$ or, equivalently, $x(t)^2 + y(t)^2 + z(t)^2 = 1$.
The initial condition can now be expressed as $$\begin{cases} \hat{r}(0) = \hat{r}_0 \\ \varphi(0) = \varphi_0 \\ \left.\frac{d \varphi(t)}{d t}\right\rvert_{t=0} = \omega_0 \end{cases} \tag{4}\label{NA4}$$ and the final condition as $$\begin{cases} \hat{r}(1) = \hat{r}_1 \\ \varphi(1) = \varphi_1 \\ \left.\frac{d \varphi(t)}{d t}\right\rvert_{t=1} = \omega_1 \end{cases} \tag{5}\label{NA5}$$ Note that $\omega_0$ denotes the initial rate of rotation around axis $\hat{r}_0$, and $\omega_1$ the final rate of rotation around axis $\hat{r}_1$, and differ from OP's usage.
There are two different approaches I can think of. The first one is to interpolate $\hat{r}(t)$ and $\varphi(t)$ separately. The second one is to use spherical quadratic interpolation (with an extra intermediate unit quaternion), or spherical cubic interpolation (with two extra intermediate unit quaternions).
Let's examine the first option.
We can do simple spherical linear interpolation for the rotation axis vector $\hat{r}(t)$. When the rotation axis $\hat{r}(t)$ varies (i.e., $\hat{r}_0 \ne \hat{r}_1$ and $\hat{r}_0 \ne -\hat{r}_1$), we have $\lVert\hat{r}_0\times\hat{r}_1\rVert\ne 0$, and $$\hat{a} = \frac{ \hat{r}_0 \times \hat{r}_1 }{\left\lVert \hat{r}_0 \times \hat{r}_1 \right\rVert}$$ If $\hat{r}_1 = -\hat{r}_0$, use $\hat{a} = \vec{r}_0$, but negate $\varphi_1$ and $\omega_1$. If $\hat{r}_0 = \hat{r}_1$, then use $\hat{a} = \hat{r}_0 = \hat{r}_1$.
If we define angle $\theta$ as $$\begin{cases} \cos(\theta) = \frac{\hat{r}_0 \cdot \hat{r}_1}{\left\lVert \hat{r}_0 \right \rVert \, \left \lVert \hat{r}_1 \right \rVert} = \hat{r}_0 \cdot \hat{r}_1 \\ \sin(\theta) = \frac{\left\lVert\hat{r}_0 \times \hat{r}_1\right\rVert}{\left\lVert\hat{r}_0 \right\rVert \, \left\lVert \hat{r}_1 \right\rVert} \end{cases} = \left\lVert\hat{r}_0 \times \hat{r}_1\right\rVert$$ then $\hat{r}(t)$ is unit vector $\hat{r}_0$ rotated around axis $\hat{a}$ by angle $t \theta$. You can use e.g. Rodrigues' rotation formula to implement this. The rotated vector is obviously also an unit vector.
That leaves the function $\varphi(t)$. Let's define it as a simple cubic function, $$\varphi(t) = \phi_0 + t \phi_1 + t^2 \phi_2 + t^3 \phi_3$$ The initial and final conditions are now $$\begin{cases} \varphi(0) = \varphi_0 \\ \varphi(1) = \varphi_1 \\ \left.\frac{d \varphi(t)}{d t}\right\rvert_{t=0} = \omega_0 \\ \left.\frac{d \varphi(t)}{d t}\right\rvert_{t=1} = \omega_1 \end{cases}$$ Solving the system of four equations gives us the coefficients $$\begin{cases} \phi_0 = \varphi_0 \\ \phi_1 = \omega_0 \\ \phi_2 = 3 (\varphi_1 - \varphi_0) - \omega_1 - 2 \omega_0 \\ \phi_3 = 2 (\varphi_0 - \varphi_1) + \omega_0 + \omega_1 \end{cases}$$ which fulfills both the initial and final rotations $\varphi_0$ and $\varphi_1$, but also the rates of rotation $\omega_0$ and $\omega_1$.
Although this approach requires the least number of additional assumptions, you will have to evaluate both $\sin$ and $\cos$ at every iteration; and although the axis of rotation $\hat{r}(t)$ traces an arc of $\theta$ of a great circle on the unit sphere, there is no guarantees the rotation itself is any kind of "minimum" or "maximum" choice anymore, because the amount of rotation around the axis varies.