No need to use a PID. You can have a much simpler hysteresis control.
Thermal capacity of water is around $$4184 J/(kg·K)$$
In my case, the heater´s power was actually $$2kW = 2kJ/s$$
This means, it heats 1kg of water for about 0.5 K per second. In my case, 2kg water usually.
The bigger problem is precise measurement of the water temperature because of thermal convection. You want a circulator in there to keep thermal differences minimal.
Once you have established good measurement, you can do a simple hysteresis control, in my case I switched the heater on with a mechanical relays for 1s to have a 0.25 K rise.
Temp reading error is going to be around 0.5 K anyways, so don't bother with too much of a regulation.
For a purely resistive load, you will be fine with a simple relays, which also does the electrical isolation for you.
If you want to go for electronic switches, an optotriac will be just fine.
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.
Best Answer
It looks to me like OUT1 is a SPDT relay contact and OUT2 is a SPST relay output.
SPDT means single pole double throw i.e. a changeover relay contact
SPST means single pole signle throw i.e. a normally open relay contact.
Both relay contacts can be presumed to be able to work with a 250VAC AC connected load at up to 3 amps but, without the exact model number and technical manual this is still a bit of an assumption.