% lolimotDemo5 Demo 5. Dynamic process of first order with 2 inputs and 2 outputs

% LoLiMoT - Nonlinear System Identification Toolbox
% Torsten Fischer, 20-October-2011
% Institute of Mechanics & Automatic Control, University of Siegen, Germany
% Copyright (c) 2012 by Prof. Dr.-Ing. Oliver Nelles


% Training
lolimotObject= lolimot;                   % Generate an empty net and data set structure

u1 = randn(300,1);
u2 = sin(linspace(0,2*pi,300))';
y1 = zeros(300,1);
y2 = zeros(300,1);
for k = 2:300
    y1(k) = atan(u1(k-1)) + atan(u2(k-1)) + 0.9*y1(k-1);
    y2(k) = u1(k-1).^2 + u1(k-1)*u2(k-1) + 0.7*y1(k-1) + 0.8*y2(k-1);
end
lolimotObject.input = [u1 u2];
lolimotObject.output = [y1 y2];

lolimotObject.xInputDelay = cell(1,2);
lolimotObject.xOutputDelay = cell(1,2);
lolimotObject.zInputDelay = cell(1,2);
lolimotObject.zOutputDelay = cell(1,2);
lolimotObject.xInputDelay{1} = [1]; lolimotObject.xInputDelay{2} = [1];        % Delayed inputs
lolimotObject.xOutputDelay{1} = [1]; lolimotObject.xOutputDelay{2} = [1];      % Delayed outputs
lolimotObject.zInputDelay{1} = [1]; lolimotObject.zInputDelay{2} = [1];        % Delayed inputs
lolimotObject.zOutputDelay{1} = [1]; lolimotObject.zOutputDelay{2} = [1];      % Delayed outputs

lolimotObject.maxNumberOfLM = 10;                          % Termination criterion for maximal number of LLMs
lolimotObject.minError = 0.02;                              % Termination criterion for minimal error
lolimotObject.kStepPrediction = inf;                        % Simulation not one-step-ahead prediction
lolimotObject.history.displayMode = true;                   % display information

% Train net
lolimotObject = lolimotObject.train;

% Find best performance/complexity trade-off
[~, idxBest] = min(lolimotObject.history.penaltyLossFunction);
lolimotObject.leafModels = lolimotObject.history.leafModelIter{idxBest};

% Generalization
u1G = randn(270,1);
u2G = sin(linspace(0,2*pi,length(u1G)))';
y1G = zeros(length(u1G),1);
y2G = zeros(length(u1G),1);
for k = 2:length(y1G)
    y1G(k) = atan(u1G(k-1)) + atan(u2G(k-1)) + 0.9*y1G(k-1);
    y2G(k) = u1G(k-1).^2 + u1G(k-1)*u2G(k-1) + 0.7*y1G(k-1) + 0.8*y2G(k-1);
end
uG = [u1G u2G];
yG = [y1G y2G];

% Simulate net
yGModel = calculateModelOutput(lolimotObject, uG, yG);
JG = calcGlobalLossFunction(lolimotObject ,yG, yGModel);

figure
subplot(1,2,1)
lolimotObject.plotModel([],1)
subplot(1,2,2)
lolimotObject.plotModel([],2)
