Search code examples
pythonnumpyscipyscipy-optimizescipy-optimize-minimize

scipy.optimize.minimize results differ between Python 2.x-3.x


Basically, I have a nonlinear constrained problem using the SLSQP solver in scipy.optimize.minimize. Unfortunately, the problem (same file, same code) is returning different results on different computers (one Windows, one Linux). The scipy version is the same (1.2.1). Here is my code:

import numpy as np
from scipy.optimize import minimize

class OptimalAcc():
    def __init__(self, v0, tg, tr, D0, sgr, l, t0, a0, b0, 
                       rho_t=0.5, rho_u=0.5, vM=15, vm=2.78, aM=2.5, am=-2.9):

        # Problem constants
        self.v0 = v0
        self.D0 = D0
        self.sgr = sgr
        self.l = l
        self.T = tg + tr
        self.D = tg / self.T
        self.t0 = t0
        self.a0 = a0
        self.b0 = b0        
        self.rho_t = rho_t
        self.rho_u = rho_u
        self.vM = vM
        self.vm = vm
        self.aM = aM
        self.am = am

    def cost_fn(self, x):

        # Acceleration profile variables
        t = x[:1]
        a = x[1:2]
        b = x[2:3]

        # Objective function
        f = self.rho_t*x[:1]+self.rho_u*(a**2*t**3/3 +
                                         a*b*t**2 +
                                         b**2*t)
        return f

    def solve(self):

        # Inequality constraints
        ineq = ({'type':'ineq',
             'fun':lambda x: np.array([self.aM - x[2],
                                       x[2]-self.am,
                                       x[0],
                                       self.vM - (self.v0 + x[2]*x[0] + 0.5*x[1]*x[0]**2),
                                       self.v0 + x[2]*x[0] + 0.5*x[1]*x[0]**2 - self.vm,
                                       np.sin(np.pi*self.D - np.pi/2)-
                                       np.sin(2*np.pi*(x[0] -((self.D0*self.T)/abs(self.sgr-2)))/self.T + 3*np.pi/2 - np.pi*self.D)])})

        # Equality constraints
        eq = ({'type':'eq',
               'fun':lambda x: np.array([x[1]*x[0] + x[2],
                                         self.v0*x[0] + 0.5*x[2]*x[0]**2 + x[1]*x[0]**3/6 - self.l])})

        # Starting points
        x0 = np.array([self.t0, self.a0, self.b0])

        # Solve optimization problem
        res = minimize(self.cost_fn, 
                       x0=x0,
                       constraints=[ineq,eq],
                       options={'disp': True})

        return res

if __name__== "__main__":
    v0 = 1
    tg = 20
    tr = 20
    D0 = 1
    sgr = 1
    l = 70
    t0 = 10
    a0 = -0.1
    b0 = 1.5

    # Create instance of optimization problem class
    obj = OptimalAcc(v0, tg, tr, D0, sgr, l, t0, a0, b0)

    # Solve problem and return optimal profile
    u_t = obj.solve().x
    print('x_1:',u_t[0])
    print('x_2:',u_t[1])
    print('x_3:',u_t[2])

The Windows machine yields:

Optimization terminated successfully.    (Exit mode 0)
            Current function value: 8.696191258640086
            Iterations: 7
            Function evaluations: 35
            Gradient evaluations: 7
x_1: 13.508645429307041
x_2: -0.06874922875473621
x_3: 0.9287089606820067

I believe these results are locally optimal and I can verify the same output with fmincon in MATLAB.

However, the Linux machine yields:

Positive directional derivative for linesearch    (Exit mode 8)
            Current function value: 14.4116342889
            Iterations: 17
            Function evaluations: 147
            Gradient evaluations: 13
x_1: 7.65875894797259
x_2: -0.241800477348664
x_3: 2.5000000000000053

Clearly, the optimizer is getting stuck in the Linux computer. What could be causing this? My only guess is that there's some precision within numpy that's throwing off the numbers.


Solution

  • As discussed in the comments, the issue is most likely not related to Windows vs. Linux, but more a Python 2 vs. Python 3 one. For example, the term

    a**2*t**3/3
    

    can look different between Python 2 and 3 as there might be only integers involved (there are more examples like this in your code).

    An easy fix might be to include

    from __future__ import division
    

    at the top of your script which would take care of the differences on how divisions are performed in Python 2 and 3.