load ll4

Ny = [	1 0 1 1
      	0 1 1 1
        0 0 1 0
        0 0 0 1 ];         % denominator order
Nu = [0 0 
      0 0
      1 0
      0 1];       % numerator orders

Nd = [0 0
      0 0
      1 0
      0 1];       % transport delays

% Create the system: type = external (input/output representation)
ll4_model=external('4 vessels system',2,4);			% parameters are name, n. of inputs, n. of ouputs
ll4_model=set(ll4_model,'T',10);					% set the sampling time
ll4_model=add_dynamics(ll4_model,Ny,Nu,Nd);			% set the dynamics

inputs = ll4data(:,5:6)*1e5;
outputs  = ll4data(:,1:4);

% Build a Dataset object to hold the data
ds = dataset([inputs outputs],'variables',{'u1' 'u2' 'y1' 'y2' 'y3' 'y4'});

% The first output is a TakSug system (3 inputs, 1 output)
o1 = taksug('first output',3,1);
o1.m = 2.2;
o1.model_code = {'crossproduct' 'inversedist' 'linear'};

o2 = taksug('second output',3,1);
o2.m = 2.2;
o2.model_code = {'crossproduct' 'inversedist' 'linear'};

o3 = taksug('third output',2,1);
o3.m = 2.2;
o3.model_code = {'crossproduct' 'inversedist' 'linear'};

o4 = taksug('fourth output',2,1);
o4.m = 2.2;
o4.model_code = {'projected' 'standard' 'linear'};

% Link the dataset to the model
add_data(ll4_model,ds);

% Link the mappings to the model
add_mapping(ll4_model,o1,o2,o3,o4);

% Specify explicitely which output of the model is to be represented by which output of which mapping
ll4_model.links = [	1 1
					2 1
					3 1
					4 1	];

% Perform the identification using method fmclust (Robert Babuska's method)
identify(ll4_model,[1 2 3 4],struct('method','fmclust','n_rules',[3 3 2 2]));

% Compute the simulation error (and plot the results)
simulation_error(ll4_model,ds(:,[2 1 4 3 6 5]),1:4,'plot');
