%% Kittiphan Techakittiroj (27 Jan 1996)

%% Input
figure(1);
k = 0:1:100;
C = [50*sin(0.1*k)] + 700;
Sold = [ zeros(11,1); 100*ones(30,1); zeros(30,1); 50*ones(10,1); zeros(20,1)] + 1000;
CSold = C.*Sold';  CSold = CSold(:);        %% Make a new input from a non-linear operation
xa = [100*sin(0.1*k+3)] + 300;  xa = xa(:);
xb = [100*sin(0.2*k+3)] + 275;  xb = xb(:);
xc = [100*sin(0.3*k+3)] + 300;  xc = xc(:);
xd = [100*sin(0.4*k+3)] + 275;  xd = xd(:);
plot(k,CSold/1000,k,xa,k,xb,k,xc,k,xd);
title('an example input');
ylabel('unit and 1000$');  xlabel('time (day)');  grid;
input = [ CSold xa xb xc xd];


%% Model
figure(2);
B0 = [   1    0    0    0    0];   A0 = [ 1];  %% delay 0
B1 = [   0    0    0    0    0];   A1 = [-1];  %% delay 1
B2 = [   0    0    0    0    0];               %% delay 2
B3 = [   0    0 -600    0    0];               %% delay 3
B4 = [   0    0    0    0 -550];               %% delay 4
B5 = [   0 -500    0    0    0];               %% delay 5
B6 = [   0    0    0    0    0];               %% delay 6
B7 = [   0    0    0 -400    0];               %% delay 7
AA = [ A0 A1];
BB = [ B0 B1 B2 B3 B4 B5 B6 B7];
Model = arx2th( AA, BB, 1, 5);
y = idsim(input, Model);
plot(k,y);
title('an example output from the simulation');
ylabel('profit ($)');  xlabel('time (day)');  grid;


%% Add noise to the output
figure(3);
y_n = y + max(y)*0.1*(randn(size(y)) - 0.5);
plot(k,y_n);
title('an example output from the simulation + noise');
ylabel('profit ($)');  xlabel('time (day)');  grid;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% ARX Method
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Identification
ident_data = [y input];
na = [1];
nb = [1 1 1 1 1];
nc = [0 5 3 7 4];
[I_A I_B] = th2arx(arx(ident_data, [na nb nc]));

%% Compare
I_A
I_B

%% Identification
ident_data = [y_n input];
na = [1];
nb = [1 1 1 1 1];
nc = [0 5 3 7 4];
[I_A I_B] = th2arx(arx(ident_data, [na nb nc]));

%% Compare
I_A
I_B


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% RARX Method
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
lam = 0.9;
%% Shift the input because arx doesn't except nc = 0
input = [ [CSold(2:length(CSold));0] xa xb xc xd];

%% Identification
ident_data = [y input];
na = [1];
nb = [1 1 1 1 1];
nc = [1 5 3 7 4];
I_Model = rarx(ident_data, [na nb nc], 'ff', lam);

%% Compare
I_Model(length(ident_data),:)

%% Identification
ident_data = [y_n input];
na = [1];
nb = [1 1 1 1 1];
nc = [1 5 3 7 4];
I_Model = rarx(ident_data, [na nb nc], 'ff', lam);

%% Compare
I_Model(length(ident_data),:)
