%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Demonstração da função ssest
% System Identification Toolbox
% Identificação do modelo de estado contínuo usando um método PEM
%--------------------------------------------------------------------------

% Limpa o espaço de trabalho e as figuras

clear all
close all
clc

%--------------------------------------------------------------------------
% Dados entrada/saída do sistema a identificar

%Este ficheiro contém  vectores com amostras de 3 variáveis: u, y e
% ts (instantes de amostragem)
load iodata1sintetic.mat;
Ts=ts(2)-ts(1);  % intervalo de amostragem (s)

% Remove a tendência linear por forma a que os dados para identificação
% sejam de média nula
y_all=detrend(ys);
u=detrend(us);

% Filtragem seguida de eliminação de tendências
af = 0.8;
Afilt = [1 -af];
Bfilt = (1-af)*[1 -1];
y = filter(Bfilt,Afilt,ys);

% Cria a estrutura de dados para identificação
data = iddata(y,u,Ts);
N=length(y);
t=(1:N)*Ts;

%--------------------------------------------------------------------------
% Estimação do modelo de estado

% Gama de ordens possíveis para escolher a melhor
nx = 1:10;
sys = ssest(data,nx,'form','canonical');
nest=5; %ordem do sistema 
%sys = ssest(data,nest,'form','canonical');

%--------------------------------------------------------------------------
% Validação do modelo - comparação da saída do modelo com os dados

compare(data,sys)

% Converte o modelo de estado para função de transferência

A=sys.A;
B=sys.B;
C=sys.C;
D=sys.D;
[num den]=ss2tf(A,B,C,D);
proc_zeros=roots(num)
proc_poles=roots(den)

pzmap(sys)

%--------------------------------------------------------------------------
% Gráficos

% Dados de entrada/saída para identificação
figure

subplot(211)
gg=plot(t,y);
set(gg,'Linewidth',1.5)
gg=xlabel('Tempo (s)');
set(gg,'Fontsize',14);
gg=ylabel('y');
set(gg,'Fontsize',14);
axis([0 N*Ts -1.2*max(abs(y)) 1.1*max(abs(y))])

gg=title('Dados para identificação');
set(gg,'Fontsize',14);

subplot(212)
gg=plot(t,u);
set(gg,'Linewidth',1.5)
gg=xlabel('Tempo (s)');
set(gg,'Fontsize',14);
gg=ylabel('u');
set(gg,'Fontsize',14);
axis([0 N*Ts -1.2*max(abs(u)) 1.1*max(abs(u))])

%--------------------------------------------------------------------------
% Resposta no tempo ao escalão do sistema identificado

Tsim=6;
x0=zeros(nest,1);

sim('identmodel');

figure

gg=plot(tt,yt);
set(gg,'LineWidth',1.5);
gg=xlabel('Tempo (s)');
set(gg,'FontSize',14);
gg=ylabel('y_m');
set(gg,'FontSize',14);
gg=title('Resposta ao escalão do modelo');
set(gg,'FontSize',14);

%--------------------------------------------------------------------------
% Eof
