Closed Loop system identification toolbox

Hi all,

Is there any MATLAB Toolbox available for Closed Loop System Identification because MATLAB's current system identification toolkit is for identification of open loop systems?

I already knew one tool named CLOSID but it is not as flexible as MATLAB's sys ID toolkit and i am facing some problems in using it.

Is there anyway of using MATLAB's system identification toolkit for CLOSED LOOP SYSTEMS?

[Taiwo Ajayi]

Please how / where can one get this CLOSID tool?

[Walter Roberson]

Publication is at https://pure.tudelft.nl/portal/en/publications/closid--a-matlab-toolbox-for-closedloop-system-identification(ac8b9f82-151b-47b0-97e1-f1d65eac0b33)/export.html

At one point the file was available from an ftp site that now no longer exists. Later it was placed on http://www.tn.tudelft.nl/mmr/downloads/software/closid/closid30.zip which now does not exist either; it looks like the file was removed by 2007.

Ah, I found a copy! http://www.pudn.com/Download/item/id/818840.html . However, this is at pudn, so it is very much a Trust At Your Own Risk.

[Taiwo Ajayi]

Thanks. I took a chance but it looks like they are missing some of the files and there are backward version compatibility issues. Now, I'm curious to know on which MATLAB version the original poster used the CLOSID tool.

Answer 1

[Rajiv Singh]

A model in System Identification Toolbox represents the equation y = Gu + He where G is the transfer function between y and u and H is the transfer function between y and e. G is called the "measured" component and H the "noise component" of the overall transfer function between y and [u, e].

In state space case, the equation is:

x(k+1) = A x(k) + B u(k) + K e(k)
  y(k) = C x(k) + D u(k) + e(k)

H in this case is:

x(k+1) = A x(k) + K e(k)
  yh(k) = C x(k) + e(k)

When you estimate a state space model using PEM, you estimate the values of one or more of A, B, C, D and K. In order to have a nontrivial noise component in your model, make sure that K is not fixed at zero. You do this by using "DisturbanceModel"/'estimate' PV pair in the estimation command or by setting the DisturbanceModel property of an existing IDSS model to 'estimate':

model = pem(data, NX, 'dist', 'estimate')
model2 = idss(drss(NX)) % or some other template you created using IDSS command
model2.DisturbanceMOdel = 'estimate'
model2 = pem(data, model2, 'focus', 'prediction')

Now, the claim is that estimating K along with A, B, C, D helps identify the plant even though it was operating under closed loop. Why this is the case is best answered in the reference I mentioned before. Here is a simple description: if you are measuring signals right at the I/O ports of a plant, you are indeed measuring the effect of changes in input at the plant's output. However the input signal contains the portion of the past output (because of feedback). This by itself is not a problem. But you are not only feeding back the previous outputs but also any disturbances that might have affected those previous output values. This disturbance is thus correlated with the system inputs. Adding a sufficiently flexible "H" helps you separate out the effect of disturbances on the plant output.

In other model structures, H can be suitably (even independently) defined. For example, in a Box Jenkins polynomial model, you have separate transfer functions for H and G whose orders your can pick independently (see IDPOLY, BJ; see also ARMAX).

There are persistence of excitation and other requirements for getting good results from closed loop data. Perhaps a better reference than I mentioned before is:

   System Identification — Theory For the User, Lennart Ljung, Section 13.4-13.5, 2nd ed, PTR Prentice Hall, Upper Saddle River, N.J., 1999.

[Jason Nicholson]

Thanks for the information. This is still relevant even though this was posted 3 years ago.

[Mohammad]

Hi Rajiv,

Could you please tell me where that disturbance component (signal) is exactly located in the closed loop system?

BTW, if I choose discrete-time option, then the state-space form will be x(t+Ts)=..., if I set Ts=1 is that equal to x[k+1]?

Bests, Mohammad

Answer 2

[Rajiv Singh]

You need data that comes from measurements of a system operating in closed loop. How are your measuring these input/output signals? I mean, are the measurements made across the actual plant to be modeled or are they made at some "outside the loop" reference points? If it is the former, open loop identification using a flexible noise component could still yield a good model of the plant. If it is the latter, you need some way of extracting out the plant model from the overall closed loop model. This reference might be useful:

U. Forssell and L. Ljung. Closed-loop identification revisited. Automatica, 35(7): 1215-1241, Jul 1999.

[hasan]

Hi Rajiv,

Thanks for your reply.

In my case, i am measuring inputs/outputs across the Actual plant. Can you please elaborate your term "Flexible noise component"? i have also modelled my system using matlab's system identification by using 5th order state space model. If i add disturbance to it then will this be good enough for closed loop identification? If yes, then how is it possible that open loop system resembles the closed loop system just by adding noise to it?

Regards,

hasan

[Siamak]

Hasan Can you tell me where I can find CLOSEID plugin? I have a closed loop system and I would like to find transfer function Thanks Siamak

Answer 3

[Anirudhh Ravi]

Hello all,

Can anyone please share the downloaded CLOSID toolbox? I am having a hard time finding it online.

Many thanks!

Answer 4

[antonio]

Hi,

I need help with code for closed loop identification. In my code there is a part "frek_est" line 57 and 58. I don't know which function is better for frequency estimation. Can someone help me with that or send example code for closed loop identification.

This is direct procedure for closed loop identification.

Thanks!!

%process
num=100;
den=[5 6 81 16];
sys=tf(num,den);%process

%regulator
Kr=0.1;
TI=0.5;
Gr=tf(Kr*[0.5 1],[0.5 0]);

Gz=feedback(Gr*sys,1); %closed loop

b=0.2; %variance for sum
Ts=0.4; %time
%model discretization
G=c2d(sys,Ts);
Grd=c2d(Gr,Ts,'Tustin'); %Tustin dicvretization
Gzd=feedback(Grd*G,1); 

Gzud=G/(1+Grd*G); %transfer function according to reference
Gzvd=1/(1+Grd*G); %transfer function output according to v
Gzvud=-Grd/(1+Grd*G); %transfer function according to u
Gzuud=1/(1+Grd*G); %transfer function according to r2

T=1/(1+Grd*G); %sensitivity function
%PRBS signal 
Nr=255;  %period

M=4; %4 repetition
Mu=1; 

%simulation system in closed loop
uprbs=idinput([Nr,1,M],'prbs'); %signal r2
yprbs1=lsim(Gzud,uprbs); %signal y
uprbs1=lsim(Gzuud,uprbs); %signal u
yprbs=yprbs1+b*randn(size(uprbs)); %noisy signal y
uprbs=uprbs1+b*randn(size(uprbs)); %noisy signal u

%open loop simulation
yp_ls=lsim(G,uprbs)+b*randn(size(uprbs)); %signal y
up_ls=uprbs; %signal u

%extraction data
u=uprbs((M-Mu)*Nr+1:end);
y=yprbs((M-Mu)*Nr+1:end);

u_ls=up_ls((M-Mu)*Nr+1:end); % u and y for ls nethode
y_ls=yp_ls((M-Mu)*Nr+1:end);

Ndata=length(u);
N=Ndata;
z=iddata(y,u,Ts);
z_ls=iddata(y_ls,u_ls,Ts);

Model=oe(z,[3,3,1]);
omega=(2*pi/N)*[0:N-1];
[G_cl,Ycl,Ucl,w, idx]=frek_est(u,y,Ts,1);
[G_ol,Yol,Uol,w_ol, idx1]=frek_est(u_ls,y_ls,Ts,1);
G_cl=etfe(z,u,y);
idx=5;

Gfresp=squeeze(freqresp(G,omega/Ts));
Tfresp=squeeze(freqresp(T,omega/Ts));

%plotting
figure
subplot(211)
semilogx(omega(idx),20*log10(abs(G_cl(idx))),'d-','LineWidth',1);
hold on
semilogx(omega(idx),20*log10(abs(G_ol(idx))),'o-','LineWidth',1);
semilogx(omega(idx),20*log10(abs(Gfresp(idx))),'r-','LineWidth',1.5);
xlabel('Normalna frekvencija','Interpreter','Latex','FontSize', 14);
ylabel('Amplituda','Interpreter','Latex','FontSize',14)
grid on
subplot(212)
semilogx(omega(idx),angle(G_cl(idx)),'d-','LineWidth',1);
hold on
semilogx(omega(idx),angle(G_ol(idx)),'o-','LineWidth',1);
semilogx(omega(idx),angle(Gfresp(idx)),'r-','LineWidth',1.5);
xlabel('Normalna frekvencija','Interpreter','Latex','FontSize', 14)
ylabel('Amplituda','Interpreter','Latex','FontSize',14)
grid on

mod_ls=oe(z_ls,[3,3,1]);

figure,
step(Model),hold on,step(mod_ls),step(G,'k'),legend('Model','Model_LS','Proces'),grid on;
podaci=[w,abs(Ucl) abs(Ycl) 20*log10(abs(G_cl)) 180/pi*angle(G_cl) abs(Uol) abs(Yol) 20*long10(abs(G_ol)) ];