📄 rrdesign.html
字号:
% 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 = [];
</pre><p>The linear and bound constraints are defined as constant coefficient matrices (A, Aeq) or vectors (b, beq, lb, ub).</p><pre class="codeinput"><span class="comment">%Inequality constraints A*x <= b</span>
A = [];
b = [];
<span class="comment">% Equality constraints Aeq*x = beq</span>
Aeq = [s.Lf 0 -s.Lr 0]; <span class="comment">% level car</span>
beq = 0;
<span class="comment">% Set lower and upper bounds</span>
lb = [10000; 100; 10000; 100];
ub = [100000; 10000; 100000; 10000];
</pre><p>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.
</p><pre class="codeinput">load <span class="string">tradOptimProblem</span>
optimtool(tradOptimProblem)
disp(<span class="string">'Hit Start to run solver in optimtool GUI'</span>)
</pre><pre class="codeoutput">Hit Start to run solver in optimtool GUI
</pre><p>Using a traditional optimization approach, we found that the optimal design was one where x = [kf, cf, kr, cr] = [13333, 2225,
10000, 1927].
</p><pre class="codeinput"><span class="comment">% Command line equivalent to using optimtool GUI</span>
<span class="comment">% Set solver options</span>
options = optimset;
options = optimset(options,<span class="string">'Display'</span> ,<span class="string">'iter'</span>);
options = optimset(options,<span class="string">'TolFun'</span> ,1e-08);
options = optimset(options,<span class="string">'LargeScale'</span> ,<span class="string">'off'</span>);
options = optimset(options,<span class="string">'PlotFcns'</span> ,{ @optimplotx @optimplotfval });
<span class="comment">% Define constraint and objective function (to handle parameters)</span>
myConst = @(x) myNonlCon(x,s);
myObjFcn = @(x) myCostFcn(x,s);
<span class="comment">% Run Optimization</span>
tic
[x,cost] = fmincon(myObjFcn,x0,A,b,Aeq,beq,lb,ub,myConst,options);
toc
<span class="comment">% Display results thus far</span>
rrplot([cost0 cost],[x0(f); x(f)]',[x0(r); x(r)]')
<span class="comment">% close optimtool GUI</span>
optimGUI = com.mathworks.toolbox.optim.OptimGUI.getOptimGUI;
optimGUI.close
</pre><pre class="codeoutput">
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.
</pre><img vspace="5" hspace="5" src="rrdesign_03.png"> <img vspace="5" hspace="5" src="rrdesign_04.png"> <p>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.
</p>
<h2>Ensuring Suspension System Reliability<a name="11"></a></h2>
<p>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.
</p>
<p>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
</p>
<p><img vspace="5" hspace="5" src="rrdesign_eq14838.png"> </p>
<p>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.
</p>
<p>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.
</p><pre class="codeinput">load <span class="string">suspnSurvivalData</span>
dfittool(miles3,[],[],<span class="string">'Damping Ratio 0.3'</span>)
dfittool(miles5,[],[],<span class="string">'Damping Ratio 0.5'</span>)
dfittool(miles7,[],[],<span class="string">'Damping Ration 0.7'</span>)
disp(<span class="string">'Change the Display Type to Survivor Function in dropdown menu'</span>)
disp(<span class="string">'Then fit Weibull distributions to the data using New Fit Button'</span>)
<span class="comment">% Here is the reproduced final figure from dfittool (generated from</span>
<span class="comment">% dfittool session). Note you can load the session I used into the</span>
<span class="comment">% distribution fitting tool under File --> Load Session --></span>
<span class="comment">% suspnSurvivalData.dfit</span>
survivalFit(miles3,miles5,miles7)
xlabel(<span class="string">'Miles Driven'</span>)
</pre><pre class="codeoutput">Change the Display Type to Survivor Function in dropdown menu
Then fit Weibull distributions to the data using New Fit Button
</pre><img vspace="5" hspace="5" src="rrdesign_05.png"> <img vspace="5" hspace="5" src="rrdesign_06.png"> <p>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.
</p><pre class="codeinput"><span class="comment">% Close dfittool</span>
dffig = dfgetset(<span class="string">'dffig'</span>);
close(dffig)
type <span class="string">myNonlConR</span>
</pre><pre class="codeoutput">
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 = [];
</pre><p>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.
</p><pre class="codeinput"><span class="comment">% We'll use the command line equiavlaent to run the solver here to show</span>
<span class="comment">% results.</span>
<span class="comment">% Define constraint and objective function (to handle parameters)</span>
myConst = @(x) myNonlConR(x,s);
myObjFcn = @(x) myCostFcn(x,s);
<span class="comment">% Run Optimization</span>
tic
[xr,costr] = fmincon(myObjFcn,x,A,b,Aeq,beq,lb,ub,myConst,options);
toc
<span class="comment">% Plot results thus far</span>
rrplot([cost0 cost costr],[x0(f); x(f); xr(f)]',[x0(r); x(r); xr(r)]')
s0 = s; <span class="comment">% save inital parameter set for later reference</span>
</pre><pre class="codeoutput">
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.
</pre><img vspace="5" hspace="5" src="rrdesign_07.png"> <img vspace="5" hspace="5" src="rrdesign_08.png"> <h2>Optimizing for Robustness<a name="14"></a></h2>
<p>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.
</p>
<p>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.
</p>
<p>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.
</p>
<p><i><b>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.</b></i></p><pre class="codeinput">nRuns = 10000;
<span class="comment">% Probability distributions (assumed)</span>
front = 40 + raylrnd(40,nRuns,1); <span class="comment">% additional mass for front passengers (kg)</span>
back = raylrnd(40,nRuns,1); <span class="comment">% additional mass for rear passengers/cargo (kg)</span>
trunk = raylrnd(10,nRuns,1); <span class="comment">% additional mass for luggage (kg)</span>
</pre><p>The total mass, center of gravity, and moment of inertia are adjusted to account for the changes in mass distribution of the
automobile.
</p><pre class="codeinput">mcMb = 1200 + front + back + trunk;
<span class="comment">% Computed as a perturbation from the existing center of mass</span>
cm = (front.*s.rf - back.*s.rr - trunk.*s.rt)./mcMb;
<span class="comment">% a positive cm indicates the CM moved forward</span>
<span class="comment">% a negative cm indicates the CM moved backwards</span>
<span class="comment">%Adjust the moment of inertia about the old center of mass</span>
mcIyy = s.Iyy + front.*s.rf.^2 + back.*s.rr.^2 + trunk.*s.rt.^2;
<span class="comment">% Use parallel axis theorem to move moment of inertia to the new CM</span>
mcIyy = mcIyy - mcMb.*cm.^2;
<span class="comment">% Adjust the distance measurements to the new CG</span>
mcLf = s.Lf - cm;
mcLr = s.Lr + cm;
mcrf = s.rf - cm;
mcrr = s.rr + cm;
mcrt = s.rt + cm;
<span class="comment">% Put mc parameters in structure to pass to optimization solver</span>
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);
</pre><p>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.
</p><pre class="codeinput"><span class="comment">% Set solver options</span>
options = optimset;
options = optimset(options,<span class="string">'Display'</span> ,<span class="string">'Iter'</span>);
options = optimset(options,<span class="string">'TolFun'</span> ,1e-08);
options = optimset(options,<span class="string">'LargeScale'</span> ,<span class="string">'off'</span>);
options = optimset(options,<span class="string">'PlotFcns'</span> ,{ @optimplotx @optimplotfval });
<span class="comment">% Define constraint and objective function (to handle parameters)</span>
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -