This is sketch of how to solve your problem. I still don't understand the boundary conditions/ initial conditions as stated. I have set $a,b,d$ to one, and have made generic boundary conditions.
Let $\Omega \subset \mathbb{R}^d$ be a bounded open domain. We are interested in solving the following coupled systems of PDE. This is pretty hairy.
$$
\begin{align*}
0 &= (\partial_{tt} - \nabla^2 + \partial_t) f - g&\quad \mathrm{for}\ (x,t)\in \Omega \times \mathbb{R}^+ \\
0 &= (\partial_t-\nabla^2 + 1) g - \partial_t f&\quad \mathrm{for} (x,t) \in \Omega \times \mathbb{R}^+ \\
f &= 0 &\quad \mathrm{for} (x,t)\in \Omega\times \{0\}\\
f_t &= 0&\quad \mathrm{for} (x,t)\in \Omega\times \{0\}\\
g &= 0& \quad \mathrm{for} (x,t)\in \Omega\times \{0\}\\
f_t &= s_1 &\quad \mathrm{for} (x,t)\in \partial\Omega \times \mathbb{R}^+\\
g &= s_2 &\quad \mathrm{for} (x,t)\in \partial\Omega \times \mathbb{R}^+
\end{align*}
$$
The fist two equations are the PDE. The second three are initial conditions, and the last two are boundary conditions.
Step #1 Homogenize the PDE's boundary conditions. We need to turn the boundary conditions into "forcing terms"
Let $f_{bc}$ be the solution to the following PDE
\begin{align*}
-\nabla^2 f_{bc} &= 0 &\quad \mathrm{for} (x,t)\in \Omega \times \mathbb{R}^+\\
\frac{\mathrm{d}}{\mathrm{d}t} f_{bc} &= s_1 & \quad \mathrm{for} (x,t)\in \partial\Omega\times \mathbb{R}^+
\end{align*}
With simple boundary conditions and simple geometry this $f_{bc}$ might just be a constant in space and vary in time.
Let $g_{bc}$ be the solution of a similar PDE
\begin{align*}
-\nabla^2 g_{bc} &= 0&\quad \mathrm{for} (x,t)\in \Omega \times \mathbb{R}^+\\
g_{bc} &= s_2& \quad \mathrm{for} (x,t)\in \partial\Omega \times \mathbb{R}^+.
\end{align*}
Now we have a set of new PDE's.
$$
\begin{align*}
-(\partial_{tt} + \partial_t )f_{bc} + g_{bc}&= (\partial_{tt} - \nabla^2 + \partial_t) (f-f_{bc}) -(g-g_{bc})&\quad \mathrm{for}\ (x,t)\in \Omega \times \mathbb{R}^+ \\
-(\partial_t +1 )g_{bc} + (f_{bc})_t&= (\partial_t-\nabla^2 + 1) (g-g_{bc}) - \partial_t (f-f_{bc})&\quad \mathrm{for} (x,t) \in \Omega \times \mathbb{R}^+ \\
f-f_{bc} &= -f_{bc} &\quad \mathrm{for} (x,t)\in \Omega\times \{0\}\\
f_t-(f_{bc})_t &= -(f_{bc})_t&\quad \mathrm{for} (x,t)\in \Omega\times \{0\}\\
g-g_{bc} &= -g_{bc}& \quad \mathrm{for} (x,t)\in \Omega\times \{0\}\\
f_t-(f_{bc})_t &= 0 &\quad \mathrm{for} (x,t)\in \partial\Omega \times \mathbb{R}^+\\
g-g_{bc} &= 0 &\quad \mathrm{for} (x,t)\in \partial\Omega \times \mathbb{R}^+
\end{align*}
$$
Notice what we have done, we have changed a PDE with with inhomogeneous boundary conditions and no "forcing functions" into one with homogeneous boundary conditions with forcing functions.
Step #2.
Denote $\{\Phi_i\}_{i=1}^\infty$ as a sequence of eigenfunctions for the associated homogeneous Poisson equation, ordered by increasing eigenvalue. We mean to say that $\int_{\Omega} \Phi_i \Phi_j \mathrm{d} x = \delta_{i,j}$ and $\int_{\Omega} (\nabla \Phi_i) \cdot (\nabla \Phi_j) \mathrm{d}x = \delta_{i,j} \lambda_i^2$. (For the cylindrical case we have $\Phi_i$ are the Bessel functions.)
Now, these eigenfunctions are complete so we can expand $(f-f_{bc})$ as
$$
(f-f_{bc})(x,t) = \sum_{i=1}^\infty \tilde{f}_i(t) \Phi_i(x)
$$
and $$(g-g_bc)(x,t) = \sum_{i=1}^\infty \tilde{g}_i(t) \Phi_i(x)$$
Now multiplying the PDE by $\Phi_j$ integrating by parts the Laplacian term we get
$$
\underbrace{\int_{\Omega} \Phi_i ((-\partial_{tt} - \partial_t )f_{bc} + g_{bc} )\, \mathrm{d} x}_{\alpha_i(t)} = (\tilde{f}_i)_{tt} + (\tilde{f}_i)_t + \lambda_i^2 \tilde{f}_i - \tilde{g}_i
$$
and
$$
\underbrace{\int_{\Omega} \Phi_i ((-\partial_t+1) g_{bc} + \partial_t f_{bc})\mathrm{d}x }_{\beta_i(t)} = (\tilde{g}_i)_t+\lambda_i^2 \tilde{g}_i + \tilde{g}_i - (\tilde{f}_i)_t
$$
The beauty is that, at this point $\alpha_i$ and $\beta_i$ are all known, and $\tilde{f}_i$ only couples with $\tilde{g}_i$.
Step #3. Solve ODE.
We write a state-space representation of the ODE's.
$$
\begin{align*}
\frac{\mathrm{d}}{\mathrm{d} t}\left\lbrack
\begin{array}{c}
\tilde{f}_i\\
(\tilde{f}_i)_t\\
\tilde{g}_i
\end{array}\right\rbrack + \left\lbrack
\begin{array}{ccc}
0&-1&0\\
\lambda_i^2& 1& -1\\
0& -1 & \lambda_i^2+1
\end{array}
\right\rbrack
\left\lbrack
\begin{array}{c}
\tilde{f}_i\\
(\tilde{f}_i)_t\\
\tilde{g}_i
\end{array}
\right\rbrack
= \left\lbrack
\begin{array}{c}
0\\
\alpha_i\\
\beta_i\\
\end{array}
\right\rbrack
\end{align*}
$$
Solve this system how ever you want. This solves $f-f_{bc}$ and $g-g_{bc}$ so to get $f$ and $g$ you just add the boundary functions. Do you see how to handle the initial conditions? I can explain any steps that I skipped or ones that are unclear.
