Discrete time - Z domain - upsampling - numerical issue

I'm working with a Z domain transfer function. The TF is a simple PI filter that I oversampled (W(z) -> Q(z^N) ) followed by a single pole low pass filter.

As required I included a ZOH.

Fc = 1e9; 
Fs = 10e9;

Tc = 1/Fc;
Ts = 1/Fs;

alpha = 1/15;     % Integral
beta  = 1;        % Proportional

R_load = 1e3;
Tau = R_load*10e-12;

N = round(Fs/Fc);  % N = 10

z = tf('z',Ts);

eta = Tau/Ts;

PI = (beta + alpha*(z^N)/(z^N - 1));
ZOH = 1/N*(1-z^(-N))/(1-z^(-1));
Z = R_load/(1 + eta*(z-1));

L = PI*ZOH*Z;

bode(L);

What I obtain is this:

But what I supposed to obtain was something like that, without any resonance peak:

In fact, if I don't upsample and I run the system at the same frequency (good approximation at low frequency):

Fc = 1e9; 
Fs = 1e9;

Tc = 1/Fc;
Ts = 1/Fs;

alpha = 1/15;     % Integral
beta  = 1;        % Proportional

R_load = 1e3;
Tau = R_load*10e-12;

N = round(Fs/Fc);  % N = 1

z = tf('z',Ts);

eta = Tau/Ts;

PI = (beta + alpha*(z^N)/(z^N - 1));
ZOH = 1/N*(1-z^(-N))/(1-z^(-1));
Z = R_load/(1 + eta*(z-1));

L = PI*ZOH*Z;

bode(L);

Where ZOH = 1/N*(1-z^(-N))/(1-z^(-1)) should be equal to 1 (at every frequency in band).

That is what I obtain:

That is what I obtain if I force ZOH = 1, and what I supposed to obtain also in the oversampled case:

I think it is a numerical issue, but how can avoid that?

[Birdman]

I suggest you to redefine the ZOH or use Tustin instead. Actually it would be better if you express this situation in a Simulink model because this way it is really hard to understand what is going on.