%% 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
wbOptimSetUp % Load in geometry
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.
if license('test','virtual_reality_toolbox')
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])
image(imread('schematic.png','BackgroundColor',[1 1 1]))
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*.
load optimtoolProblem
% 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.
load wbPSOpt
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)
% -----------------------
% Start with default options
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
% suspension system. The first approach illustrated traditional gradient
% base techniques that work well for well-defined problems. We then found
% a condition where the gradient based method did not find a solution,
% since we do not have a completely defined problem and were unlucky in the
% choice of start point. The second approach, pattern search, was used to
% illustrate non-gradient based methods that work well for problems that
% are ill-defined or do not have well define gradients or objective
% functions.