Main Content

Vectorized Surrogate Optimization for Custom Parallel Simulation

This example shows how to use the surrogateopt UseVectorized option to perform custom parallel optimization. You can use this technique when you cannot use the UseParallel option successfully. For example, the UseParallel option might not apply to a Simulink® simulation that requires parsim for parallel evaluation. Optimizing a vectorized parallel simulation involves considerable overhead, so this technique is most useful for time-consuming simulations.

The parallel strategy in this example is to break up the optimization into chunks of size N, where N is the number of parallel workers. The example prepares N sets of parameters in a Simulink.SimulationInput vector, and then calls parsim on the vector. When all N simulations are complete, surrogateopt updates the surrogate and evaluates another N sets of parameters.

Model System

This example attempts to fit the Lorenz dynamical system to uniform circular motion over a short time interval. The Lorenz system and its uniform circular approximation are described in the example Fit an Ordinary Differential Equation (ODE).

The Lorenz_system.slx Simulink model implements the Lorenz ODE system. This model is included when you run this example using the live script.

The fitlorenzfn helper function at the end of this example calculates points from uniform circular motion. Set circular motion parameters from the example Fit an Ordinary Differential Equation (ODE) that match the Lorenz dynamics reasonably well.

x = zeros(8,1);
x(1) = 1.2814;
x(2) = -1.4930;
x(3) = 24.9763;
x(4) = 14.1870;
x(5) = 0.0545;
x(6:8) = [13.8061;1.5475;25.3616];

This system does not take much time to simulate, so the parallel time for the optimization is not less than the time to optimize serially. The purpose of this example is to show how to create a vectorized parallel simulation, not to provide a specific example that runs better in parallel.

Objective Function

The objective function is to minimize the sum of squares of the difference between the Lorenz system and the uniform circular motion over a set of times from 0 through 1/10. For times xdata, the objective function is

objective = sum((fitlorenzfn(x,xdata) - lorenz(xdata)).^2) - (F(1) + F(end))/2

Here, lorenz(xdata) represents the 3-D evolution of the Lorenz system at times xdata, and F represents the vector of squared distances between corresponding points in the circular and Lorenz systems. The objective subtracts half of the values at the endpoints to best approximate an integral.

Consider the uniform circular motion as the curve to match, and modify the Lorenz parameters in the simulation to minimize the objective function.

Calculate Lorenz System for Specific Parameters

Calculate and plot the Lorenz system for Lorenz's original parameters.

model = 'Lorenz_system';
open_system(model);
in = Simulink.SimulationInput(model);
% params [X0,Y0,Z0,Sigma,Beta,Rho]
params = [10,20,10,10,   8/3, 28]; % The original parameters Sigma, Beta, Rho
in = setparams(in,model,params);
out = sim(in);

yout = out.yout;
h = figure;
plot3(yout{1}.Values.Data,yout{2}.Values.Data,yout{3}.Values.Data,'bx');
view([-30 -70])

Figure contains an axes object. The axes contains a line object which displays its values using only markers.

Calculate Uniform Circular Motion

Calculate the uniform circular motion for the x parameters given previously over the time interval in the Lorenz calculation, and plot the result along with the Lorenz plot.

tlist = yout{1}.Values.Time;
M = fitlorenzfn(x,tlist);
hold on
plot3(M(:,1),M(:,2),M(:,3),'kx')
hold off

Figure contains an axes object. The axes object contains 2 objects of type line. One or more of the lines displays its values using only markers

The objfun helper function at the end of this example calculates the sum of squares difference between the Lorenz system and the uniform circular motion. The objective is to minimize this sum of squares.

ssq = objfun(in,params,M,model)
ssq = 26.9975

Fit Lorenz System in Parallel

To optimize the fit, use surrogateopt to modify the parameters of the Simulink model. The parobjfun helper function at the end of this example accepts a matrix of parameters, where each row of the matrix represents one set of parameters. The function calls the setparams helper function at the end of this example to set parameters for a Simulink.SimulationInput vector. The parobjfun function then calls parsim to evaluate the model on the parameters in parallel.

Open a parallel pool and specify N as the number of workers in the pool.

pool = gcp('nocreate'); % Check whether a parallel pool exists
if isempty(pool) % If not, create one
    pool = parpool;
end
N = pool.NumWorkers
N = 6

Set the BatchUpdateInterval option to N and set the UseVectorized option to true. These settings cause surrogateopt to pass N points at a time to the objective function. Set the initial point to the parameters used earlier, because they give a reasonably good fit to the uniform circular motion. Set the MaxFunctionEvaluations option to 600, which is an integer multiple of the 6 workers on the computer used for this example.

options = optimoptions('surrogateopt','BatchUpdateInterval',N,...
    'UseVectorized',true,'MaxFunctionEvaluations',600,...
    'InitialPoints',params);

Set bounds of 20% above and below the current parameters.

lb = 0.8*params;
ub = 1.2*params;

For added speed, set the simulation to use fast restart.

set_param(model,'FastRestart','on');

Create a vector of N simulation inputs for the objective function.

simIn(1:N) = Simulink.SimulationInput(model);

For reproducibility, set the random stream.

rng(1)

Optimize the objective in a vectorized parallel manner by calling parobjfun.

tic
[fittedparams,fval] = surrogateopt(@(params)parobjfun(simIn,params,M,model),lb,ub,options)

Figure Optimization Plot Function contains an axes object. The axes object with title Best Function Value: 23.073, xlabel Iteration, ylabel Function value contains an object of type scatter. This object represents Best function value.

surrogateopt stopped because it exceeded the function evaluation limit set by 
'options.MaxFunctionEvaluations'.
fittedparams = 1×6

   10.3789   19.9803   10.0073    9.4956    2.8269   27.9263

fval = 23.0730
paralleltime = toc
paralleltime = 397.1877

The objective function value improves (decreases). Display the original and improved values.

disp([ssq,fval])
   26.9975   23.0730

Plot the fitted points.

figure(h)
hold on
in = setparams(in,model,fittedparams);
out = sim(in);
yout = out.yout;
plot3(yout{1}.Values.Data,yout{2}.Values.Data,yout{3}.Values.Data,'rx');
legend('Unfitted Lorenz','Uniform Motion','Fitted Lorenz')
hold off

Figure contains an axes object. The axes object contains 3 objects of type line. One or more of the lines displays its values using only markers These objects represent Unfitted Lorenz, Uniform Motion, Fitted Lorenz.

To close the model, you must first disable fast restart.

set_param(model,'FastRestart','off');
close_system(model)

Conclusion

When you cannot use the UseParallel option successfully, you can optimize a simulation in parallel by setting the surrogateopt UseVectorized option to true and the BatchUpdateInterval option to a multiple of the number of parallel workers. This process speeds up the parallel optimization, but involves overhead, so is best suited for time-consuming simulations.

Helper Functions

The following code creates the fitlorenzfn helper function.

function f = fitlorenzfn(x,xdata)
theta = x(1:2);
R = x(3);
V = x(4);
t0 = x(5);
delta = x(6:8);
f = zeros(length(xdata),3);
f(:,3) = R*sin(theta(1))*sin(V*(xdata - t0)) + delta(3);
f(:,1) = R*cos(V*(xdata - t0))*cos(theta(2)) ...
    - R*sin(V*(xdata - t0))*cos(theta(1))*sin(theta(2)) + delta(1);
f(:,2) = R*sin(V*(xdata - t0))*cos(theta(1))*cos(theta(2)) ...
    - R*cos(V*(xdata - t0))*sin(theta(2)) + delta(2);
end

The following code creates the objfun helper function.

function f = objfun(in,params,M,model)
in = setparams(in,model,params);
out = sim(in);
yout = out.yout;
vals = [yout{1}.Values.Data,yout{2}.Values.Data,yout{3}.Values.Data];
f = sum((M - vals).^2,2);
f = sum(f) - (f(1) + f(end))/2;
end

The following code creates the parobjfun helper function.

function f = parobjfun(simIn,params,M,model)
N = size(params,1);
f = zeros(N,1);
for i = 1:N
    simIn(i) = setparams(simIn(i),model,params(i,:));
end
simOut = parsim(simIn,'ShowProgress','off'); % Suppress output
for i = 1:N
    yout = simOut(i).yout;
    vals = [yout{1}.Values.Data,yout{2}.Values.Data,yout{3}.Values.Data];
    g = sum((M - vals).^2,2);
    f(i) = sum(g) - (g(1) + g(end))/2;
end
end

The following code creates the setparams helper function.

function pp = setparams(in,model,params)
% parameters [X0,Y0,Z0,Sigma,Beta,Rho]
pp = in.setVariable('X0',params(1),'Workspace',model);
pp = pp.setVariable('Y0',params(2),'Workspace',model);
pp = pp.setVariable('Z0',params(3),'Workspace',model);
pp = pp.setVariable('Sigma',params(4),'Workspace',model);
pp = pp.setVariable('Beta',params(5),'Workspace',model);
pp = pp.setVariable('Rho',params(6),'Workspace',model);
end

See Also

| (Simulink)

Related Topics