% FILENAME: OPT1.M
% DATE: OCT 1993 
% VERSION: 2.0
% REVISED: JAN 1995
% COPYRIGHT 1994-5 by T Marlin and T. Marlin    
% This program will call an the Nelder-Meade simplex optimization 
% routine in order to find the best values for the tuning constants 
% Kc Ti Td.
%**********************************************************************

Y0 = Y0i;


while 1
  clc
  disp('*************************************************************')
  disp('*              WELCOME TO THE OPTIMIZER                     *')
  disp('*                                                           *')
  disp('*************************************************************')
  disp('')
  disp('SELECT THE APPROPRIATE MENU ITEM')
  disp('')
  disp('MODIFY...')
  disp('                                             PRESENT VALUES')
  fprintf('1) Time to run simulations for Optimizations      %4.2f\n',time)
  fprintf('2) Time step                                      %4.3f\n',deltat)
  
  fprintf('3) Specify New Set Point                          %4.2f\n',SP)
  
  fprintf('4) Specify Magnitude of Disturbance               %4.2f\n',D)
  
  fprintf('5) Enter Tolerance for Tuning Constants           %4.3g\n',toll)
  fprintf('6) Enter Tolerance for Objective Function         %4.3g\n',toll2)
  if isempty(junk2)
      junk2=junk22;
  end
  junk22=junk2;
  if junk2 == 1
   disp('7) Specify Tuning Paramaters                     Kc Ti Td ')
  elseif junk2 == 2 
   disp('7) Specify Tuning Paramaters                      Kc Ti  ')
  elseif junk2 == 3
   disp('7) Specify Tuning Paramaters                      Kc Td ')
  end
  disp('8) Optimize for a Disturbance')
  disp('9) Optimize for a Set Point change')
  disp('10) Plot Results of Optimization ')
  disp('11) Save Plot to META-file')
  disp('12) Return to Main Menu')      
  disp('')
  Opz = input('Enter the desired selection:   ');
  
  if isempty(Opz), break, end
      
  Opz = round(Opz);
  if (Opz<1) | (Opz>12)
      disp('')
      disp('Not a valid selection')
      disp('')
      disp('Press ENTER to continue')
      pause;
  end    
  
  if Opz == 1
      junk = input('Enter value for end time of integration:  ') ;
      if junk <= 0
       disp('')
       disp(' Entry not accepted, total time must be positive')
       disp('')
       disp(' Press any key to continue')
       pause
       junk=time;
      end

      if isempty(junk)
          junk=time;
      end
      time=junk;
      if (time/deltat) > lim_it
          disp('')
          disp('Exceeded matrix limitations, time recalculated') 
          disp('')
          disp('Press Enter to continue') 
          time= round(lim_it*deltat);
          pause;
      end
      mvtest=2; cn_t=1; scn_t=1; maxit=time/deltat;

  elseif Opz == 2
      junk = input('Enter value for time step:  ') ;
      if junk <= 0
       disp('')
       disp(' Entry not accepted, time step must be positive')
       disp('')
       disp(' Press any key to continue')
       pause
       junk=deltat;
      end

      if isempty(junk)
         junk=deltat;
      end
      deltat=junk;
      if (time/deltat) > lim_it
          disp('')
          disp('Exceeded matrix limitations, time step recalculated')
          disp('')
          disp('Press Enter to continue')
          deltat=round(time/lim_it);
          pause;
      end
      mvtest=2; cn_t=1; scn_t=1; maxit=time/deltat;

  elseif Opz == 3
      junk = input('Enter a value for the SET POINT:  ');
      if isempty(junk)
          junk=SP;
      end
      SP=junk;
      junk= input('Enter the time at which to implement the set point:  ');
      if junk < 0
       disp('')
       disp(' Must be positive')
       disp('')
       disp(' Press any key to continue')
       pause
       junk=tset;
      end
      if junk < theta(1)
       clc
       disp('')
       disp(' Due to program limitations, T set-point must be > dead time')
       disp('')
       disp(' Press any key to continue ...')
       pause
       junk=tset;
      end

       if isempty(junk)
           junk=0.05*time;
       end
       tset=junk;
       betaa=tset; mvtest=2; cn_t=1; scn_t=1;

   elseif Opz == 4      
       junk = input('Enter the magnitude of the disturbance:  ');
       if isempty(junk)
           junk=D;
       end
       D=junk;
       junk = input('Enter the time of the disturbance:   ');
      if junk < 0
       disp('')
       disp(' Must be positive')
       disp('')
       disp(' Press any key to continue')
       pause
       junk=tdist;
      end
      if junk < theta(1)
       clc
       disp('')
       disp(' Due to program limitations, Time of disturbance must be > dead time')
       disp('')
       disp(' Press any key to continue ...')
       pause
       junk=tset;
      end

       if isempty(junk)
           junk=0.05*time;
       end
       tdist=junk;
       betaa=tdist; mvtest=2; cn_t=1; scn_t=1;

   elseif Opz == 5   
        junk = input('Tolerance =  ');
        if junk < 0
         disp('')
         disp(' Tolerance must be positive')
         disp('')
         disp(' Press any key to continue')
         pause
         junk=toll;
        end

        if isempty(junk)
           junk=toll;
        end
        toll=junk;

   elseif Opz == 6
        junk = input('Tolerance =  ');
        if junk < 0
         disp('')
         disp(' Tolerance must be positive')
         disp('')
         disp(' Press any key to continue')
         pause
         junk=toll2; 
        end

        if isempty(junk)
           junk=toll2; 
        end   
        toll2=junk;

   elseif Opz == 7
        disp('')
        disp('1) Kc Ti Td')
        disp('2) Kc Ti (default)')
        disp('3) Kc Td')
        
        junk2 = input('Select 1 2 or 3 : ');
        if isempty(junk2)
           junk2=junk22;
        end
        junk22=junk2;

        if junk2 == 1
          Y0=zeros(1,3);
          Y0=[kc ti td];
        elseif junk2 == 2 
          Y0=zeros(1,2);
          Y0 =[kc ti];  
        elseif junk2 == 3  
          Y0=zeros(1,2);
          Y0=[kc td];
        end
        Y0i=Y0;

   elseif Opz == 8 | Opz == 9
     o_func=0;
     clc
     disp('')
     disp(' Optimization procedure will require plenty of time !!! ')
     disp('')
     junk1=input(' Begin Optimization ...[y/n]','s');
     if junk1 == 'n', break;  end
      
     if isempty(junk1), break;  end 
     
     pp=0;
     mm= -1;
     if Opz == 8
       change=1; betaa=tdist;
     else
       change=2; betaa=tset;
     end
     p=1;
     mvbound1
     useopt=0;
     clc
     disp('')
     disp('')
     disp('')
     disp('')
     disp('')
     disp('')
     disp(' Optimization in progress.... Please wait several moments. ')
     disp('')
     disp('')
     disp(' Numerical calculations in progress will be momentarily displayed.... ')

     [Best, option_s] = fmins('integ4a',Y0,[1 toll toll2]);    % Calls Fmins

      x_ = Best;  ob_j = option_s(8); 
      optimal                          % Calls file to display results.
   
   elseif Opz == 10 | Opz == 11
       if Opz == 10
          save_graph=0;
       else
          save_graph=1; 
       end   
       if op_ts == 0 
         disp('')
         disp('Nothing to plot, run an optimization before using this option') 
         disp('')
         disp('Press Enter to continue')
         pause;
       else 
         grph2
       end
   elseif Opz == 12    
       break         % quit menu
   end
end