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

 

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