Infeasible problems occurred when using 'sort' operator

[Zhi Shang]

Dear Professor,

I am a beginner of yalmip. 

My problem is to minmize the objective with constraint that the neighboring variables are in correct order, so the 'sort' operator is used as following code. However, the solution information shows 'Infeasible problem'.

Could you please tell me the reason of error or give me some advices?

Thanks in advance.

Best regards

========================================

[MM,JJ]=size(Force);

Alpha=0.3;

lb=ones(MM,1)*min(Y(:,Cpts_idx));                                           

ub=ones(MM,1)*max(Y(:,Cpts_idx));

P = sdpvar(MM,JJ);                                                          

F=[];

obj=0;

for i=1:MM

    for j=1:JJ

        Ec=0;

        Eg=0;

        Eall=0;

        nbors=nbor_idx{i,j}';                          % nbors is corresponded to the row index of neighboring points which are sorted in correct order

        [~,loc]=sort(P(nbors,j));          

        F=[F,loc==nbors];                              % the constraint is to ensure correct order of optimized values

        F=[F,lb(i,j)<=P(i,j)<=ub(i,j)];

        

        Ec==abs(P(i,j)-Y(i,Cpts_idx(j)));

        Eg=-Force(i,j)*(P(i,j)-Y(i,Cpts_idx(j)));

        Eall=Alpha*Ec+(1-Alpha)*Eg;

        obj=obj+Eall;

    end

end

sol = optimize(F,obj);

[Johan Löfberg]

As you haven't given us the data, it is impossible to answer anything. Typical mistake is the P is square (n=m) and thus defined as symmetric

Are you sure you need the sort operator (which is extremely complicated to implement, leads to MILP constraints). Don't you have something like some indicies idx, and you want some P(idx,j) to be sorted? If so, you should simply write diff(P(idx,j))>=0 which is a simple linear constraint

Having absolute values in equality constraints is also very bad, as it leads to MILP models. You should simply write obj = obj + alpha*abs(P...) + (1-alpha)*(-Force...). Since you minimize, YALMIP will be able to form a linear model of the objective involving the convex abs operator (for alpha>=0)

[Zhi Shang]

Dear Professor,

Thanks a lot for your reply, the diff operator is suitable for the problem.

I modify the code as you advised, but I don't understand the symmetric problem about P, as P is not square and all variables defined by P are totally different.

Now another error arise when I change the row number of 'Force' (MM) from 200 to 300,  the solution informations show 'Unbounded objective function 'or 'infeasible problem' . No errors occured when MM=200. 

The data is attached. Thanks in advance.

Best regards.

=====================================================

tic

clear 

clc

close all

load Force.mat

load Y.mat

load nbor_idx.mat

CPts=2;

NN=3;

[MM,JJ]=size(Force);

Alpha=0.3;

solution=[];

for var=1:NN-1                                              

    yalmip('clear')

    P = sdpvar(MM,CPts);  

    E=sdpvar(MM,CPts);

    F=[];

    obj=0;

    for i=1:MM

        for j=1:CPts

            idx=(var-1)*(CPts+1)+j+1; 

            idx2=(var-1)*CPts+j;

            lb=min(Y(:,idx));                                           

            ub=max(Y(:,idx));

            nbors=nbor_idx{i,idx2};                           

            F=[F,implies(diff(P(nbors,j))>=1e-5,E(i,j)==0)];

            F=[F,1<=E<=10];

            F=[F,lb<=P(i,j)<=ub];

            obj=obj+Alpha*abs(P(i,j)-Y(i,idx))+(1-Alpha)*(-Force(i,idx2)*(P(i,j)-Y(i,idx))+E(i,j)*1e10);   %E is infinite if P(nbors,j) is in incorrect order        

        end         

    end

    sol = optimize(F,obj);

    if sol.problem == 0

         solution = [solution value(P)]

    else

         display('Something went wrong!');

         sol.info

         yalmiperror(sol.problem)

         break;

    end

end

toc

[Johan Löfberg]

Since for your particular data Force is non-square, P will be non-square and you will not have any problems with symmetry. Never hurts to use the 'full' flag though in case you run on a square problem later

However, your model is bad

To begin with, this is nonsense, as E never can be zero according to second constraint

F=[F,implies(diff(P(nbors,j))>=1e-5,E(i,j)==0)];

F=[F,1<=E<=10];

Secondly, you have data Force in the objective, and some elements of that matrix is inf, meaning most solvers will crash and/or get confused about the objective 

Third, having the weight 1e10 on some parameters is a recipie for numerical disaster, as the objective will be huge and the solver easily can get lost and think the problem simply is infeasible

[Johan Löfberg]

and this comment

%E is infinite if P(nbors,j) is in incorrect order

is not consistent with your model. What you (try to) obtain in your constraint is that E is zero when it is in correct order, nothing else. The model does not exclude the case that it is in correct order, and E is zero (well, your bound constraints are inconsistent so E can never be zero, but that was another issue)

[Zhi Shang]

Thanks! the infinite element in Force will be dealed later.

what about the following code?As I still need to consider E is infinite if  P(nbors,j) is in incorrect order,or can you give me some advices to handle the issue?

 F=[F,implies(diff(P(nbors,j))>=1e-5,E(i,j)==0)];
 F=[F,implies(diff(P(nbors,j))<=-1e-5,E(i,j)==10)];
 F=[F,0<=E<=10];

[Johan Löfberg]

You cannot have inf in the data defining the model (in Force). The problem is not well-defined to the solver. For instance, minimize inf*x subject to x>=0 does not have any solution really, as inf*0 is not a number. You should do P(isinf(Force)) = Y(isinf(force)), or something like that, not sure what you're trying to encode 

And E being infinite is not in any sense encoded by your model. You are now trying to say if P is an increasing sequence, then E is 0, and if P is a decreasing sequence, E is 10. All other cases leaves E undefined (but at 0 due to optimality as you penalize it). 

[Johan Löfberg]

Aren't you simply trying to encode a slack penalty? No reason to work with implications etc

 F=[F,E(i,j)>=0, diff(P(nbors,j)>=1e-5-E(i,j)];

and then add sum(E) or E(:)'*E(:) (weighted with something more numerically sensible than 1e10) to the objective.

[Zhi Shang]

So how to define E for other cases? for example, also define E is 10 if P is in any order that not consisent with incresing sequence? 

[Johan Löfberg]

If you absolutely must have that combinatorial form of the slack, instead of a simple cheap classical linear slack, I would use a binary variable E (which thus is either 0 or 1, scale to 10 if you wish...)

 F=[F,diff(P(nbors,j)>=1e-5-E(i,j)*(max(ub))];

If they are all sorted, E can be 0. If any element violates sorting, E will be forced to be 1. Since the difference between elements cannot be larger than max(ub), the big-M constant can be taken as such