MATLAB Examples

Electron Problem Optimization

This demo was adapted from a 2009 digest article: Improving Optimization Performance with Parallel Computing

This demo walks through setting up and running an electrostatics optimization problem using parallel computing to reduce solution time.


Problem Description

This example was inspired by Dolan, Moré, and Munson, “Benchmarking Optimization Software with COPS 3.0”, Argonne National Laboratory Technical Report. ANL/MCS-TM-273, February 2004.

This example illustrates how to formulate and solve an optimization problem in MATLAB. Consider N electrons in a conducting body shown below. The electrons arrange themselves to minimize their potential energy subject to the constraint of lying inside the conducting body. At the minimum total potential energy, all the electrons lie on the boundary of the body. Because the electrons are indistinguishable, there is no unique minimum for this problem (permuting the electrons in one solution gives another valid solution).

Elapsed time is 0.008745 seconds.

The optimization goal is to minimize the total potential energy of the electrons subject to the constraint that the electrons remain within the conducting body. The objective function, potential energy, is the sum of the inverses of the distances between each electron pair (i,j = 1, 2, 3,… E):

$$energy=\sum_{i<j}{1\over \vert x_i\vert - \vert x_j \vert}$$

The objective function is defined in sumInvDist

type sumInvDist
electronProblem.objective = @sumInvDist;
function energy = sumInvDist(x)

E = length(x);
energy = 0;
for i = 1:3:E-5
    for j = i+3:3:E-2
        dist   = x(i:i+2) - x(j:j+2);
        if norm(dist) == 0
            energy = 1e9;
            energy = energy + 1/sqrt(dist'*dist);

The constraints that define the boundary of the conducting body are

$$z\leq -\vert x\vert - \vert y \vert$$

$$x^2+y^2+(z+1)^2 \leq 1$$

As written, the first inequality is a nonsmooth nonlinear constraint because of the absolute values on x and y. Absolute values can be linearized to simplify the optimization problem. This constraint can be written in linear form as a set of four constraints for each electron, i.

$$x_{i,3} - x_{i,1} - x_{i,2} \leq 0$$

$$x_{i,3} - x_{i,1} + x_{i,2} \leq 0$$

$$x_{i,3} + x_{i,1} - x_{i,2} \leq 0$$

$$x_{i,3} + x_{i,1} + x_{i,2} \leq 0$$

The indices 1, 2, and 3 refer to the x, y, and z coordinates, respectively.

This problem can be solved with the nonlinear constrained solver fmincon in Optimization Toolbox.

Linear inequality equations (Aineq*x = bineq) are defined for E electrons.

E = 25; % number of electrons
B = [ 1  1  1;
     -1  1  1;
      1 -1  1;
     -1 -1  1];
electronProblem.Aineq = sparse(4*E,3*E);
for i = 1:E
    electronProblem.Aineq(4*i-3:4*i, 3*i-2:3*i) = B;
electronProblem.bineq = sparse(4*E,1);

Nonlinear constraints are defined in a nlConstraint function

type nlConstraint
electronProblem.nonlcon = @nlConstraint;
function [c,ceq] = nlConstraint(x)
        E = length(x)/3;
        % inequality constraint
        c = zeros(E,1);
        i = 1:E;
        c = (x(3*i-2).^2 + x(3*i-1).^2 + (x(3*i)+1).^2 - 1).';
        % equality constraint (none)
        ceq = [];

The starting condition will set all electrons to initially lie in the center of the conducting body

% Generate a pseudo-random set of electrons using sobol sequence
x0 = net(sobolset(3),E)-0.5;
% Normalize to 1/2 line within a 1/2 radius sphere
x0(:,3) = x0(:,3) - 2;
x0 = reshape(x0',E*3,1)/2;
electronProblem.x0 = x0;

The problem is defined. We'll now define the solver and options to use.

electronProblem.options = optimset('Algorithm','interior-point',...
electronProblem.solver = 'fmincon';

save electronProblem.mat electronProblem E

Run the Solver fmincon

Solve the problem and show the solution and solution time.

x = fmincon(electronProblem);
Elapsed time is 50.652165 seconds.

Run in Parallel

Run the problem on 4 cores

electronProblem.options.UseParallel = 'Always';
x = fmincon(electronProblem);
matlabpool close
Error using ==> matlabpool at 125
Could not find a scheduler when using the configuration 'speedy-4core'.
To learn how to modify your configuration to uniquely identify your scheduler execute the command:
  <a href="matlab:docsearch('Programming with User Configurations')">docsearch('Programming with User Configurations')</a>

Error in ==> electronProblemOptimization at 106