function control = GPC(cdata, sdata, xdata, ydata, A, B, d, N2, rho, Nu)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% function control = GPC(cdata, sdata, xdata, ydata, A, B, d, N2, rho, Nu)
%%
%% Creators Kittiphan Techakittiroj
%%          Thaveesak Krisornsuthasinee
%%          Yujira Jirapinyo
%% Created  20 April 1996 (at ESCI, Case Western Reserve University)
%% Modified 21 April 1996 (at ESCI, Case Western Reserve University)
%% Requires dtp_sp.m
%% MATLAB   version 4.2c (it should work with the lower version also because
%%          there is no special command)
%%
%% Calculate the control signal by using GPC method
%%
%% Output: control    control signal
%% Input:  cdata      vector of setpoints  [ w(t+1); w(t+2); ...; w(t+N2)]
%%         sdata      vector of difference in inputs  [ du(t-1); du(t-2); ...; du(t-N2)]
%%         xdata      vector of inputs  [ u(t-1); u(t-2); ...; u(t-nb-d)]
%%         ydata      vector of measuring outputs  [ y(t-1); y(t-2); ...; y(t-na)]
%%         A, B       system parameters
%%         d          delay of the system
%%         N2         number of looking ahead
%%         rho        weight on changing in input signal
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  if nargin<10
    Nu = N2;
    end;

  ff = [];  GG = [];
  for k=1:N2,
    [E, F] = dpt_sp(k, ([A 0] - [0 A]));
    G  = conv(E,B);
    GG = [[zeros(1, N2-k) G(1:k)]' GG];
    G  = G(k+1:length(G));
    if (length(G)>0)
      f = G*sdata(1:length(G))' + F*ydata';
      else
        f = F*ydata';
      end;
    ff = [ff f];
    end;

  GG = GG(:,1:Nu);
  G  = GG'*GG;
  G  = inv(G + rho*eye(size(G)))*GG';
  control = xdata(1) + G(1,:)*(cdata - ff)';
