It's not a complete answer but I hope that it could be of some help.
You could rewrite the first system as
$$
\begin{cases}
P(n) = K_P E(n) \\
I(n) = I(n-1) + \frac{K_P}{T_I} E(n) \Delta t \\
D(n) = K_P T_D \frac{E(n) - E(n-1)}{\Delta t}
\end{cases}
$$
Where \$E(n) = G(n) - target(n)\$ and \$\Delta t\$ is your sampling interval. Note that \$T_D\$ and \$T_I\$ are not defined as gains. \$K_I = \frac{K_P}{T_I}\$ and \$K_D = K_P T_I\$ are respectively the integral gain and the derivative gain.
Now you can rewrite the system as a single function of the error.
$$
PID(n) = P(n) + I(n) + D(n)
$$
$$
I(n-1) = PID(n-1) - P(n-1) - D(n-1) \\
= PID(n-1) - K_P E(n-1) - K_P T_D \frac{E(n-1) - E(n-2)}{\Delta t}
$$
$$
PID(n) = K_P E(n) + PID(n-1) - K_P E(n-1) - K_P T_D \frac{E(n-1) - E(n-2)}{\Delta t} + \frac{K_P}{T_I} E(n) \Delta t + K_P T_D \frac{E(n) - E(n-1)}{\Delta t} \\
= PID(n-1) + K_P \left(\left(1 + \frac{\Delta t}{T_I} + \frac{T_D}{\Delta t} \right)E(n) - \left(1 + 2\frac{T_D}{\Delta t} \right)E(n-1) + \frac{T_D}{\Delta t} E(n-2) \right)
$$
The second one is a bit more complex to rewrite as a single equation but you can do it in a similar way. The result should be
$$
R(n) = K_1 R(n-1) - (\gamma K_0 + K_2) R(n-2) + (1+\gamma) (PID(n) - K_1 PID(n-1) + K_2 PID(n-2))
$$
Now you only need to substitute the equation of the PID in order to obtain the equation of the regulator as function of the error.
The reason why transfer functions work so well for linear time-invariant (LTI) systems (and don't for non-linear systems) is that they are the Laplace transform (or, in discrete time) the Z-transform of the system's impulse response, which completely characterizes the behavior of such systems. I.e., the impulse response, or, equivalently, the transfer function is all you need to know to compute the response of the system to any input signal. The reason for this is that any input signal can be written as an integral (or a sum in discrete time) of scaled and shifted delta impulses, and due to linearity and time invariance, the response to this integral/sum equals the integral/sum of scaled and shifted impulse responses:
$$x(t)=\int_{-\infty}^{\infty}x(\tau)\delta(t-\tau)d\tau\quad\Longrightarrow\quad
y(t)=\int_{-\infty}^{\infty}x(\tau)h(t-\tau)d\tau\tag{1}$$
where \$x(t)\$ is the input signal, \$y(t)\$ is the system's response, and \$h(t)\$ is the impulse response. Equation (1) is the convolution integral, which completely describes the system's behavior. In the transform domain (i.e., Laplace transform, Fourier transform, or Z-transform), convolution becomes multiplication.
For non-linear systems you can of course compute or measure the system's response to an impulse, but this function does not tell you anything about the system's response to other input signals, or even to a scaled impulse. This is why the impulse response and its transform do not have any significance for non-linear systems.
Best Answer
It is basic block diagram algebra.
First, you write out the algebraic equations and solve for the unknowns \$x1\$ and \$x2\$.
$$\text{x1}=\frac{0.05 (\text{Ua}-0.1 \text{x2})}{0.01 s+1}$$ $$\text{x2}=\frac{-\text{Mx}+\text{x1}-2 Y}{0.5 s+1}$$
Next you write the output equation as \$Y= \frac{1}{s} x2\$ and solve for \$Y\$ in terms of the two inputs \$Ua\$ and \$Mx\$. This will give you the transfer function model from the two inputs to the output \$Y\$.
Then obtain a state-space realization.
Since \$Mx\$ is a disturbance, you want to remove the column of the input matrix that corresponds to that before computing the controllability matrix.