___ /\_ \ _____ __ __ __\//\ \ ___ ___ ___ __ /\ '__`\ /'__`\ /\ \/\ \ \ \ \ / __`\ /' __` __`\ /'__`\ \ \ \L\ \/\ \L\.\_\ \ \_/ | \_\ \_/\ \L\ \__/\ \/\ \/\ \/\ __/ \ \ ,__/\ \__/.\_\\ \___/ /\____\ \____/\_\ \_\ \_\ \_\ \____\ \ \ \/ \/__/\/_/ \/__/ \/____/\/___/\/_/\/_/\/_/\/_/\/____/ \ \_\ \/_/

# SUPERBall Bot Strain Gauge System ID Work

OrderReductionTesting

```clc
%clear all
close all
```

## System ID

Here I take a bandwidth limited white-noise input, excite a system with it, use arx to build a model of the input/output.

```%load strain
Ts = 1/1000; % Time-step
na = 3; %order denominator model
nb = 3; %order numerator model
nc = 2; %order numerator noise model (BJ and ARMAX)
nd = 2; %order denominator noise model (BJ)

output1 = output(8100:8.1*10^4,:);

U = output1(:,3);
U1 = output1(:,2);
Y = output1(:,1);
N = length(U);

% U = U(8100:(8.2*10^4));
% U1 = U1(8100:(8.2*10^4));
% Y = Y(8100:(8.2*10^4));

% Remove nans if they exist.
unan=isnan(U);
ynan=isnan(Y);
for i=1:length(U)
if (unan(i))
U(i) = U(i-1);
end
if (ynan(i))
Y(i) = Y(i-1);
end
end
unan=isnan(U);
ynan=isnan(Y);
figure;
plot([unan ynan]);

% Linearizing data, if needed.
K = 0;
V = (U -0*K*(U.^3).*(Y));
Y = (2*3.14/2048)*Y;

plot([V,Y])
``` Estimate models using various system ID methods.

```highorder = 100;

data=iddata(Y,U1,Ts);
data1=iddata(Y,U,Ts);

data = [data, data1];
data = data(:,1,:) %keep only one of the outputs MISO
nk = delayest(data1);

modelARX = arx(data, [highorder [highorder highorder] [nk nk]]);
modelARX3 = arx(data, [3 [3 3] [nk nk]]);
modelARMAX = armax(data,[na [nb nb] nc [nk nk]]);
modelOE = oe(data, [[nb nb] [na na] [nk nk]]);
modelBJ = bj(data, [[nb nb] nc nd [na na] [nk nk]]);
```
```Time domain data set with 72901 samples.
Sampling interval: 0.001

Outputs      Unit (if specified)
y1

Inputs       Unit (if specified)
u1
u2

```

Convert to statespace and returns a decomposition of the model into the observable and unobservable subspaces.

```a = idss(modelARX);
[A,B,C,T,k] = obsvf(a.A,a.B,a.C);
sys = ss(A,B,C,a.D,modelARX.Ts);

a = idss(modelARX3);
[A,B,C,T,k] = obsvf(a.A,a.B,a.C);
sysARX3 = ss(A,B,C,a.D,modelARX3.Ts);

a = idss(modelARMAX);
[A,B,C,T,k] = obsvf(a.A,a.B,a.C);
sysARMAX = ss(A,B,C,a.D,modelARMAX.Ts);

a = idss(modelOE);
[A,B,C,T,k] = obsvf(a.A,a.B,a.C);
sysOE = ss(A,B,C,a.D,modelOE.Ts);

a = idss(modelBJ);
[A,B,C,T,k] = obsvf(a.A,a.B,a.C);
sysBJ = ss(A,B,C,a.D,modelBJ.Ts);
```

## Order Reduction

Computes the Hankel singular values hsv of the dynamic system sys. In state coordinates that equalize the input-to-state and state-to-output energy transfers, the Hankel singular values measure the contribution of each state to the input/output behavior. Hankel singular values are to model order what singular values are to matrix rank.

We then use balred to reduce the model order through use of a balanced realization.

```figure;
hsvd(sys); % Check out what orders contain the most amount of energy.
sys3 = balred(sys, 3)
sys4 = balred(sys, 4)
sys5 = balred(sys, 5)

figure;
bode(sys3);
title('3rd Order BR');
figure;
bode(sys4);
title('4th Order BR');
figure;
bode(sys5);
title('5th Order BR');
figure;
bode(sys);
title('100th Order ARX');
figure;
bode(sysARX3);
title('3rd Oder ARX');
figure;
bode(sysOE);
title('OE Model');
figure;
bode(sysBJ);
title('BJ Model');
figure;
bode(sysARMAX);
title('ARMAX Model');
```
```
a =
x1        x2        x3
x1     0.953    0.2052   0.01946
x2  -0.09065     1.098    0.2412
x3   0.06411   -0.1237    0.8559

b =
u1          u2
x1  -1.856e-06   -0.002165
x2   0.0001254   0.0001829
x3  -7.085e-05   0.0002413

c =
x1        x2        x3
y1   -0.1186  -0.01792   -0.1334

d =
u1          u2
y1  -6.173e-05  -0.0008695

Sampling time (seconds): 0.001
Discrete-time state-space model.

a =
x1         x2         x3         x4
x1     0.7566    -0.1398     0.5856    -0.2578
x2     0.2515     0.9966    -0.1714     0.1034
x3     0.0448    0.06568     0.7407     0.1338
x4  -0.007778   -0.01769    0.06589     0.9607

b =
u1          u2
x1  -0.0001264   -0.005856
x2    4.28e-05   0.0006873
x3   6.846e-05    0.002724
x4  -2.002e-05    -0.00096

c =
x1       x2       x3       x4
y1  0.04525  -0.1192  -0.1168  0.01011

d =
u1         u2
y1  1.887e-05  0.0004616

Sampling time (seconds): 0.001
Discrete-time state-space model.

a =
x1          x2          x3          x4          x5
x1      0.6903     -0.5091     -0.3041      0.4811     -0.1035
x2      0.1544      0.7294     -0.1825      0.8653  -0.0009059
x3      0.2542      0.1104       1.092     -0.2477     0.06413
x4    -0.05073     0.06509     0.07993      0.4394    -0.03276
x5     0.02531     -0.1257     -0.0902      0.3841      0.9944

b =
u1          u2
x1  -7.601e-05   0.0007801
x2  -0.0001636   -0.004855
x3   3.768e-05    0.001123
x4   0.0001206    0.003997
x5  -7.311e-05   -0.001823

c =
x1        x2        x3        x4        x5
y1  -0.04945  -0.03785    0.1087    0.0277    0.0246

d =
u1          u2
y1  -4.194e-06  -0.0001146

Sampling time (seconds): 0.001
Discrete-time state-space model.
```         ## Simulation Of Models To Real Data

Using lsim to look at the differences between the actual data, and what our models will predict the output to be.

```Ts = 1/1000;

V = [U1, U];
T = Ts*(0:length(V)-1);

y3 = lsim(sys3, V);
y4 = lsim(sys4, V);
y5 = lsim(sys5, V);
ysys = lsim(sys, V);
yARX3 = lsim(sysARX3, V);
yARMAX = lsim(sysARMAX, V);
yOE = lsim(sysOE, V);
yBJ = lsim(sysBJ, V);

figure
plot([Y y3 y4 y5 ysys])
legend('Real','3rd Order BR','4th Order BR','5th Order BR', 'BR');
title('Simulation of Order-Reduction and Balanced Realization');

plot([Y yARX3 yARMAX yOE yBJ])
legend('Real','ARX','ARMAX','OE', 'BJ');
title('Simulation w/o Balanced Realization');

figure
plot([Y-y3 Y-y4 Y-y5 Y-ysys])
title('Error');
legend('3rd Order BR','4th Order BR','5th Order BR', 'BR');

figure
plot([Y-yARX3 Y-yARMAX Y-yOE Y-yBJ])
title('Error');
legend('ARX','ARMAX','OE', 'BJ');
```   ## Non-Linear Model

I'm plotting a non-linear ARX model of the system here with 5 regressors and a wavelet network with 57 units which was modeled to the same iddata object. First we'll look at the non-linearity of the wavenet and then simulate the output strain with both the non-linear model and a linearized verion of the nonlinear autoregressive approximation using linapp. Linapp does /not/ use a tangent linearization approach and should work better for larger ranges of input 'u'.

```figure;
plot(nlarx1);
figure;
a = linapp(nlarx1,V);
a = idss(a);
a = ss(a.A, a.B, a.C, a.D, 1/1000);
plot([Y sim(nlarx1, V) lsim(a, V)]);
legend('strain', 'non-linear model', 'linearized nlarx');
```    