[SciPy-User] Problem using linprog

classic Classic list List threaded Threaded
4 messages Options
Reply | Threaded
Open this post in threaded view
|

[SciPy-User] Problem using linprog

Montgomery-Smith, Stephen
I am trying to solve a linear programming problem.  The constraint is of
the form A.x <= 0.  But linprog gives an answer that doesn't satisfy the
constraint.

The attached program gives A.x as

[-2.32109228  2.32017594  4.71436317  3.6433767  -4.26629574  2.32384597
 -1.96166184 -4.96206197]

which definitely doesn't satisfy the constraint.  Is this a bug, or some
subtle floating point error?

Program follows (also as attachment):

from scipy.optimize import linprog
import numpy as np

A = [[0.5919650431077654, -0.5271408402306996, 0.6096719792636803,
1.2379670854947114, 0.2656040423387233, -0.972363043155988],
[-0.5914974900295467, -0.5266568950860249, 0.6105433925177587,
1.258297461476007, -0.285688537323182, 0.9726089241528251],
[-0.593015674004932, 0.5280764198909397, 0.6078385518701857,
-1.1964319796886902, -0.2223431679788034, -0.9740888117098865],
[0.5935986604093653, 0.5285277328950352, 0.6068764832493029,
-1.1752312553140132, 0.19916734259906424, 0.976063912714949],
[0.593015674004932, -0.5280764198909397, -0.6078385518701857,
-1.1964319796886902, -0.2223431679788034, -0.9740888117098865],
[-0.5935986604093653, -0.5285277328950352, -0.6068764832493029,
-1.1752312553140132, 0.19916734259906424, 0.976063912714949],
[-0.5919650431077654, 0.5271408402306996, -0.6096719792636803,
1.2379670854947114, 0.2656040423387233, -0.972363043155988],
[0.5914974900295467, 0.5266568950860249, -0.6105433925177587,
1.258297461476007, -0.285688537323182, 0.9726089241528251]]
e = [0, 0, 0, 0, 0, -1]
bounds = [(None, None), (None, None), (None, None), (None, None), (None,
None), (0, 1)]
b = [0]*len(A)
result = linprog(e, A_ub = A, b_ub = b, bounds = bounds)
print np.matmul(A, result.x)


_______________________________________________
SciPy-User mailing list
[hidden email]
https://mail.python.org/mailman/listinfo/scipy-user

test.py (1K) Download Attachment
Reply | Threaded
Open this post in threaded view
|

Re: Problem using linprog

Chris Waigl
Hi Stephen,

The problem appears to be singular around the solution. A very quick exploration shows me that if you replace your upper bound b by a very small epsilon > 0, you get a stable result. 

For example: 
b = np.zeros(8) + 0.001
 fun: -0.11764011575264395
 message: 'Optimization terminated successfully.'
 nit: 6
 slack: array([0.        , 0.        , 0.40742577, 0.        , 0.40742577, 0.        , 0.        , 0.        , 0.88235988])
 status: 0
 success: True
 x: array([0.        , 0.        , 0.        , 0.0834722 , 0.41811509, 0.11764012])
And for print(np.dot(A, result.x)) I get [ 0.001 0.001 -0.00307426 0.001 -0.00307426 0.001 0.001 0.001 ]

In the objective function, y_2 = y_4 = -3.0742577 * epsilon, and the other 6 values also converge towards zero when epsilon -> 0 . 

If I read your problem correctly, your objective function is simply (-1) times x_5, the last element of x. The approach above would converge towards the trivial solution, x = 0, but your solution above minimizes f(x) by maximizing x_5 at 1. If we pick out an x_5, then the problem collapses to a new problems to find [x_0 ... x_4] so that A[:, 0:7] * [x_0 ... x_4]' < b, where b is (-1) * the last column of your A. But the objective function is now indeterminate, so there is nothing to optimize.  

HTH,

Chris

On Mon, Nov 26, 2018 at 11:43 AM Montgomery-Smith, Stephen <[hidden email]> wrote:
I am trying to solve a linear programming problem.  The constraint is of
the form A.x <= 0.  But linprog gives an answer that doesn't satisfy the
constraint.

The attached program gives A.x as

[-2.32109228  2.32017594  4.71436317  3.6433767  -4.26629574  2.32384597
 -1.96166184 -4.96206197]

which definitely doesn't satisfy the constraint.  Is this a bug, or some
subtle floating point error?

Program follows (also as attachment):

from scipy.optimize import linprog
import numpy as np

A = [[0.5919650431077654, -0.5271408402306996, 0.6096719792636803,
1.2379670854947114, 0.2656040423387233, -0.972363043155988],
[-0.5914974900295467, -0.5266568950860249, 0.6105433925177587,
1.258297461476007, -0.285688537323182, 0.9726089241528251],
[-0.593015674004932, 0.5280764198909397, 0.6078385518701857,
-1.1964319796886902, -0.2223431679788034, -0.9740888117098865],
[0.5935986604093653, 0.5285277328950352, 0.6068764832493029,
-1.1752312553140132, 0.19916734259906424, 0.976063912714949],
[0.593015674004932, -0.5280764198909397, -0.6078385518701857,
-1.1964319796886902, -0.2223431679788034, -0.9740888117098865],
[-0.5935986604093653, -0.5285277328950352, -0.6068764832493029,
-1.1752312553140132, 0.19916734259906424, 0.976063912714949],
[-0.5919650431077654, 0.5271408402306996, -0.6096719792636803,
1.2379670854947114, 0.2656040423387233, -0.972363043155988],
[0.5914974900295467, 0.5266568950860249, -0.6105433925177587,
1.258297461476007, -0.285688537323182, 0.9726089241528251]]
e = [0, 0, 0, 0, 0, -1]
bounds = [(None, None), (None, None), (None, None), (None, None), (None,
None), (0, 1)]
b = [0]*len(A)
result = linprog(e, A_ub = A, b_ub = b, bounds = bounds)
print np.matmul(A, result.x)

_______________________________________________
SciPy-User mailing list
[hidden email]
https://mail.python.org/mailman/listinfo/scipy-user


--

_______________________________________________
SciPy-User mailing list
[hidden email]
https://mail.python.org/mailman/listinfo/scipy-user
Reply | Threaded
Open this post in threaded view
|

Re: Problem using linprog

Montgomery-Smith, Stephen
Thank you.  That helps a little.  What is the definition of singular in
this context?

You definitely understood my problem.  The correct answer should have
been 0, as you correctly surmised.

I am trying to find a test to see if a union of half planes captures all
of Euclidean n-space.

On 11/26/18 5:40 PM, Chris F Waigl wrote:

> Hi Stephen,
>
> The problem appears to be singular around the solution. A very quick
> exploration shows me that if you replace your upper bound b by a very
> small epsilon > 0, you get a stable result. 
>
> For example: 
> b = np.zeros(8) + 0.001
>
>  fun: -0.11764011575264395
>  message: 'Optimization terminated successfully.'
>  nit: 6
>  slack: array([0.        , 0.        , 0.40742577, 0.        , 0.40742577, 0.        , 0.        , 0.        , 0.88235988])
>  status: 0
>  success: True
>  x: array([0.        , 0.        , 0.        , 0.0834722 , 0.41811509, 0.11764012])
>
> And for print(np.dot(A, result.x)) I get [ 0.001 0.001 -0.00307426 0.001
> -0.00307426 0.001 0.001 0.001 ]
>
> In the objective function, y_2 = y_4 = -3.0742577 * epsilon, and the
> other 6 values also converge towards zero when epsilon -> 0 . 
>
> If I read your problem correctly, your objective function is simply (-1)
> times x_5, the last element of x. The approach above would converge
> towards the trivial solution, x = 0, but your solution above minimizes
> f(x) by maximizing x_5 at 1. If we pick out an x_5, then the problem
> collapses to a new problems to find [x_0 ... x_4] so that A[:, 0:7] *
> [x_0 ... x_4]' < b, where b is (-1) * the last column of your A. But the
> objective function is now indeterminate, so there is nothing to optimize.  
>
> HTH,
>
> Chris
>
> On Mon, Nov 26, 2018 at 11:43 AM Montgomery-Smith, Stephen
> <[hidden email] <mailto:[hidden email]>> wrote:
>
>     I am trying to solve a linear programming problem.  The constraint is of
>     the form A.x <= 0.  But linprog gives an answer that doesn't satisfy the
>     constraint.
>
>     The attached program gives A.x as
>
>     [-2.32109228  2.32017594  4.71436317  3.6433767  -4.26629574  2.32384597
>      -1.96166184 -4.96206197]
>
>     which definitely doesn't satisfy the constraint.  Is this a bug, or some
>     subtle floating point error?
>
>     Program follows (also as attachment):
>
>     from scipy.optimize import linprog
>     import numpy as np
>
>     A = [[0.5919650431077654, -0.5271408402306996, 0.6096719792636803,
>     1.2379670854947114, 0.2656040423387233, -0.972363043155988],
>     [-0.5914974900295467, -0.5266568950860249, 0.6105433925177587,
>     1.258297461476007, -0.285688537323182, 0.9726089241528251],
>     [-0.593015674004932, 0.5280764198909397, 0.6078385518701857,
>     -1.1964319796886902, -0.2223431679788034, -0.9740888117098865],
>     [0.5935986604093653, 0.5285277328950352, 0.6068764832493029,
>     -1.1752312553140132, 0.19916734259906424, 0.976063912714949],
>     [0.593015674004932, -0.5280764198909397, -0.6078385518701857,
>     -1.1964319796886902, -0.2223431679788034, -0.9740888117098865],
>     [-0.5935986604093653, -0.5285277328950352, -0.6068764832493029,
>     -1.1752312553140132, 0.19916734259906424, 0.976063912714949],
>     [-0.5919650431077654, 0.5271408402306996, -0.6096719792636803,
>     1.2379670854947114, 0.2656040423387233, -0.972363043155988],
>     [0.5914974900295467, 0.5266568950860249, -0.6105433925177587,
>     1.258297461476007, -0.285688537323182, 0.9726089241528251]]
>     e = [0, 0, 0, 0, 0, -1]
>     bounds = [(None, None), (None, None), (None, None), (None, None), (None,
>     None), (0, 1)]
>     b = [0]*len(A)
>     result = linprog(e, A_ub = A, b_ub = b, bounds = bounds)
>     print np.matmul(A, result.x)
>
>     _______________________________________________
>     SciPy-User mailing list
>     [hidden email] <mailto:[hidden email]>
>     https://mail.python.org/mailman/listinfo/scipy-user
>
>
>
> --
> Chris Waigl . [hidden email] <mailto:[hidden email]> .
> [hidden email] <mailto:[hidden email]>
> http://eggcorns.lascribe.net . http://chryss.eu
>
> _______________________________________________
> SciPy-User mailing list
> [hidden email]
> https://mail.python.org/mailman/listinfo/scipy-user
>
_______________________________________________
SciPy-User mailing list
[hidden email]
https://mail.python.org/mailman/listinfo/scipy-user
Reply | Threaded
Open this post in threaded view
|

Re: Problem using linprog

Chris Waigl
Hi Stephen,

Sorry, I'm merely a user of numerical mathematics, and last time I studied the simplex algorithm was 25 years ago. What I do see is that your objective function depends only on x_5, the last element of the x vector, and that you are trying to maximise x_5, but constrain it to the interval [0, 1] where a) the expected maximum is the lower bound of your interval and b) something weird is going on around x_5 = 1, where I guess the simplex algorithm breaks down. If you relax the bound, say from -0.5 to 1.5, the algorithm finds your expected result easily.

Chris

On Mon, Nov 26, 2018 at 6:36 PM Montgomery-Smith, Stephen <[hidden email]> wrote:
Thank you.  That helps a little.  What is the definition of singular in
this context?

You definitely understood my problem.  The correct answer should have
been 0, as you correctly surmised.

I am trying to find a test to see if a union of half planes captures all
of Euclidean n-space.

On 11/26/18 5:40 PM, Chris F Waigl wrote:
> Hi Stephen,
>
> The problem appears to be singular around the solution. A very quick
> exploration shows me that if you replace your upper bound b by a very
> small epsilon > 0, you get a stable result. 
>
> For example: 
> b = np.zeros(8) + 0.001
>
>  fun: -0.11764011575264395
>  message: 'Optimization terminated successfully.'
>  nit: 6
>  slack: array([0.        , 0.        , 0.40742577, 0.        , 0.40742577, 0.        , 0.        , 0.        , 0.88235988])
>  status: 0
>  success: True
>  x: array([0.        , 0.        , 0.        , 0.0834722 , 0.41811509, 0.11764012])
>
> And for print(np.dot(A, result.x)) I get [ 0.001 0.001 -0.00307426 0.001
> -0.00307426 0.001 0.001 0.001 ]
>
> In the objective function, y_2 = y_4 = -3.0742577 * epsilon, and the
> other 6 values also converge towards zero when epsilon -> 0 . 
>
> If I read your problem correctly, your objective function is simply (-1)
> times x_5, the last element of x. The approach above would converge
> towards the trivial solution, x = 0, but your solution above minimizes
> f(x) by maximizing x_5 at 1. If we pick out an x_5, then the problem
> collapses to a new problems to find [x_0 ... x_4] so that A[:, 0:7] *
> [x_0 ... x_4]' < b, where b is (-1) * the last column of your A. But the
> objective function is now indeterminate, so there is nothing to optimize.  
>
> HTH,
>
> Chris
>
> On Mon, Nov 26, 2018 at 11:43 AM Montgomery-Smith, Stephen
> <[hidden email] <mailto:[hidden email]>> wrote:
>
>     I am trying to solve a linear programming problem.  The constraint is of
>     the form A.x <= 0.  But linprog gives an answer that doesn't satisfy the
>     constraint.
>
>     The attached program gives A.x as
>
>     [-2.32109228  2.32017594  4.71436317  3.6433767  -4.26629574  2.32384597
>      -1.96166184 -4.96206197]
>
>     which definitely doesn't satisfy the constraint.  Is this a bug, or some
>     subtle floating point error?
>
>     Program follows (also as attachment):
>
>     from scipy.optimize import linprog
>     import numpy as np
>
>     A = [[0.5919650431077654, -0.5271408402306996, 0.6096719792636803,
>     1.2379670854947114, 0.2656040423387233, -0.972363043155988],
>     [-0.5914974900295467, -0.5266568950860249, 0.6105433925177587,
>     1.258297461476007, -0.285688537323182, 0.9726089241528251],
>     [-0.593015674004932, 0.5280764198909397, 0.6078385518701857,
>     -1.1964319796886902, -0.2223431679788034, -0.9740888117098865],
>     [0.5935986604093653, 0.5285277328950352, 0.6068764832493029,
>     -1.1752312553140132, 0.19916734259906424, 0.976063912714949],
>     [0.593015674004932, -0.5280764198909397, -0.6078385518701857,
>     -1.1964319796886902, -0.2223431679788034, -0.9740888117098865],
>     [-0.5935986604093653, -0.5285277328950352, -0.6068764832493029,
>     -1.1752312553140132, 0.19916734259906424, 0.976063912714949],
>     [-0.5919650431077654, 0.5271408402306996, -0.6096719792636803,
>     1.2379670854947114, 0.2656040423387233, -0.972363043155988],
>     [0.5914974900295467, 0.5266568950860249, -0.6105433925177587,
>     1.258297461476007, -0.285688537323182, 0.9726089241528251]]
>     e = [0, 0, 0, 0, 0, -1]
>     bounds = [(None, None), (None, None), (None, None), (None, None), (None,
>     None), (0, 1)]
>     b = [0]*len(A)
>     result = linprog(e, A_ub = A, b_ub = b, bounds = bounds)
>     print np.matmul(A, result.x)
>
>     _______________________________________________
>     SciPy-User mailing list
>     [hidden email] <mailto:[hidden email]>
>     https://mail.python.org/mailman/listinfo/scipy-user
>
>
>
> --
> Chris Waigl . [hidden email] <mailto:[hidden email]> .
> [hidden email] <mailto:[hidden email]>
> http://eggcorns.lascribe.net . http://chryss.eu
>
> _______________________________________________
> SciPy-User mailing list
> [hidden email]
> https://mail.python.org/mailman/listinfo/scipy-user
>
_______________________________________________
SciPy-User mailing list
[hidden email]
https://mail.python.org/mailman/listinfo/scipy-user


--

_______________________________________________
SciPy-User mailing list
[hidden email]
https://mail.python.org/mailman/listinfo/scipy-user