% mimo_dat 
% Training- und Testdaten des CSTR, nicht isotherm
% wichtig sind:
% (Ca_in ist nicht messbar, soll deshalb auch nicht benutzt werden)
% (delta_Q soll auch nicht (als Stellgroesse) benutzt werden)
% delta_F, Tn (Eingangstemperatur) 
% und Ausgangsgroesse Cbn
% ausserdem werden spaeter noch t und f_bias gebraucht
%
% Parâmetros da equação de Arrehnius, i.e., 
%   k(i)= k0(i)*exp(EA(i)/T);
%
%	Kr=[k1  EA1  
%           k2  EA2  
%           k3  EA3 ]
%p.ex.:   
     Kr=[
        1.287e12  -9758.3   4.2;
        1.287e12  -9758.3  -11.0;
        9.043e09  -8560    -41.85;
        ];
%
    Konstante=[
         0.9342*3.01;    % rho*Cp
         4032.0*0.215;   % kw*Ar
         0.1;            % 1/(mk*Cpk)
         10;             % Vr             [m3]
    ];
%

rausch=0; % (1) valores MEDIDOS (com ruído) (0) valores simulados
Cain   =    5.1;
% Tin    = 104.38;
Tin=107.5;
op=2; 
fbias   = [ 8, 14.2, 19 ,  25.50 ]  ; f_bias=fbias(op);
%Qk_bias= -4500;
Qk_bias=-1113.5;

%x0=fsolve('cstr_x',[1.235;0.9;134.14;128.95],[],[],Kr,Konstante, [ Cain; f_bias; Qk_bias; Tin]); % vrs. 5.2
options = optimset('LargeScale','off','MaxFunEvals',1000000); %vrs 5.3

x0 = fsolve('cstr_x',[1.235;0.9;134.14;128.95],options,Kr,Konstante,[Cain, f_bias, Qk_bias,Tin]);% vrs 5.3
Cb_0=x0(2);
 
%  PRBS : Sinais já acrescidos dos BIAS, colocar ZERO nos blocos BIAS em mimo_cst.mdl
tend = 8; t = [0:1/60:tend]; let = length(t); % in Stunden
% Valores dos Patamares min e máx. do sinal APRBS
%       [ F    Q    Ca_in  T_in ]
%minp = [ 25   30   40     40   ];
%maxp = [ 200  250  500    500  ];
% Para as constantes de tempo do processo (valores aproximados):
minp =  [ 15   18   25     25   ];  % Calculado em relação ao bias
maxp =  [ 310  320  410    410  ];  % Calculado em relação ao bias

% Amplitude do sinal
amp  =  [ 17 3750 .1     .1   ];

% VARIÁVEIS MANIPULADAS
delta_F = amp(1)*prbs(let,minp(1),maxp(1),0.8,1)'+f_bias;
for i = 1:5, delta_F(10:let)=(delta_F(10:let)+delta_F(9:let-1))/2; end

delta_Q = amp(2)*prbs(let,minp(2),maxp(2),0.8,1)'+Qk_bias;
for i = 1:5, delta_Q(10:let)=(delta_Q(10:let)+delta_Q(9:let-1))/2; end
%delta_Q=zeros(1,let);

% DISTÚRBIOS
Ca_in = amp(3)*prbs(let,minp(3),maxp(3),0.8,1)'+Cain;
for i = 1:5, Ca_in(10:let)=(Ca_in(10:let)+Ca_in(9:let-1))/2; end
T_in = amp(4)*prbs(let,minp(4),maxp(4),0.8,1)'+Tin;
for i = 1:5, T_in(10:let)=(T_in(10:let)+T_in(9:let-1))/2; end

% Manter variáveis  de entrada constantes
%T_in = Tin * ones(size(T_in));
%Ca_in = Cain * ones(size(Ca_in));
%delta_Q = Qk_bias * ones(size(delta_Q));
%delta_F = f_bias * ones(size(delta_F));

% Degrau em F e Q partindo do estado estacionário
%delta_F(1:10)		= f_bias;
%delta_F(11:25)	= f_bias + 10;
%delta_F(26:90)	= f_bias + 6;
%delta_F(91:180)	= f_bias + 2;
%delta_F(181:250)	= f_bias - 4;
%delta_F(251:330)	= f_bias - 8;
%delta_F(331:410)	= f_bias - 12;
%delta_F(410:end)	= f_bias - 6;
%delta_Q(1:70)		= Qk_bias;
%delta_Q(71:170)	= Qk_bias + 0.4 * Qk_bias;
%delta_Q(171:270)	= Qk_bias - 0.4 * Qk_bias;
%delta_Q(271:400)	= Qk_bias - 0.8 * Qk_bias;
%delta_Q(401:end)	= Qk_bias + 0.5 * Qk_bias;

% Simulação no simulink
delta_f_Q = [delta_F' delta_Q'];

ts = sim('mimo_cst',tend); % vrs 5.3
Cbs = [Cb  Cb_medido]; % Cb_medido ist mit Rauschen 
Temp= [T T_medido];


if rausch, Cb=Cb_medido; T=T_medido; end

Cbn= spline(ts,Cb,t');
Tn = spline(ts,T, t');

Ca_in		= spline(ts,Stellgroesse(:,1),t');
delta_F	= spline(ts,Stellgroesse(:,2),t');
delta_Q	= spline(ts,Stellgroesse(:,3),t');
T_in		= spline(ts,Stellgroesse(:,4),t');

data = [Cbn Tn delta_F delta_Q Ca_in T_in];
t =t';
dat.xy = {[t data(:,1)]  [t data(:,2)]  [t data(:,3)]  [t data(:,4)]  [t data(:,5)]  [t data(:,6)]};
dat.iy={[1] [2] [3] [4] [5] [6]};
dat.title = {'Cbn','T','F','Q','Ca_in','T_in'};
%mpcplot(dat);
save test data t
%save mimodat data t
%clear all
