# 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)`
```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```
```Starting parallel pool (parpool) using the 'local' profile ... Connected to the parallel pool (number of workers: 6). ```
`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(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)``` ```surrogateopt stopped because it exceeded the function evaluation limit set by 'options.MaxFunctionEvaluations'. ```
```fittedparams = 1×6 10.5627 19.8962 9.8420 8.9616 2.5723 27.9687 ```
```fval = 23.6361 ```
`paralleltime = toc`
```paralleltime = 457.9271 ```

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

`disp([ssq,fval])`
``` 26.9975 23.6361 ```

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