% Template MATLAB code for prediction from a channel as a 2D line plot using PLOT function.
% Sprungantwort eines PT1-Systems

% Anpassen der Randbedingungen an die eigene Applikation
readChannelID = [167884];  % Anpassen an eigenen Kanal
fieldID1      = [1 2];     % für die ersten beiden Felder (Temperatur, Feuchte)
NumberOfData  = 50;        % die letzten 50 Datenpunkte holen

est_time = 5; % letzten 5 Minuten zur Identifikation (Learning)
pred_time= 5; % 5 Minuten Vorhersage der Zukunft (Prädiktion)

% Channel Read API Key, if your channel is private, then enter the read API Key between the '' below: 
readAPIKey = '';

%% Read Data %
[data, time] = thingSpeakRead(readChannelID, 'Field', fieldID1, 'NumPoints', NumberOfData, 'ReadKey', readAPIKey);

% Auswahl der Zeitachse
i       = find(time > time(end)-minutes(est_time)); % Alle Daten aus Lernfenster
dat     = data(i,:);                                
t       = minutes(time(i)-time(i(1)));            % Zeitachse Minuten
t_est   = [t; (t(end)+.1:t(end)+pred_time)'];     % mit Prädiktionshorizont
time_est= time(i(1))+minutes(t_est);              % Zeitachse mit Datum und Stunde

modell   = @(p,t) p(1).*(1-exp(p(2).*t))+p(3);    % Sprungantwort PT1
Ordnung  = 1;                                     % vs. Polynom 1. Ordnung

% Für alle gewünschten Felder 
for feldindex = 1:length(fieldID1)                % Modell für jedes Feld
  
  mess     = data(i,feldindex);                   % Messdaten
  
  % "Grey-Box" Sprungantwort PT1
  p0       = [1; -1; mess(1)];                    % Startparameter
  popt_PT1 = fminsearch(@(p) norm(mess - modell(p,t)),p0); % Kalibrierung
  sse_PT1  = norm(mess-modell(popt_PT1,t));       % Fehler 
  
  % Black-Box - Polynom 1. Ordnung
  popt_P = polyfit(t,mess,Ordnung);               % Kalibrierung
  sse_P  = norm(mess-polyval(popt_P,t));          % Fehler
  
  if sse_PT1 < sse_P                              % bessere Modell gewinnt
      disp('PT1')
      Vorhersage(:,feldindex) = modell(popt_PT1,t_est);
  else
      Vorhersage(:,feldindex) = polyval(popt_P,t_est);
      disp('Polynom')
  end  
end

plot(time, data,'o',time_est,Vorhersage);