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

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

%
%% 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

% 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(100)

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

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

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

disp([ssq,fval])

% 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

% Copyright 2020 The MathWorks, Inc.