Get parallel simulations in objective function of optimization Tool (GA,PSO)

Hey Guys

i have an optimization problem using the constrainted PSO:

https://de.mathworks.com/matlabcentral/fileexchange/25986-constrained-particle-swarm-optimization

The syntax and general working of the constrainted PSO is the same as GA.

In my objective function, i have to generate a simOut with data from a simulink simulation. One calculation needs about 7 hours of time, since the simulation needs about 3 seconds.To increase the algorithm speed, i want the simulations to run parallel. But it doesnt seem to work if i use "parpool('local')". So far ich changed the settings to "UseVectorized" so have the whole population in one objective function. My code looks like this:

parpool('local')
    options = psooptimset('Display','iter','OutputFcns',@myoutputFunc,'UseParallel',true,...
                          'Generations',200,'PopulationSize',40,'Vectorized','on');
    
    [xopt,fval,exitflag,output,population,scores] = pso(@objective,9,[],[],[],[],lb,ub,@constraint,options);

    function obj = objective(x)     
    % WANT TO: check the number of workers , in my case 6 
            % WANT TO: set parallel simulations equal to work number
           for p=1:length(x(:,1)) %WANT TO: set parallel working parfor loop to execute the parallel simulations
     
           SimuGeometrics = calculateGeometrics(x(p,:));            %SimuGeometrics is a structure of lengths and angle of different body that are used for simulink
           %start simulation
           simOut = sim('Mechanism.slx','CaptureErrors','on','SrcWorkspace','current');
            %Result
            pos1 = [simOut.xx.data, simOut.yy.data, simOut.zz.data];
            
            obj(p,1) = calculateDeviation(pos1);
            
            disp(['x: ',num2str(x(p,:)) ,' obj: ',num2str(obj(p,1))]);
           end
    end
 

I know that there are some functions like parfor or parsim but i dont know how to include them into my Code.

Of cause, my WANTO TO's are only ideas how i think it can work, if for example setting a parfor loop is nonsense, you can suggest another solution

Edit:

@Abdolkarim Mohammadi-Balani

i saw your Answer from :https://de.mathworks.com/matlabcentral/answers/595402-how-to-use-the-parallel-computing-toolbox-of-matlab-to-optimize-the-parameters-of-a-simulink-model-w#answer_496378

it seems like you described a possible solution for my problem. I will try out what you described today and will share my solution if i succeed, but it would help me alot if u could share your expertise on this.

I hope somebody can help me with this problem.

Answer 1

[Abdolkarim Mohammadi]

Hi Marcus.

I haven't used that File Exchange submission for PSO, but the general workflow in such situations is to:

1. Pass the whole population to the objective function. This is done by setting the UseVectorized option to true.
2. When the whole population is received by the objective function, it must create a SimulationInput object. You can learn about creating it here. Your SimulationInput object is an array, where each element represents a chromosome in the population.
3. Then you pass the single SimulationInput object to parsim().
4. Next, simulations are run and the SimulationOutput object is returned.
5. Finally, you calculate the objective function values from the SimulationOutput object and return them to the solver.

The objective funciton's code might look something like this:

function f = ObjFun(x)

% get population size 
PopulationSize = size(x, 1);

% create the SimluationInput object 
for i1 = 1:PopulationSize 
    SimIn(i1) = Simulink.SimulationInput('MyModel');
    SimIn(i1) = SimIn(i1).setVariable('MyVar1', x(i1,1))
    SimIn(i1) = SimIn(i1).setVariable('MyVar2', x(i1,2))
    % do it for all of the variables
end

% run the simulations and get the SimulationOutput object 
SimOut = parsim(SimIn); 

% initialize the arrays containing the objective function values 
f = zeros (PopulationSize, 1);

% fill f with objective function values 
for i1 = 1:PopulationSize 
    f(i1) = SimOut(i1).MyOutput; 
end

end

Read about setVariable function here.

[Marcus Heßeling]

Hey @Abdolkarim Mohammadi-Balani,

thank you very much for your help, it worked perfectly. In your other post you were talking about speeding up the simulation by using RapidAccelerator or FastRestart. RA doesnt work, since it doesnt provide any output data, but the Accelerator does work. The FastRestart doesnt work, and im curious:

As far as i understand it, in the objective i build 40 SimIn Models for every particle. Every SimIn Model is only performed once, so the FastRestart wouldnt be of any use in this case right ? Since i only change the variables in my Models, a FastRestart could be possible? is there a way to implement it into my Code, for example buy only make SimIn Models equal to my number of workers and then reuse the SimIn until i reach the 40 particles ? The problem i see is in the unequal numbers of iterations (40 / 6 workers)

For everyone else who want to know my made code, its provided in the following. Its very similar to Abdolkarims Version ( Particle in PSO is the same as Chromosomes in GA)

function obj = objective(x)
 
PopulationSize = size(x, 1);
startAngle=1;
 
    % create the SimluationInput object for every particle
    for i1 = 1:PopulationSize 
        SimuGeometrics = calculateGeometrics(x(i1,:)); 
        % Set Model
        SimIn(i1) = Simulink.SimulationInput('Mechanism');
 
        %Set Model Parameters
 
        %SimIn(i1) = SimIn(i1).setModelParameter('SimulationMode','accelerator');
        %SimIn(i1) = SimIn(i1).setModelParameter('UseFastRestart', 'on');
 
        %Set Variables
        SimIn(i1) = SimIn(i1).setVariable('SimuGeometrics', SimuGeometrics);
        SimIn(i1) = SimIn(i1).setVariable('startAngle', startAngle);
    end
    
    %start simulation
     SimOut = parsim(SimIn);
     
     
    %load wanted Model 
    [perfectModel] = perfectModelLoading;
    
    
    %hand over the objectives for every particle and display results     
    obj =zeros(PopulationSize,1);  
    
    for i2 = 1:PopulationSize
            
        [disVec]= calculateDeviation(perfectModel,SimOut(i2));    
 
        obj(i2,1) = mean(disVec) + max(disVec);
    
    disp(['Particle: ',num2str(i2),' x: ',num2str(x(i2,:)) ,' obj: ',num2str(obj(i2,1))]);
    end
 
end

[Abdolkarim Mohammadi]

1. When you use parsim(), you do not need to do anything further. Simulink performs one experiment for each element of the SimIn array, and outputs the results in SimOut, which has the same size as SimIn.

2. Fast restart is very useful here. Again, you do not need to handle it separately, it is an option in parsim(), as mentioned in parsim() documentation. You just need to modify parsim() like this:

SimOut = parsim (SimIn, 'UseFastRestart', 'on');

You can read about fast restart here. Simulink passes different stages, in order to execute a simulation. It first compiles the model, then executes it, and finally terminates. Model compilation can take a long time in complex model. When you enable fast restart, Simulink compiles the model once and executes it multiple times. This can save a lot of time for parsim(). Please read the documentation thoroughly.

[Marcus Heßeling]

Thank you @Abdolkarim Mohammadi-Balani

When i turn on FastRestart, my model runs way faster. Unfortunatly , in FR mode the output value remains the same and i cant find the reason for the non- change of the objective. I already checked the link you provided about FastRestart and also read the limitations here and the part of "Why Did the Simulation Output Stay the Same?" in the Tuning documentation. As Far as i know i dont make any structural changes as i only change the length parameters of my solids in each simulation. To make it more applicable for you, i attached the matlab script and the simulink model in which i test the FastRestart mode. (Run Speedtests script)

Hopefully theres an easy way to see why my model has no compatibility with FastRestart. Thanks for your useful tips man, i appreciate it a lot !

[Abdolkarim Mohammadi]

I ran your code, but it threw an error. I really don't know why fast restart makes the outputs not to change for different particles. I use this approach a lot, and it works as expected. Maybe you should not use accelerator with fast restart.

[Abdolkarim Mohammadi]

By the way, you can use the metaheuristic algorithm I recently proposed, called Golden Eagle Optimizer (GEO), instead of PSO. The source code has similar syntax compared to MATLAB's own ga(). In addition, a graphical user interface is also developed.

https://www.mathworks.com/matlabcentral/fileexchange/84430-golden-eagle-optimizer-toolbox

[Marcus Heßeling]

Hey @Abdolkarim Mohammadi-Balani

Wow thanks man, this algorithm is very interesting, i will definitely give this algorithm a try.

I see that the Toolbox was developed in R2020a and my Model was made wir R2020b, so i attached the file again for R2020a. Just use this in the command window if the Library doesnt work:

set_param(gcs,'EnableLBRepository','on')

If it still threws an error or you can't find a reason for failed Restart in the Model, then we can quit the solving of the FRestart model(PS. My model runs in Normal mode since one particle only has 2 sec Simulation time. Accelerator even increase the Sim() time. The compiling takes the majority of the whole simulation prozess , thats why i wanted to Use FastRestart)

Nevertheless i want to thank you again for you help, your optimization tips saved me a great deal of time and trys. Hopefully your GEO gets a good reputation.

[Abdolkarim Mohammadi]

I reviewed your code and I found out that you didn't use the proper methodology for multiple simulation runs. You wrote:

for i0 = ...
    simOut = sim('Opt_9Facets_Model_Parallel.slx','CaptureErrors','on','SrcWorkspace','current',...
                        'FastRestart','on');
end

Which is the worst way to run multiple simulations. I told you that you should use SimulationInput objects, not calling sim() using model name. Using sim() with model name as the first input and fast restart enabled, prevents Simulink engine from updating simulation parameters in future runs. That's why you are getting the same objective function value in each run. The first run with fast restart LOCKS the model, and future runs do not read parameters from workspace (either base workspace or model workspace).

[Irfan Khan]

@Abdolkarim Mohammadi I used the objective function code and simulation is very fast now. But since i am new matlab user i made some mistakes. I fitness value graphs is abnormal. and the values in varriable solution2 are not giving stabilized output.

The background is i am tuning a PID controller in simulink using the genetic algorithim optimization tool. My cost function is "to workspace" (ITAE). In serial simulation the old objective function was working fine but in parallel simulation objective function, the output is not correct becuase i am not sure what to right in this line

cost(i2,1) = SimOut.ITAE(end);

Parallel Simulation objective function which is running but not tunning

function cost= tunning(k)

% get population size

PopulationSize = size(k, 1);

% create the SimluationInput object

for i1 = 1:PopulationSize

SimIn(i1) = Simulink.SimulationInput('test');

SimIn(i1) = SimIn(i1).setVariable('k(1)', k(i1,1));

SimIn(i1) = SimIn(i1).setVariable('k(2)', k(i1,2));

SimIn(i1) = SimIn(i1).setVariable('k(3)', k(i1,3));

% do it for all of the variables

end

% run the simulations and get the SimulationOutput object

SimOut = sim("test.slx",'CaptureErrors','on','SrcWorkspace','current');

%cost= ITAE;

% initialize the arrays containing the objective function values

cost= zeros (PopulationSize, 1);

% fill cost with objective function values

for i2 = 1:PopulationSize

cost(i2,1) = SimOut.ITAE(end);

end

end

Serial Objective Function which is slow but give right result

function cost = tunning_serial(k)

assignin('base','k',k);

sim('test.slx');

cost= ITAE(length(ITAE));

end

My Simulink Model using PID controller

Live Editor Optimization using Genetic Algorithim

% Set nondefault solver options

options2 = optimoptions("ga","UseVectorized",true,"CreationFcn",...

"gacreationuniform","SelectionFcn","selectiontournament","MutationFcn",...

"mutationadaptfeasible","Display","iter","PlotFcn","gaplotbestf");

% Solve

[solution2,objectiveValue] = ga(@tunning,nVariables,[],[],[],[],...

zeros(nVariables,1),repmat(5,nVariables,1),[],[],options2);

Abnormal Fitnessvalue graphs

I guess this line Simout.ITAE(end) is the problem.

I have attached my objective function files