1%% -*- texinfo -*-
2%% Robust control of a mass-damper-spring system.
3##
4%% Type @code{which MDSSystem} to locate,
5##
6%% @code{edit MDSSystem} to open and simply
7##
8%% @code{MDSSystem} to run the example file.
9
10% ===============================================================================
11% Robust Control of a Mass-Damper-Spring System     Lukas Reichlin    August 2011
12% ===============================================================================
13% Reference: Gu, D.W., Petkov, P.Hr. and Konstantinov, M.M.
14%            Robust Control Design with Matlab, Springer 2005
15% ===============================================================================
16
17% Tabula Rasa
18clear all, close all, clc
19
20% ===============================================================================
21% System Model
22% ===============================================================================
23%                +---------------+
24%                | d_m   0    0  |
25%          +-----|  0   d_c   0  |<----+
26%      u_m |     |  0    0   d_k |     | y_m
27%      u_c |     +---------------+     | y_c
28%      u_k |                           | y_k
29%          |     +---------------+     |
30%          +---->|               |-----+
31%                |     G_nom     |
32%        u ----->|               |-----> y
33%                +---------------+
34
35% Nominal Values
36m_nom = 3;   % mass
37c_nom = 1;   % damping coefficient
38k_nom = 2;   % spring stiffness
39
40% Perturbations
41p_m = 0.4;   % 40% uncertainty in the mass
42p_c = 0.2;   % 20% uncertainty in the damping coefficient
43p_k = 0.3;   % 30% uncertainty in the spring stiffness
44
45% State-Space Representation
46A =   [            0,            1
47        -k_nom/m_nom, -c_nom/m_nom ];
48
49B1 =  [            0,            0,            0
50                -p_m,   -p_c/m_nom,   -p_k/m_nom ];
51
52B2 =  [            0
53             1/m_nom ];
54
55C1 =  [ -k_nom/m_nom, -c_nom/m_nom
56                   0,        c_nom
57               k_nom,            0 ];
58
59C2 =  [            1,            0 ];
60
61D11 = [         -p_m,   -p_c/m_nom,   -p_k/m_nom
62                   0,            0,            0
63                   0,            0,            0 ];
64
65D12 = [      1/m_nom
66                   0
67                   0 ];
68
69D21 = [            0,            0,            0 ];
70
71D22 = [            0 ];
72
73inname = {'u_m', 'u_c', 'u_k', 'u'};   % input names
74outname = {'y_m', 'y_c', 'y_k', 'y'};  % output names
75
76G_nom = ss (A, [B1, B2], [C1; C2], [D11, D12; D21, D22], ...
77            'inputname', inname, 'outputname', outname);
78
79G = G_nom('y', 'u');                   % extract output y and input u
80
81
82% ===============================================================================
83% Frequency Analysis of Uncertain System
84% ===============================================================================
85
86% Uncertainties: -1 <= delta_m, delta_c, delta_k <= 1
87[delta_m, delta_c, delta_k] = ndgrid ([-1, 0, 1], [-1, 0, 1], [-1, 0, 1]);
88
89% Bode Plots of Perturbed Plants
90w = logspace (-1, 1, 100);             % frequency vector
91Delta = arrayfun (@(m, c, k) diag ([m, c, k]), delta_m(:), delta_c(:), delta_k(:), 'uniformoutput', false);
92G_per = cellfun (@lft, Delta, {G_nom}, 'uniformoutput', false);
93
94figure (1)
95bode (G_per{:}, w)
96legend off
97
98
99% ===============================================================================
100% Mixed Sensitivity H-infinity Controller Design (S over KS Method)
101% ===============================================================================
102%                                    +-------+
103%             +--------------------->|  W_p  |----------> e_p
104%             |                      +-------+
105%             |                      +-------+
106%             |                +---->|  W_u  |----------> e_u
107%             |                |     +-------+
108%             |                |    +---------+
109%             |                |  ->|         |->
110%  r   +    e |   +-------+  u |    |  G_nom  |
111% ----->(+)---+-->|   K   |----+--->|         |----+----> y
112%        ^ -      +-------+         +---------+    |
113%        |                                         |
114%        +-----------------------------------------+
115
116% Weighting Functions
117s = tf ('s');                          % transfer function variable
118W_p = 0.95 * (s^2 + 1.8*s + 10) / (s^2 + 8.0*s + 0.01);  % performance weighting
119W_u = 10^-2;                           % control weighting
120
121% Synthesis
122K_mix = mixsyn (G, W_p, W_u);          % mixed-sensitivity H-infinity synthesis
123
124% Interconnections
125L_mix = G * K_mix;                     % open loop
126T_mix = feedback (L_mix);              % closed loop
127
128
129% ===============================================================================
130% H-infinity Loop-Shaping Design (Normalized Coprime Factor Perturbations)
131% ===============================================================================
132
133% Settings
134W1 = 8 * (2*s + 1) / (0.9*s);          % precompensator
135W2 = 1;                                % postcompensator
136factor = 1.1;                          % suboptimal controller
137
138% Synthesis
139K_ncf = ncfsyn (G, W1, W2, factor);    % positive feedback controller
140
141% Interconnections
142K_ncf = -K_ncf;                        % negative feedback controller
143L_ncf = G * K_ncf;                     % open loop
144T_ncf = feedback (L_ncf);              % closed loop
145
146
147% ===============================================================================
148% Plot Results
149% ===============================================================================
150
151% Bode Plot
152figure (2)
153bode (K_mix, K_ncf)
154
155% Step Response
156figure (3)
157step (T_mix, T_ncf, 10)                % step response for 10 seconds
158
159% ===============================================================================
160