Search code examples
pythonnumpysympyminimization

Using minimise function ('SLSQP' method) in sympy with free and fixed parameters


I am still a beginner in python, so I am sorry if this is too trivial. I want to calculate the minimum value of a function which has 12 variables in total. Of these 12 variables, 10 are fixed at a given value and the remaining 2 is left free to compute the minimum. Here is an example of my code.

import numpy as np
from sympy import *
from scipy.optimize import minimize
init_printing(use_unicode=True)

X_1,X_2,Y_1,Y_2,X_c1,X_c2,Y_c1,Y_c2,a_1,a_2,b_1,b_2,t_1,t_2,psi_1,psi_2= symbols('X_1 X_2 Y_1 Y_2 X_c1 X_c2 Y_c1 Y_c2 a_1 a_2 b_1 b_2 t_1 t_2 psi_1 psi_2')

    X_1=X_c1 + (a_1 * cos(t_1) * cos(psi_1)) - ((b_1) * sin(t_1)* sin(psi_1))
    X_2=X_c2 + (a_2 * cos(t_2) * cos(psi_2)) - ((b_2) * sin(t_2)* sin(psi_2))
    Y_1=Y_c1 + (a_1 * cos(t_1) * sin(psi_1)) + ((b_1) * sin(t_1)* cos(psi_1))
    Y_2=Y_c2 + (a_2 * cos(t_2) * sin(psi_2)) + ((b_2) * sin(t_2)* sin(psi_2))

param=(t_1,t_2,X_c1,X_c2,Y_c1,Y_c2,a_1,a_2,b_1,b_2,psi_1,psi_2) #12 parameters, 10 are fixed and 2 are free.
free_param=(t_1,t_2) #These are my two free parameters
D=((X_2-X_1)**2 + (Y_2-Y_1)**2)**0.5 #Expression to be minimised
distance=lambdify(param, D, modules='numpy')

Following piece of code has been based on this link: Want to do multi-variation minimize with sympy

#Build Jacobian:
jac_D=[D.diff(x) for x in param]
jac_distance=[lambdify(param, jf, modules='numpy') for jf in jac_D]

def vector_distance(zz):
""" Helper for receiving vector parameters """
return distance(zz[0], zz[1], zz[2], zz[3], zz[4], zz[5], zz[6], zz[7], zz[8], zz[9], zz[10], zz[11])

def jac_vector_distance(zz):
""" Jacobian Helper for receiving vector parameters """
return np.array([jfn(zz[0], zz[1], zz[2], zz[3], zz[4], zz[5], zz[6], zz[7], zz[8], zz[9], zz[10], zz[11]) for jfn in jac_distance])

zz0 = np.array([np.pi/2, np.p1/2]) #Guess values for t_1 and t_2

Now I want to fix the values of the other 10 variables. I thought of using constrains. (I want X_c1=150, X_c2=2.03 and so on as shown below)

cons=({'type': 'eq',
   'fun' : lambda x: np.array([X_c1-150])},
  {'type': 'eq',
   'fun' : lambda x:np.array([X_c2-2.03)]}, 
  {'type': 'eq',
    'fun': lambda x:np.array([Y_c1-152])},
  {'type': 'eq',
   'fun' : lambda x: np.array([Y_c2-2.31])},
  {'type': 'eq',
   'fun' : lambda x:np.array([a_1-5])}, 
  {'type': 'eq',
    'fun': lambda x:np.array([a_2-3])},
  {'type': 'eq',
   'fun' : lambda x: np.array([b_1-9])},
  {'type': 'eq',
   'fun' : lambda x:np.array([b_2-4])}, 
  {'type': 'eq',
    'fun': lambda x:np.array([psi_1-np.pi/2])},
  {'type': 'eq',
   'fun' : lambda x: np.array([psi_2-np.pi/4])},
  )

 bnds=((0,np.2pi), (0,np.2pi)) # My free parameters can take values between 0 and 2pi.
 rslts = minimize(vector_distance, zz0, method='SLSQP', jac=jac_vector_distance, constraints=cons, bounds=bnds)

This returns the following error:

    ---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
TypeError: can't convert expression to float

During handling of the above exception, another exception occurred:

SystemError                               Traceback (most recent call last)
SystemError: <built-in function hasattr> returned a result with an error set

During handling of the above exception, another exception occurred:

SystemError                               Traceback (most recent call last)
SystemError: <built-in function hasattr> returned a result with an error set

During handling of the above exception, another exception occurred:

SystemError                               Traceback (most recent call last)
SystemError: <built-in function hasattr> returned a result with an error set

During handling of the above exception, another exception occurred:

SystemError                               Traceback (most recent call last)
SystemError: <built-in function hasattr> returned a result with an error set

During handling of the above exception, another exception occurred:

SystemError                               Traceback (most recent call last)
SystemError: <built-in function hasattr> returned a result with an error set

During handling of the above exception, another exception occurred:

SystemError                               Traceback (most recent call last)
SystemError: <built-in function hasattr> returned a result with an error set

During handling of the above exception, another exception occurred:

SystemError                               Traceback (most recent call last)
SystemError: <built-in function hasattr> returned a result with an error set

During handling of the above exception, another exception occurred:

SystemError                               Traceback (most recent call last)
SystemError: <built-in function hasattr> returned a result with an error set

During handling of the above exception, another exception occurred:

SystemError                               Traceback (most recent call last)
<ipython-input-18-fc64da7d0cae> in <module>()
----> 1 rslts = minimize(vector_distance, zz0, method='SLSQP', jac=jac_vector_distance, constraints=cons)

/users/vishnu/anaconda3/lib/python3.5/site-packages/scipy/optimize/_minimize.py in minimize(fun, x0, args, method, jac, hess, hessp, bounds, constraints, tol, callback, options)
    453     elif meth == 'slsqp':
    454         return _minimize_slsqp(fun, x0, args, jac, bounds,
--> 455                                constraints, callback=callback, **options)
    456     elif meth == 'dogleg':
    457         return _minimize_dogleg(fun, x0, args, jac, hess,

/users/vishnu/anaconda3/lib/python3.5/site-packages/scipy/optimize/slsqp.py in _minimize_slsqp(func, x0, args, jac, bounds, constraints, maxiter, ftol, iprint, disp, eps, callback, **unknown_options)
    404 
    405         # Call SLSQP
--> 406         slsqp(m, meq, x, xl, xu, fx, c, g, a, acc, majiter, mode, w, jw)
    407 
    408         # call callback if major iteration has incremented

/users/vishnu/anaconda3/lib/python3.5/site-packages/sympy/core/expr.py in __float__(self)
    219         # to fail, and if it is we still need to check that it evalf'ed to
    220         # a number.
--> 221         result = self.evalf()
    222         if result.is_Number:
    223             return float(result)

/users/vishnu/anaconda3/lib/python3.5/site-packages/sympy/core/evalf.py in evalf(self, n, subs, maxn, chop, strict, quad, verbose)
   1359 
   1360         """
-> 1361         from sympy import Float, Number
   1362         n = n if n is not None else 15
   1363 

/users/vishnu/anaconda3/lib/python3.5/importlib/_bootstrap.py in _handle_fromlist(module, fromlist, import_)

SystemError: <built-in function hasattr> returned a result with an error set

Solution

  • It seems that you are minimizing distance between two ellipse. You don't need sympy to do this. Here is an example:

    from math import sin, cos, hypot, pi
    from scipy import optimize
    import numpy as np
    
    def ellipse(xc, yc, a, b, psi):
        a_cos_p = a * cos(psi)
        a_sin_p = a * sin(psi)
        b_cos_p = b * cos(psi)
        b_sin_p = b * sin(psi)
        def f(t):
            cos_t = cos(t)
            sin_t = sin(t)
            x = xc + cos_t * a_cos_p - sin_t * b_sin_p
            y = yc + cos_t * a_sin_p + sin_t * b_cos_p
            return x, y
        return f
    
    def min_dist_between_ellipses(el1, el2):
        def dist(pars):
            t1, t2 = pars.tolist()
            x1, y1 = el1(t1)
            x2, y2 = el2(t2)
            return hypot(x1 - x2, y1 - y2)
    
        r = optimize.minimize(dist, (0, 0))
        return r.x.tolist(), dist(r.x)
    
    xc1 = 150
    xc2 = 2.03
    yc1 = 152
    yc2 = 2.31
    a1 = 5
    a2 = 3
    b1 = 9
    b2 = 4
    psi1 = pi / 2
    psi2 = pi / 4
    
    elpars1 = xc1, yc1, a1, b1, psi1
    elpars2 = xc2, yc2, a2, b2, psi2
    el1 = ellipse(*elpars1)
    el2 = ellipse(*elpars2)
    
    print((min_dist_between_ellipses(el1, el2)))
    
    x1, y1 = np.array([el1(t) for t in np.linspace(0, 2*np.pi, 100)]).T
    x2, y2 = np.array([el2(t) for t in np.linspace(0, 2*np.pi, 100)]).T
    print(np.hypot(x1[:, None] - x2[None, :], y1[:, None] - y2[None, :]).min())
    

    outputs:

    ([2.098535986219504, 0.03199718973020122], 200.25805791197473)
    200.259630185