# Parfor

In this post we'll explore how to use Matlab's parfor with a CasADi nonlinear program.

The nonlinear program (NLP) of interest is the following: \begin{align} \displaystyle \underset{x,y} {\text{minimize}}\quad &\displaystyle (1-x)^2+(y-x^2)^2 \newline \text{subject to} \, \quad & x^2+y^2 \leq r \end{align} Note that $r$ is a free parameter. Our goal is to collect each solution of the NLP as we loop over $r$. We could imagine performing such task when tracing a pareto front of a multi-objective optimization, for which the parameter would be the a weight to combine those objectives.

# Problem construction inside the loop

The first way to work with parfor and CasADi is to do the problem construction entirely in the loop:

N = 250;

fsol = zeros(1,N);
rs = linspace(1,3,N);

parfor i=1:N
r = rs(i);

v = [x;y];
f = (1-x)^2+(y-x^2)^2;
g = x^2+y^2;
nlp = struct('x', v, 'f', f, 'g', g);

% Create IPOPT solver object

% Solve the NLP
res = solver('x0' , [2.5 3.0],...      % solution guess
'lbx', -inf,...           % lower bound on x
'ubx',  inf,...           % upper bound on x
'lbg', -inf,...           % lower bound on g
'ubg',  r);               % upper bound on g

% Store the solution
fsol(i) = full(res.f);
end

plot(rs,fsol,'o')

We obtain a nice plot:

# Problem construction outside the loop

Constructing CasADi expressions and Functions incurs an initialization cost. Therefore, we always advise to construct the solver fully outside of the loop, and simply make numerical calls to it in the loop.

This advice collided with our limited support for parfor so far. CasADi 3.5 comes with improved support. The loop body may now contain numerical calls to CasADi Functions defined outside that loop:

x = casadi.SX.sym('x');

v = [x;y];
f = (1-x)^2+(y-x^2)^2;
g = x^2+y^2;
nlp = struct('x', v, 'f', f, 'g', g);

% Create IPOPT solver object

tic
parfor i=1:N
% Solve the NLP
res = solver('x0' , [2.5 3.0],...      % solution guess
'lbx', -inf,...           % lower bound on x
'ubx',  inf,...           % upper bound on x
'lbg', -inf,...           % lower bound on g
'ubg',  rs(i));           % upper bound on g

% Store the solution
fsol(i) = full(res.f);
end
toc

plot(rs,fsol,'o')

Note 1: for this example the parameter entered trivially through 'ubg'. In general, you may need to work with a parametric NLP:

p = casadi.SX.sym('p');
g = x^2+y^2-p;
nlp = struct('x', v, 'p', p, 'f', f, 'g', g);

res = solver('x0' , [2.5 3.0],...      % solution guess
...
'p',  rs(i));

Note 2: constructing CasADi Functions inside the loop with CasADi symbols declared outside the loop is forbidden, but inefficient anyway.

# Implementation details and Python

We did not actually code anything parfor-specific in CasADi. The reason that CasADi 3.5 supports the construct above now is because we implemented serialization/deserialization of CasADi Functions.

Parfor uses this under the hood to transfer problems/results to/from different threads. Such mechanism is often encounter for parallelization, for example in Python's multiprocessing module. In fact, the above code can easily be rewritten for Python:

from multiprocessing import Pool

N = 250
rs = np.linspace(1,3,N)

x = SX.sym('x')
y = SX.sym('y')

v = vertcat(x,y)
f = (1-x)**2+(y-x**2)**2
g = x**2+y**2
nlp = {'x': v, 'f': f, 'g': g}

# Create IPOPT solver object

def optimize(r):
res = solver(x0=[2.5,3.0],      # solution guess
lbx=-inf,          # lower bound on x
ubx=inf,           # upper bound on x
lbg=-inf,          # lower bound on g
ubg=r)         # upper bound on g

return float(res["f"])

b = Pool(2)
fsol = b.map(optimize, rs)