octave_packages/control-2.3.52/Madievski.m

   1 %% -*- texinfo -*-
   2 %% Frequency-weighted controller reduction.
   3
   4 % ===============================================================================
   5 % Frequency Weighted Controller Reduction       Lukas Reichlin      December 2011
   6 % ===============================================================================
   7 % Reference: Madievski, A.G. and Anderson, B.D.O.
   8 %            Sampled-Data Controller Reduction Procedure
   9 %            IEEE Transactions of Automatic Control
  10 %            Vol. 40, No. 11, November 1995
  11 % ===============================================================================
  12
  13 % Tabula Rasa
  14 clear all, close all, clc
  15
  16 % Plant
  17 Ap1 = [  0.0         1.0
  18          0.0         0.0     ];
  19
  20 Ap2 = [ -0.015       0.765
  21         -0.765      -0.015   ];
  22
  23 Ap3 = [ -0.028       1.410
  24         -1.410      -0.028   ];
  25
  26 Ap4 = [ -0.04        1.85
  27         -1.85       -0.04    ];
  28
  29 Ap = blkdiag (Ap1, Ap2, Ap3, Ap4);
  30
  31 Bp = [   0.026
  32         -0.251
  33          0.033
  34         -0.886
  35         -4.017
  36          0.145
  37          3.604
  38          0.280   ];
  39
  40 Cp = [  -0.996      -0.105       0.261       0.009      -0.001      -0.043       0.002      -0.026   ];
  41
  42 Dp = [   0.0     ];
  43
  44 P = ss (Ap, Bp, Cp, Dp);
  45
  46 % Controller
  47 Ac = [  -0.4077      0.9741      0.1073      0.0131      0.0023     -0.0186     -0.0003     -0.0098
  48         -0.0977     -0.1750      0.0215     -0.0896     -0.0260      0.0057      0.0109     -0.0105
  49          0.0011      0.0218     -0.0148      0.7769      0.0034     -0.0013     -0.0014      0.0011
  50         -0.0361     -0.5853     -0.7701     -0.3341     -0.0915      0.0334      0.0378     -0.0290
  51         -0.1716     -2.6546     -0.0210     -1.4467     -0.4428      1.5611      0.1715     -0.1318
  52         -0.0020      0.0950      0.0029      0.0523     -1.3950     -0.0338     -0.0062      0.0045
  53          0.1607      2.3824      0.0170      1.2979      0.3721     -0.1353     -0.1938      1.9685
  54         -0.0006      0.1837      0.0048      0.1010      0.0289     -0.0111     -1.8619     -0.0311  ];
  55
  56 Bc = [  -0.4105
  57         -0.0868
  58         -0.0004
  59          0.0036
  60          0.0081
  61         -0.0085
  62         -0.0004
  63         -0.0132  ];
  64
  65 Cc = [  -0.0447     -0.6611     -0.0047     -0.3601     -0.1033      0.0375      0.0427     -0.0329  ];
  66
  67 Dc = [   0.0     ];
  68
  69 K = ss (Ac, Bc, Cc, Dc);
  70
  71 % Controller Reduction
  72 Kr4 = spaconred (P, K, 4, 'feedback', '-')
  73 Kr2 = spaconred (P, K, 2, 'feedback', '-')
  74
  75 % Open Loop
  76 L = P * K;
  77 Lr4 = P * Kr4;
  78 Lr2 = P * Kr2;
  79
  80 % Closed Loop
  81 T = feedback (L);
  82 Tr4 = feedback (Lr4);
  83 Tr2 = feedback (Lr2);
  84
  85 % Frequency Range
  86 w = {1e-2, 1e1};
  87
  88 % Bode Plot of Controller
  89 [mag, pha, w] = bode (K, w);
  90 [magr4, phar4, wr4] = bode (Kr4, w);
  91 [magr2, phar2, wr2] = bode (Kr2, w);
  92
  93 mag = 20 * log10 (mag);
  94 magr4 = 20 * log10 (magr4);
  95 magr2 = 20 * log10 (magr2);
  96
  97 figure (1)
  98 subplot (2, 1, 1)
  99 semilogx (w, mag, wr4, magr4, wr2, magr2)
 100 axis ('tight')
 101 ylim (__axis_margin__ (ylim))
 102 grid ('on')
 103 title ('Bode Diagrams of K and Kr')
 104 ylabel ('Magnitude [dB]')
 105
 106 subplot (2, 1, 2)
 107 semilogx (w, pha, wr4, phar4, wr2, phar2)
 108 axis ('tight')
 109 ylim (__axis_margin__ (ylim))
 110 grid ('on')
 111 xlabel ('Frequency [rad/s]')
 112 ylabel ('Phase [deg]')
 113 legend ('K (8 states)', 'Kr (4 states)', 'Kr (2 states)', 'location', 'southwest')
 114
 115 % Step Response of Closed Loop
 116 [y, t] = step (T, 100);
 117 [yr4, tr4] = step (Tr4, 100);
 118 [yr2, tr2] = step (Tr2, 100);
 119
 120 figure (2)
 121 plot (t, y, tr4, yr4, tr2, yr2)
 122 grid ('on')
 123 title ('Step Response of Closed Loop')
 124 xlabel ('Time [s]')
 125 ylabel ('Output [-]')
 126 legend ('K (8 states)', 'Kr (4 states)', 'Kr (2 states)', 'Location', 'SouthEast')