octave_packages/control-2.3.52/optiPID.m

   1 %% -*- texinfo -*-
   2 %% Numerical optimization of a PID controller using an objective function.
   3 %% The objective function is located in the file @command{optiPIDfun}.
   4 %% Type @code{which optiPID} to locate, @code{edit optiPID} to open
   5 %% and simply @code{optiPID} to run the example file.
   6
   7 % ===============================================================================
   8 % optiPID                          Lukas Reichlin                       July 2009
   9 % ===============================================================================
  10 % Numerical Optimization of an A/H PID Controller
  11 % Reference: Guzzella, L. (2007) Analysis and Synthesis of SISO Control Systems.
  12 %            vdf Hochschulverlag, Zurich
  13 % Required OCTAVE Packages: control, optim (and its dependencies)
  14 % Required MATLAB Toolboxes: Control, Optimization
  15 % ===============================================================================
  16
  17 % Tabula Rasa
  18 clear all, close all, clc;
  19
  20 % Global Variables
  21 global P t dt mu_1 mu_2 mu_3
  22
  23 % Plant
  24 numP = [1];
  25 denP = conv ([1, 1, 1], [1, 4, 6, 4, 1]);
  26 P = tf (numP, denP);
  27
  28 % Relative Weighting Factors: PLAY AROUND WITH THESE!
  29 mu_1 = 1;               % Minimize ITAE Criterion
  30 mu_2 = 10;              % Minimize Max Overshoot
  31 mu_3 = 20;              % Minimize Sensitivity Ms
  32
  33 % Simulation Settings: PLANT-DEPENDENT!
  34 t_sim = 30;             % Simulation Time [s]
  35 dt = 0.05;              % Sampling Time [s]
  36 t = 0 : dt : t_sim;     % Time Vector [s]
  37
  38 % A/H PID Controller: Ms = 2.0
  39 [gamma, phi, w_gamma, w_phi] = margin (P);
  40
  41 ku = gamma;
  42 Tu = 2*pi / w_gamma;
  43 kappa = inv (dcgain (P));
  44
  45 disp ('optiPID: Astrom/Hagglund PID controller parameters:');
  46 kp_AH = ku * 0.72 * exp ( -1.60 * kappa  +  1.20 * kappa^2 )
  47 Ti_AH = Tu * 0.59 * exp ( -1.30 * kappa  +  0.38 * kappa^2 )
  48 Td_AH = Tu * 0.15 * exp ( -1.40 * kappa  +  0.56 * kappa^2 )
  49
  50 C_AH = optiPIDctrl (kp_AH, Ti_AH, Td_AH);
  51
  52 % Initial Values
  53 C_par_0 = [kp_AH; Ti_AH; Td_AH];
  54
  55 % Optimization
  56 if (exist ('fminsearch'))
  57   warning ('optiPID: optimization starts, please be patient ...');
  58 else
  59   error ('optiPID: please load/install optim package to proceed');
  60 end
  61
  62 C_par_opt = fminsearch (@optiPIDfun, C_par_0);
  63
  64 % Optimized Controller
  65 disp ('optiPID: optimized PID controller parameters:');
  66 kp_opt = C_par_opt(1)
  67 Ti_opt = C_par_opt(2)
  68 Td_opt = C_par_opt(3)
  69
  70 C_opt = optiPIDctrl (kp_opt, Ti_opt, Td_opt);
  71
  72 % Open Loop
  73 L_AH = P * C_AH;
  74 L_opt = P * C_opt;
  75
  76 % Closed Loop
  77 T_AH = feedback (L_AH, 1);
  78 T_opt = feedback (L_opt, 1);
  79
  80 % A Posteriori Stability Check
  81 disp ('optiPID: closed-loop stability check:');
  82 st_AH = isstable (T_AH)
  83 st_opt = isstable (T_opt)
  84
  85 % Stability Margins
  86 disp ('optiPID: gain margin gamma [-] and phase margin phi [deg]:');
  87 [gamma_AH, phi_AH] = margin (L_AH)
  88 [gamma_opt, phi_opt] = margin (L_opt)
  89
  90 % Plot Step Response
  91 [y_AH, t_AH] = step (T_AH, t);
  92 [y_opt, t_opt] = step (T_opt, t);
  93
  94 figure (1)
  95 plot (t_AH, y_AH, 'b', t_opt, y_opt, 'r')
  96 grid ('on')
  97 title ('Step Response')
  98 xlabel ('Time [s]')
  99 ylabel ('Output [-]')
 100 legend ('A/H', 'Optimized', 'Location', 'SouthEast')
 101
 102 % ===============================================================================