## Create Efficient Optimization Problems

When a linear problem has integer constraints, `solve`

calls
`intlinprog`

to obtain the solution. For
suggestions on obtaining a faster solution or more integer-feasible points, see Tuning Integer Linear Programming.

Before you start solving the problem, sometimes you can improve the formulation of your problem constraints or objective. Usually, the software can create expressions for the objective function or constraints more quickly in a vectorized fashion rather than in a loop. This speed difference is especially large when an optimization expression is subject to automatic differentiation; see Automatic Differentiation in Optimization Toolbox.

Usually, it is easier to create an efficient loop by writing it as a separate
function, as described in Create for Loop for Static Analysis. In this case,
create the optimization expression containing the function by using `fcn2optimexpr`

, as explained in Static Analysis of Optimization Expressions.

Suppose that your objective function is

$$\sum _{i=1}^{30}{\displaystyle \sum _{j=1}^{30}{\displaystyle \sum _{k=1}^{10}{x}_{i,j,k}{b}_{k}{c}_{i,j}}}},$$

where `x`

, `b`

and `c`

are
optimization variables. Two general ways to formulate this objective function are as
follows:

Use a

`for`

loop.x = optimvar('x',30,30,10); b = optimvar('b',10); c = optimvar('c',30,30); tic expr = optimexpr; for i = 1:30 for j = 1:30 for k = 1:10 expr = expr + x(i,j,k)*b(k)*c(i,j); end end end toc

Elapsed time is 307.459465 seconds.

Here,

`expr`

contains the objective function expression. To have the software evaluate this expression efficiently, write it as a separate function:function expr = loopme(x,b,c) expr = 0; ni = size(x,1); nj = size(x,2); nk = size(x,3); for i = 1:ni for j = 1:nj for k = 1:nk expr = expr + x(i,j,k)*b(k)*c(i,j); end end end end

Include the expression using

`fcn2optimexpr`

.tic expr = fcn2optimexpr(@loopme,x,b,c,OutputSize=[1,1]); toc

Elapsed time is 23.888763 seconds.

Include the objective function in the problem.

`problem = optimproblem("Objective",expr);`

Use a vectorized statement. Vectorized statements generally run faster than a

`for`

loop that has not been modified for static analysis. You can create a vectorized statement in several ways.Expand

`b`

and`c`

. To enable term-wise multiplication, create constants that are the same size as`x`

.tic bigb = reshape(b,1,1,10); bigb = repmat(bigb,30,30,1); bigc = repmat(c,1,1,10); expr = sum(sum(sum(x.*bigb.*bigc))); toc

Elapsed time is 0.013631 seconds.

Loop once over

`b`

.tic expr = optimexpr; for k = 1:10 expr = expr + sum(sum(x(:,:,k).*c))*b(k); end toc

Elapsed time is 0.044985 seconds.

Create an expression by looping over

`b`

and then summing terms after the loop.tic expr = optimexpr(30,30,10); for k = 1:10 expr(:,:,k) = x(:,:,k).*c*b(k); end expr = sum(expr(:)); toc

Elapsed time is 0.039518 seconds.

Observe the speed difference between a vectorized and nonvectorized implementation of the example Constrained Electrostatic Nonlinear Optimization Using Optimization Variables. This example was timed using automatic differentiation in R2020b.

N = 30; x = optimvar('x',N,'LowerBound',-1,'UpperBound',1); y = optimvar('y',N,'LowerBound',-1,'UpperBound',1); z = optimvar('z',N,'LowerBound',-2,'UpperBound',0); elecprob = optimproblem; elecprob.Constraints.spherec = (x.^2 + y.^2 + (z+1).^2) <= 1; elecprob.Constraints.plane1 = z <= -x-y; elecprob.Constraints.plane2 = z <= -x+y; elecprob.Constraints.plane3 = z <= x-y; elecprob.Constraints.plane4 = z <= x+y; rng default % For reproducibility x0 = randn(N,3); for ii=1:N x0(ii,:) = x0(ii,:)/norm(x0(ii,:))/2; x0(ii,3) = x0(ii,3) - 1; end init.x = x0(:,1); init.y = x0(:,2); init.z = x0(:,3); opts = optimoptions('fmincon','Display','off'); tic energy = optimexpr(1); for ii = 1:(N-1) jj = (ii+1):N; % Vectorized tempe = (x(ii) - x(jj)).^2 + (y(ii) - y(jj)).^2 + (z(ii) - z(jj)).^2; energy = energy + sum(tempe.^(-1/2)); end elecprob.Objective = energy; disp('Vectorized computation time:') [sol,fval,exitflag,output] = solve(elecprob,init,'Options',opts); toc

Vectorized computation time: Elapsed time is 1.838136 seconds.

tic energy2 = optimexpr(1); % For nonvectorized comparison for ii = 1:(N-1) for jjj = (ii+1):N; % Not vectorized energy2 = energy2 + ((x(ii) - x(jjj))^2 + (y(ii) - y(jjj))^2 + (z(ii) - z(jjj))^2)^(-1/2); end end elecprob.Objective = energy2; disp('Non-vectorized computation time:') [sol,fval,exitflag,output] = solve(elecprob,init,'Options',opts); toc

Non-vectorized computation time: Elapsed time is 204.615210 seconds.

The vectorized version is about 100 times faster than the nonvectorized version.

Compare with the static analysis version. The code for the `myenergy`

function appears at the end of this example.

elecprob.Objective = fcn2optimexpr(@myenergy,N,x,y,z);

Solve the problem by calling `solve`

. Time the solution.

disp('Static analysis computation time:') tic [sol,fval,exitflag,output] = solve(elecprob,init,'Options',opts); toc

Static analysis computation time: Elapsed time is 1.293513 seconds.

Static analysis provides the lowest computation time.

This code creates the `myenergy`

function.

function energy = myenergy(N, x, y, z) energy = 0; for ii = 1:(N-1) for jj = (ii+1):N energy = energy + ((x(ii) - x(jj))^2 + (y(ii) - y(jj))^2 + (z(ii) - z(jj))^2)^(-1/2); end end end