Closed-loop system identification

Hi All

Just curious if anyone has any experience on this. Supposed I have a system which is perturbed through "disturbance", normally a random signal, and I can measure u and y. Is it theoretically possible to identify G1(s) or other equivalent model (linear or nonlinear)? what should I do?

I have tried several methods available in system identification toolbox and consider directly u and y as if it is from an open-loop system, but so far the results are not good. Just think that it may have something to do with the fact that the system operates in feedback (?)

Any help is appreciated

 Thanks you

Answer 1

[Aquatris]

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;

   G1(s) = T(s)/H(s)

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.

% mass-spring parameters
k = 10; % [N/m]
m = 2;  % [kg]

% system equation of motion
%    xd = Ax + Bu
%    y  = Cx + Du

A = [0    1;
     -k/m 0];
B = [0
     1/m];
C = [1 0];
D = 0;
G1 = ss(A,B,C,D); % marginally stable open loop
                   % input is force, output is position

% controller 
Kp = 10; Ki = 8; Kd = 5; % parameters to be tuned
s = tf('s');
G2 = Kp + Ki/s + Kd*s;   

% closed loop systems
H = feedback(1,G1*G2); % output for u 
T = feedback(G1,G2);   % output for y

% define time
dt = 0.001;
t = 0:dt:15-dt;
L = length(t);
f = 1/dt*(0:(L/2))/L; % frequency for fft

% impulse excitation signal design
e = zeros(length(t),1);
e(11:15) = 1;

% fft of excitation signal
Y = fft(e);
P2 = (Y/L);
P1 = P2(1:L/2+1,:);
P1(2:end-1,:) = 2*P1(2:end-1,:);
P1_input = P1;

% simulate response of the system
[u,t2] = lsim(H,e,t) ; % u measurement
[y,t2] = lsim(T,e,t) ; % y measurement

% fft of y signal
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; % FRD of T

% fft of u signal
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; % FRD of H

% calcualte FRD for G1
G1_frd = T_frd./H_frd;

%compare real system to calcualted FRD
[mag,phase] = bode(G1,f*2*pi);
loglog(f,squeeze(mag),f ,abs(G1_frd),'r--')
xlabel('frequency [Hz]')
ylabel('Magnitude [abs]')

[Husni Rois Ali]

Dear Aquatris

Thank you so much for your super clear explanation and also for the example. Following this, I have another question: If the assumption of being able to add signal e(s) to u(s) is removed, can we still be able to identify G1(s)? In other words the only signal I have is u and y, both due to the "disturbance". If not, what's the reason behind this?

Sorry if my question is a bit strange, but it comes from practical point of view. In what I'm doing now, G1 and G2 are two large power systems, u and y can be voltage and current at a node, and the "disturbance" can be, for instance, lots of customers who are constantly switching on and off their loads at home. So what often practically available is only u and y, and we cannot mess with this.

Look forward for your follow-up

Thank you

Ali

[Aquatris]

IMHO, I do not think it is impossible but since you have no control over u, it is highly unlikely that you can get an accurate model using whatever u and y measurement you have. In principle, whether it is in feedback or not (this is different than adding an external signal);

 y(s) = G1(s)*u(s).

However, if u(s) does not have values in a wide range of frequencies, the model you identify with your measurements will be a local fit. An extreme case is, you measure u(t) = 1 and y(t) = 3. Then using these measurements, a model of G1(s) can be obtained as G1(s) = 3. However, this is valid only if u(t) is constant and probably will give a wrong estimation of y if a sine wave is applied.

I would suggest you use as much measurements for u and y as possible and hope u actually excites the G1 dynamics for identification purposes.

[Husni Rois Ali]

Thank you so much for your help