function [x,nf,time] = rrDesignParallel(nRuns)
% Copyright 2010 The MathWorks, Inc.
%% 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;
%%
% 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' ,'None');
options = optimset(options,'TolFun' ,1e-08);
options = optimset(options,'Algorithm' ,'active-set');
% Define constraint and objective function (to handle parameters)
myConst = @(x) myNonlConRR(x,s);
myObjFcn = @(x) myCostFcnRR(x,s);
%% Run optimization
x0 = [s.kf0 s.cf0 s.kr0 s.cr0];
% 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];
tic
[x,costRR,exitflag,output] = fmincon(myObjFcn,x0,A,b,Aeq,beq,lb,ub,myConst,options);
time = toc;
nf = output.funcCount;