Possible bug in coneprog.m function

I am running a trajectory optimisation problem.

The formulation is convexified via slack variables and linearisation and I was testing coneprog as solver.

I noticed some issues with the function that I would like to report and verify if it is my mistake or a bug.

My code is completely deterministic. There is no randomisation and the execution should be 100% reproducible.

Nevertheless, running the same code multiple times with the same input and options leads to different results.

Unlinke fmincon, with coneprog it is impossible to access internal calculation and to understand what is going on, therefore I do not have many insights.

My guess is that a random noise is applied to the initial guess.

In theory this should not be a problem since a convex problem has only one minimum.

In practice, the solver sometimes finds an optimal solution and other times declares the problem unfeasible.

I can solve the same exact problem with an NLP formulation using fmincon, so I know for sure the proble is solvable.

I tried to change the initial conditions and the problem persists.

Sometimes, different minima are found in different calls of the solver.

This is not due to the presence of a flat in the solution space because the fval differ by 1 to 2 orders of magnitude.

Another issue I am having is connected to the iterative procedure I have implemented.

At first I obtain an anlytical initial guess that is numerically propagated.

The constraints are linearised around this reference trajectory and the coneprog is called to solve the problem.

Once the optimal control is obtained, the solution is again propagated and if the error is not below the specified tolerance the process is repeated iteratively.

For some initial conditions the solver works smoothly, but for other it stops at the 5-th or 6-th iteration (declaring the problem unfeasible) when it is really close to the finalise the solution.

Again, it is important to highlight that I am 100% confident that the problem is feasible.

It is important to highlight that some degree of reproducibility is still present.

Even if the solver behaves differently, the number of possible behaviours is limited to just a few outcomes, depending on the initial conditions.

In general, for every problem there are 2 or 3 possible outcomes.

Moreover, some outcomes may occur more frequently than others.

I have spent days trying to find the issue in my code, but now, after testing with fmincon, I am starting to be confident that the problem is in coneprog.

I have also noticed that the solver behaves differently depending on the optimalityTolerance.

This was later confirmed when I found out that the problem has already been reported on this forum here.

An other bug I noticed in some solutions is that both Primal and Dual Feasibility report a value, but if I check manually the constraint satisfaction, the feasibility is off by several order of magnitudes.

I am 100% confident that the check I performed is correct because it was performed with a breakpoint placed right after the call to coneprog as shown below.

I have also stored the inputs of coneprog right before the function was called.

In my opinion, if the problem is coneprog as I think it is, the bug is related to the update step of the current guess.

Have you experienced similar problem with coneprog, or are you aware of some limitations it may have?

Did you already encountered such problem and understood where it is coming from?

The main motivations why I think coneprog is the problem are the different outputs produced from a deterministic code, and the fact that, when the problem converges, the solution is correct.

As a final comment, I think it would be really useful to have the freedom of specifying an initial guess.

If I know already a solution that is both an interior point and very close to the optimal one, this would cut down considerably the computational time.

----------------------------------------------------------------------------------------------------------------------------

Below you can find an example of the solver converging to an unfeasible point after some iterations.

The target value of poc is 1e-10.

The formulation of the optimsiation problem is based on a direct transcription to discretise a continuous trajectory into discrete nodes.

The solution vector x is built as:

x = [ux1 uy1 uz1 u1 ... uxN uyN uzN uN b1 b2]

where:

• [uxi uyi uzi] is the control vector at the i-th node.
• ui is the magnitude of the control vector at the i-th node.
• [b1 b2] are some end conditions that must satisfy a linear constraint

The only cone constraints present are the one that impose:

• sqrt(uxi^2 + uyi^2 + uzi^2) <= 1
• sqrt(uxi^2 + uyi^2 + uzi^2) <= ui

The linear objective function minimises the control effort defined as:

• 

In the example below, the trajectory is discretised into 852 nodes.

This is clear from the initial Fval of every iteration that corresponds to the number of nodes, as the initial guess is probably all ones.

The last call of coneprog declares the problem unfeasible and it is the only iteration where Fval increases after the initial guess.

The problem with the constraint satisfication can be understood from the values of the Infeasibilities.

To help the problem converges I defined a lower bound to uxi, uyi, and uzi of -1, and a lower bound to ui of 0.

Similarly, I defined an upper bound to uxi, uyi, uzi, and ui of 1.

An Fval of 1.118894e+03 in the last line of the example below, clearly goes against the upper bound constraint, as the maximum value of Fval allowed by the upper bound is 852.

Instead the sovler is analysing a solution with Fval of approximately 1119, but it thinks that the constraints are satisfied, which is impossible.

Maybe my interpretation of Primal and Dual Feasibility is wrong, but my intuition is suggesting me that coneprog may have some problems.

Iter           Fval  Primal Infeas    Dual Infeas    Duality Gap    Time
   0   8.520000e+02   3.351509e-01   2.500000e-01   8.338790e-02    0.43
   1   2.903535e+02   5.653968e-02   4.217479e-02   1.406747e-02    0.46
   2   5.760140e+01   1.155766e-02   8.621236e-03   2.875627e-03    0.50
   3   1.826360e+01   3.567321e-03   2.660981e-03   8.875744e-04    0.54
   4   1.587390e+01   2.115580e-03   1.578080e-03   5.263711e-04    0.58
   5   1.496351e+01   1.979895e-03   1.476868e-03   4.926117e-04    0.61
   6   1.477297e+01   1.930240e-03   1.439829e-03   4.802573e-04    0.65
   7   1.407791e+01   1.740965e-03   1.298643e-03   4.331644e-04    0.69
   8   1.195150e+01   1.139652e-03   8.501032e-04   2.835533e-04    0.73
   9   1.118074e+01   9.249616e-04   6.899590e-04   2.301369e-04    0.76
  10   9.263514e+00   3.932808e-04   2.933610e-04   9.785102e-05    0.80
  11   8.435984e+00   1.662211e-04   1.239898e-04   4.135699e-05    0.84
  12   8.137082e+00   8.259246e-05   6.160841e-05   2.054958e-05    0.89
  13   7.834045e+00   4.811473e-06   3.589036e-06   1.197128e-06    0.93
  14   7.817072e+00   4.211285e-07   3.141435e-07   1.047797e-07    0.97
  15   7.815469e+00   1.297491e-08   9.741467e-09   3.228250e-09    1.02
  16   7.815429e+00   2.773558e-09   3.724366e-09   6.900807e-10    1.05
Optimal solution found.

poc =

   4.7669e-05

Iter           Fval  Primal Infeas    Dual Infeas    Duality Gap    Time
   0   8.520000e+02   4.148309e-01   2.500000e-01   8.339003e-02    0.24
   1   2.936365e+02   7.220152e-02   4.351262e-02   1.451408e-02    0.28
   2   3.227398e+01   9.133979e-03   5.504640e-03   1.836128e-03    0.32
   3   2.999204e+01   5.626122e-03   3.390612e-03   1.130973e-03    0.36
   4   2.891739e+01   5.422358e-03   3.267812e-03   1.090012e-03    0.40
   5   2.734221e+01   4.962369e-03   2.990597e-03   9.975440e-04    0.44
   6   2.355320e+01   3.776002e-03   2.275627e-03   7.590585e-04    0.48
   7   1.765904e+01   1.858229e-03   1.119871e-03   3.735444e-04    0.53
   8   1.493406e+01   9.914144e-04   5.974810e-04   1.992958e-04    0.57
   9   1.321215e+01   4.284890e-04   2.582311e-04   8.613560e-05    0.61
  10   1.247209e+01   1.862607e-04   1.122510e-04   3.744245e-05    0.65
  11   1.203916e+01   4.726695e-05   2.848568e-05   9.501685e-06    0.68
  12   1.190390e+01   4.898257e-06   2.951962e-06   9.846561e-07    0.73
  13   1.188906e+01   1.900382e-07   1.145391e-07   3.820181e-08    0.77
  14   1.188847e+01   9.042867e-09   5.942560e-09   1.817813e-09    0.81
Optimal solution found.

poc =

   1.7472e-08

Iter           Fval  Primal Infeas    Dual Infeas    Duality Gap    Time
   0   8.520000e+02   3.740554e-01   2.500000e-01   8.338857e-02    0.23
   1   2.549188e+02   4.172979e-02   2.789011e-02   9.302867e-03    0.26
   2   9.392526e+01   1.561744e-02   1.043792e-02   3.481612e-03    0.30
   3   2.667086e+01   4.410425e-03   2.947708e-03   9.832206e-04    0.34
   4   2.407913e+01   2.549053e-03   1.703660e-03   5.682632e-04    0.39
   5   2.355160e+01   2.499060e-03   1.670247e-03   5.571180e-04    0.43
   6   2.324912e+01   2.449195e-03   1.636919e-03   5.460015e-04    0.47
   7   2.249923e+01   2.308274e-03   1.542735e-03   5.145858e-04    0.50
   8   2.040632e+01   1.851989e-03   1.237777e-03   4.128657e-04    0.54
   9   1.970151e+01   1.692941e-03   1.131477e-03   3.774090e-04    0.59
  10   1.744126e+01   1.126199e-03   7.526951e-04   2.510647e-04    0.62
  11   1.620420e+01   8.023922e-04   5.362790e-04   1.788781e-04    0.66
  12   1.474391e+01   3.945190e-04   2.636768e-04   8.795053e-05    0.70
  13   1.427557e+01   2.644519e-04   1.767465e-04   5.895455e-05    0.73
  14   1.370019e+01   1.011430e-04   6.759892e-05   2.254791e-05    0.77
  15   1.348178e+01   3.984925e-05   2.663325e-05   8.883635e-06    0.81
  16   1.336507e+01   7.537767e-06   5.037870e-06   1.680402e-06    0.85
  17   1.334002e+01   6.477134e-07   4.329037e-07   1.443954e-07    0.88
  18   1.333772e+01   1.179424e-08   7.999240e-09   2.629302e-09    0.93
  19   1.333768e+01   1.381986e-10   8.290199e-09   3.080877e-11    0.99
Optimal solution found.

poc =

   5.3035e-10

Iter           Fval  Primal Infeas    Dual Infeas    Duality Gap    Time
   0   8.520000e+02   3.472650e-01   2.500000e-01   8.338749e-02    0.23
   1   2.309201e+02   2.520373e-02   1.814445e-02   6.052081e-03    0.27
   2   9.405942e+01   1.047608e-02   7.541850e-03   2.515584e-03    0.31
   3   2.721726e+01   3.116075e-03   2.243298e-03   7.482520e-04    0.35
   4   2.416710e+01   1.805118e-03   1.299525e-03   4.334564e-04    0.40
   5   2.375956e+01   1.781739e-03   1.282694e-03   4.278426e-04    0.43
   6   2.358202e+01   1.764206e-03   1.270072e-03   4.236324e-04    0.47
   7   2.304355e+01   1.703461e-03   1.226341e-03   4.090461e-04    0.52
   8   2.181656e+01   1.527846e-03   1.099914e-03   3.668761e-04    0.56
   9   2.124914e+01   1.448597e-03   1.042861e-03   3.478463e-04    0.60
  10   1.948615e+01   1.170737e-03   8.428266e-04   2.811248e-04    0.64
  11   1.783054e+01   8.717316e-04   6.275695e-04   2.093258e-04    0.68
  12   1.708765e+01   7.336925e-04   5.281936e-04   1.761789e-04    0.72
  13   1.550449e+01   3.894897e-04   2.803981e-04   9.352676e-05    0.76
  14   1.514484e+01   3.155406e-04   2.271613e-04   7.576962e-05    0.80
  15   1.441642e+01   1.496860e-04   1.077607e-04   3.594357e-05    0.85
  16   1.411738e+01   7.908921e-05   5.693722e-05   1.899141e-05    0.90
  17   1.383348e+01   1.063875e-05   7.658959e-06   2.554645e-06    0.95
  18   1.380209e+01   2.769706e-06   1.993947e-06   6.650794e-07    0.99
  19   1.379178e+01   1.701668e-07   1.225354e-07   4.086155e-08    1.03
  20   1.379114e+01   8.779631e-09   6.671732e-09   2.108221e-09    1.07
Optimal solution found.

poc =

   1.8216e-10

Iter           Fval  Primal Infeas    Dual Infeas    Duality Gap    Time
   0   8.520000e+02   3.384156e-01   2.500000e-01   8.338731e-02    0.24
   1   2.283493e+02   2.314639e-02   1.709909e-02   5.703387e-03    0.29
   2   9.475565e+01   9.793019e-03   7.234463e-03   2.413050e-03    0.33
   3   2.755188e+01   2.925713e-03   2.161331e-03   7.209104e-04    0.37
   4   2.438422e+01   1.697792e-03   1.254221e-03   4.183445e-04    0.41
   5   2.398865e+01   1.677057e-03   1.238903e-03   4.132353e-04    0.45
   6   2.383017e+01   1.662889e-03   1.228437e-03   4.097444e-04    0.49
   7   2.333142e+01   1.612151e-03   1.190955e-03   3.972422e-04    0.53
   8   2.217540e+01   1.460099e-03   1.078629e-03   3.597758e-04    0.56
   9   2.173362e+01   1.404707e-03   1.037709e-03   3.461270e-04    0.61
  10   2.035602e+01   1.210646e-03   8.943487e-04   2.983093e-04    0.65
  11   1.865148e+01   9.312689e-04   6.879625e-04   2.294694e-04    0.68
  12   1.727485e+01   7.016610e-04   5.183427e-04   1.728928e-04    0.72
  13   1.622825e+01   5.067634e-04   3.743647e-04   1.248691e-04    0.76
  14   1.539423e+01   3.289696e-04   2.430219e-04   8.105977e-05    0.80
  15   1.498430e+01   2.416941e-04   1.785483e-04   5.955466e-05    0.83
  16   1.446995e+01   1.222391e-04   9.030248e-05   3.012032e-05    0.87
  17   1.418069e+01   5.599477e-05   4.136539e-05   1.379739e-05    0.91
  18   1.397083e+01   5.693174e-06   4.205758e-06   1.402827e-06    0.95
  19   1.394958e+01   5.174181e-07   3.822538e-07   1.274944e-07    0.99
  20   1.394766e+01   4.663196e-08   3.460641e-08   1.149035e-08    1.03
  21   1.394748e+01   2.032849e-09   2.212183e-09   5.009042e-10    1.07
Optimal solution found.

poc =

   1.0958e-10

Iter           Fval  Primal Infeas    Dual Infeas    Duality Gap    Time
   0   8.520000e+02   3.341137e-01   2.500000e-01   8.338698e-02    0.24
   1   9.059917e+02   1.639868e-01   1.227028e-01   4.092728e-02    0.28
   2   1.017406e+03   7.906474e-02   5.916006e-02   1.973271e-02    0.31
   3   1.124331e+03   4.121461e-03   3.083876e-03   1.028620e-03    0.35
   4   1.118894e+03   6.591392e-06   4.931998e-06   1.645058e-06    0.39
   5   1.118894e+03   9.882157e-09   7.394312e-09   2.466356e-09    0.43
   6   1.118894e+03   1.481584e-11   1.109217e-11   3.698267e-12    0.46
Problem is infeasible.

This was tested on R2025a and R2023a.

[Cris LaPierre]

Please report potential bugs here: https://supportcases.mathworks.com/mwsupport/s/casetypeselection?language=en_US&c__caseParameter=productusage&s_tid=srlp_product

[chicken vector]

Reported.

[Torsten]

What eactly do you think is a bug in "coneprog" ? From what you report I would infer that "fmincon" is better suited to solve your problem than "coneprog". But from my understanding giving bad results is not a bug . What are the exitflags you get for the several cases ?

[chicken vector]

You are correct highlighting the difference, but I do not agree in saying that fmincon is better suited.

If the problem is convex and a global minimum exists, coneprog MUST find it.

This is an intrinsic property of a convex solver.

If it does not find it, then either the optimisation problem is not well poned, or the solver might have some limitations.

What makes me call these limitations a bug is the fact that, in this case, a deterministic code does not always produce the same output.

But I agree that the difference is subtle in this case and bug may sound too strong.

exitflag is either 1 for convergence or -2 for unefasible.

[chicken vector]

Also I forgot to mention that Feasibility and a-posteriori check of the constraints do not agree.

If that is the case, that would definetly be a bug.

[Torsten]

Also I forgot to mention that Feasibility and a-posteriori check of the constraints do not agree.

If that is the case, that would definetly be a bug.

Thus you got exitflag = 1 although the constraints were not satisfied ?

[chicken vector]

In this case each one of the SOC constraint was not satisfied by a large margin.

[Torsten]

In this case each one of the SOC constraint was not satisfied by a large margin.

That's surprising if your inputs for the coneprog constraints in "socConstraint" were really as intended.

[chicken vector]

I placed the break point right after the call to coneprog so I am sure I am using the same matrices.

I also tested using different methods to be sure:

1. I saved the results right after coneprog execution
2. I run check function inside the code right after coneprog execution to be sure the break-point was not the issue
3. I executed the same check inside the coneprog function file both before and after the interior point function shown below.

[x, fval, exitflag, output, lambda] = optim.coneprog.interiorPointMethod( ...

f,socConstraints,Aineq,bineq,Aeq,beq,lb,ub,options,coneprogStartTime);

[chicken vector]

At this point I managed to improve how the constraints are imposed.

The solver behaves somewhat better but the problem reported are still there, just less frequent.

I have a case that, when the solver behaves correctly, converges in 3 iterations.

Sometimes it does, other times the problem mentioned above (unsatisfied constraint with null feasibility) occur.

Most of the times occurs at the third iteration, other times at the second.

It seems like there are 3 possible outcomes to this optimisation, and when the error occurs, the solution is always the same.

I am starting to think that the source of error rests in a part of coneprog that use some randomisation.

Maybe it is used to add a perturbation to the initial guess.

Nevertheless, this perturbation, if present, is causing the problem, probably affecting the precision required for the initial guess to be an interior-point.

I think that the problem could easily be solved by allowing the user to specify an initial guess.

[chicken vector]

Here the example I was mentioning above.

The variable n is equal to 1704, meaning that every cone constraint is not satisfied.

Running this code multiple times result in different outcomes.

sometimes the third iteration converges, other times this error occurs at the 3rd iteration, more rarely at 2nd.

Iter           Fval  Primal Infeas    Dual Infeas    Duality Gap    Time
   0   0.000000e+00   1.974758e-02   3.190621e-01   5.861005e-04    0.29
   1   4.667885e-01   2.187830e-03   3.534882e-02   6.493395e-05    0.35
   2   5.452636e+00   6.393111e-04   1.032936e-02   1.897451e-05    0.45
   3   8.599331e+00   3.316953e-04   5.359208e-03   9.844587e-06    0.54
   4   8.327029e+00   2.916956e-04   4.712933e-03   8.657413e-06    0.63
   5   8.422668e+00   2.187673e-04   3.534629e-03   6.492930e-06    0.72
   6   8.691548e+00   1.513222e-04   2.444916e-03   4.491183e-06    0.80
   7   8.901929e+00   5.648062e-05   9.125587e-04   1.676323e-06    0.89
   8   8.978177e+00   2.166699e-05   3.500740e-04   6.430678e-07    0.97
   9   9.005291e+00   7.488899e-06   1.209983e-04   2.222676e-07    1.06
  10   9.017622e+00   6.383928e-07   1.031453e-05   1.894725e-08    1.14
  11   9.018050e+00   4.782065e-07   7.726395e-06   1.419299e-08    1.23
  12   9.018637e+00   1.818530e-09   2.939137e-08   5.397327e-11    1.32
  13   9.018638e+00   6.939495e-11   2.811608e-09   8.092324e-14    1.40
Optimal solution found.
Iter           Fval  Primal Infeas    Dual Infeas    Duality Gap    Time
   0   0.000000e+00   1.981009e-02   3.190444e-01   5.860946e-04    0.25
   1   4.761665e-01   2.276325e-03   3.666055e-02   6.734659e-05    0.31
   2   5.438280e+00   6.481829e-04   1.043908e-02   1.917692e-05    0.41
   3   8.687641e+00   3.367964e-04   5.424155e-03   9.964343e-06    0.50
   4   8.419633e+00   2.968237e-04   4.780388e-03   8.781724e-06    0.60
   5   8.541256e+00   2.260515e-04   3.640592e-03   6.687883e-06    0.69
   6   8.787613e+00   1.419160e-04   2.285577e-03   4.198678e-06    0.82
   7   9.027472e+00   4.371756e-05   7.040778e-04   1.293413e-06    0.92
   8   9.100130e+00   1.042201e-05   1.678480e-04   3.083420e-07    1.01
   9   9.112941e+00   2.079457e-06   3.348996e-05   6.152212e-08    1.11
  10   9.114953e+00   1.296767e-06   2.088462e-05   3.836571e-08    1.20
  11   9.116343e+00   7.275913e-08   1.171798e-06   2.152627e-09    1.29
  12   9.116374e+00   1.096068e-10   1.774742e-09   3.242816e-12    1.39
Optimal solution found.
Iter           Fval  Primal Infeas    Dual Infeas    Duality Gap    Time
   0   0.000000e+00   1.980905e-02   3.191624e-01   2.932378e-04    0.24
   1   5.846916e+01   1.793433e-02   2.889571e-01   2.654860e-04    0.30
   2   1.617623e+02   1.210288e-02   1.950010e-01   1.791617e-04    0.39
   3   1.417412e+02   7.042323e-03   1.134656e-01   1.042491e-04    0.48
   4   1.525567e+02   5.522407e-03   8.897676e-02   8.174943e-05    0.57
   5   1.714116e+02   4.637908e-03   7.472575e-02   6.865600e-05    0.65
   6   2.058631e+02   6.542215e-04   1.054079e-02   9.684588e-06    0.74
   7   2.056674e+02   1.059580e-04   1.707190e-03   1.568519e-06    0.83
   8   2.058230e+02   2.484433e-05   4.002907e-04   3.677762e-07    0.91
   9   2.058922e+02   2.545499e-06   4.101293e-05   3.768157e-08    1.01
  10   2.059006e+02   4.204285e-07   6.773894e-06   6.223670e-09    1.10
  11   2.059022e+02   2.200719e-08   3.545490e-07   3.257439e-10    1.18
  12   2.059023e+02   4.168059e-11   7.147533e-10   4.883665e-13    1.27
Optimal solution found.
K>> c = socConstraint;
n = numel(c);
idx = false(1,n);
for k = 1 : n
  idx(k) = norm(c(k).A*x-c(k).b) > c(k).d'*x-c(k).gamma;
end

sum(idx)


ans =

        1704

[Torsten]

What do you get if you compute the "mean constraint violation" ?

summe = 0;
count = 0;
for k = 1:n
    value = norm(c(k).A*x-c(k).b) -(c(k).d'*x-c(k).gamma );
    if value > 0
      summe = summe + value^2;
      count = count + 1;      
    end
end
mean_constraint_violation = sqrt(summe)/count

[chicken vector]

Below the cone constraint violation and the respective coneprog iter output.

Note that I removed the bounds constraint. I thought they were helping instead it behaves better without them.

They were already redundant with the cone constraint, so not strictly necessary.

The peculiarity is that the violations are mostly integers (3 and 4).

I do not have any explanation for that, but the control should be bounded between -1 and 1, instead, for this run, uyi = 3 for most nodes.

Iter           Fval  Primal Infeas    Dual Infeas    Duality Gap    Time
   0   0.000000e+00   1.980905e-02   3.191624e-01   2.932378e-04    0.26
   1   3.738905e+01   1.820873e-02   2.933783e-01   2.695480e-04    0.32
   2   5.919189e+01   6.927068e-03   1.116086e-01   1.025430e-04    0.41
   3   9.100675e+01   1.814142e-03   2.922937e-02   2.685516e-05    0.50
   4   9.579543e+01   4.494799e-04   7.241998e-03   6.653752e-06    0.59
   5   9.670716e+01   1.039742e-04   1.675227e-03   1.539153e-06    0.68
   6   9.691739e+01   2.598008e-05   4.185898e-04   3.845890e-07    0.77
   7   9.696804e+01   5.054429e-06   8.143671e-05   7.482185e-08    0.86
   8   9.697853e+01   1.730429e-07   2.788060e-06   2.561593e-09    0.95
   9   9.697884e+01   6.657302e-09   1.071695e-07   9.843973e-11    1.04
  10   9.697885e+01   5.081288e-11   7.906141e-10   1.477319e-13    1.13
Optimal solution found.

[Torsten]

I do not have any explanation for that, but the control should be bounded between -1 and 1, instead, for this run, uyi = 3 for most nodes.

And you defined -1 and 1 as bound constraints for the respective variables in lb and ub ?

I think without code that demonstrates this behaviour we are not able to help.

[chicken vector]

Not in this run, but in previous ones yes. However the cone constraint impose a stricter feasible domain so they are not necessary.

I understand the code is required but I cannot legally share it.

I tried to ask the forum first, becasue maybe someone else already noticed this behaviour.

If necessary, I can share my screen in case someone from the staff wants to dig deeper.

Nevertheless, thanks for your interest and time.

[Torsten]

Not in this run, but in previous ones yes. However the cone constraint impose a stricter feasible domain so they are not necessary.

Note that the constraints

• sqrt(uxi^2 + uyi^2 + uzi^2) <= 1
• sqrt(uxi^2 + uyi^2 + uzi^2) <= ui

don't ensure that ui <= 1 (if it is that what you mean).

[chicken vector]

It is true, but the objective function does.

Regardless, the first constraint is not satisfied if uyi = 3.

[Torsten]

It is true, but the objective function does.

Not necessarily. ui could appear in other linear equality/inequality constraints (Aeq, beq, A, b) that would make a value ui <= 1 impossible.

Is it correct that in reality, you want to implement the quadratic equality constraint

ui^2 = uxi^2 + uyi^2 + uzi^2

ui >= 0

?

[chicken vector]

In general the following constraint

sqrt(uxi^2 + uyi^2 + uzi^2) <= 1

imposes that every component is, for sure, smaller than 1.

Then, by imposing

sqrt(uxi^2 + uyi^2 + uzi^2) <= ui

while minimising ui, should cause one of the two:

• the problem is feasible and a solution is found where uxi, uyi, uzi and ui are always smaller than 1
• the problem is unfeasible

Instead what is happening in the example above is

uyi > 1

which is clearly against the constraint

sqrt(uxi^2 + uyi^2 + uzi^2) <= 1

but coneprog still thinks the solution is feasible.

The only issue that comes to my head is the LinearSolver option that is set to 'auto'.

Maybe the selection process of coneprog is not fully deterministic.

I will run some checks when I will find the time, today is a busy day unfortunately.

[Torsten]

In general the following constraint

sqrt(uxi^2 + uyi^2 + uzi^2) <= 1

imposes that every component is, for sure, smaller than 1.

Then, by imposing

sqrt(uxi^2 + uyi^2 + uzi^2) <= ui

while minimising ui, should cause one of the two:

• the problem is feasible and a solution is found where uxi, uyi, uzi and ui are always smaller than 1
• the problem is unfeasible

I agree that uxi, uyi and uzi should be between -1 and 1. If these are the only constraints, I agree that ui should be equal to sqrt(uxi^2 + uyi^2 + uzi^2), thus <= 1. If there are additional constraints on ui formulated via A, b, Aeq, beq, lb, ub, then ui could become > 1.

[chicken vector]

The bug is probably in the compiling step.

If you run the code below multiple times (try at least 10 times) you should see different results:

load data.mat
x = coneprog(f, socConstraint, A, b, Aeq, beq, [], [], opts);

Iter           Fval  Primal Infeas    Dual Infeas    Duality Gap    Time
   0   0.000000e+00            nan   3.316314e-01            nan    0.43
   1   0.000000e+00            nan   3.316314e-01            nan    0.58
Search direction is small and current iterate is not within the specified constraint and/or optimality tolerance.

If instead you run this code once you will always see the same result:

load data.mat
n = 3410;
N = 10;
x = NaN(n,N);
for j = 1 : N
    x(:,j) = coneprog(f, socConstraint, A, b, Aeq, beq, [], [], opts);
end

Iter           Fval  Primal Infeas    Dual Infeas    Duality Gap    Time
   0   0.000000e+00            nan   3.316314e-01            nan    0.33
   1   0.000000e+00            nan   3.316314e-01            nan    0.43
Search direction is small and current iterate is not within the specified constraint and/or optimality tolerance.
Iter           Fval  Primal Infeas    Dual Infeas    Duality Gap    Time
   0   0.000000e+00            nan   3.316314e-01            nan    0.35
   1   0.000000e+00            nan   3.316314e-01            nan    0.45
Search direction is small and current iterate is not within the specified constraint and/or optimality tolerance.
Iter           Fval  Primal Infeas    Dual Infeas    Duality Gap    Time
   0   0.000000e+00            nan   3.316314e-01            nan    0.33
   1   0.000000e+00            nan   3.316314e-01            nan    0.43
Search direction is small and current iterate is not within the specified constraint and/or optimality tolerance.
Iter           Fval  Primal Infeas    Dual Infeas    Duality Gap    Time
   0   0.000000e+00            nan   3.316314e-01            nan    0.35
   1   0.000000e+00            nan   3.316314e-01            nan    0.45
Search direction is small and current iterate is not within the specified constraint and/or optimality tolerance.
Iter           Fval  Primal Infeas    Dual Infeas    Duality Gap    Time
   0   0.000000e+00            nan   3.316314e-01            nan    0.33
   1   0.000000e+00            nan   3.316314e-01            nan    0.43
Search direction is small and current iterate is not within the specified constraint and/or optimality tolerance.
Iter           Fval  Primal Infeas    Dual Infeas    Duality Gap    Time
   0   0.000000e+00            nan   3.316314e-01            nan    0.33
   1   0.000000e+00            nan   3.316314e-01            nan    0.43
Search direction is small and current iterate is not within the specified constraint and/or optimality tolerance.
Iter           Fval  Primal Infeas    Dual Infeas    Duality Gap    Time
   0   0.000000e+00            nan   3.316314e-01            nan    0.32
   1   0.000000e+00            nan   3.316314e-01            nan    0.43
Search direction is small and current iterate is not within the specified constraint and/or optimality tolerance.
Iter           Fval  Primal Infeas    Dual Infeas    Duality Gap    Time
   0   0.000000e+00            nan   3.316314e-01            nan    0.32
   1   0.000000e+00            nan   3.316314e-01            nan    0.42
Search direction is small and current iterate is not within the specified constraint and/or optimality tolerance.
Iter           Fval  Primal Infeas    Dual Infeas    Duality Gap    Time
   0   0.000000e+00            nan   3.316314e-01            nan    0.33
   1   0.000000e+00            nan   3.316314e-01            nan    0.44
Search direction is small and current iterate is not within the specified constraint and/or optimality tolerance.
Iter           Fval  Primal Infeas    Dual Infeas    Duality Gap    Time
   0   0.000000e+00            nan   3.316314e-01            nan    0.35
   1   0.000000e+00            nan   3.316314e-01            nan    0.46
Search direction is small and current iterate is not within the specified constraint and/or optimality tolerance.

Also, another important detail is that when I run the optimisation uploading data form a .mat file i obtain completely difference results from when I run them inside my function.

Moreover, when using the .mat file, convergence is never achieved, while inside my function it works.

Another difference is that sometimes the solve stops after two iterations as the Feasibility is a NaN.

This never occurred in the original function, even if the same data was used.

[Torsten]

I get nan values repeatedly. Which MATLAB version do you use ?

load data.mat
[x,~,exitflag,~] = coneprog(f, socConstraint, A, b, Aeq, beq, [], [], opts);

Iter           Fval  Primal Infeas    Dual Infeas    Duality Gap    Time
   0   0.000000e+00            nan   3.316314e-01            nan    0.38
   1   0.000000e+00            nan   3.316314e-01            nan    0.55
Search direction is small and current iterate is not within the specified constraint and/or optimality tolerance.

exitflag

exitflag = -7

[chicken vector]

I am using R2025a, but I previously tested it also in R2023a.

This is what happens when I run it multiple times.

I saved the data with the breakpoint inside my code, but I never get the NaNs inside my original code. only when I run them as a script from the .mat file.

Do you get the NaNs 10 out of 10 times?

If yes, are you using Windows 11?

tbd
Iter           Fval  Primal Infeas    Dual Infeas    Duality Gap    Time
   0   0.000000e+00            nan   3.316314e-01            nan    0.47
   1   0.000000e+00            nan   3.316314e-01            nan    0.57
Search direction is small and current iterate is not within the specified constraint and/or optimality tolerance.
>> tbd
Iter           Fval  Primal Infeas    Dual Infeas    Duality Gap    Time
   0   0.000000e+00            nan   3.316314e-01            nan    0.33
   1   0.000000e+00            nan   3.316314e-01            nan    0.45
Search direction is small and current iterate is not within the specified constraint and/or optimality tolerance.
>> tbd
Iter           Fval  Primal Infeas    Dual Infeas    Duality Gap    Time
   0   0.000000e+00            nan   3.316314e-01            nan    0.30
   1   0.000000e+00            nan   3.316314e-01            nan    0.44
Search direction is small and current iterate is not within the specified constraint and/or optimality tolerance.
>> tbd
Iter           Fval  Primal Infeas    Dual Infeas    Duality Gap    Time
   0   0.000000e+00   1.315934e-01   3.316314e-01   4.054770e-04    0.27
   1   1.045507e+02   1.127653e-01   2.841824e-01   3.474146e-04    0.33
   2   1.984316e+05   5.411436e-04   1.363748e-03   6.828660e-06    0.43
   3   5.334492e+06   2.020939e-05   5.092812e-05   2.547574e-07    0.51
   4   6.046077e+06   1.793287e-05   4.492304e-05   2.247182e-07    0.60
   5   7.813494e+06   1.387049e-05   3.411377e-05   1.706471e-07    0.69
   6   9.797389e+06   1.141189e-05   2.437971e-05   1.219545e-07    0.77
   7   9.119914e+06   1.200243e-05   2.305982e-05   1.153520e-07    0.85
   8   9.246453e+06   1.174153e-05   2.264574e-05   1.132806e-07    0.92
   9   9.578867e+06   1.150282e-05   8.767908e-06   4.385965e-08    1.01
  10   9.955200e+06   1.130209e-05   7.036273e-06   3.519751e-08    1.09
  11   1.103410e+07   1.011822e-05   1.539801e-06   7.702536e-09    1.17
  12   1.123519e+07   9.951630e-06   1.165619e-08   5.830776e-11    1.25
  13   1.123606e+07   9.951021e-06   1.576483e-09   7.885931e-12    1.33
  14   1.123616e+07   9.952947e-06   1.696436e-10   8.485478e-13    1.41
  15   1.123617e+07   9.952944e-06   8.131822e-12   4.065292e-14    1.49
  16   1.123617e+07   9.952275e-06   8.011148e-12   4.016886e-14    1.57
  17   1.123617e+07   9.952275e-06   8.010980e-12   4.013633e-14    1.65
  18   1.123617e+07   9.952275e-06   8.010979e-12   4.008761e-14    1.73
  19   1.123617e+07   9.952275e-06   8.010979e-12   4.008756e-14    1.82
Search direction is small and current iterate is not within the specified constraint and/or optimality tolerance.

[chicken vector]

Tested in Matlab 2023a on the same pc and I get the same result.

Tested on the laptop of my co-worker on R2025a, she never gets the NaNs, but always different solutions.

She has the same computer model as me.

[Torsten]

but I never get the NaNs inside my original code. only when I run them as a script from the .mat file.

That's what I did. I ran it here on MATLAB answers using R2025a (green arrow RUN) as a script with the data from the .mat-file. And the result is always the same (independent of the numer of times I start it).

load data.mat
[x,~,exitflag,~] = coneprog(f, socConstraint, A, b, Aeq, beq, [], [], opts);

Iter           Fval  Primal Infeas    Dual Infeas    Duality Gap    Time
   0   0.000000e+00            nan   3.316314e-01            nan    0.40
   1   0.000000e+00            nan   3.316314e-01            nan    0.55
Search direction is small and current iterate is not within the specified constraint and/or optimality tolerance.

exitflag

exitflag = -7

[x,~,exitflag,~] = coneprog(f, socConstraint, A, b, Aeq, beq, [], [], opts);

Iter           Fval  Primal Infeas    Dual Infeas    Duality Gap    Time
   0   0.000000e+00            nan   3.316314e-01            nan    0.32
   1   0.000000e+00            nan   3.316314e-01            nan    0.43
Search direction is small and current iterate is not within the specified constraint and/or optimality tolerance.

exitflag

exitflag = -7

[x,~,exitflag,~] = coneprog(f, socConstraint, A, b, Aeq, beq, [], [], opts);

Iter           Fval  Primal Infeas    Dual Infeas    Duality Gap    Time
   0   0.000000e+00            nan   3.316314e-01            nan    0.32
   1   0.000000e+00            nan   3.316314e-01            nan    0.43
Search direction is small and current iterate is not within the specified constraint and/or optimality tolerance.

exitflag

exitflag = -7

[x,~,exitflag,~] = coneprog(f, socConstraint, A, b, Aeq, beq, [], [], opts);

Iter           Fval  Primal Infeas    Dual Infeas    Duality Gap    Time
   0   0.000000e+00            nan   3.316314e-01            nan    0.32
   1   0.000000e+00            nan   3.316314e-01            nan    0.43
Search direction is small and current iterate is not within the specified constraint and/or optimality tolerance.

exitflag

exitflag = -7

[x,~,exitflag,~] = coneprog(f, socConstraint, A, b, Aeq, beq, [], [], opts);

Iter           Fval  Primal Infeas    Dual Infeas    Duality Gap    Time
   0   0.000000e+00            nan   3.316314e-01            nan    0.33
   1   0.000000e+00            nan   3.316314e-01            nan    0.44
Search direction is small and current iterate is not within the specified constraint and/or optimality tolerance.

exitflag

exitflag = -7

[x,~,exitflag,~] = coneprog(f, socConstraint, A, b, Aeq, beq, [], [], opts);

Iter           Fval  Primal Infeas    Dual Infeas    Duality Gap    Time
   0   0.000000e+00            nan   3.316314e-01            nan    0.32
   1   0.000000e+00            nan   3.316314e-01            nan    0.42
Search direction is small and current iterate is not within the specified constraint and/or optimality tolerance.

exitflag

exitflag = -7

[x,~,exitflag,~] = coneprog(f, socConstraint, A, b, Aeq, beq, [], [], opts);

Iter           Fval  Primal Infeas    Dual Infeas    Duality Gap    Time
   0   0.000000e+00            nan   3.316314e-01            nan    0.32
   1   0.000000e+00            nan   3.316314e-01            nan    0.42
Search direction is small and current iterate is not within the specified constraint and/or optimality tolerance.

exitflag

exitflag = -7

[x,~,exitflag,~] = coneprog(f, socConstraint, A, b, Aeq, beq, [], [], opts);

Iter           Fval  Primal Infeas    Dual Infeas    Duality Gap    Time
   0   0.000000e+00            nan   3.316314e-01            nan    0.36
   1   0.000000e+00            nan   3.316314e-01            nan    0.47
Search direction is small and current iterate is not within the specified constraint and/or optimality tolerance.

exitflag

exitflag = -7

[x,~,exitflag,~] = coneprog(f, socConstraint, A, b, Aeq, beq, [], [], opts);

Iter           Fval  Primal Infeas    Dual Infeas    Duality Gap    Time
   0   0.000000e+00            nan   3.316314e-01            nan    0.33
   1   0.000000e+00            nan   3.316314e-01            nan    0.43
Search direction is small and current iterate is not within the specified constraint and/or optimality tolerance.

exitflag

exitflag = -7

[x,~,exitflag,~] = coneprog(f, socConstraint, A, b, Aeq, beq, [], [], opts);

Iter           Fval  Primal Infeas    Dual Infeas    Duality Gap    Time
   0   0.000000e+00            nan   3.316314e-01            nan    0.32
   1   0.000000e+00            nan   3.316314e-01            nan    0.45
Search direction is small and current iterate is not within the specified constraint and/or optimality tolerance.

exitflag

exitflag = -7

[chicken vector]

I think the problem occurs during the compiling step and, probably, that is the reason why the bug does not show up in Matlab Answers.

To demonstrate this I first run the inputs in data.mat on Matlab Answers and, as expected, NaNs are returned.

Then I load the solution that my code finds and check its feasibility which is always below the ConstraintTolerance.

This is to show that an optimal solution exists and those NaN should not appear.

I didn't compute the OptimalityTolerance so I cannot demonstrate that is the optimal solution (this is an act of faith), but it is a feasible solution nonetheless.

It follows that coneprog should not return NaNs.

% Load coneprog input:
load data.mat;

% Run optimisation (expected NaN on Matlab Answers):
[x,~,exitflag] = coneprog(f, socConstraint, A, b, Aeq, beq, [], [], opts);

Iter           Fval  Primal Infeas    Dual Infeas    Duality Gap    Time
   0   0.000000e+00            nan   3.316165e-01            nan    0.42
   1   0.000000e+00            nan   3.316165e-01            nan    0.58
Search direction is small and current iterate is not within the specified constraint and/or optimality tolerance.

% NaN solution outputs x = zeros and exitflag = -7:
[all(x==0) exitflag]

ans = 1×2

1 -7

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

% Load solution computed inside code when it converges:
load solution.mat;

% Check linear inequality constraint satisfaction:
linineq = all(A*x - b <= opts.ConstraintTolerance);

% Check linear equality constraint satisfaction:
lineq = all(abs(Aeq*x - beq) <= opts.ConstraintTolerance);

% Check SOC constraint satisfaction:
N = numel(socConstraint);
socs = false(1,N);
for j = 1 : N
    socs(j) = norm(socConstraint(j).A*x - socConstraint(j).b) - socConstraint(j).d'*x + socConstraint(j).gamma <= opts.ConstraintTolerance;
end
soc = all(socs);

% Overall constraint satisfication:
all([linineq lineq soc])

ans = logical
   1

Answer 1

[Torsten]

This works for dense matrices/vectors, but in sparse mode, it gives exitflag = -7 as above. Maybe it's worth asking MATLAB support how this can be avoided.

n = 100;
index = 1;
bsoc = zeros(4*n+2,1);
f = zeros(4*n+2,1);
for i = 1:n
  Asoc = zeros(4*n+2);
  dsoc = zeros(4*n+2,1);
  Asoc(index,index)     = 1.0;
  Asoc(index+1,index+1) = 1.0;
  Asoc(index+2,index+2) = 1.0;
  gamma = -1.0;
  socConstraints(2*i-1) = secondordercone(Asoc,bsoc,dsoc,gamma);
  %socConstraints(2*i-1) = secondordercone(sparse(Asoc),sparse(bsoc),sparse(dsoc),gamma);
  dsoc(index+3) = 1.0;
  gamma = 0;
  socConstraints(2*i) = secondordercone(Asoc,bsoc,dsoc,gamma);
  %socConstraints(2*i) = secondordercone(sparse(Asoc),sparse(bsoc),sparse(dsoc),gamma);
  f(index+3) = 1;
  index = index + 4;
end
options = optimoptions("coneprog",Display='iter');
[x,~,exitflag,~] = coneprog(f, socConstraints,[],[],[],[],[],[],options)

Iter           Fval  Primal Infeas    Dual Infeas    Duality Gap    Time
   0   0.000000e+00   5.000000e-01   3.333333e-01   4.975124e-03    9.04
   1   0.000000e+00   7.496256e-04   4.997504e-04   7.458961e-06   14.05
   2  -4.125786e-17   1.123877e-06   7.492513e-07   1.118286e-08   18.97
   3  -6.119795e-20   1.684974e-09   1.123316e-09   1.676591e-11   24.10
Optimal solution found.

x = 402×1

1.0e-20 * 0 0 0 -0.0770 0 0 0 -0.0159 0 0 0 -0.0770 0 0 0

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

exitflag = 1

[chicken vector]

Thank you for your time @Torsten, you have been very helpful and, most of all, you most probably manage to find the origin of the problem.

The bug is indeed connected to sparsity and I'm guessing it is occuring when compiling.

If you post the answer I would be very glad to accept it.

% Load coneprog input:
load data.mat;

% Run optimisation with sparse SOC constraints:
[~,~,sparseExitflag] = coneprog(f, socConstraint, A, b, Aeq, beq, [], [], opts);

Iter           Fval  Primal Infeas    Dual Infeas    Duality Gap    Time
   0   0.000000e+00            nan   3.316165e-01            nan    0.54
   1   0.000000e+00            nan   3.316165e-01            nan    0.75
Search direction is small and current iterate is not within the specified constraint and/or optimality tolerance.

% Make SOC constraints matrices full:
N = numel(socConstraint);
for j = 1 : N
    socConstraint(j).A = full(socConstraint(j).A);
    socConstraint(j).b = full(socConstraint(j).b);
    socConstraint(j).d = full(socConstraint(j).d);
    socConstraint(j).gamma = full(socConstraint(j).gamma);
end

% Run optimisation with full SOC constraints:
[~,~,fullExitflag] = coneprog(f, socConstraint, A, b, Aeq, beq, [], [], opts);

Iter           Fval  Primal Infeas    Dual Infeas    Duality Gap    Time
   0   0.000000e+00   1.308390e-01   3.316165e-01   5.860879e-04    0.47
   1   1.056167e+00   6.315378e-02   1.600657e-01   2.828948e-04    0.54
   2   3.308565e+01   5.091152e-02   1.290373e-01   2.280561e-04    0.65
   3   6.991760e+01   2.113367e-02   5.356412e-02   9.466741e-05    0.79
   4   8.093810e+01   8.019225e-03   2.032504e-02   3.592179e-05    0.90
   5   8.723166e+01   1.832851e-03   4.645433e-03   8.210181e-06    1.01
   6   8.832996e+01   5.419241e-04   1.373528e-03   2.427527e-06    1.12
   7   8.864836e+01   1.870717e-04   4.741406e-04   8.379800e-07    1.24
   8   8.878390e+01   3.622492e-05   9.181348e-05   1.622680e-07    1.36
   9   8.881244e+01   3.349583e-06   8.489652e-06   1.500432e-08    1.47
  10   8.881521e+01   1.743829e-07   4.419871e-07   7.811406e-10    1.58
  11   8.881536e+01   3.017681e-10   8.281729e-10   1.186783e-12    1.69
Optimal solution found.

% Output:
[sparseExitflag fullExitflag]

ans = 1×2

-7 1

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

[Torsten]

@chicken vector

If you get an answer from support, it might be helpful for other users if you post it here. But I think it is already a known problem that switching to full matrices can help solving convergence issues. There are some allusions on the documentation side for coneprog:

https://uk.mathworks.com/help/optim/ug/coneprog.html

"Switching to full data instead of sparse can help the solver to complete successfully."

Maybe directly contacting the coneprog developers is an option:

https://cvxopt.org/userguide/coneprog.html

[chicken vector]

They are investigating the issue.

When I will get an answer I will post it here.

Thanks again!