Designing for Reliability and Robustness
by Stuart Kozola
Code used in Jan 2008 Digest Article of same title. http://http://www.mathworks.com/company/newsletters/digest/2008/jan/
No design is free from uncertainty or natural variation. For example, how will the design be used? How will it respond to environmental factors or to changes in manufacturing or operational processes? These kinds of uncertainty compound the challenge of creating designs that are reliable and robust –designs that perform as expected over time and are insensitive to changes in manufacturing, operational, or environmental factors.
Using an automotive suspension system as an example, this article describes tools and techniques in MATLAB®, Statistics Toolbox™, and Optimization Toolbox™ software that let you extend a traditional design optimization approach to account for uncertainty in your design, improving quality and reducing prototype testing and overall development effort.
We begin by designing a suspension system that minimizes the forces experienced by front- and rear-seat passengers when the automobile travels over a bump in the road. We then modify the design to account for suspension system reliability; we want to ensure that the suspension system will perform well for at least 100,000 miles. We conclude our analysis by verifying that the design is resilient to, or unaffected by, changes in cargo and passenger mass.
Contents
% verify installed products are available
checkLicense
Performing Traditional Design Optimization
Our Simulink suspension system model has two inputs – the bump starting and ending height – and eight adjustable parameters. We can modify the four parameters that define the front and rear suspension system stiffness and damping rate: kf, kr, cf, cr. The remaining parameters are defined by applying passenger and cargo loading to the vehicle and are not considered to be design parameters.
clc, close all, clear all % Model parameters % Constants s.g = 9.81; % gravity (m/s^2) % Car Geometry s.Lf = 0.9; % front hub displacement from body CG (m) s.Lr = 1.2; % rear hub displacement from body CG (m) s.rf = 0.5*s.Lf; % location of front passengers (m) s.rr = 0.9*s.Lr; % location of rear passengers (m) s.rt = 1.1*s.Lr; % location of trunk (m) s.Mb = 1300; % nominally loaded car mass (kg) (1200 empty) s.Iyy = 2100; % body moment of inertia about y-axis (kgm^2) % Initial suspension design s.kf0 = 19600; % front suspension stiffness (N/m) s.cf0 = 2200; % front suspension damping (N/(m/s)) s.kr0 = 14700; % rear suspension stiffness (N/m) s.cr0 = 2000; % rear suspension damping (N/(m/s)) % Bump height s.hb = 0.1; % Show simulink model and dialog open('mldemo_suspnslow.mdl') open_system('mldemo_suspnslow/Suspension Model')
The model outputs are angular acceleration about the center of gravity (thetadotddot) and vertical acceleration (zdotdot). The Figure below illustrates the model response for our initial design to a simulated bump in the road.
% Initial Design (Run Simulink Model with Initial Design)
f = [1 3];
r = [2 4];
x0 = [s.kf0 s.cf0 s.kr0 s.cr0];
cost0 = myCostFcn_slow(x0,s);
Our goal is to set the parameters kf, kr, cf, and cr to minimize the discomfort that front- and rear-seat passengers experience as a result of traveling over a bump in the road. We use acceleration as a proxy for passenger discomfort. The design optimization problem can be summarized as follows:
Objective: Minimize peak and total acceleration Design variables: Front/rear spring/shock absorber design (kr, kr, cf, cr)
Design constraints: Car is level when at rest.
Suspension system maintains a natural frequency of vibration below 2 Hz.
Damping ratio remains between 0.3 and 0.5.This problem is nonlinear in both response and design constraints. To solve it, a nonlinear optimization solver is required. The Optimization Toolbox solver fmincon is designed specifically for this type of problem.
We begin by casting our optimization problem into the form accepted by fmincon. The table below summarizes the problem formulation that fmincon accepts and the definition of the suspension problem in MATLAB syntax.
suspnHtmlTable
| fmincon Standard Form | Suspension Problem (MATLAB m-code) | |||||||
|---|---|---|---|---|---|---|---|---|
| Objective | min f(x) | myCostFcn(x,simParms) | ||||||
| Design variables | x | x = [kf, cf, kr, cr] | ||||||
| Nonlinear constraints |
|
myNonlCon(x,simParms) | ||||||
| Linear constraints |
|
|
||||||
| Bound constraints | lb <= x <= ub |
|
% close original simulink model close_system type myCostFcn
function mycost = myCostFcn(x,simParms)
%% Unpack simParms structure for Simulink Use
struct2var(simParms)
%% Extract suspension variables
kf = x(1); cf = x(2);
kr = x(3); cr = x(4);
%% Run Simulink model
% Initial Conditions for Simulink Model
theta0 = 0; % initial pitch (rad)
thetadot0 = 0; % initial pitch rate (rad/s)
Z0 = -0.5*Mb*g/(kf+kr); % initial equilibrium position, assumes full car (m)
Zdot0 = 0; % initial bounce rate (m/s)
simTime = [0 8];
sim('mldemo_suspnfast.mdl',simTime,...
simset('SrcWorkspace', 'current', 'DstWorkspace', 'current'))
%% Compute cost
totalAccel = (Zdotdot + rf * thetadotdot).^2 + ...
(Zdotdot - rr * thetadotdot).^2;
%% Return final value
mycost = 0.5*max(totalAccel)/56 + 0.5 * sum(totalAccel)/904;
assignin('base','totAccel',totalAccel);
Nonlinear constraints are defined in the M-file function myNonlCon that returns values for c(x) and ceq(x).
type myNonlCon
function [c,ceq] = myNonlCon(x,simParms)
%% mynonlcon.m Nonlinear constraints for fmincon.
% Extract suspension variables
kf = x(1);
cf = x(2);
kr = x(3);
cr = x(4);
% Extract parameters needed for constraints
struct2var(simParms)
%% Define Desired Damping Coefficient
%cd = 0.40;
% Can tolerate a +/-10% change from desired damping coefficient
cupper = 0.5; % upper limit for damping coefficients
clower = 0.3; % lower limit for damping coefficients
%% Compute Constraints
Mf = Mb*Lr/(Lf+Lr)/2;
Mr = Mb*Lf/(Lf+Lr)/2;
% Inequality constraints c <= 0
c = [sqrt(kf/Mf)/(2*pi)-2;... % fn <= 2 Hz for front
sqrt(kr/Mr)/(2*pi)-2;... % fn <= 2 Hz for rear
cf/(2*sqrt(kf*Mf))-cupper;... % damping ratio for front
clower-cf/(2*sqrt(kf*Mf));...
cr/(2*sqrt(kr*Mf))-cupper;... % damping ratio for rear
clower-cr/(2*sqrt(kr*Mf))];
ceq = [];
The linear and bound constraints are defined as constant coefficient matrices (A, Aeq) or vectors (b, beq, lb, ub).
%Inequality constraints A*x <= b A = []; b = []; % Equality constraints Aeq*x = beq Aeq = [s.Lf 0 -s.Lr 0]; % level car beq = 0; % Set lower and upper bounds lb = [10000; 100; 10000; 100]; ub = [100000; 10000; 100000; 10000];
We defined and solved our problem using the Optimization Tool graphical user interface (optimtool), which simplifies the tasks of defining an optimization problem, choosing an appropriate solver, setting solver options, and running the solver.
load tradOptimProblem optimtool(tradOptimProblem) disp('Hit Start to run solver in optimtool GUI')
Hit Start to run solver in optimtool GUI
Using a traditional optimization approach, we found that the optimal design was one where x = [kf, cf, kr, cr] = [13333, 2225, 10000, 1927].
% Command line equivalent to using optimtool GUI % Set solver options options = optimset; options = optimset(options,'Display' ,'iter'); options = optimset(options,'TolFun' ,1e-08); options = optimset(options,'LargeScale' ,'off'); options = optimset(options,'PlotFcns' ,{ @optimplotx @optimplotfval }); % Define constraint and objective function (to handle parameters) myConst = @(x) myNonlCon(x,s); myObjFcn = @(x) myCostFcn(x,s); % Run Optimization tic [x,cost] = fmincon(myObjFcn,x0,A,b,Aeq,beq,lb,ub,myConst,options); toc % Display results thus far rrplot([cost0 cost],[x0(f); x(f)]',[x0(r); x(r)]') % close optimtool GUI optimGUI = com.mathworks.toolbox.optim.OptimGUI.getOptimGUI; optimGUI.close
Max Line search Directional First-order
Iter F-count f(x) constraint steplength derivative optimality Procedure
0 5 0.997217 0
1 10 0.997217 3.638e-012 1 -2.85e-008 0.000127
2 15 0.910198 3.638e-012 1 -0.0764 9.17e-005 Hessian modified
3 20 0.903431 4.19e-006 1 -0.0067 7.35e-005
4 25 0.903434 1.233e-012 1 3.66e-006 7.43e-005 Hessian modified
5 30 0.903431 1.141e-013 1 -3.55e-006 7.38e-005 Hessian modified
6 35 0.898039 9.111e-006 1 -0.00536 7.81e-005 Hessian modified
7 40 0.898045 9.391e-012 1 5.54e-006 7.79e-005 Hessian modified
8 45 0.898043 3.638e-012 1 -1.24e-006 7.77e-005 Hessian modified
9 50 0.45899 0.008591 1 -0.368 0.00116 Hessian modified
10 55 0.462055 0 1 0.0031 0.000109
11 60 0.462055 0 1 -1.5e-011 0.000109 Hessian modified
Optimization terminated: first-order optimality measure less than options.TolFun
and maximum constraint violation is less than options.TolCon.
Active inequalities (to within options.TolCon = 1e-006):
lower upper ineqlin ineqnonlin
3 3
5
Elapsed time is 4.942308 seconds.
The top figure above shows a standard Optimization Toolbox solution progress plot. The top plot shows the current value of the design variables for the current solver iteration, which at iteration 11 is the final solution. The bottom plot shows the objective function value (passenger discomfort relative to the initial design) across solver iterations. This plot shows that our initial design (iteration 0) had a discomfort level of 1, while the optimal design, found after 11 iterations, has a discomfort level of 0.46 – a reduction of 54% from our initial design.
Ensuring Suspension System Reliability
Our optimal design satisfies the design constraints, but is it a reliable design? Will it perform as expected over a given time period? We want to ensure that our suspension design will perform as intended for at least 100,000 miles.
To estimate the suspension's reliability, we use historical maintenance data for similar suspension system designs (below). The horizontal axis represents the driving time, reported as miles. The vertical axis shows how many suspension systems had degraded performance requiring repair or servicing. The different sets of data apply to suspension systems with different damping ratios. The damping ratio is defined as
where c is the damping coefficient (cf or cr), k is the spring stiffness (kf or kr), and M is the amount of mass supported by the front or rear suspension. The damping ratio is a measure of the stiffness of the suspension system.
We fit Weibull distributions to the historical data using the Distribution Fitting Tool (dfittool). Each fit provides a probability model that we can use to predict our suspension system reliability as a function of miles driven. Collectively, the three Weibull fits let us predict how the damping ratio affects the suspension system reliability as a function of miles driven. For example, the optimal design found previously has a damping ratio for the front and rear suspension of 0.5. Using the plots in dfittool, we can expect that after 100,000 miles of operation, our design will have 88% of the original designs operating without the need for repair. Conversely, 12% of the original designs will require repair before 100,000 miles.
load suspnSurvivalData dfittool(miles3,[],[],'Damping Ratio 0.3') dfittool(miles5,[],[],'Damping Ratio 0.5') dfittool(miles7,[],[],'Damping Ration 0.7') disp('Change the Display Type to Survivor Function in dropdown menu') disp('Then fit Weibull distributions to the data using New Fit Button') % Here is the reproduced final figure from dfittool (generated from % dfittool session). Note you can load the session I used into the % distribution fitting tool under File --> Load Session --> % suspnSurvivalData.dfit survivalFit(miles3,miles5,miles7) xlabel('Miles Driven')
Change the Display Type to Survivor Function in dropdown menu Then fit Weibull distributions to the data using New Fit Button
We want to improve our design so that it has a 90% survival rate at 100,000 miles of operation. We add this reliability constraint to our traditional optimization problem by adding a nonlinear constraint to myNonlCon.
% Close dfittool dffig = dfgetset('dffig'); close(dffig) type myNonlConR
function [c,ceq] = myNonlConR(x,simParms)
%% mynonlcon.m Nonlinear constraints for fmincon.
% Extract suspension variables
kf = x(1);
cf = x(2);
kr = x(3);
cr = x(4);
% Extract parameters needed for constraints
struct2var(simParms)
%% Define Reliability Limit and Desired Damping Coefficient Range
Plimit = 0.90; % max probability of shock absorber failure for lifetime
A = @(dampRatio) -1.0129e+005.*dampRatio.^2 -28805.*dampRatio + 2.1831e+005;
B = @(dampRatio) 1.6865.*dampRatio.^2 -1.8534.*dampRatio + 4.1507;
Ps = @(miles,dampRatio) 1 - wblcdf(miles,A(dampRatio),B(dampRatio));
mileage = 100000;
% Can tolerate a +/-10% change from desired damping coefficient
cupper = 0.5; % upper limit for damping coefficients
clower = 0.3; % lower limit for damping coefficients
%% Compute Constraints
Mf = Mb*Lr/(Lf+Lr)/2;
Mr = Mb*Lf/(Lf+Lr)/2;
cdf = cf/(2*sqrt(kf*Mf));
cdr = cr/(2*sqrt(kr*Mf));
% Inequality constraints c <= 0
c = [sqrt(kf/Mf)/(2*pi)-2;... % fn <= 2 Hz for front
sqrt(kr/Mr)/(2*pi)-2;... % fn <= 2 Hz for rear
clower-cdf;... % damping ratio for front
clower-cdr;... % damping ratio for rear
Plimit- Ps(mileage,cdf);... % front shock absorber reliability constraint
Plimit- Ps(mileage,cdr)]; % rear shock absorber reliability constraint
ceq = [];
We solve the optimization problem as before using optimtool. The results, summarized in the figure below, show that including the reliability constraint changed the design values for cf and cr and resulted in a slightly higher discomfort level. The reliability-based design still performs better than the initial design.
% We'll use the command line equiavlaent to run the solver here to show % results. % Define constraint and objective function (to handle parameters) myConst = @(x) myNonlConR(x,s); myObjFcn = @(x) myCostFcn(x,s); % Run Optimization tic [xr,costr] = fmincon(myObjFcn,x,A,b,Aeq,beq,lb,ub,myConst,options); toc % Plot results thus far rrplot([cost0 cost costr],[x0(f); x(f); xr(f)]',[x0(r); x(r); xr(r)]') s0 = s; % save inital parameter set for later reference
Max Line search Directional First-order
Iter F-count f(x) constraint steplength derivative optimality Procedure
0 5 0.462055 0.01374 Infeasible start point
1 10 0.481347 0.001134 1 0.0206 3.53e+003
2 15 0.483375 9.086e-006 1 0.00204 3.15
3 20 0.483391 5.922e-010 1 1.66e-005 0.000418
4 25 0.483391 1.11e-016 1 -1.46e-009 4.39e-005 Hessian modified
Optimization terminated: magnitude of directional derivative in search
direction less than 2*options.TolFun and maximum constraint violation
is less than options.TolCon.
Active inequalities (to within options.TolCon = 1e-006):
lower upper ineqlin ineqnonlin
5
6
Elapsed time is 1.667676 seconds.
Optimizing for Robustness
Our design is now reliable--it meets our life and design goals—but it may not be robust. Operation of the suspension-system design is affected by changes in the mass distribution of the passengers or cargo loads. To be robust, the design must be insensitive to changes in mass distribution.
To account for a distribution of mass loadings in our design, we use Monte Carlo simulation, repeatedly running the Simulink model for a wide range of mass loadings at a given design. The Monte Carlo simulation will result in the model output having a distribution of values for a given set of design variables. The objective of our optimization problem is therefore to minimize the mean and standard deviation of passenger discomfort.
We replace the single simulation call in myCostFcn with the Monte Carlo simulation and optimize for the set of design values that minimize the mean and standard deviation of total passenger discomfort. We’ll assume that the mass distribution of passengers and trunk loads follow Rayleigh distributions and randomly sample the distributions to define conditions to test our design performance.
NOTE: This published file used nRuns = 10000, which should take about 2-3 hours on a machine with 2.2GHz Intel Processor and 2 GB of RAM to run. After publishing this file nRuns was changed to 100, which you can run in short order and view results and will have results that are slightly different then the article.
nRuns = 10000; % Probability distributions (assumed) front = 40 + raylrnd(40,nRuns,1); % additional mass for front passengers (kg) back = raylrnd(40,nRuns,1); % additional mass for rear passengers/cargo (kg) trunk = raylrnd(10,nRuns,1); % additional mass for luggage (kg)
The total mass, center of gravity, and moment of inertia are adjusted to account for the changes in mass distribution of the automobile.
mcMb = 1200 + front + back + trunk; % Computed as a perturbation from the existing center of mass cm = (front.*s.rf - back.*s.rr - trunk.*s.rt)./mcMb; % a positive cm indicates the CM moved forward % a negative cm indicates the CM moved backwards %Adjust the moment of inertia about the old center of mass mcIyy = s.Iyy + front.*s.rf.^2 + back.*s.rr.^2 + trunk.*s.rt.^2; % Use parallel axis theorem to move moment of inertia to the new CM mcIyy = mcIyy - mcMb.*cm.^2; % Adjust the distance measurements to the new CG mcLf = s.Lf - cm; mcLr = s.Lr + cm; mcrf = s.rf - cm; mcrr = s.rr + cm; mcrt = s.rt + cm; % Put mc parameters in structure to pass to optimization solver s.mcMb = mcMb; s.mcIyy= mcIyy; s.mcLf = mcLf; s.mcLr = mcLr; s.mcrf = mcrf; s.mcrr = mcrr; s.mcrt = mcrt; s.nRuns = length(s.mcMb);
The optimization problem, including the reliability constraints, is solved as before in optimtool. The results are shown below. The robust design has an average acceleration level that is higher than that in the reliability-based design, resulting in a design with higher damping coefficient values. The discomfort measure reported in the figure below is an average value for the robust design case.
% Set solver options options = optimset; options = optimset(options,'Display' ,'Iter'); options = optimset(options,'TolFun' ,1e-08); options = optimset(options,'LargeScale' ,'off'); options = optimset(options,'PlotFcns' ,{ @optimplotx @optimplotfval }); % Define constraint and objective function (to handle parameters) myConst = @(x) myNonlConRR(x,s); myObjFcn = @(x) myCostFcnRR(x,s); % Run optimization tic [xRR,costRR] = fmincon(myObjFcn,xr,A,b,Aeq,beq,lb,ub,... myConst,options); toc % Plot results with mean solution for robust case rrplot([cost0 cost costr costRR],[x0(f); x(f); xr(f); xRR(f)]',... [x0(r); x(r); xr(r); xRR(r)]')
Max Line search Directional First-order
Iter F-count f(x) constraint steplength derivative optimality Procedure
0 5 0.489457 0.001926 Infeasible start point
1 10 0.492796 2.59e-005 1 0.00338 32.5
2 15 0.492842 4.811e-009 1 4.66e-005 0.0123
3 20 0.492842 4.441e-016 1 6.03e-009 4.47e-005 Hessian modified
Optimization terminated: magnitude of directional derivative in search
direction less than 2*options.TolFun and maximum constraint violation
is less than options.TolCon.
Active inequalities (to within options.TolCon = 1e-006):
lower upper ineqlin ineqnonlin
41703
51703
Elapsed time is 8830.742192 seconds.
The figure below displays a scatter-matrix plot that summarizes variability seen in discomfort as a result of different mass loadings for the robust design case.
The diagonals show a histogram of the variables listed on the axis. The plots above and below the diagonal are useful for quickly finding trends across variables. The histograms for front, back, and trunk represent the distribution of the inputs to our simulation. The histogram for discomfort shows that it is concentrated around the value 0.47 and is approximately normally distributed. The plots below the diagonal do not show a strong trend in discomfort with trunk loading levels, indicating that our design is robust to changes in this parameter. There is a definite trend associated with the front loading on discomfort. Front loading appears to be approximately linear with a min of 0.43 and a max of 0.52. A trend between back loading and discomfort can also be seen, but from this plot it is difficult determine if it is linear. From this plot, it is difficult to determine whether our design is robust with respect to back loading.
figure
cost = totAccel/904 * 0.5 + pkAccel/56 * 0.5;
Xdata = [front back trunk cost'];
varNames = {'Front (kg)','Back (kg)','Trunk (kg)',{'Normalized', 'Discomfort'}};
[h,ax,bax] = plotmatrix(Xdata);
for ii = 1:length(varNames)
xlabel(ax(end,ii),varNames{ii})
ylabel(ax(ii,1),varNames{ii})
end
Using a cumulative probability plot of the discomfort results (below), we estimate that 90% of the time, passengers will experience less than 50% of the discomfort they would have experienced with the initial design. We can also see that our new design maintains a normalized level of discomfort below 0.52 nearly 100% of the time. We therefore conclude that our optimized design overall is robust to expected variation in loadings and will perform better than our initial design.
% final figure generated from dfittool cumprobPlot(cost) xlabel('Normalized Discomfort')
Design Trade-Offs
This article showed how MATLAB, Statistics Toolbox, and Optimization Toolbox can be used to capture uncertainty within a simulation-based design problem in order to find an optimal suspension design that is reliable and robust.
We began by showing how to reformulate a design problem as an optimization problem that resulted in a design that performed better than the initial design. We then modified the optimization problem to include a reliability constraint. The results showed that a trade-off in performance was required to meet our reliability goals.
We completed our analysis by including the uncertainty that we would expect to see in the mass loading of the automobile. The results showed that we derived a different design if we accounted for uncertainty in operation and quantified the expected variability in performance. The final design traded performance to maintain reliability and robustness.
Products used
Resources
Using Statistics to Analyze Uncertainty in System Models
Improving an Engine Cooling Fan Using Design for Six Sigma Techniques