interp2

interp2 overloads interp2, with additional flags for a creating mixed-integer sos2-based approximation

Syntax

y = interp2(xdata,ydata,zdata,x,y,'method')

Example (nonlinear optimization)

Create a nonlinear function over a grid and find the global minima of a spline interpolated version of the nonlinear function, using the global solver BMIBNB (naturally silly to do as we simply can use the data to find the global minimum)

[X,Y] = meshgrid(-3:.5:3,-3:.5:3);
Z = 8*(X-.5).^2 + 100*sin(5*X).^2+((Y).^2-2*X.^2).^2+15*(Y+3).^2;

sdpvar x y
z = interp2(X,Y,Z,x,y,'spline');
optimize([-3 <= [x y] <= 3],z,sdpsettings('solver','bmibnb'));
mesh(X,Y,Z);hold on;plot3(value(x),value(y),value(z),'*')

Of course, in this case we could have worked with the nonlinear function directly.

sdpvar x y
z = 8*(x-.5)^2 + 100*sin(5*x)^2+((y)^2-2*x.^2).^2+15*(y+3).^2;
optimize([-3 <= [x y] <= 3],z,sdpsettings('solver','bmibnb'));
plot3(value(x),value(y),value(z),'r*')

Example (linear optimization)

With the flags ‘sos2’ or ‘milp’ a sos2 representation of the piecewise affine function interpolating the supplied data points is created, and thus requires a sos2-capable mixed-integer linear solver.

Create a more challenging nonlinear function than above

N = 40;
M = 40;
[xi,yi] = meshgrid(-1:2/(M-1):1,-1:2/(N-1):1);
zi = (xi).^4 + (yi).^4 + xi.*yi + sin(xi.*yi*4*pi).*(1 + xi + yi);
clf;mesh(xi,yi,zi)

We can try to find a minimizer by overloading an interpolation (here using splines) of the data-points. This will lead to a nonlinear program

sdpvar x y z
z = interp2(xi,yi,zi,x,y,'spline');
optimize([],z)
value([x y z])

Almost surely, a local nonlinear solver will fail to find a good solution. The same happens when we simply use the nonlinear function directly

z =  (x).^4 + (y).^4 + x.*y + sin(x.*y*4*pi).*(1 + x + y);
optimize([-1 <= [x y] <= 1],z)
value([x y z])

A global solver finds the globally optimal solution easily

z =  (x).^4 + (y).^4 + x.*y + sin(x.*y*4*pi).*(1 + x + y);
optimize([-1 <= [x y] <= 1],z,sdpsettings('solver','bmibnb'))
value([x y z])

An alternative way to attack this with a hope of a global solution is to approxmate the objective using a sos2 representation and thus solve the problem using a mixed-integer linear solver.

z = interp2(xi,yi,zi,x,y,'sos2');
optimize([],z)
value([x y z])
clf;mesh(xi,yi,zi);hold on
l = line([value(x) value(x)],[value(y) value(y)],[5 value(z)]);
set(l,'LineWidth',4)

THe precision in this solver is limited by the granularity of our grid. Let us start our local solver in the computed point and try to improve it

z =  (x).^4 + (y).^4 + x.*y + sin(x.*y*4*pi).*(1 + x + y);
optimize([-1 <= [x y] <= 1],z,sdpsettings('usex0',1))
value([x y z])
l = line([value(x) value(x)],[value(y) value(y)],[5 value(z)]);
set(l,'LineWidth',4,'Color','red');

Implementation

In the nonlinear case when the operator is implemented using the callback framework, the operator does not compute derivatives, which can lead to slow computations (and in general, computing a single function value is expensive). All bounding in the global solver is done using sampling heuristics.