Entry 3266

clc
close all
clear all

%%INIZIALIZZAZIONE DELLE VARIABILI E DELLE STRUTTURE DATI UTILIZZATE DAL MODELLO:

% 1 - Caricamento valori del ramo ascendente del major loop ricavati sperimentalmente; 
p2f =fopen('ramo_ascendente_1.txt','r');
fu = fscanf(p2f, '%f',[1 13]);
fclose(p2f);

% 2 - Definizione vettori alfa e beta per la memorizzazione della "storia"
% del sistema;
alfa = zeros(1,30);
beta = zeros(1,30);

% 3 - Calcolo dell'uscita all'istante zero (condizione di reset della
% memoria);
f(1) = fu(indice(u(1)));    %indice(): restituisce il corretto indice del vettore in base al valore dell'ingresso;

% 4 - Inizializzazione dei valori estremi che delimitano il "piano di piano
% di Preisach";

alfa(1) = 0.180;
beta(1) = -0.180;

% 5 - Inizilizzazione variabili di controllo del modello (numero di
% segmenti orizzontali, valore precedente dell'ingresso, orientamento
% ultimo segmento della L(t));

n = 1;
up = u(1);
L = 0;

%6 - Inizializzazione della sequenza dei valori di input da dare al
%sistema;

u = [-0.180 -0.150 -0.120 -0.100 -0.080 -0.050 -0.010 0.050 0.080 0.100 0.120 0.150 0.180 0.150 0.120 0.100 0.080 0.050 -0.010 -0.050 -0.080 -0.100 -0.120 -0.150 -0.180];


%% MODELLO DI PREISACH NON LINEARE:

for i=2:length(u)
    calcola = 1;
    
    % caso in cui l'ingresso è pari al valore di saturazione positiva
    % (negativa):
    if( u(i) == 0.180 )
        alfa(n+1) = u(i);
        n = 1;
        f(i) = fu(indice(u(i)));
        calcola = 0;
    elseif( u(i) == -0.180 )
        f(i) = fu(indice(u(i)));
        calcola = 0;
        n = 1;
    end    
        
    
while (calcola)
   
  % Distinzione dei casi in base all'orientamento dell'ultimo segmento della L(T):  
  switch (L)

   case 0 % L(t) orizzontale 

       
   % u(k)>u(k-1):
   
      if ( u(i) > up )       
        if ( u(i) < alfa(n) )
           alfa(n+1) = u(i); %aggiungi un nuovo alfa
           calcola = 0;
         elseif ( u(i) >= alfa(n) )  
           alfa(n) = u(i);  %aggiorna alfa
             if ( n>1 )
                 n = 1;
             end
           calcola = 0;
        end
      end
      

   % u(k)<u(k-1):
   
      if ( u(i) < up )
         if ( u(i) > beta(n) )  %aggiungi un nuovo beta
             beta(n+1) = u(i);
             n = n + 1;
             L = 1;           %passa al caso L(t) verticale per il calcolo effettivo
         elseif ( u(i) <= beta(n) )
             beta(n) = u(i);    %aggiorna beta 
             n = n + 1;
             L = 1;           %passa al caso L(t) verticale per il calcolo effettivo
         end
      end
 

  case 1 % L(t) verticale

   % u(k)>u(k-1):
   
      if ( u(i) > up )
        if ( u(i) < alfa(n) ) 
           alfa(n+1) = u(i);   %aggiungi un nuovo alfa
           L = 0;               %passa al caso L(t) orizzontale per il calcolo effettivo
        elseif( u(i) >= alfa(n) )
           alfa(n) = u(i);     %aggiorna alfa
           n = n - 1;
           L = 0;               %passa al caso L(t) orizzontale per il calcolo effettivo
        end
        
   % u(k)<u(k-1):
   
      elseif ( u(i) < up )
         if ( u(i) > beta(n-1) )
              beta(n) = u(i);    %aggiorna beta
              calcola = 0;
         elseif ( u(i) < beta (n-1) )
             beta(n-1) = u(i);   %aggiorna beta
             n = n - 1;
             calcola = 0;
         end
    end
  end
end
P = 0;

% implementazione formula (4.3) pag. 57

if n >= 2
for k=2:n
    P_1 = preisach_dib(alfa(k),beta(k-1),u(i));   %preisach_dib(): funzione che prendendo in ingresso alfa, beta
    if (beta(k) == u(i))                          %e il valore corrente dell'ingresso restituisce il valore di P(alfa,beta,u);
        P_2 = 0;
    else
    P_2 = preisach_dib(alfa(k),beta(k),u(i));
    end
    Pk = P_1 - P_2;
    P = Pk + P;
end
    f(i) = fu(indice(u(i))) + P;
else
f(i) = fu(indice(u(i)));
end
up = u(i);

end