Mercurial > octave

changeset 30316:c08c73d96985

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
qp.m: Workaround apparent bug in glpk which causes infeasible solution error (bug #38353) * qp.m: Check that first component of first return value from glpk call used to construct initial guess at solution is not the same as the the second return value from glpk.
author Rik <rik@octave.org>
date Sun, 21 Nov 2021 20:35:10 -0800
parents 8a80c89b021e
children 7fc817342393
files scripts/optimization/qp.m
diffstat 1 files changed, 17 insertions(+), 10 deletions(-) [+]
line wrap: on
line diff
--- a/scripts/optimization/qp.m	Thu Dec 24 22:38:46 2020 -0500
+++ b/scripts/optimization/qp.m	Sun Nov 21 20:35:10 2021 -0800
@@ -393,18 +393,25 @@
           lb = [-Inf(n-n_eq,1); zeros(n_in,1)];
           ub = [];
           ctype = repmat ("L", n_in, 1);
-          [P, dummy, status] = glpk (ctmp, Atmp, btmp, lb, ub, ctype);
-          if ((status == 0)
-              && all (abs (P(n-n_eq+1:end)) < rtol * (1 + norm (btmp))))
-            ## We found a feasible starting point
-            if (n_eq > 0)
-              x0 = xbar + Z*P(1:n-n_eq);
+          [P, FMIN, status] = glpk (ctmp, Atmp, btmp, lb, ub, ctype);
+          ## FIXME: Test based only on rtol occasionally fails (Bug #38353).
+          ## This seems to be a problem in glpk in which return value XOPT(1) is the
+          ## same as FMIN.  Workaround this by explicit test
+          if (status != 0)
+            info = 6;  # The problem is infeasible
+          else
+            if (all (abs (P(n-n_eq+2:end)) < rtol * (1 + norm (btmp)))
+                && (P(n-n_eq+1) < rtol * (1 + norm (btmp))
+                    || P(n-n_eq+1) == FMIN))
+              ## We found a feasible starting point
+              if (n_eq > 0)
+                x0 = xbar + Z*P(1:n-n_eq);
+              else
+                x0 = P(1:n);
+              endif
             else
-              x0 = P(1:n);
+              info = 6;  # The problem is infeasible
             endif
-          else
-            ## The problem is infeasible
-            info = 6;
           endif
         endif
       else