Documentation |
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_{∞} design techniques to design a controller for a nominal quarter-car model. Then, you use μ 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.
On this page… |
---|
Quarter-Car Suspension Model
This example uses the quarter-car model of the following illustration to design active suspension control laws.
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.
$$\begin{array}{l}{\dot{x}}_{1}={x}_{2},\\ {\dot{x}}_{2}=-\frac{1}{{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}=\frac{1}{{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.
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'))
ans = -0.0000 +56.2731i -0.0000 -56.2731i
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'))
ans = 0.0000 +22.9734i 0.0000 -22.9734i
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)=\frac{1}{\frac{1}{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],'SampleStateDim',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_{∞} 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.
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}, which is shaped by a weighting function W_{road}.
Sensor noise on both measurements. This noise is modeled as normalized signals, d_{2} and d_{3}, which are shaped by weighting functions W_{d2} and W_{d3}.
You can reinterpret the control objectives as a disturbance rejection goal. 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_{∞} 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_{∞} 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, (βW_{sd},(1 – β)W_{ab}). These weights use a blending parameter, β, 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 (β = 0.01), balancing comfort and handling (β = 0.5), and emphasizing handling (β = 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, β. 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_{∞} controller for each value of the blending parameter, β. 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_{∞} norms, gamma.
gamma
gamma = 0.9410 0.6724 0.8877
The three H_{∞} controllers achieve closed-loop H_{∞} 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)=\{\begin{array}{c}a\left(1-\mathrm{cos}8\pi t\right),\text{\hspace{1em}}0\le t\le 0.25\\ 0,\text{\hspace{1em}}\text{otherwise}\text{.}\text{\hspace{1em}}\text{\hspace{1em}}\text{\hspace{1em}}\text{\hspace{1em}}\text{\hspace{1em}}\end{array}$$
This signal corresponds to a road bump of height 5 cm.
Create a vector that represents the road disturbance.
t = 0:0.005:1; roaddist = zeros(size(t)); roaddist(1:51) = 0.025*(1-cos(8*pi*t(1:51)));
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)') 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)') legend('Open-loop','Comfort','Balanced','Suspension','Road Disturbance','location','SouthEast')
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_{∞} 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 μ-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 β = 0.5.
Examine the resulting μ-synthesis controller.
size(Krob)
State-space model with 1 outputs, 2 inputs, and 11 states.
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)') legend('Open-loop','Robust design','location','SouthEast')
These responses are similar to those obtained with the balanced H_{∞} 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_{∞} and the robust controllers, to compare the results.
Compute an uncertain closed-loop model with the balanced H_{∞} 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 robustperf 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);
Krob is a model array containing a reduced-order controller of every order from 1 up to the original 11 states.
Compute the robust performance margin for each reduced controller.
CLP = lft(qcaric(:,:,2),Kred); ropt = robustperfOptions('Sensitivity','off','Display','off','Mussv','a'); PM = robustperf(CLP,ropt);
Compare the robust performance of the reduced- and full-order controllers.
plot(StateOrders,[PM.LowerBound],'b-o',... StateOrders,repmat(1/RPmuval,[1 NS]),'r'); title('Robust performance as a function of controller order') legend('Reduced order','Full order','Location','SouthEast')
There is no significant difference in robust performance between the 8th- and 11th-order controllers. Therefore, 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.