www.gusucode.com > robust 案例源码程序 matlab代码 > robust/ActiveSuspensionControlDesignExample.m
%% Active Suspension Control Design
% This example shows how to use robust control techniques to design an active
% suspension system for a quarter car body and wheel assembly model. In
% this example, you use ${H_\infty }$ design techniques to design
% a controller for a nominal quarter-car model. Then, you use ${\mu}$
% synthesis to design a robust controller that accounts for uncertainty
% in the model.
%
% Conventional _passive suspensions_ use a spring and damper between the
% car body and wheel assembly. The spring-damper characteristics are adjusted
% to emphasize one of several conflicting objectives such as passenger comfort,
% road holding, and suspension deflection. _Active suspensions_ use a feedback-controller
% hydraulic actuator between the chassis and wheel assembly, which allows
% the designer to better balance these objectives.
%% Quarter-Car Suspension Model
% This example uses the quarter-car model of the following illustration
% to design active suspension control laws.
%
%
% <<../hinfdk1.png>>
%
% The mass $m_{b}$ represents the car chassis (body) and the mass $m_{w}$
% represents the wheel assembly. The spring $k_{s}$ and damper $b_{s}$
% represent the passive spring and shock absorber placed between the car
% body and the wheel assembly. The spring $k_{t}$ models the compressibility
% of the pneumatic tire. The variables $x_{b}$, $x_{w}$, and $r$ are the
% car body travel, the wheel travel, and the road disturbance, respectively.
% The force $f_{s}$, which is applied between the body and wheel assembly,
% is controlled by feedback. This force represents the active component
% of the suspension system.
%
% The following state-space equations describe the quarter-car dynamics.
%
% $$\renewcommand{\arraystretch}{1.5}
% \begin{array}{rcl}
% {{\dot x}_1} &=& {x_2},\\
% {{\dot x}_2} &=& - {1 \over {{m_b}}}\left[ {{k_s}\left( {{x_1} - {x_3}} \right) + {b_s}\left( {{x_2} - {x_4}} \right) - {f_s}} \right],\\
% {{\dot x}_3} &=& {x_4}, \\
% {{\dot x}_4} &=& {1 \over {{m_w}}}\left[ {{k_s}\left( {{x_1} - {x_3}} \right) + {b_s}\left( {{x_2} - {x_4}} \right) - {k_t}\left( {{x_3} - r} \right) - {f_s}} \right].
% \end{array}$$
%
% The state variables in the system are defined as ${x_1}: = {x_b}$,
% ${x_2}: = {\dot x_b}$, ${x_3}: = {x_w}$, and ${x_4}: = {\dot x_w}$.
%%
% Define the physical parameters of the system. For this example, use the
% following values.
mb = 300; % kg
mw = 60; % kg
bs = 1000; % N/m/s
ks = 16000 ; % N/m
kt = 190000; % N/m
%%
% Use these equations and parameter values to construct a state-space model,
% |qcar|, of the quarter-car suspension system.
A = [ 0 1 0 0; ...
[-ks -bs ks bs]/mb ; ...
0 0 0 1; ...
[ks bs -ks-kt -bs]/mw];
B = [0 0; 0 10000/mb ; 0 0; [kt -10000]/mw];
C = [1 0 0 0; 1 0 -1 0; A(2,:)];
D = [0 0; 0 0; B(2,:)];
qcar = ss(A,B,C,D);
qcar.StateName = {'body travel xb (m)';'body vel (m/s)';...
'wheel travel xw (m)';'wheel vel (m/s)'};
qcar.InputName = {'r';'fs'};
qcar.OutputName = {'xb';'sd';'ab'};
%%
% The model inputs are the road disturbance, $r$, and actuator force,
% $f_{s}$. The model outputs are the car body travel, $x_{b}$, suspension
% deflection ${s_d} = {x_b} - {x_w}$, and car body acceleration
% ${a_b} = {{\ddot x}_s}$.
%
%%
% The transfer function from actuator to body travel and acceleration has
% an imaginary-axis zero. Examine the zero of this transfer function.
tzero(qcar({'xb','ab'},'fs'))
%%
% The natural frequency of this zero, 56.27 rad/s, is called the _tire-hop
% frequency_.
%%
% The transfer function from the actuator to suspension deflection also
% has an imaginary-axis zero. Examine this zero.
zero(qcar('sd','fs'))
%%
% The natural frequency of this zero, 22.97 rad/s, is called the _rattlespace
% frequency_.
%%
% Plot the frequency response of the quarter-car model from inputs, $(r,f_{s})$,
% to outputs, $(a_{b},s_{d})$. Both zeros are visible on the Bode plot.
bodemag(qcar({'ab','sd'},'r'),'b',qcar({'ab','sd'},'fs'),'r',{1 100});
legend('Road disturbance (r)','Actuator force (fs)','Location','southwest')
title({'Gain from road dist (r) and actuator force (fs)';...
'to body accel (ab) and suspension travel (sd)'})
%%
% Road disturbances influence the motion of the car and suspension:
%
% * Small body acceleration influences passenger comfort.
% * Small suspension travel contributes to good road handling. Further,
% limits on the actuator displacement constrain the allowable travel.
%%
% Because of the imaginary axis zeros, feedback control cannot improve the
% response from road disturbance ($r$) to body acceleration ($a_{b}$) at
% the tire-hop frequency. Similarly, feedback control cannot improve the
% response from $r$ to suspension deflection ($s_{d}$) at the rattlespace
% frequency. Moreover, there is an inherent trade-off between passenger
% comfort and suspension deflection. Any reduction of body travel at low
% frequency results in an increase of suspension deflection. This trade-off
% arises because of the relationship $x_{w} = x_{b} - s_{d}$ and the fact
% that $x_{w}$ roughly follows $r$ at low frequency (less than 5 rad/s).
%% Hydraulic Actuator Model
% The hydraulic actuator used for active suspension control is connected
% between the body mass, $m_{b}$, and the wheel assembly mass, $m_{w}$.
% The nominal actuator dynamics can be represented by the first-order transfer
% function:
%
% $$ActNom\left( s \right) = {1 \over {{1 \over {60}}s + 1}}.$$
%
%
% The maximum displacement is 0.05 m.
%
%%
% The nominal actuator model approximates the physical actuator dynamics.
% You can model variations between the actuator model and the physical device
% as a family of actuator models. You can also use this approach to model
% variations between the passive quarter-car model and actual vehicle dynamics.
% The resulting family of models comprises a nominal model with a frequency-dependent
% amount of uncertainty.
%
% Create an uncertain model that represents this family of models.
ActNom = tf(1,[1/60 1]);
Wunc = makeweight(0.40,15,3);
unc = ultidyn('unc',[1 1],'SampleStateDimension',5);
Act = ActNom*(1 + Wunc*unc);
Act.InputName = 'u';
Act.OutputName = 'fs';
%%
% At low frequency, below 3 rad/s, the model can vary up to 40% from its
% nominal value. Around 3 rad/s, the percentage variation starts to increase.
% The uncertainty crosses 100% at 15 rad/s, and reaches 2000% at approximately
% 1000 rad/sec. The weighting function, |Wunc|, reflects this profile and
% is used to modulate the amount of uncertainty as a function of frequency.
% The result |Act| is an uncertain state-space model of the actuator.
%%
% Examine the uncertain actuator model by plotting the frequency response
% of 20 randomly sampled models from |Act|.
bodeplot(Act,'b',Act.NominalValue,'r+',logspace(-1,3,120))
title('Nominal and 20 random actuator models')
%%
% The plus (+) marker denotes the nominal actuator model. The blue solid
% lines represent the randomly sampled models.
%
%% Design Objectives for H-Infinity Synthesis
% To use ${H_\infty }$ synthesis algorithms, you must express your design
% objectives as a single cost function to be minimized. For the quarter-car
% model, the main control objectives are formulated in terms of passenger
% comfort and road handling. These objectives relate to body acceleration,
% $a_{b},$ and suspension travel, $s_{d}$. Other factors that influence the
% control design include:
%
% * Characteristics of the road disturbance
% * Quality of the sensor measurements for feedback
% * Limits on the available control force
%
% Use weights to model external disturbances and quantify the design objectives,
% as shown in the following diagram.
%
% <<../hincontroldes82.png>>
%
% The feedback controller uses the measurements $y_{1}$ and $y_{2}$ of the
% suspension travel, $s_{d}$, and body acceleration, $a_{b}$, to compute
% the control signal, $u$. This control signal drives the hydraulic actuator.
% There are three external sources of disturbance:
%
% * The road disturbance, $r$, which is modeled as a normalized signal
% $d_{1}$ shaped by a weighting function $W_{road}$.
%
% * Sensor noise on both measurements. This noise is modeled as normalized
% signals $d_{2}$ and $d_{3}$ shaped by weighting functions
% $W_{d2}$ and $W_{d3}$.
%
% You can reinterpret the control objectives as a disturbance rejection
% goal. In this interpretation, the goal is to minimize the impact of the
% disturbances, $d_{1}$, $d_{2}$, and $d_{3}$, on a weighted combination of
% suspension travel $(s_{d})$, body acceleration $(a_{b})$, and control
% effort $(u)$. You can consider the ${H_\infty }$ norm (peak gain) as the
% measure of the effect of the disturbances. Then, you can meet the
% requirements by designing a controller that minimizes the ${H_\infty }$
% norm from the disturbance inputs, $d_{1}$, $d_{2}$, and $d_{3}$, to the
% error signals, $e_{1}$, $e_{2}$, and $e_{3}$.
%
%%
% Create the weighting functions that model the design objectives.
Wroad = ss(0.07);
Wroad.u = 'd1';
Wroad.y = 'r';
Wact = 8*tf([1 50],[1 500]);
Wact.u = 'u';
Wact.y = 'e1';
Wd2 = ss(0.01);
Wd2.u = 'd2';
Wd2.y = 'Wd2';
Wd3 = ss(0.5);
Wd3.u = 'd3';
Wd3.y = 'Wd3';
%%
% The constant weight |Wroad = 0.07| models broadband road deflections of
% magnitude 7 cm. |Wact| is a highpass filter. This filter penalizes high-frequency
% content in the control signal, and thus limits control bandwidth. |Wd2|
% and |Wd3| model broadband sensor noise of intensity 0.01 and 0.5, respectively.
% In a more realistic design, |Wd2| and |Wd3| would be frequency dependent
% to model the noise spectrum of the displacement and acceleration sensors.
% The inputs and outputs of all weighting functions are named to facilitate
% interconnection. The notation |u| and |y| are shorthand for the |InputName|
% and |OutputName| properties, respectively.
%%
% Specify target functions for the closed-loop response of the system from
% the road disturbance, $r$, to the suspension deflection, $s_{d}$, and
% body acceleration, $a_{b}$.
HandlingTarget = 0.04 * tf([1/8 1],[1/80 1]);
ComfortTarget = 0.4 * tf([1/0.45 1],[1/150 1]);
Targets = [HandlingTarget; ComfortTarget];
%%
% Because of the actuator uncertainty and imaginary-axis zeros, the targets
% attenuate disturbances only below 10 rad/s. These targets represent the
% goals of passenger comfort (small car body acceleration) and adequate
% road handling (small suspension deflection).
%%
% Plot the closed-loop targets and compare to the open-loop response.
bodemag(qcar({'sd','ab'},'r')*Wroad,'b',Targets,'r--',{1,1000})
grid, title('Response to road disturbance')
legend('Open-loop','Closed-loop target')
%
%%
% The corresponding performance weights $W_{sd}$ and $W_{ab}$ are the reciprocals
% of the comfort and handling targets. To investigate the trade-off between
% passenger comfort and road handling, construct three sets of weights,
% $(\beta {W_{sd}},(1 - \beta ){W_{ab}})$. These weights use
% a blending parameter, $\beta$, to modulate the trade-off.
beta = reshape([0.01 0.5 0.99],[1 1 3]);
Wsd = beta/HandlingTarget;
Wsd.u = 'sd';
Wsd.y = 'e3';
Wab = (1-beta)/ComfortTarget;
Wab.u = 'ab';
Wab.y = 'e2';
%%
% |Wsd| and |Wab| are arrays of weighting functions that correspond to three
% different trade-offs: emphasizing comfort ($\beta$ = 0.01),
% balancing comfort and handling ($\beta$ = 0.5), and emphasizing
% handling ($\beta$ = 0.99).
%%
% Connect the quarter-car plant model, actuator model, and weighting functions
% to construct the block diagram of the plant model weighted by the objectives.
sdmeas = sumblk('y1 = sd + Wd2');
abmeas = sumblk('y2 = ab + Wd3');
ICinputs = {'d1';'d2';'d3';'u'};
ICoutputs = {'e1';'e2';'e3';'y1';'y2'};
qcaric = connect(qcar(2:3,:),Act,Wroad,Wact,Wab,Wsd,Wd2,Wd3,...
sdmeas,abmeas,ICinputs,ICoutputs);
%%
% |qcaric| is an array of three models, one for each value of the blending
% parameter, $\beta$. Also, the models in |qcaric| are uncertain,
% because they contain the uncertain actuator model Act.
%% Nominal H-Infinity Synthesis
% Use |hinfsyn| to compute an ${H_\infty }$ controller for each
% value of the blending parameter, $\beta$. |hinfsyn| ignores
% the uncertainty in the plant models and synthesizes a controller for the
% nominal value of each model.
ncont = 1;
nmeas = 2;
K = ss(zeros(ncont,nmeas,3));
gamma = zeros(3,1);
for i=1:3
[K(:,:,i),~,gamma(i)] = hinfsyn(qcaric(:,:,i),nmeas,ncont);
end
%%
% The weighted plant model has one control input (|ncont|), the hydraulic
% actuator force. The model also has two measurement outputs (|nmeas|),
% which give the suspension deflection and body acceleration.
%%
% Examine the resulting closed-loop ${H_\infty }$ norms, |gamma|.
gamma
%%
% The three ${H_\infty }$ controllers achieve closed-loop ${H_\infty }$
% norms of 0.94 (emphasizing comfort), 0.67 (balancing comfort and handling),
% and 0.89 (emphasizing handling).
%%
% Construct closed-loop models of the quarter-car plant with the synthesized
% controller, corresponding to each of the three blending parameter values.
% Compare the frequency response from the road disturbance to $x_{b}$, $s_{d}$,
% and $a_{b}$ for the passive and active suspensions.
K.u = {'sd','ab'}; K.y = 'u';
CL = connect(qcar,Act.Nominal,K,'r',{'xb';'sd';'ab'});
clf
bodemag(qcar(:,'r'),'b', CL(:,:,1),'r-.', ...
CL(:,:,2),'m-.', CL(:,:,3),'k:',{1,140})
grid
legend('Open-loop','Comfort','Balanced','Handling','location','SouthEast')
title('Body travel, suspension deflection, and body acceleration due to road')
%%
% The solid blue line corresponds to the open-loop response. The other lines
% are the closed-loop frequency responses for the different comfort and
% handling blends. All three controllers reduce suspension deflection and
% body acceleration below the rattlespace frequency (23 rad/s).
%% Time-Domain Evaluation
% To further evaluate the three designs, perform time-domain simulations
% using the following road disturbance signal $r(t)$:
%
% $$r\left( t \right) = \left\{ {\matrix{
% {a\left( {1 - \cos 8\pi t} \right),\quad 0 \le t \le 0.25} \cr
% {0,\quad {\rm{otherwise}}{\rm{.}}\quad \quad \quad \quad \quad } \cr
% } } \right.$$
%
%%
% This signal corresponds to a road bump of height 5 cm.
%%
% Create a vector that represents the road disturbance.
t = 0:0.0025:1;
roaddist = zeros(size(t));
roaddist(1:101) = 0.025*(1-cos(8*pi*t(1:101)));
%%
% Build the closed-loop model using the synthesized controller.
SIMK = connect(qcar,Act.Nominal,K,'r',{'xb';'sd';'ab';'fs'});
%%
% |SIMK| is a model array containing three closed-loop models, one for each
% of the three blending parameter values. Each model in the array represents
% a closed loop that is built from the original quarter-car plant model,
% the nominal actuator model, and the corresponding synthesized controller.
%%
% Simulate and plot the time-domain response of the closed-loop models to
% the road disturbance signal.
p1 = lsim(qcar(:,1),roaddist,t);
y1 = lsim(SIMK(1:4,1,1),roaddist,t);
y2 = lsim(SIMK(1:4,1,2),roaddist,t);
y3 = lsim(SIMK(1:4,1,3),roaddist,t);
clf
subplot(221)
plot(t,p1(:,1),'b',t,y1(:,1),'r.',t,y2(:,1),'m.',t,y3(:,1),'k.',t,roaddist,'g')
title('Body travel')
ylabel('x_b (m)')
subplot(222)
plot(t,p1(:,3),'b',t,y1(:,3),'r.',t,y2(:,3),'m.',t,y3(:,3),'k.',t,roaddist,'g')
title('Body acceleration')
ylabel('a_b (m/s^2)')
% configure legend
h = legend('Open-loop','Comfort','Balanced','Suspension','Road Dist.',...
'Location','southeast');
h.FontSize = 7;
h.Box = 'off';
h.Position = [0.68, 0.58, 0.26, 0.2];
subplot(223)
plot(t,p1(:,2),'b',t,y1(:,2),'r.',t,y2(:,2),'m.',t,y3(:,2),'k.',t,roaddist,'g')
title('Suspension deflection')
xlabel('Time (s)')
ylabel('s_d (m)')
subplot(224)
plot(t,zeros(size(t)),'b',t,y1(:,4),'r.',t,y2(:,4),'m.',t,y3(:,4),'k.',t,roaddist,'g')
title('Control force')
xlabel('Time (s)')
ylabel('f_s (N)')
%%
% The simulations show that the body acceleration is smallest for the controller
% emphasizing passenger comfort. Body acceleration is largest for the controller
% emphasizing suspension deflection. The balanced design achieves a good
% tradeoff between body acceleration and suspension deflection.
%% Robust μ Design
% So far you designed ${H_\infty }$ controllers that meet the performance
% objectives for the nominal actuator model. However, this model is only
% an approximation of the true actuator. To make sure that controller performance
% is maintained even with model error and uncertainty, you must design the
% model to have _robust performance_. In this part of the example, you use
% $\mu$-synthesis to design a controller that achieves robust
% performance for the entire family of actuator models that takes uncertainty
% into account.
%%
% Use D-K iteration to synthesize a controller for the quarter-car model
% with actuator uncertainty.
[Krob,~,RPmuval] = dksyn(qcaric(:,:,2),nmeas,ncont);
%%
% The model |qcaric(:,:,2)| is the weighted quarter-car model for the uncertain
% model that corresponds to the blending variable $\beta$ = 0.5.
%%
% Examine the resulting $\mu$-synthesis controller.
size(Krob)
%%
% Build the closed-loop model using the robust controller, |Krob|.
Krob.u = {'sd','ab'};
Krob.y = 'u';
SIMKrob = connect(qcar,Act.Nominal,Krob,'r',{'xb';'sd';'ab';'fs'});
%%
% Simulate and plot the nominal time-domain response to a road bump with
% the robust controller.
p1 = lsim(qcar(:,1),roaddist,t);
y1 = lsim(SIMKrob(1:4,1),roaddist,t);
clf
subplot(221)
plot(t,p1(:,1),'b',t,y1(:,1),'r',t,roaddist,'g')
title('Body travel'), ylabel('x_b (m)')
subplot(222)
plot(t,p1(:,3),'b',t,y1(:,3),'r')
title('Body acceleration'), ylabel('a_b (m/s^2)')
subplot(223)
plot(t,p1(:,2),'b',t,y1(:,2),'r')
title('Suspension deflection'), xlabel('Time (s)'), ylabel('s_d (m)')
subplot(224)
plot(t,zeros(size(t)),'b',t,y1(:,4),'r')
title('Control force'), xlabel('Time (s)'), ylabel('f_s (N)')
h = legend('Open-loop','Robust design','Location','southeast');
h.Box = 'off';
h.FontSize = 8;
%%
% These responses are similar to those obtained with the balanced ${H_\infty }$
% controller.
%%
% Examine the effect of the robust controller on variability caused by model
% uncertainty. To do so, simulate the response to a road bump for 120 actuator
% models randomly sampled from the uncertain model, |Act|. Perform this
% simulation for both the ${H_\infty }$ and the robust controllers,
% to compare the results.
%
% Compute an uncertain closed-loop model with the balanced ${H_\infty }$
% controller, |K|. Sample this model, simulate the sampled models, and plot
% the results.
CLU = connect(qcar,Act,K(:,:,2),'r',{'xb','sd','ab'});
nsamp = 120;
rng('default');
figure(1)
clf
lsim(usample(CLU,nsamp),'b',CLU.Nominal,'r.',roaddist,t)
title('Nominal "balanced" design')
%%
% Compute an uncertain closed-loop model with the balanced robust controller,
% |Krob|. Sample this model, simulate the sampled models, and plot the results.
CLU = connect(qcar,Act,Krob,'r',{'xb','sd','ab'});
figure(2)
clf
lsim(usample(CLU,nsamp),'b',CLU.Nominal,'r.',roaddist,t)
title('Robust "balanced" design')
%%
% The robust controller reduces variability caused by model uncertainty,
% and delivers more consistent performance.
%% Controller Simplification
% The robust controller |Krob| has eleven states. It is often the case that
% controllers synthesized with |dksyn| have high order. You can use the
% model reduction functions to find a lower-order controller that achieves
% the same level of robust performance. Use |reduce| to generate approximations
% of various orders. Then, use |robgain| to compute the robust performance
% margin for each reduced-order approximation.
%%
% Create an array of reduced-order controllers.
NS = order(Krob);
StateOrders = 1:NS;
Kred = reduce(Krob,StateOrders);
%%
% |Kred| is a model array containing a reduced-order controller of every
% order from 1 up to the original 11 states.
%%
% Next use |robgain| to compute the robust performance margin
% for each reduced-order approximation. The performance goals are met when
% the closed-loop gain is less than $\gamma=1$. The robust performance
% margin measures how much uncertainty can be sustained without
% degrading performance (exceeding $\gamma=1$). A margin of 1
% or more indicates that we can sustain 100% of the specified uncertainty.
% Compute robust performance margin for each reduced controller
gamma = 1;
CLP = lft(qcaric(:,:,2),Kred);
for k=1:NS
PM(k) = robgain(CLP(:,:,k),gamma);
end
% Compare robust performance of reduced- and full-order controllers
PMfull = PM(end).LowerBound;
plot(StateOrders,[PM.LowerBound],'b-o',...
StateOrders,repmat(PMfull,[1 NS]),'r');
title('Robust performance margin as a function of controller order')
legend('Reduced order','Full order','location','SouthEast')
%%
% The robust performance margin is well below 1 for controllers of
% order 7 and lower. However
% there is no significant difference in performance margin between the 8th-
% and 11th-order controllers, so you can safely replace |Krob| by its
% 8th-order approximation.
%%
Krob8 = Kred(:,:,8);
%%
% You now have a simplified controller, |Krob8|, that provides robust control
% with a balance between passenger comfort and handling.