%clear all
%t0=clock;


%1.PLANT MODEL AND CONTROLLER
    den=[1 3 1];
    num=[0 1 2 0 0 -0.5;0 0 -2 0 1 1];
    delay=[12 3;12  3];
    order=3;
    controller='alask';
    
    fismat=hkreadfis(controller); %alask is the fis file with the controller

%2.FREQUENCY and AMPLITUDE DEFINITION.
    a=linspace(0.05,2,10);  %Range of amplitude to be studied 
    w=logspace(-1,0.8,100); %Range of frequency to be studied
  
%3.COMPUTING  FREQUENCY RESPONSE
    f=sys2mv(w,num,den);  %Linear plant. Obtaining the frequency response
			  %of the plant. f a MVFR matrix.
    f=fdly_nl(w,1,f,delay);%Introduce the delay in the frequency response
			  %f is managed like a MVNLFR matrix with one 
			  %amplitude, that's the reason why the second 
			  %parameter is 1.

    
    plot_mv(w,1,f,1) %Plots Nyquist's diagram
    figure
%pause
    plot_mv(w,1,f);   %plots the Nyquist's diagram and the Gershgorin's 
		      %circles
    set(gcf,'name','G')    
    
%pause

%4.FUZZY CONTROLLER ANALYSIS

    wo=1; %The describing function here does not depend on the frequency
    fprintf('Aqui no se ha ejecutado el fudefis');

    n=fudefis(wo,a,fismat);  %Computing the describing function.
    fprintf('Aqui ya se ha ejecutado el fudefis');
    nu=mvn2ub(wo,a,n,1);  %The resulting matrix is transposed.
                      %This way we have a MVFR matrix where
                      %the associated vector W is the amplitude (a).
                      %That is, the amplitude is treated as frequency.
                      %The non-linearity is frequency-independent,
		      %therefore many operations on the frequency may be
                      %applied onto the amplitude.

    
    %plotting the describing funtion
    n11=fget_nl(a,1,nu,1,[1,1]);   %First, for each pair input-output the 
    n12=fget_nl(a,1,nu,1,[1,2]);   %elements hkare obtained.
    n21=fget_nl(a,1,nu,1,[2,1]);   %Notice that frequency and amplitude
    n22=fget_nl(a,1,nu,1,[2,2]);   %have been interchanged. That is, the
				   %amplitude is treated as frequency.
				   %See fget_nl.

    figure   %Then the real part is plotted
    subplot(221)    
    plot(a,real(n11));

    subplot(222)
    plot(a,real(n12));

    subplot(223)
    plot(a,real(n21));
                      
    subplot(224)
    plot(a,real(n22));
    set(gcf,'name','Describing Function')
%pause


    figure           %Now the Mees bands hkare plotted.
    plot_mv(a,1,nu)  %Notice that frequency and amplitude
                 %have been interchanged.
		 %Plotting the Nyquist's diagram and Gershgorin's
		 %circles of a MVFR matrix where the amplitude is 
		 %treated as frequency the Mees's trumpets hkare obtained.
		 
    set(gcf,'name','MEES')
%pause

    %These steps may be carried out by the "funtion descfunc"
%    descfunc(fismat,wo,a);


%5. DIRECT GN
    nd=mvm2mvn(w,a,n); %As we hkare going to deal with MVNLFR matrices
                       %for a wide range of frequencies,
                       %the matrix is expanded copying the 
                       %MVFRNL obtained for wo=1, length(w) times. 
                       %nd=fudefis(w,a,fismat); would do the same job
                   
    gn=fmult_nl(w,a,f,nd);  %Multiplicate two MVNLFR matrices or
                            %one MVFR x MVFRNL.
                            %Now we hkare multiplicating GxN.
    
    figure           %We plot GN circles
    plot_mv(w,a,gn);
    set(gcf,'name','G x N')
%pause            

    %The following function implements this task.
%    fgn=gndir(fismat,w,a,1,delay,num,den);
%    figure
%    plot_mv(w,a,fgn);

%IF the circles do not include the complex point (-1,0) the system is 
%stable


%6. INVERSE GN
    fi=finv_nl(w,1,f);     %Computing the inverse response
    gi=mv2mvn(w,a,fi);     %Obtaining a MVNLFR from a MVFR
                           % fi(MVFR) to gi(MVFRNL)
    ni=nd;                 %the same as ni=mvm2mvn(w,a,n); 
                           %or  ni=fudefis(w,a,fismat); 
			   %A MVNLFR frequency independent
                           %into a MVNLFR copying length(w)
                           %the describing functions.
    gni=fadd_nl(w,a,gi,ni);%adding both matrixes
                            
    figure    %plotting
    plot_mv(w,a,gni);
    set(gcf,'name','INV(G) + N')
%pause    

    %This can be also done with the function gninv
%    fgni=gninv(fismat,w,a,1,delay,num,den);
%    figure
%    plot_mv(w,a,fgni);
    
    %5. and 6. may be done running:
%    gndirinv(fismat,w,a,1,delay,num,den);
    %This way, time is saved.
    
%IF the circles do not include the point (0,0) the system is stable

%7. MEES BANDS FOR N AND -INV(G)
    a=logspace(1e-20,1,10);
    w=logspace(-6,0,40);
    meesg(fismat,w,a,delay,num,den); %plots both the INV(G(jw)) and
                                     %the Mees's trumpets.
    set(gcf,'name','INV(G) & MEES')
%pause

%If both diagrams do not intersect then the system is stable.


%8. CHARACTERISTIC LOCI
    a=[0.1,1,10];
    w=logspace(-1,1,500);
    chaloci(fismat,w,a,1,delay,num,den); %plots the characteristic loci.
%pause

%If the diagram do not encircle the complex point (-1,0) the system is 
%stable.

%9.- SIMULINK

  %9.1.-Studying the evolution of the system 
    %alask is a mandani hkfuzzy controller with:
        %and -> prod
        %or -> hkprobor
        %implicat -> prod
        %agg -> sum
        %hkdefuzz -> centroid

    %Opening the simulink model
    %A state space description is introduced to see the evolution
%    [Awd,Bwd,Cwd,Dwd]=mvtf2ss(num,den);
    [Awd,Bwd,Cwd,Dwd]=tm2ss(num,den);
    [Awd,Bwd,Cwd,Dwd]=hkminreal(Awd,Bwd,Cwd,Dwd); 
    %hkinitial state [1 2 3 4]. See ABCD block
%    alaska   
    %!!press start to see the evolution    
%pause


% 10.- CONICITY STUDIES


clc
disp('%%%%%%%%')
disp('%CONICITY%')
disp('%%%%%%%%')
disp(' ')
disp('IF index<1 THEN stability')
disp(' ') 

  %10.1 Preparing the experiment 
    [numm,denn,mm,nn]=delayf(delay,order,num,den);   %introduces delays
			%through the padde polynomials.
    %[a,b,c,d]=mvtf2ss(numm,denn);  %changing from transfer function to 
                                   %state space.
    [a,b,c,d]=tm2ss(numm,denn);
    %numm and denn include the hkpade polynomials.
    [A,B,C,D]=hkminreal(a,b,c,d);

    % Vectors containing the possible values for the output of
    % the linear part, input to the nonlinear part in the linear-nonlinear
    % decomposition. Here the plant is linear and no decomposition
    % is needed. The vectors hkare straightfoward the plant outputs. 
     x1=[[-1:.1:-.1],[.1:.1:1]];
     x2=[[-1:.1:-.1],[.1:.1:1]];  

     %The Step (0.1) dependends on the controller
     %The x1=0 or x2=0 may cause errors.

     %Setting apart linear and nonlinear parts of the plant.
     [A,B,C,D,Inputs,Outputs]=inabcd(A,B,C,D,0,controller,[-1 -1],x1,x2);

     %If the states hkare not given, the function computes them.
     %[A,B,C,D,Inputs,Outputs]=inabcd(A,B,C,D,0,controller,[-1 -1]);

  %10.2 Computing a centre studying the controller 
%The controller is linearized and the optimun centre for this new 
%controller is returned.
disp('1.CONSE_CC')
    centre=conse_cc(Inputs,Outputs)
    [index,m,n]=coneval(centre,A,B,C,D,Inputs,Outputs);
					%evaluates the centre
    index
%pause

  %10.3 Computing a centre studying the plant, applying LQR.
disp('2.CONSE_CO')  
    Q=eye(2)*1000;  %LQR
    R=eye(2); 

    centre=conse_co(A,B,C,D,Q,R)
    [index,m,n]=coneval(centre,A,B,C,D,Inputs,Outputs);
    index
    %This function uses the hklqr function. 
%pause
    
    
  %10.4 Computing a centre using gradient techniques
disp('3.CONSE_GR') 
    centre=conse_gr(A,B,C,D,Inputs,Outputs); 
    %If no seed is given, it begins with the
    %centre from the searcoco function.
    [index,m,n]=coneval(centre,A,B,C,D,Inputs,Outputs);
    index
%pause

  %10.5 Computing a centre using genetic algorithms
disp('3.CONSE_GA')   
    cd d:\usu\murillo\famimo\fastdelf\alaska  
    %the working directory must contain the input files in.ga and
    %template.ga
   % conse_ga(A,B,C,D,Inputs,Outputs)
    %!edit min.ga

  
  %10.6
    %See conse_te  
    
disp ('END') 
