Vectorized Surrogate Optimization for Custom Parallel Simulation

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])

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

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)

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

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