rrdesign.m

%% 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.

%%

% 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
%%
% The design objective is defined as an M-file function myCostFcn that
% accepts two inputs: the design vector x and simParms.  x
% contains our design variables for the suspension system. simParms is a
% structure that passes in the remaining defining parameters of the
% Simulink model (Mb, Lf, Lr, and Iyy).  myCostFcn runs the suspension model
% defined by x and simParms and returns a measure of passenger discomfort,
% calculated as the weighted average of the peak and total acceleration, as
% shown below.  Passenger discomfort is normalized so that our
% initial design has a discomfort level of 1.

% close original simulink model
close_system

type myCostFcn

%%
% Nonlinear constraints are defined in the M-file function myNonlCon that
% returns values for c(x) and ceq(x).  
type myNonlCon

%%
% 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')
%%
% 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
%%
% 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
%
%
% $$\eta = {c\over{2\sqrt{k M}}}$$
% 
% 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')

%%
% 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
%% 
% 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
%% Optimizing for Robustness
% Our design is now reliable--it meets our life and design goalsbut 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.  Well 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 = 100;%00;
% 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)]')

%%
% 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
% <http://www.mathworks.com/products/matlab MATLAB>
%
% <http://www.mathworks.com/products/simulink Simulink>
%
% <http://www.mathworks.com/products/optimization Optimization Toolbox>
%
% <http://www.mathworks.com/products/statistics Statistics Toolbox> 

%% Resources
% <http://www.mathworks.com/company/newsletters/digest/2007/may/uncertainity.html Using Statistics to Analyze Uncertainty in System Models>
%
% <http://www.mathworks.com/company/newsletters/digest/2006/july/cooling.html Improving an Engine Cooling Fan Using Design for Six Sigma Techniques>
%
% <http://www.mathworks.com/company/newsletters/digest/nov04/modeling.html Modeling Survival Data with Statistics Toolbox>