TL;DR: No, that problem is not hard.
Synopsis: After remapping of $\Bbb Z_p$ by an involution $x\to\overline x$, the function $(x,y)\to\overline{-f(\overline x,\overline y)}$ is mostly associative. We massage it into an Abelian finite group $(\mathcal S,\boxplus)$. This makes $\overline{r_n}$ a linear function of $\overline a$ and $\overline b$ with known scalar coefficients, allowing to efficiently solve the question's problem.
Define $\delta$ and $\hat z$ in $\Bbb Z_p$ with $\delta=(p+1)/2$ and $\hat z=(-3\,\delta^2-z)\bmod p$. The reason will become apparent later, for now consider $\hat z$ as an arbitrary fixed public element of $\Bbb Z_p$.
Define $l$ as the Legendre symbol $l=\displaystyle\biggl(\frac{\hat z}p\biggr)$, and define $m=p-l$.
When $l=+1$, define $\omega$ as a particular solution of $\omega^2=a$ in $\Bbb Z_p$, e.g. the odd one in range $[1,p)$.
Define the set $\mathcal S$ as:
$$\mathcal S=\begin{cases}
\{\infty\}\cup\Bbb Z_p&\text{when }l=-1\\
\{\infty\}\cup\Bbb Z_p-\{0\}&\text{when }l=0\\
\{\infty\}\cup\Bbb Z_p-\{\omega,p-\omega\}&\text{when }l=+1
\end{cases}$$
Define internal law $\boxplus$ in $\mathcal S$ as:
$$x\boxplus y=\begin{cases}
y&\text{when }x=\infty\\
x&\text{when }x\ne\infty\text{ and }y=\infty\\
\infty&\text{when }x\ne\infty\text{ and }y=-x\\
(x+y)^{-1}(x\,y+\hat z)&\text{otherwise}
\end{cases}$$
$(\mathcal S,\boxplus)$ is¹,²,⁶ a finite Abelian group of $m$ elements with neutral $\infty$. With the convention $-\infty=\infty$, the opposite $-x$ of $x$ in group $(\mathcal S,\boxplus)$ is² computed as in $(\Bbb Z_p,+)$ when $x\ne\infty$.
Define scalar multiplication $\boxtimes: \Bbb Z\times\mathcal S\mapsto \mathcal S$ as:
$$k\boxtimes x=\begin{cases}
\infty&\text{when }k=0\\
(-k)\boxtimes(-x)&\text{when }k<0\\
\bigl((k-1)\boxtimes x\bigr)\boxplus x&\text{otherwise}
\end{cases}$$
and for all $x$ in $\mathcal S$ and integers $k$, $k'$, it holds²:
$$\begin{align}
(k\boxtimes x)\,\boxplus\,(k'\boxtimes x)&=(k+k')\boxtimes x&\text{and}\\
k\boxtimes (k'\boxtimes x)&=(k\,k')\boxtimes x&\text{and}\\
k\boxtimes x&=(k\bmod m)\boxtimes x
\end{align}$$
For $x$ in $\{\infty\}\cup\Bbb Z_p$ define $\overline x=\begin{cases}
\infty&\text{when }x=\infty\\
\delta-x&\text{otherwise}\\
\end{cases}$.
Define $\hat{\mathcal S}$ as the set of $\overline x$ for all $x$ in $\mathcal S$. It holds $\hat{\mathcal S}=\mathcal S\iff l=-1$.
A function $f:\hat{\mathcal S}\times\hat{\mathcal S}\mapsto\hat{\mathcal S}$ compatible⁴ with the question's $f$ can²,⁵ now be defined as:
$$f(x,y)=\overline{-(\overline x\boxplus\overline y)}$$
The «restricted commutativity»³ property that $f\bigl(f(x,y),f(y',z)\bigr)=f\bigl(f(x,y'),f(y,z)\bigr)$ is a consequence of the associativity and commutativity of $\boxplus$ , and the fact for all $x$ in $\mathcal S$ it holds $\overline{(\overline x)}=x$.
Define:
$$\begin{align}
\hat s_0&=\overline a\\
\hat s_1&=\overline b\\
\hat s_n&=-(\hat s_{n-2}\boxplus\hat s_{n-1})\text{ when }n>1\\
\hat r_n&=-(\hat s_n\boxplus\hat s_0)\\
\end{align}$$
and for all $n$ it holds²: $\hat s_n=\overline{s_n}$ and $\hat r_n=\overline{r_n}$.
From this we can efficiently compute² integers $u_n$ and $v_n$ in $\Bbb Z_m$ such that:
$$\overline{r_n}=(u_n\boxtimes\overline a)\,\boxplus\,(v_n\boxtimes\overline b)$$
This allows² to solve for $a$ given $b$ and $r_n$, knowing parameters $p$, $z$, $n$. That's significantly easier than a discrete logarithm problem. When $\gcd(u_n,m)=1$, the unique solution is $$a=\overline{(u_n^{-1}\bmod m)\boxtimes(\overline{r_n}\,\boxplus\,(-v_n\boxtimes\overline b))}$$
Notes:
¹ : In particular, $\boxplus$ is² associative!
² : Proof is left as an exercise to the reader.
³ : See document linked in the question.
⁴ : When $l\ne-1$ one of the input of $f$ can be excluded from $S$; assimilate it to $\infty$.
⁵ : $\delta$ and $\hat z$ are chosen such that $\overline{(x+y-1)^{-1}(-x\,y+x+y+z)}=-(\overline x\boxplus\overline y)$.
⁶ : This group has been identified there, if not named.
