Assuming you can add excitation signal to u, and with large enough excitation signal magnitude, you can identify G1 without knowing G2. Here is how (assuming SISO but can be extended to MIMO);
You design an excitation signal e (sine sweep or PRBS) and add it to u. The transfer function from e to u and e to y is;
u(s)/e(s) = 1/(1+G1(s)*G2(s)) = H(s)
y(s)/e(s) = G1(s)/(1+G1(s)*G2(s)) = T(s)
Frequency response data for these two transfer functions, H(s) and T(s) can be calculated using Fourier transform of the signals since they are measured. From here you can simply do;
This would give you an FRD. Then, you can fit a model to the FRD.
This process requires careful selection of the excitation signal to allow necessary inversions. Also e(s) should be designed such that the system response is dominated by e(s) not the disturbance.
Here is an example code. G1 is a spring-mass system. G2 is PID controller. There are no disturbances or noise in measurements in the code, but you can add and play around to see their effects.
k = 10;
m = 2;
A = [0 1;
-k/m 0];
B = [0
1/m];
C = [1 0];
D = 0;
G1 = ss(A,B,C,D);
Kp = 10; Ki = 8; Kd = 5;
s = tf('s');
G2 = Kp + Ki/s + Kd*s;
H = feedback(1,G1*G2);
T = feedback(G1,G2);
dt = 0.001;
t = 0:dt:15-dt;
L = length(t);
f = 1/dt*(0:(L/2))/L;
e = zeros(length(t),1);
e(11:15) = 1;
Y = fft(e);
P2 = (Y/L);
P1 = P2(1:L/2+1,:);
P1(2:end-1,:) = 2*P1(2:end-1,:);
P1_input = P1;
[u,t2] = lsim(H,e,t) ;
[y,t2] = lsim(T,e,t) ;
Y = fft(y);
P2 = (Y/L);
P1 = P2(1:L/2+1,:);
P1(2:end-1,:) = 2*P1(2:end-1,:);
P1_y = P1;
T_frd = P1_y./P1_input;
Y = fft(u);
P2 = (Y/L);
P1 = P2(1:L/2+1,:);
P1(2:end-1,:) = 2*P1(2:end-1,:);
P1_u = P1;
H_frd = P1_u./P1_input;
G1_frd = T_frd./H_frd;
[mag,phase] = bode(G1,f*2*pi);
loglog(f,squeeze(mag),f ,abs(G1_frd),'r--')
xlabel('frequency [Hz]')
ylabel('Magnitude [abs]')