Hassan Bevrani
Professor, University of Kurdistan
Fall 2023
H. Bevrani University of Kurdistan 1
Contents
1. Mass-Spring-Damper Control System
2. State-Space Model
3. Nominal Parameters and Uncertainties
4. Mixed-Sensitivity Problem
5. H Control Design
6. Simulation Results
H. Bevrani University of Kurdistan 2
Reference
1. M. Hirata, Practical Robust Control, CORONA Press , 2017 (In Japanese).
H. Bevrani University of Kurdistan 3
4th-order Mass-Spring-Damper System with Uncertainty
System and Dynamic Model:
H. Bevrani University of Kurdistan 4
Continue
H. Bevrani University of Kurdistan 5
Modeling Using MATLAB Codes (from MATLAB Program 6)
The nominal parameters and uncertainty are Parameters Nominal Value
assumed as:
Uncertainty: 20% perturbation in !!
%% 4th order Mass-Spring-Damper System
%% Parameter definition
m1 = 0.8; m2 = 0.2; k1 = 100;
k2 = ureal('k2',300,'percent',20);
c1 = 1; c2 = 0.3; Ks = 100; Bode Diagram
%% Define the M, K, C matrices of the motion equation 20
M = [ m1, 0 ; 0, m2 ];
0
C = [ c1+c2, -c2 ; -c2, c2 ];
Magnitude (dB)
K = [ k1+k2, -k2 ; -k2, k2 ]; -20
F = [ Ks ; 0 ];
%% State space realization -40
iM = inv(M); -60
Ap = [ zeros(2,2), eye(2,2) ; -iM*K, -iM*C ]; 0
Bp = [ zeros(2,1) ; iM*F ]; -90
Cp = [ 0 1 0 0 ];
Phase (deg)
Dp = 0; -180
%% Bode plot of the plant -270
P = ss(Ap,Bp,Cp,Dp);
bode(P,{1e0,1e2}); % Bode plot -360
10 0 10 1 10 2
Frequency (rad/s)
H. Bevrani University of Kurdistan 6
Uncertainties and its Cover (from MATLAB Program 7)
Considering the multiplicative uncertainties, one can give the following weighting
transfer function to cover all uncertainties.
Bode Diagram
%% Uncertainties (Deltam) and its Cover function (Wm) 20
m
Wm
%% Multiplicative uncertainty model
0
w = logspace(0,3,100); % Definition of frequency vector
P_g = ufrd(P,w); % Frequency response calculation
-20
%% Computing multiplicative uncertainty:
Dm_g = (P_g - P_g.nominal)/P_g.nominal;
Magnitude (dB)
-40
%% Definition of weight Wm
s = tf('s'); -60
Wm = 3*s^2/(s^2+18*s+45^2);
-80
%% Gain plot
bodemag(Dm_g,'--',Wm,'r-',w);
legend('\Deltam','Wm'); -100
-120
0 1 2 3
10 10 10 10
Frequency (rad/s)
H. Bevrani University of Kurdistan 7
H Control Design Using Mixed-Sensitivity Problem
Performance Requirement: To have a small enough sensitivity function at low
frequency, we can select the following performance weighting function:
Mixed-Sensitivity Problem: Now, we can formulate the
control design problem via a mixed-sensitivity problem
(MSP), which can be considered in form of input-side Input-side mixed-sensitivity problem
MSP or output-side MSP.
Output-side mixed-sensitivity problem
H. Bevrani University of Kurdistan 8
Framework Preparation
Using the input-side MSP the following configuration can be reached. To satisfy the
full rankness of the D12 in the given generalized plant, a new input with gain of
" = 5×10"# is added.
H. Bevrani University of Kurdistan 9
Building the Generalized Plant
MATLAB codes for plotting the uncertainty and performance weighting functions
and building the generalized plant. Bode Diagram
30
%% Determine of weighting functions Ws
s = tf('s'); Wt
20
Ws = 15/(s + 1.5e-2); % Ws
Wm = 3*s^2/(s^2+18*s+45^2);
10
Wt = Wm; % Wt
Weps = 5e-4; % Weps
figure(3); 0
w=logspace(0,3,100);
Magnitude (dB)
bodemag(Ws,Wt,'--',w); -10
legend('Ws','Wt');
%% Building of generalized plant -20
Pn = P.nominal;
systemnames = 'Pn Ws Wt Weps'; -30
inputvar = '[w1; w2; u]';
outputvar = '[Ws; Wt; Pn+Weps]'; -40
input_to_Pn = '[w1 - u]';
input_to_Ws = '[w1 - u]'; -50
input_to_Wt = '[ u ]';
input_to_Weps = '[ w2 ]'; -60
G = sysic; 10
0
10
1
10
2
10
3
Frequency (rad/s)
H. Bevrani University of Kurdistan 10
H Control Design (MATLAB Codes)
%% Hinf Control Design using Mixed sensitivity problem
clear all; close all;
rng('default'); % Initializing random numbers
%% Parameter definition
m1 = 0.8; m2 = 0.2; k1 = 100; c1 = 1; c2 = 0.3; Ks = 100; k2 = ureal('k2',300,'percent',20);
%% Define the M, K, C matrices of the motion equation
M = [ m1, 0 ; 0, m2 ]; C = [ c1+c2, -c2 ; -c2, c2 ]; K = [ k1+k2, -k2 ; -k2, k2 ]; F = [ Ks ; 0 ];
%% State space realization
iM = inv(M);
Ap = [ zeros(2,2), eye(2,2) ; -iM*K, -iM*C ]; Bp = [ zeros(2,1) ; iM*F ]; Cp = [ 0 1 0 0 ]; Dp = 0;
P = ss(Ap,Bp,Cp,Dp);
%% Multiplicative uncertainty model
w = logspace(0,3,100);
P_g = ufrd(P,w); % Frequency response calculation
Dm_g = (P_g - P_g.nominal)/P_g.nominal;
s = tf('s');
Wt = 3*s^2/(s^2+18*s+45^2); %Uncertainty weighting function Wt
Ws = 15/(s + 1.5e-2); % Performance weighting function
Weps = 5e-4; % Weps
%% Building of generalized plant
Pn = P.nominal;
systemnames = 'Pn Ws Wt Weps';
inputvar = '[w1; w2; u]';
outputvar = '[Ws; Wt; Pn+Weps]';
input_to_Pn = '[w1 - u]';
input_to_Ws = '[w1 - u]';
input_to_Wt = '[ u ]';
input_to_Weps = '[ w2 ]';
G = sysic;
%% H-inf Controller Design
[K,clp,gamma_min,hinf_info] = hinfsyn(G,1,1,'display','on')
H. Bevrani University of Kurdistan 11
Synthesis Results
%% H-inf Controller Design
[K,clp,gamma_min,hinf_info] = hinfsyn(G,1,1,'display','on')
Test bounds: 0.55 <= gamma <= 1.51
gamma X>=0 Y>=0 rho(XY)<1 p/f
9.10e-01 -6.6e-16 -1.2e-13 1.708e+00 # f
1.17e+00 0.0e+00 1.5e-16 3.938e-01 p
1.03e+00 0.0e+00 1.2e-15 6.800e-01 p
9.69e-01 0.0e+00 1.0e-15 9.902e-01 p
9.39e-01 1.6e-19 -2.1e-14 1.260e+00 # f
9.54e-01 -8.6e-16 -9.6e-14 1.110e+00 # f
9.61e-01 -3.3e-16 -1.6e-14 1.047e+00 # f
Limiting gains...
9.71e-01 0.0e+00 2.4e-15 9.760e-01 p
9.71e-01 0.0e+00 0.0e+00 9.761e-01 p
Best performance (actual): 0.971
H. Bevrani University of Kurdistan 12
Controller Characteristics
Bode Diagram
40
K
P.nominal
30
%% Controller gain plot
figure(1); 20
w = logspace(0,2,100); % Definition of frequency vector
bodemag(K,P.nominal,'--',w); 10
legend('K','P.nominal');
Magnitude (dB)
0
-10
-20
-30
-40
-50
-60
0 1 2
10 10 10
Frequency (rad/s)
H. Bevrani University of Kurdistan 13
Frequency Response Evaluation
%% Evaluation of closed-loop characteristics
L = P*K;
T = feedback(L,1); % T = L/(1+L)
Bode Diagram
S = feedback(1,L); % S = 1/(1+L) 30
M = feedback(P,K); % M = P/(1+L) T.nomial
figure(2); 1/Wt
bodemag(T.nominal,'-',1/Wt,':',S.nominal,'--',1/Ws,'-.',w); S.nominal
1/Ws
20
legend('T.nomial','1/Wt','S.nominal','1/Ws');
ylim([-30 30]);
>> pole(T) 10
Magnitude (dB)
ans =
0
1.0e+03 *
-2.4272 + 0.0000i
-0.0522 + 0.1324i -10
-0.0522 - 0.1324i
-0.1257 + 0.0545i
-0.1257 - 0.0545i -20
-0.0011 + 0.0436i
-0.0011 - 0.0436i
-0.0180 + 0.0180i -30
10 0 10 1 10 2
-0.0180 - 0.0180i Frequency (rad/s)
-0.0005 + 0.0099i
-0.0005 - 0.0099i
H. Bevrani University of Kurdistan 14
Closed-Loop Time Response Evaluation
figure(3); step(T,2);
figure(4); subplot(211),impulse(M,2); subplot(212),impulse(Pn,2);
H. Bevrani University of Kurdistan 15
Project: Report 7
Consider your dynamic system :
Using uncertainty and performance weighing functions (W2 and W1)
in the previous project steps, Design an H∞ controller.
Deadline: The day before next Meeting
Please only use this email address: bevranih18@gmail.com
H. Bevrani University of Kurdistan 16
Thank You!
H. Bevrani University of Kurdistan 17
