Contents

Example: Example from Reynders

close all; clear; clc;

LTI model

% LTI system matrices
h = 1; % Sample time
n = 4; % number of states
A22 = [0.5321 0.8349; -0.8349 0.5231];
A33 = [-0.9165 0.1313; -0.1313 -0.9615];
B2 = [0.0012; 0.0007];
C2 = [-0.0267 2.6723];
C3 = [-0.0792 7.9236];
D = -0.008;
K2 = [0.0060; 0.0009];
K3 = [0.0242; -0.1218];
A = [A22 K2*C3; zeros(2) A33];
B = [B2; zeros(2,1)];
K = [K2; K3];
C = [C2 C3];

% open-loop system
OL = ss(A,[B K],C,[D eye(1)],1);

Open-loop identification experiment

Simulation of the model in open loop

% input signals
N = 10000; % number of samples
t = (0:N-1)';   % time samples
u = 1000*randn(N,1); % excitation signal

% noise
e = randn(N,1); % noise signal

% simulation
y0 = lsim(OL,[u zeros(N,1)],t);
y = lsim(OL,[u e],t);
disp('Signal to noise ratio (SNR) (open-loop)')
snr(y,y0)

% Defining a number of constants
p = 15;     % past window size
f = 15;     % future window size
[u,Du,y,Dy] = sigscale(u,y); % signal scaling

% PBSID-opt
[S,x] = dordvarx(u,y,f,p,'tikh','gcv');
figure, semilogy(S,'x');
title('Singular values')
x = dmodx(x,n);
[Ai,Bi,Ci,Di,Ki] = dx2abcdk(x,u,y,f,p);
Dat = iddata(y',u',h);
Mi = abcdk2idss(Dat,Ai,Bi,Ci,Di,Ki);
Mi

% Variance-accounted-for (by Kalman filter)
yest = predict(Mi,Dat);
x0 = findstates(Mi,Dat);
disp('VAF of identified system')
vaf(y,yest.y)

% PBSID-opt (greybox)
[S,x,xd,xs] = dordvarx_grey(u,y,f,p,'tikh','gcv');
figure, semilogy(S,'x');
title('Singular values');
[x,R1,R2,R3] = dmodx_grey(x,xd,xs,n,2,4);
figure, plot(R1,'*');
title('Canonical values between full state and deterministic only');
figure, plot(R2,'*');
title('Canonical values between full state and stochastic only');
figure, plot(R3,'*');
title('Canonical values between deterministic and stochastic only');
[Ap,Bp,Cp,Dp,Kp] = dx2abcdk(x,u,y,f,p);
Dat = iddata(y',u',h);
Mp = abcdk2idss(Dat,Ap,Bp,Cp,Dp,Kp);
Mp

% Variance-accounted-for (by Kalman filter)
yest = predict(Mp,Dat);
x0 = findstates(Mp,Dat);
disp('VAF of identified system')
vaf(y,yest.y)
Signal to noise ratio (SNR) (open-loop)

ans =

   17.4062

State-space model:  x(t+Ts) = A x(t) + B u(t) + K e(t)
                       y(t) = C x(t) + D u(t) + e(t)
 
A = 
                        x1           x2           x3           x4
           x1      0.51518     -0.83698     0.010729    -0.050514
           x2      0.83289      0.53982    0.0089433     -0.02801
           x3    0.0012153  -0.00039965     -0.97023     -0.26433
           x4   -0.0016516   0.00088268      0.13543     -0.96658
 
 
B = 
                        u1
           x1     -0.02301
           x2    -0.014251
           x3    0.0014799
           x4   0.00079533
 
 
C = 
                        x1           x2           x3           x4
           y1      -6.3747        3.427       2.8118      0.72174
 
 
D = 
                        u1
           y1     -0.42905
 
 
K = 
                        y1
           x1    -0.023525
           x2    -0.012381
           x3    -0.079619
           x4    -0.049436
 
 
x(0) = 
                          
           x1   0.00060748
           x2  -0.00096233
           x3     0.017172
           x4     0.022878
 
Estimated using PBSIDopt from data set Dat
Loss function 0.0065966 and FPE 0.00663486
Sampling interval: 1                      
                                          
VAF of identified system

ans =

   99.3403

Warning: Y is not full rank. 
State-space model:  x(t+Ts) = A x(t) + B u(t) + K e(t)
                       y(t) = C x(t) + D u(t) + e(t)
 
A = 
                        x1           x2           x3           x4
           x1      0.51663     -0.83786    -0.043114     -0.05224
           x2      0.83268      0.53981    -0.019834    -0.037556
           x3     0.020949   -0.0094696     -0.86688     -0.23178
           x4      0.01322   -0.0068769      0.20412     -0.83112
 
 
B = 
                        u1
           x1    -0.021548
           x2    -0.014234
           x3     0.016234
           x4    0.0091957
 
 
C = 
                        x1           x2           x3           x4
           y1      -6.3855       3.4277        1.589      0.17497
 
 
D = 
                        u1
           y1     -0.42804
 
 
K = 
                        y1
           x1    -0.027462
           x2     -0.01823
           x3    -0.042363
           x4   -0.0092331
 
 
x(0) = 
                          
           x1    0.0030623
           x2  -0.00065717
           x3     0.011354
           x4      0.10385
 
Estimated using PBSIDopt from data set Dat
Loss function 0.0175901 and FPE 0.0176921 
Sampling interval: 1                      
                                          
VAF of identified system

ans =

   98.2408

Identification results

% Plot eigenvalues
figure
hold on
title('poles of identified LTI systems')
[cx,cy] = pol2cart(linspace(0,2*pi),ones(1,100));
plot(cx,cy,'k');
plot(real(eig(OL.a)),imag(eig(OL.a)),'k+','LineWidth',0.1,'MarkerEdgeColor','k', 'MarkerFaceColor','k', 'MarkerSize',10);
plot(real(eig(Mi.a)),imag(eig(Mi.a)),'rx');
plot(real(eig(Mp.a)),imag(eig(Mp.a)),'gx');
axis([-1 1 -1 1]);
axis square
legend('STABBND','TRUE','PBSID-opt','GREY','Location','East');
hold off

% Bodediagram (open loop)
OLi = ss(Mi);
OLp = ss(Mp);
OLi = Dy*OLi*inv([Du 0; 0 1]);
OLp = Dy*OLp*inv([Du 0; 0 1]);
figure, bodemag(OL(1,1),'k',OLi(1,1),'r',OLp(1,1),'g');
figure, bodemag(OL(1,2),'k',OLi(1,2),'r',OLp(1,2),'g');