coneprog: Infeasibilities are satisfied in last step from iterative display but exitflag is still -2 (NoFeasiblePointFound)

I am trying to solve a trajectory optimization problem with coneprog under the problem-based approach and am having convergence issues. The last iteration of the iterative display appears to show that the problem is feasibly solved:

Iter           Fval  Primal Infeas    Dual Infeas    Duality Gap    Time
1   0.000000e+00   2.989297e-01   2.661135e-01   2.396722e-03    0.01
2   1.671582e+00   1.371375e-01   1.220827e-01   1.099524e-03    0.02
3   3.247714e+00   4.683911e-02   4.169716e-02   3.755409e-04    0.02
4   5.466823e+00   2.962828e-02   2.637572e-02   2.375500e-04    0.02
5   6.748515e+00   1.679716e-02   1.495319e-02   1.346743e-04    0.02
6   3.367365e+00   1.788542e-03   1.592198e-03   1.433995e-05    0.03
7   5.484853e+00   4.005449e-06   3.565735e-06   3.211443e-08    0.03
8   5.484859e+00   6.005173e-09   5.346126e-09   4.817329e-11    0.04
9   5.484864e+00   1.200136e-11   1.091778e-11   1.141043e-13    0.04

However, I am still getting a -2 exit flag and the message:

Problem is infeasible.

Is there anyway to extract anymore information about why this solution is infeasible? The functions where the solver operations occur are .p files so I can't find what's going wrong on the inside. Also, when a solution is returned infeasible all output fields are empty except the error message. I am using optimality and constraint tolerances of 1e-8 for reference, but even when I drop them down to 1e-1 just for the sake of trying to get a feasible solution returned it still does not work. Does this mean maybe that the solution converges but upper/lower bounds are violated? Any advice would be much appreciated.

[Torsten]

Can't you provide a feasible solution by "precomputing" your constraints ?

[Matt J]

I suggest setting MaxIterations=9 which, judging from the iterative display that you posted, should force the algorithm to terminate at the 9th iteration with exitflag=0. Then the output arguments from coneprog() should be non-trivial and you can better inspect whether the solution x up to that point is truly feasible or infeasible.

[Nick]

@Torsten Maybe I am misunderstanding what you mean, but I don't see the benefit of precomputing the constraints with my initial guess because it will surely be infeasible. And I don't have access to the output from infeasible runs of coneprog.

@Matt J I tried this and get the exitflag=0, but all the outputs are empty just like if the exitflag were to be -2. I thought maybe that is because I am using the problem-based approach with the solve function, but I put a breakpoint inside coneprog right after the .p function file call at line 98 and all of its outputs are also empty.

[Torsten]

@Nicholas Nurre

I mean: Can't you use "fsolve" to determine a feasible point before you call "coneprog" ?

[Matt J]

Yes, or use fmincon to find the nearest feasible point to your initial point. If fmincon comes up empty, that will add to the evidence that your problem is infeasible. Regardless, with fmincon, you will definitely be able to track the iteration history using the OutputFcn option.

[Nick]

Ok I understand now. I'll try this. Thanks for the advice @Torsten @Matt J. I still find it weird the solver won't output any other information even when the maximum number of iterations has been reached

[Bruno Luong]

https://www.mathworks.com/help/optim/ug/coneprog.html#mw_04f9ab9a-aa5e-4961-bb6e-ef0f94c0602e

"exitflag — Reason coneprog stopped -2: No feasible point was found."

What else information do you need?

If you are looking for a fish that living in mars and MATLAB tells you it doesn't exist, do you want to know it because you ask for a "fish" or because you ask for "living in mars"? The right answer is neither, but both combined is impossible.

[Torsten]

The fact that coneprog cannot find a feasible point does not necessarily mean that no feasible point exists ...

[Bruno Luong]

Indeed but OP asks why the solver does not tell more details when it says it fails to detect feasible point. What the solver can tell more????

[Torsten]

I would have expected that coneprog at least reports the actual x-value after "MaxIterations" iterations.

[Bruno Luong]

If that is the case it's easy, just make a wrapper around the p function and capture the last call input. I still don't understand how he would use it. But can we just let him answer the question?

[Nick]

As for my original question I understand that my problem is infeasible. I am now wondering why when even if a problem is retruned infeasible (or even if the max iterations is hit) the solver returns empty outputs. I was hoping to maybe gain insight as to why my problem is infeasible (if it could even help) by looking at the outputs returned by the solver like how other solvers like fmincon will return outputs no matter the exitflag. This behavior is reproducible if you want to try yourself with the example used by @Alan Weiss in his answer (url: https://www.mathworks.com/help/optim/ug/discretized-optimal-trajectory-cone-programming.html ). If you adjust the Maxiterations to something like 1, all output values are empty and no information would be gained. It does not seem that wrapping the p function (interiorPointMethod inside of coneprog) will help because the p function is only called once and the iterates are performed inside of it and it still returns empty outputs. This is in R2021b by the way. Maybe this behavior is by desisn though. I don't know

[Alan Weiss]

The solver is designed to return empty outputs for infeasible problems. Otherwise, people are tempted to use the outputs as solutions, which they are not.

If you want to investigate why your problem is infeasible, you can try relaxing constraints until you get a feasible solution. There are several ways of relaxing constraints. You can give a linear constraint more room, such as writing A*x <= b + 10, or some other positive added constant other than 10. For cone constraints you can include them one at a time until your problem becomes infeasible. Or you can relax them similarly to the linear constraint: \norm(A*x - b) <= d*x - gamma + 10, or some other positive added constant instead of 10. Sometimes it is an interaction between bounds and other constraints that makes a problem infeasible.

You mentioned that you are doing a trajectory optimization problem. If it is similar to the one in the documentation, then I suggest that you model your problem on the problem in the documentation, which works, and change things one at a time until you reach an issue.

Good luck,

Alan Weiss

MATLAB mathematical toolbox documentation

[Nick]

Ok I see. Thank you for such quick and helpful responses!

[Torsten]

@Bruno Luong

Can you show how to "make a wrapper around the p function and capture the last call input" ? I honestly don't know how this can be done in the case of "coneprog".

Answer 1

[Alan Weiss]

I think that you have a misunderstanding about what the iterative display shows. For example, consider the documentation example about trajectory optimization Discretized Optimal Trajectory, Problem-Based. If you insert iterative display in an infeasible problem there, you get the following display:

p0 = [0 0 0];
pF = [5 10 3];
myprob = setupproblem(1);
opts = optimoptions("coneprog",Display="iter");
[solm,fvalm,eflagm,outputm] = solve(myprob,Options=opts);

Solving problem using coneprog.
Iter           Fval  Primal Infeas    Dual Infeas    Duality Gap    Time
   1   9.800000e-01   3.681614e-01   3.311258e-01   3.148330e-03    0.10
   2   2.051010e+00   1.708611e-01   1.536731e-01   1.461117e-03    0.12
   3   1.886027e+00   8.505728e-02   7.650085e-02   7.273667e-04    0.12
   4   2.587887e+00   1.176815e-02   1.058432e-02   1.006352e-04    0.12
   5   2.287761e+00   2.885400e-05   2.595136e-05   2.467444e-07    0.13
   6   2.286884e+00   2.348651e-07   3.892967e-08   3.701098e-10    0.13
Problem is infeasible.

function trajectoryprob = setupproblem(T,air)
    if nargin == 1
        air = false;
    end
    N = 50;
    g = [0 0 -9.81];
    p0 = [0 0 0];
    pF = [5 10 3];
    Amax = 25;
    t = T/N;
    p = optimvar("p",N,3);
    v = optimvar("v",N,3);
    a = optimvar("a",N-1,3,"LowerBound",-Amax,"UpperBound",Amax);
    trajectoryprob = optimproblem;
    s = optimvar("s",N-1,"LowerBound",0,"UpperBound",3*Amax);
    trajectoryprob.Objective = sum(s)*t;
    scons = optimconstr(N-1);
    for i = 1:(N-1)
        scons(i) = norm(a(i,:)) <= s(i);
    end
    acons = optimconstr(N-1);
    for i = 1:(N-1)
        acons(i) = norm(a(i,:)) <= Amax;
    end
    vcons = optimconstr(N+1,3);
    vcons(1,:) = v(1,:) == [0 0 0];
    if air
        vcons(2:N,:) = v(2:N,:) == v(1:(N-1),:)*exp(-t) + t*(a(1:(N-1),:) + repmat(g,N-1,1));
    else
        vcons(2:N,:) = v(2:N,:) == v(1:(N-1),:) + t*(a(1:(N-1),:) + repmat(g,N-1,1));
    end
    vcons(N+1,:) = v(N,:) == [0 0 0];
    pcons = optimconstr(N+1,3);
    pcons(1,:) = p(1,:) == p0;
    if air
        pcons(2:N,:) = p(2:N,:) == p(1:(N-1),:) + t*(1+exp(-t))/2*v(1:(N-1),:) + t^2/2*(a(1:(N-1),:) + repmat(g,N-1,1));
    else
        pcons(2:N,:) = p(2:N,:) == p(1:(N-1),:) + t*v(1:(N-1),:) + t^2/2*(a(1:(N-1),:) + repmat(g,N-1,1));
    end
    pcons((N+1),:) = p(N,:) == pF;
    trajectoryprob.Constraints.acons = acons;
    trajectoryprob.Constraints.scons = scons;
    trajectoryprob.Constraints.vcons = vcons;
    trajectoryprob.Constraints.pcons = pcons;
end

The iterative display looks a lot like yours. This problem is definitely infeasible, it is not a question of twiddling some options to make it feasible. I know that this problem is infeasible because setupproblem(1) specifies too small a time for the problem to be solved with the given limits on applied acceleration (see the example for details). Perhaps something like that is going on with your example.

Alan Weiss

MATLAB mathematical toolbox documentation

[Torsten]

No chance to get "solm" after the 6th iteration ?

[Nick]

Ok yes I had a false understanding of what the iterative display is showing. Thank you for this example. As @Torsten mentioned, do you know if we should be getting any output at iterations that are not converged? Like if we were to set the maximum iterations to something smaller?