rrdesign.html

全面系统地给出鲁棒设计的源程序及运算实例

myConst = @(x) myNonlConRR(x,s);
myObjFcn = @(x) myCostFcnRR(x,s);

<span class="comment">% Run optimization</span>
tic
[xRR,costRR] = fmincon(myObjFcn,xr,A,b,Aeq,beq,lb,ub,<span class="keyword">...</span>
    myConst,options);
toc

<span class="comment">% Plot results with mean solution for robust case</span>
rrplot([cost0 cost costr costRR],[x0(f); x(f); xr(f); xRR(f)]',<span class="keyword">...</span>
       [x0(r); x(r); xr(r); xRR(r)]')
</pre><pre class="codeoutput">
                                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.
</pre><img vspace="5" hspace="5" src="rrdesign_09.png"> <img vspace="5" hspace="5" src="rrdesign_10.png"> <p>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.
         </p>
         <p>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.
         </p><pre class="codeinput">figure
cost = totAccel/904 * 0.5 + pkAccel/56 * 0.5;
Xdata = [front back trunk cost'];
varNames = {<span class="string">'Front (kg)'</span>,<span class="string">'Back (kg)'</span>,<span class="string">'Trunk (kg)'</span>,{<span class="string">'Normalized'</span>, <span class="string">'Discomfort'</span>}};
[h,ax,bax] = plotmatrix(Xdata);
<span class="keyword">for</span> ii = 1:length(varNames)
    xlabel(ax(end,ii),varNames{ii})
    ylabel(ax(ii,1),varNames{ii})
<span class="keyword">end</span>
</pre><img vspace="5" hspace="5" src="rrdesign_11.png"> <p>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.
         </p><pre class="codeinput"><span class="comment">% final figure generated from dfittool</span>
cumprobPlot(cost)
xlabel(<span class="string">'Normalized Discomfort'</span>)
</pre><img vspace="5" hspace="5" src="rrdesign_12.png"> <h2>Design Trade-Offs<a name="19"></a></h2>
         <p>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.
         </p>
         <p>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.
         </p>
         <p>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.
         </p>
         <h2>Products used<a name="20"></a></h2>
         <p><a href="http://www.mathworks.com/products/matlab">MATLAB</a></p>
         <p><a href="http://www.mathworks.com/products/simulink">Simulink</a></p>
         <p><a href="http://www.mathworks.com/products/optimization">Optimization Toolbox</a></p>
         <p><a href="http://www.mathworks.com/products/statistics">Statistics Toolbox</a></p>
         <h2>Resources<a name="21"></a></h2>
         <p><a href="http://www.mathworks.com/company/newsletters/digest/2007/may/uncertainity.html">Using Statistics to Analyze Uncertainty in System Models</a></p>
         <p><a href="http://www.mathworks.com/company/newsletters/digest/2006/july/cooling.html">Improving an Engine Cooling Fan Using Design for Six Sigma Techniques</a></p>
         <p><a href="http://www.mathworks.com/company/newsletters/digest/nov04/modeling.html">Modeling Survival Data with Statistics Toolbox</a></p>
         <p class="footer"><br>
            Published with MATLAB&reg; 7.5<br></p>
      </div>
      <!--
##### SOURCE BEGIN #####
%% 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.
%
fprintf('%s',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