Losing Accuracy with Piece-Wise Linear Term on Objective Function

[Lun Yang]

Dear Prof. Löfberg,

I am trying to solve a optimization problem with very simply quadratic objective function plus a piece-wise linear term for variable 1 only. For simplicity, I am setting the piece-wise linear term to be y = 0, which shouldn't affect the solution at all. 

Please see the code below. Interestingly, for number of pieces equals to 50, I am getting the exact solution of [0.2 0.2 0.2 0.2 0.2], which is correct. But the when I set numberOfPieces to 500, I am getting [0.196 0.201 0.201 0.201 0.201], which is a 2% difference for variable 1.

Could you please help me take a look, if my implementation is correct? Is there a way to retain accuracy when the numberOfPieces is getting larger? I am using cplex, btw. 

Thank you very muck.

Sincerely,
Lun Yang

n = 5;

S = eye(n);

w = sdpvar(n,1);

F = [sum(w)==1, 1>=w>=0];

Obj = w'*S*w;

% Adding Piece-Wise Linear Term on Objective Function for Variable 1, flat line at y = 0 actually

numberOfPieces = 50;                             % Please change this to 500 and rerun

xs = linspace(0,1,numberOfPieces+1);   % x points

ys = zeros(1,numberOfPieces+1);           % y points

slopes = zeros(1,numberOfPieces);        % slopes

Piecewise = binvar(numberOfPieces,1);  % Binary variables to indicate which piece w stays

F = [F, sum(Piecewise) == 1];                  % Only 1 piece is allowed

FreeVariables = sdpvar(1,1);                   % Value of the Piece-Wise Linear Term 

for i = 1:numberOfPieces

    LinearTerm = 0*w(1);                           % Should be all 0....

    PWLConstraint = implies(Piecewise(i),[FreeVariables == LinearTerm,xs(i)<=w(1)<=xs(i+1)]);

    F = [F, PWLConstraint];

end

F = [F, -1e5<=FreeVariables<=1e5];     % Artificial bound to stabilize solver

Obj = Obj + FreeVariables;

% Solve

optimize(F,Obj,sdpsettings('solver','cplex'));

value(w)

[Johan Löfberg]

cannot reproduce (cplex 12.8 win64)

CPXPARAM_MIP_Display                             1

Tried aggregator 1 time.

MIQP Presolve eliminated 14 rows and 0 columns.

Reduced MIQP has 2000 rows, 506 columns, and 4501 nonzeros.

Reduced MIQP has 500 binaries, 0 generals, 0 SOSs, and 0 indicators.

Reduced MIQP objective Q matrix has 5 nonzeros.

Presolve time = 0.00 sec. (2.76 ticks)

Probing time = 0.01 sec. (6.57 ticks)

Cover probing fixed 1 vars, tightened 0 bounds.

Tried aggregator 1 time.

MIQP Presolve eliminated 1000 rows and 0 columns.

Reduced MIQP has 1000 rows, 506 columns, and 2501 nonzeros.

Reduced MIQP has 500 binaries, 0 generals, 0 SOSs, and 0 indicators.

Reduced MIQP objective Q matrix has 5 nonzeros.

Presolve time = 0.00 sec. (2.07 ticks)

Probing time = 0.02 sec. (5.31 ticks)

MIP emphasis: balance optimality and feasibility.

MIP search method: dynamic search.

Parallel mode: deterministic, using up to 4 threads.

Node log . . .

Best integer =   1.000000e+00  Node =       0  Best node =   0.000000e+00

Best integer =   2.490050e-01  Node =       0  Best node =   2.000000e-01

Best integer =   2.004050e-01  Node =       0  Best node =   2.000000e-01

Best integer =   2.003200e-01  Node =       0  Best node =   2.000000e-01

Best integer =   2.000450e-01  Node =       0  Best node =   2.000000e-01

Best integer =   2.000000e-01  Node =       0  Best node =   2.000000e-01

>> numberOfPieces 

numberOfPieces =

   500

>> value(w)

ans =

     2.000000000000000e-01

     2.000000000000000e-01

     2.000000000000000e-01

     2.000000000000000e-01

     2.000000000000000e-01

[Johan Löfberg]

btw, if you want to do pwa models, you should probably use sos2 constructions

convenient quick approach is to use interp1

https://yalmip.github.io/command/interp1/

[Lun Yang]

Dear Prof. Löfberg,

Thank you very much for the interp1 suggestion. I have a follow up question, what should I do if I want to loop the call when x y points are changing . When I am trying to use optimizer to speed up the the call, I am getting the following error message, please see the code and error below,

n = 5;

S = eye(n);

w = sdpvar(n,1);

F = [sum(w)==1, 1>=w>=0];

Obj = w'*S*w;

numberOfPieces = 500;

xs = sdpvar(numberOfPieces,1);        % x points

ys = sdpvar(numberOfPieces,1);        % y points

FreeVariables = interp1(xs,ys,w(1),'milp');

Obj = Obj + FreeVariables;

savedModel = optimizer(F,Obj,sdpsettings('solver','cplex'),{xs,ys},w);

%%%%%%%%%%%%%%%%%%%%%% Error Message %%%%%%%%%%%%%%%%%%%%%

Error using sdpvar/model (line 66)

Failed when trying to create a model for the "interp1_internal" operator

Error in expandrecursive (line 19)

    [properties,F_graph,arguments,fcn] = model(variable,method,options,allExtStruct(ext_index),w);

Error in expandmodel>expand (line 416)

            [F_expand,failure,cause] =

            expandrecursive(recover(variables(i)),F_expand,extendedvariables,monomtable,variabletype,where,level+1,options,method,[],'convex',allExtStruct,w);

            

Error in expandmodel (line 248)

    [F_expand,failure,cause] =

    expand(index_in_extended,variables,h,F_expand,extendedvariables,monomtable,variabletype,'objective',0,options,method,[],allExtStructs,w);

    

Error in compileinterfacedata (line 116)

    [F,failure,cause,operators] = expandmodel(F,h,options);

Error in export (line 135)

    [interfacedata,recoverdata,solver,diagnostic,F] = compileinterfacedata(F,[],logdetStruct,h,options,findallsolvers);

Error in optimizer (line 187)

    [aux1,aux2,aux3,model] = export((x == repmat(pi,nIn*mIn,1))+Constraints,Objective,options,[],[],0);

And I have to comment out line 65 of sanity checks in interp1 to get this point. Could you please help take a look? Thank you very much. 

Sincerely,

Lun Yang

[Johan Löfberg]

interp1 won't work in optimizer. it's a way too complex operator

[Lun Yang]

Thank you