Code covered by the BSD License

# Optimization of a Double Wishbone Suspension System

31 Oct 2006 (Updated )

Demo files from the Webinar "Introduction to Optimization with MATLAB(R) Products" Oct. 26, 2006

walkthrough.m
%% Optimization of a Double Wishbone Suspension System
% This demo shows how to use MATLAB, Optimization Toolbox, and Genetic
% Algorithm and Direct Search Toolbox to optimize the design of a double
% wishbone suspension system.
%
% _Note: You will need to have the following products installed in order to
% run this demo: MATLAB, Simulink, Optimization Toolbox, Genetic Algorithm
% and Direct Search Toolbox, and SimMechanics. Optional: Virtual Reality
% Toolbox._
%
%% Introduction to the Problem
% We wish to optimize the response of a double wishbone suspension system.
% The response we are evaluating is the camber angle vs. travel distance as
% shown in the figure below.  Our current design is shown in green.  We
% want to achieve the profile shown in blue.
%
% Camber angle is the angle of the tire relative to the perpendicular axis
% to the road.  A negative camber angle is beneficial to handling.
%
% Travel distance is the amount of verticle motion of the car as a result
% of traveling over a bump or pothole in the road.
%
verifyInstalled                         % verify products installed
wbPlotFun(idealProfile)                 % Plot ideal profile
wbPlotFun(x,time,'initial')             % Simulate and Plot current profile

%% Double Wishbone Model
% A model of the double wishbone suspension system was created using
% Simulink with SimMechanics.  We will use this model to simulate
% the performance of our suspension system.
simName = 'DoubleWishboneVR.mdl';
else
simName = 'DoubleWishbone.mdl';
end
open(simName)
sim(simName)
%% Input Geometry Definition
% The wishbone geometry is defined according to diagram below.  We have 15
% inputs that we specify, denoted as X1 - X15, that determine the double
% wishbone connection points and rod lengths.  These are defined as:
%
%   Upper Arm Length                               X1  =  AG
%   Upper Arm Connection Point B         (X2, X3, X4)  =  (xB, yB, zB)
%   Upper Arm Connection Point C         (X5, X6, X7)  =  (xC, yC, zC)
%   Lower Arm Length                               X8  =  DG
%   Upper Arm Connection Point E      (X9, X10, X11)   =  (xE, yE, zE)
%   Upper Arm Connection Point F      (X12, X13, X14)  =  (xF, yF, zF)
%   Connecting Ling Length                        X15  =  AD
%
% All dimensions are input to the model in inches.

x   %Current geometric definition
figure('Position',[1 700 700*317/545 700])
axis equal, axis tight, axis off

%% Problem Formulation
% Our objective is to minimize the difference between the 'Actual' profile
% and the 'Ideal' profile by changing the geometry parameters X1 - X15.
% Our problem is formulated as a constrained minimization problem:
%
% Objective function:
%
% $$min f(x)$$
%
% where
%
% $$f(x)= norm[(CamberAngle_{actual} - CamberAngle_{ideal}) +$$
% $$(TravelDistance_{actual}-TravelDistance_{ideal})]$$
%
% Subject to (constraints):
%
% $$A \cdot x <= 0$$
%
% $$lb <= x <= ub$$

%% Objective Function
% Our objective function is defined to return a single value, which is a
% measure of how close the 'Actual' profile is to the 'Ideal' profile.  Our
% objective function is defined in the function |wbObjFun|.
type wbObjFun

%% Constraint Definition
% The constraint specified for the model are (refer to previous figure):
%
%  BAC Angle (degrees)               25 <=   BAC   <=  35  -|
%  EDF Angle (degrees)               15 <=   EDF   <=  30   |- A
%  Point B/C rotation (degrees)              BC    <=  10   |
%  Point E/F rotation (degrees)              EF    <=   5  -|
%  Upper Arm Length Limits           6  <=   X1    <=  16
%  Point B X-Axis Limits            10  <=   X2    <=  16
%  Point C X-Axis Limits            10  <=   X5    <=  16
%  Lower Arm Length Limits           8  <=   X8    <=  18
%  Point E X-Axis Limits             6  <=   X9    <=  14
%  Point F X-Axis Limits            12  <=   X12   <=  20
%  Ling Length                               X15   <=  18
%                                  |__________|__________|
%                                        |          |
%                                       lb         ub
%
% Coefficient matrix A and upper and lower bounds are defined in
% |wbOptimSetup|.

%% Solve the Problem Using Optimization Toolbox
% The problem can be solved using the |fmincon| solver in Optimization
% Toolbox.  To use the |optimtool| GUI to set up an run the problem, load
% the optimization problem definition in |optimtoolProblem.mat|, then type
% |optimtool| at the command line and once the GUI is open, select  _File
% --> Import Problem --> optimtoolProblem._ The problem is now defined in
% the GUI, select *Start*.

% optimtool % uncomment this line to run interactively in the GUI
%
% Command line equivalent
% -----------------------
% Start with the default options (may want to comment out this section if
% running interactively.
options = optimset;
% Modify options setting
% Note: You can use the defaults as well, this will speed up the solution
options = optimset(options,'Display' ,'iter');
options = optimset(options,'TolFun' ,0.1);
options = optimset(options,'LargeScale' ,'off');
% Request plots, add my custom plot to the mix
options = optimset(options,'OutputFcn' ,{ @updatePlot });
options = optimset(options,'PlotFcns' ,{  @optimplotx @optimplotfval });
tic
[x_x0] = fmincon(@(x) wbObjFun(x,time,idealProfile),x0,Aineq,bineq,...
[],[],lb,ub,[],options);
toc
x_x0

%% Solve the Problem Using a Different Start Point
% At this point, Id like take a look at our problem from a different
% angle.  As you may know, gradient based solvers, such as the one
% we used here, tend to fail if the objective function does not have smooth
% derivatives or the problem is ill-defined.  Let's change the start point
% of this problem, and see if |fmincon| can find a solution.
tic
try
[x_x0c] = fmincon(@(x) wbObjFun(x,time,idealProfile),x0c,Aineq,bineq,...
[],[],lb,ub,[],options);
catch
disp('We caught the following error')
nowork = lasterror;
nowork.message
end
toc

%%
% The error message is from our Simulink model and basically tells us we
% put in a bad definition for geometry that our model could not handle.
%
% Just to be clear, the problem is not with the |fmincon| solver, but with
% my knowledge of the design space in which we are looking to explore.
% |fmincon| did what I asked it to do, it just found a point that my model
% could not handle.  To use |fmincon| we need to add constraints to our
% problem definition or redefine our objective function to prevent this
% situation from happening.
%
% Im pointing this out because in practice you may often encounter problems
% where you just dont have enough information to completely define your
% problem or you may have an objective function that does not have well
% defined gradients, or even holes in the design space. This often renders
% traditional gradient based solvers, such as the one used here, ineffective.
% This is where algorithms in Genetic Algorithm and Direct Search Toolbox
% are useful.
%% Solve the Problem Using Genetic Algorithm and Direct Search Toolbox
% We can use pattern search to solve this problem without having to redefine
% our objective function, constraints, or the Simulink model.  Let's define
% our problem as before and solve.  We loaded in the problem formulation in
% the structure |wbPSOpt| in the previous step.  Let's see if pattern
% search can find a solution.
%
% _Note: You can set up your problem in the |psearchtool| GUI analogously
% to the |optimtool| GUI.  File --> Import Options --> wbPSOpt, set up the
% problem definition, and then click on *Start*._

% psearchtool % uncomment this line to run interactively
%
% Command line equivalent (comment out this section if using psearchtool)
% -----------------------
options = psoptimset;
% Modify some parameters
% Note: You can use the defaults as well, this will speed up the solution
options = psoptimset(options,'TolFun' ,0.06);
options = psoptimset(options,'MeshContraction' ,0.25);
options = psoptimset(options,'MaxMeshSize' ,1 );
options = psoptimset(options,'PollingOrder' ,'success');
options = psoptimset(options,'SearchMethod' ,@GPSPositiveBasisNp1);
options = psoptimset(options,'Display' ,'iter');
options = psoptimset(options,'Cache' ,'on');
options = psoptimset(options,'CacheTol' ,2.2204e-010);
% Request plots, add my custom plot to the mix
options = psoptimset(options,'OutputFcns' ,{ { @updatePlot } });
options = psoptimset(options,'PlotFcns' ,{  @psplotbestf @psplotbestx });
%%Run PATTERNSEARCH
tic
[x_x0c_ps] = patternsearch(@(x) wbObjFun(x,time,idealProfile),x0c,...
Aineq,bineq,[],[],lb,ub,[],options);
toc
x_x0c_ps
%%
% Pattern search found a solution.  We set the function tolerance for
% both |fmincon| and |patternsearch| to limit the total number of
% iterations required to satisfy our convergence criteria.  Decreasing this
% value would allow the solution to converge more tightly on the 'Ideal'
% profile, but would require longer solution times.

%% Summary
% We used two approaches to optimize the design of a double wishbone
% functions.