(Slight) Constraint violation without warnings.

[Luca]

Dear Johan,

First of all thanks for the great work you are doing with YALMIP. 

I have encountered an annoying problem in my code, and I do not know exactly why this is happening (without warnings from the solver or YALMIP). 

I have the following optimization problem,which can be thought as a sizing and control (over an interval 1:T) of an energy generation unit, in which a constraint is (slightly) violated all the time. 

to summarize briefly:

- I use interp2 to estimate the variable z based on the values of variables x and y, with a predetermined data set. 

- Beside the boundaries, the constraints can be expressed as  D(t)<=z(t)<=y for every t=1:1:T.

- D(t)<=z(t) is always slightly violated. 

I have the feeling that this has to do with the numerical precision of the solver, but the constraint violation is in the order of 10^-1. (e.g. 15.5<=15.4), which seems quite high.

Also, I do not get any warnings or errors from YALMIP or Gurobi. Yet, the constraint violation is there when i check the constraints with check().

Here attached a simple example code. 

Do you have any suggestions on how to get rid or reduce this problem?

Thank you.

Best regards

[Johan Löfberg]

Reproduced here. Not solver issue as exactly same solution with cplex.

Do you know if the solution really is violating the constraint (it could be that the computation inside the interp2 operator when being asked by check is wrong). When I look at the plot, it really looks like the points are on the surface, and not violating it with up to -0.14 as indicated by check

[Johan Löfberg]

Yes, this is an issue with the post-evaluation of the interp2 operator

When YALMIP sets up the interp2 model for sos2 models, it does this by triangulating the squares represented by the grid, and then doing a linear model in each triangle.

However, when interp2 is executed for a particular value of the argument such as doing value which is used by check, it is computed directly using the interp2 operator i.e. it does interp2(X,Y,Z,value(x),value(y)). This operator is however not based on triangulation, but it does bilinear interpolation over every square, and that does not yield the same interpolant (but a worse)

[X,Y] = meshgrid(0:1,0:1);

Z = X + Y;

Z(2,2) = 5;

mesh(X,Y,Z)

sdpvar x y z

optimize([z == interp2(X,Y,Z,x,y,'milp'), x == .5, y == .4])

value(z) % As expected, linear interpolation over a trangle should be exact

value(interp2(X,Y,Z,x,y,'milp')) % This post-evaluation does bilinear

value(interp2(X,Y,Z,value(x),value(y),'milp')) % this is what is done

I'll have to think about how I can combat this

[Johan Löfberg]

A fix is to change line 8 in interp2_internal to

            varargout{1} = griddata(varargin{2},varargin{3},varargin{4},varargin{1}(1),varargin{1}(2));

            return

effectively, when we use the YALMIP-specific sos2/graph/lp arguments in interp2, the interp2 operator corresponds to the griddata command

[Luca]

Dear Johan,

Thanks for the answer and the analysis!

I have changed the lines in the interp2_internal operator and, as you can see from the first 2 figures, it seems to work. 

Concerning my example code, it also works except for one point (see figure D_z). Is it the result of an interpolation error due to the Sol_z point (see fig z_3D)  is in the mesh square where there is a discontinuity in the function, therefore the interpolation error is larger?

If this is the case, it is not a code problem anymore, but only a matter of having a finer mesh grid of data, am I correct?

Thank you.

Best regards

[Johan Löfberg]

ok, missed than one point was off. it's a point on the ridge, and it's a discreprency in  how griddata/interp2 does the interpolation, vs a linear interpolation over triangles which you get using sos2 yalmip. The problem is solved as expected, it's just the post-evaluation which is wrong. 

To really get what the solver works with, instead of doing z = interp2..., introduce a new variable z and add the constraint z == interp2..., and then use z as the variable to look at afterwords

[Luca]

Yes, it works as expected. Thank you a lot for your time.