function [A,B]=adaptmulti(rlsdata)
%Adaptiver based on RLS method in the multivariable case
%rlsdata is data from the GPC function

global ny nu yrls urls w A B ypast beta P K

%declare variables
ypred=zeros(1,3);

%w is the setpoint trajectory
%A and B is modelparamters
%ypast is the past mesured values
%yrls and urls is vectors the past vector is made from
%beta is the parameter estimate in a column vector, P is the covar. matrix
%and K is the gain vector

%Save the past data
yrls=[rlsdata(1:ny);yrls]; %update the past measurement vector
urls=[rlsdata(ny+1:ny+nu);urls]; %update the past control signals vector
i=size(yrls,1); %count how many data the rls algorithm has

%Choose values for the degree of polynomnials and the deadtime
n=2;  %Degree of the A polynominal
m=4;  %Degree of the B polynominal
d=1;  %Deadtime

%Test if there is enough data to run RLS
if i>=(m+d+ny)*ny
   %Constants
   c=1000; %for initializing of the P matrix
   gama=0.9;%forgetting factor
   alfa=0.0045; %variance to the error, found from PRBS-data
  
  %Initializing of the system
  beta_init=zeros(n*ny+m*nu,ny); %initialize the beta-vector
  I=eye(m*nu+n*ny);%identity matrix
  Pinit=c*I; %initialize the P matrix
   %vector with past data from measurement and ctrl signal
   ksi=[yrls(1:ny*n);urls(1+d:nu*m+d)];
   %Compute model prediction
   if i==(m+d+ny)*ny
      ypred=ksi'*beta_init; %if the loop runs the first time, uses the initlized beta
   else
      ypred=ksi'*beta; %if the loop runs 2nd time or more, uses the beta from last run
   end
      
   %Update gain vector
   if i==(m+d+ny)*ny
      K=(Pinit*ksi)./(alfa*gama+ksi'*Pinit*ksi); 
   else
      K=(P*ksi)./(alfa*gama+ksi'*P*ksi);
   end
   
   %Update parameter estimates
   if i==(m+d+ny)*ny
      beta=beta_init+K*((yrls(1:ny))'-ypred);
      %update modelparameter A. The first coeff is always the identity matrix
      A=[eye(ny);-beta(1:ny*n,1:ny)]';
      %updata B. -beta comes from different definition of A and B in RLS and GPC
      B=[zeros(d*nu,ny);-beta(n*ny+1:n*ny+m*nu,1:ny)]';
   else
      beta=beta+K*((yrls(1:ny))'-ypred);
      A=[eye(ny);-beta(1:ny*n,1:ny)]';
      B=[zeros(d*nu,ny);-beta(n*ny+1:n*ny+m*nu,1:ny)]';
   end
         
   %Compute covariance matrix for next iteration
   if i==(m+d+ny)*ny
      P=(1/gama)*(I-K*ksi')*Pinit;
   else
      P=(1/gama)*(I-K*ksi')*P;
   end
else %else if i>m+d+1
   %if not enough data, the ARX-model found from PRBS-data is used
   %A given from the PRBS data
A=[ 1.0000         0         0   -1.6900    0.2388    0.0595    0.6971  -0.2364   -0.0571
         0    1.0000         0    0.0251   -1.9607    0.0451   -0.0256   0.9633   -0.0429
         0         0    1.0000    0.1511    0.0849   -1.6938   -0.1500  -0.0797    0.6990];
      
B=[0         0         0    0.0179    0.0055   -0.0967    0.0032  -0.0023    0.1440   -0.0029   -0.0018   -0.0362    0.0011   -0.0004  -0.0148
   0         0         0    0.0008    0.0001    0.0037    0.0004  -0.0005   -0.0005   -0.0005    0.0001   -0.0054    0.0007    0.0002   0.0010
   0         0         0    0.0009    0.0009   -0.0235    0.0049  0.0002     0.0353    0.0003   -0.0007   -0.0132    0.0040    0.0000   0.0005];
   
end %end if i>m+d+1

