
function [x_Path, Q_Path, statVec] =  DTSHS(xl, xu, ci)
    % xl    = [17 20 18 18 21];
    % xu    = [22 23 23 22 24]; 
    % ci = [ 6 6  4 6 30];
    % Simulation Time steps
    T = 200;
    % Number of Rooms
    h = 5;
    
    
    % Room Area 
    Area = [10 10 10 20 20];
    
    
    ci = ci ./ Area;
    
    
    x_Path = zeros(T,h);
    Q_Path = zeros(T, h);
    S.X = 6*ones(1,h);
    S.Q = zeros(1,h);
    statVec = zeros(1, h);
    TempMet = zeros(1, h);




    % =======================================================================
    %                          GENERATE THE DATA
    % =======================================================================




    % Model State , Randomly initialized
    % Rooms temperture
    S.X =  random('unif',15,24,[1,h]) ;
    % Heater status , ON = 1, OFF = 0
    S.Q =  random('unid',2,[1,h]) - 1 ;

    %Plot inilizations

    x_Path(1,:) = S.X(:);
    Q_Path(1,:) = S.Q(:);



    %parameter calculations

    % Sigmoidal  function paramters
    % Both Rooms have a steep thershold sigmodal function
    D = [100  100  100 100 100];
    
    % frist room has a medium funtion and high for 2nd room
    Alpha_Low = ((3/4)*xl+ (1/4)*xu);
    Alpha_Med = 0.5* (xl + xu);
    Alpha_High = ((1/4)*xl+ (3/4)*xu);
    Alpha = [ Alpha_Med(1) Alpha_High(2)  Alpha_Low(3) Alpha_Med(4) Alpha_Med(5)];



    % Room Temperture Model Parameters
    %Calculate Gamma parameter matrix
    % transfer rate to ith room's ambient
    bi = [0.038 0.02 0.04 0.05 0.055];

    %ambient temperture
    xa = 6;

    
    Gamma = [0 0 0 0 0];
    %initlize Gmma
    %     for j = 1:2
    %         if (S.Q(j) == 0)
    %             Gamma(j) = bi(j) * xa;
    %         else
    %             Gamma(j) = (bi(j) * xa) + ci(j);
    %         end
    %     end



    % transfer rate from room i to room j

    aij= [0 0.1 0 0.1 0; 0.1 0 0.11 0.05 0.06 ;0 0.11 0 0 0.1; 0.1 0.05 0 0 0.13; 0 0.06 0.1 0.13 0];
    %Calculate Segma parameter matrix
    Segma= aij + diag(-(bi + sum(aij)));
    % Noise Variance
    v_2 = 0.25 ;
    %calculate the normal distribuation Variance Vector
    Var = v_2 * ones(1,h);
    %Model Simulation for [0 N]  time interval
    for i = 1:T

        %Calculate new temp
  

        %Upadate Gamma
        Gamma = (bi .* xa) + ci .* S.Q;
        S.X = Tx_( S.X, Gamma, Segma, Var, h);

        %get new heaters status
        S.Q = Tq_( S.X, D, Alpha , h );

     

        %plot code
        x_Path(1+i,:) = S.X(:);


        Q_Path(1+i,:) = S.Q(:);
        
        %statistical summary to measure the precentage of constrain
        %satisfaction
        
%         TempMet1 = S.X - xl;
%         TempMet1(TempMet1 <0) = 0;
%         TempMet1(TempMet1 >0) = 1;
%         
%         TempMet2 = xu - S.X ;
%         TempMet2(TempMet2 <0) = 0;
%         TempMet2(TempMet2 >0) = 1;
%         
%         TempMet = TempMet1 & TempMet2;
        
            for k = 1:h
               if (S.X(k) > xl(k) &&  S.X(k) < xu(k))
                   
                    TempMet(k) = 1;
               else
                   TempMet(k) = 0;
               end
            end
        
        
        statVec = statVec + TempMet;

    end
    
    
    statVec = statVec ./ T;
    
    
    
    
    
    
    
    

   

end
    
