Code covered by the BSD License  

Highlights from
Speeding Up Optimization Problems with Parallel Computing

image thumbnail

Speeding Up Optimization Problems with Parallel Computing

by

Stuart Kozola (view profile)

 

Files from the webinar: Speeding up optimization problems with parallel computing

Editor's Notes:

This file was selected as MATLAB Central Pick of the Week

optimMotorDesign.m
%% DC Motor Design Optimization Demo
% This script shows how to perform an optimization problem using two levels
% of parallelism.  Level 1: Mulitobjective Pareto frount mapping, Level 2:
% parallel finite differences in optimization solvers.
clear all; close all; clc
% Copyright 2010 The MathWorks, Inc.
%% Motor Design Variables
D   = 6.13;                 % rotor diameter (cm)
L   = 8.648;                % rotor axial length (cm)
Lwa = 14.12;                % armature wire length (m)
Lwf = 309.45;               % field wire length (m)
ds  = 0.612;                % slot depth (cm)

d = [D L Lwa Lwf ds];
dscale = [6 8 14 309 0.6];

%% Coupling Variables b, shared with control problem
n   = 122;                   % rotational speed (rev/s)
V   = 40;                   % design voltage (V)
Pmin= 3.94;                 % minimum required power (kW)
Tmin= 5.12e-3;         	% minimum required torque (kNm)

b = [n V Pmin Tmin];
bscale = [122 40 3.94 5.12e-3];

%% Parameter Vector a (constants)
dcu = 8.94e3;               % copper density at 25C (kg/m^3)
dfe = 7.98e3;               % iron density at 25C (kg/m^3)
fi  = 0.7;                  % pole arc to pole pitch ratio
p   = 2;                    % number of poles
S   = 27;                   % number of slots (teeth on rotor)
rho = 1.8e-8;               % resistivity (ohm-m) of copper at 25C

a = [dcu,dfe,fi,p,S,rho];

%% Solve for optimal motor design
% Solve for d and update b and v parameters to start with feasible design
op = optimset('Disp','iter','Algorithm','SQP');
f  = @(x) motorDesign(x,a,b);
dopt = fmincon(f,d./dscale,[],[],[],[],zeros(size(d)),[],f,op);

%% Initial Control Design
pid = [0 0.2 0];                % [Kp Ki Kd] PID Gains
u   = [0.02 20];                % [Kf NL]    friction coefficient, speed
[w,cd,ceqd,v] = motorDesign(dopt,a,b);
motorController(pid,u,v)

%% Optimization Problem Definition
alpha = 0.5;                    % weighting for weight/performance index
z0 = [dopt pid b./bscale];      % intial starting point
LB = zeros(size(z0));           % lower bound (all positive)
LB(11) = eps;

% Objective and Constraint function
sys = @(z) motorSystem(z,a,u,alpha);

%% Solve using |fmincon|
% Baseline
op = optimset('useParallel','never','Disp','iter','Algorithm','SQP',...
              'TolCon',1e-4);
tic
znew = fmincon(sys,z0,[],[],[],[],LB,[],sys,op) ;  
toc

[wJ,c,ceq] = motorSystem(znew,a,u,alpha);
close all
plot(wJ(:,1)*10,wJ(:,2),'o')
xlabel('DC Motor Weight (kg)')
ylabel('Control Performance Index')

%% Solve using |fmincon| using Parallel Gradient Estimation
% Time using 4 workers
op = optimset('useParallel','always','Disp','iter','Algorithm','SQP',...
              'TolCon',1e-4);
files = {'motorController.m','motorDesign.m','motorSystem.m'};          
matlabpool('open','speedy-4core',4,'FileDependencies',files)
tic
znew = fmincon(sys,z0,[],[],[],[],LB,[],sys,op);   
toc
matlabpool close

%% Map Pareto Front by changing alpha on single machine
alpha = [0 0.1 0.3 0.5 0.8 0.9 0.92 0.95 0.97 1];

op = optimset('useParallel','never','Disp','iter','Algorithm','SQP',...
              'TolCon',1e-4);
%%
tic
for j = 1:length(alpha)
    sys = @(z) motorSystem(z,a,u,alpha(j));
    znew = fmincon(sys,z0,[],[],[],[],LB,[],sys,op);
    [wJ(j,:),~,~,~] = motorSystem(znew,a,u,alpha(j));
end
toc
plot(wJ(:,1)*10,wJ(:,2),'-o')
xlabel('DC Motor Weight (kg)')
ylabel('Control Performance Index')
%% Map Pareto Front by changing alpha on 10 machines
% Run each problem in parallel on a single core

op = optimset('useParallel','never','Disp','none','Algorithm','SQP',...
              'TolCon',1e-4);

matlabpool('open','speedy-4core',4,'FileDependencies',files)
tic
parfor j = 1:length(alpha)
    sys = @(z) motorSystem(z,a,u,alpha(j));
    znew = fmincon(sys,z0,[],[],[],[],LB,[],sys,op);
    [wJ(j,:),~,~,~] = motorSystem(znew,a,u,alpha(j));
end
toc
matlabpool close
plot(wJ(:,1)*10,wJ(:,2),'-o')
xlabel('DC Motor Weight (kg)')
ylabel('Control Performance Index')
%% Solve using Jobs/Tasks and Parallel Gradient Estimation
nTasks = length(alpha);
nLabs = 4;

% Find job manager
jm = findResource('scheduler','type','jobmanager', ...
                  'name','speedyJM','LookupURL','speedy-00-ah');
for j = 1:nTasks
    % Create matlabpool job
    job(j) = createMatlabPoolJob(jm); %#ok<*SAGROW>
    job(j).MaximumNumberOfWorkers = nLabs;
    job(j).FileDependencies = {'motorSystem.m','motorDesign.m',...
                               'motorController.m','motorOptimizationTask.m'};
    % Add tasks to job
    createTask(job(j),@motorOptimizationTask,1,{z0,a,u,alpha(j),LB});
end

%% Submit Job
tic
for j = 1:nTasks
    submit(job(j));
end
% and wait
for j = 1:nTasks
    waitForState(job(j),'finished');
    results{j} = getAllOutputArguments(job(j));
end
toc

%% Plot Result
for j = 1:nTasks
    [wJ(j,:),~,~,~] = motorSystem(results{j}{1},a,u,alpha(j));
end
plot(wJ(:,1)*10,wJ(:,2),'-o')
xlabel('DC Motor Weight (kg)')
ylabel('Control Performance Index')

%% Clean Up Job
destroy(job);

Contact us