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?
3 个评论
Taiwo Ajayi
2018-12-20
Please how / where can one get this CLOSID tool?
Walter Roberson
2018-12-20
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
2018-12-25
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.
采纳的回答
Rajiv Singh
2011-7-19
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.
2 个评论
Jason Nicholson
2014-6-15
Thanks for the information. This is still relevant even though this was posted 3 years ago.
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
更多回答（3 个）
Rajiv Singh
2011-7-17
2 个投票
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.
Anirudhh Ravi
2021-11-13
0 个投票
Hello all,
Can anyone please share the downloaded CLOSID toolbox? I am having a hard time finding it online.
Many thanks!
antonio
2022-7-9
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)) ];
类别
在 帮助中心 和 File Exchange 中查找有关 Data Preparation Basics 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
