Plotting optimum as a function of parameter in the objective

Question

I am trying to minimize a 2d function using scipy.optimize. Specifically I want to plot the minimum value of the function fun as a function of the parameter wjk. The problem is that I cannot pass wjk as a function argument as I am not optimizing over it. How to obtain the optimized value as a function of a parameter in the objective function?

import numpy as np
import matplotlib.pyplot as mp
from scipy.stats import lognorm
from scipy.optimize import minimize
from numdifftools import Jacobian, Hessian

def fun(y):
    li = 1e4
    wj = 0.1
    wk = 0.4
    wjk = 0.2
    pji = 1 - lognorm.cdf(y[0], 10)
    pki = 1 - lognorm.cdf(y[1], 10)
    return (li*(wj * pji + wk * pki - wjk * pji * pki ) + y[0] + y[1])

def fun_der(x):
    return Jacobian(lambda x: fun(x))(x).ravel()

def fun_hess(x):
    return Hessian(lambda x: fun(x))(x)

def main():

    y0 = [100.0, 100.0]
    b = (0, np.inf)
    bounds = (b, b)
    y=minimize(fun, y0, bounds=bounds, method='SLSQP', jac=fun_der, hess=fun_hess).__getitem__('x')
    print(y, fun(y))


main()```

Answer 1

There are several way to do this. For instance, you can define a new function fun_wjk(y) = fun(y, wjk) each time that you want to do a minimization with a different value of wjk.

import numpy as np
import matplotlib.pyplot as mp
from scipy.stats import lognorm
from scipy.optimize import minimize
from numdifftools import Jacobian, Hessian

def fun(y, wjk_value):
    li = 1e4
    wj = 0.1
    wk = 0.4
    wjk = wjk_value
    pji = 1 - lognorm.cdf(y[0], 10)
    pki = 1 - lognorm.cdf(y[1], 10)
    return (li*(wj * pji + wk * pki - wjk * pji * pki ) + y[0] + y[1])

def fun_der(x, wjk_value):
    return Jacobian(lambda x: fun(x, wjk_value))(x).ravel()

def fun_hess(x, wjk_value):
    return Hessian(lambda x: fun(x, wjk_value))(x)

def main(wjk):

    y0 = [100.0, 100.0]
    b = (0, np.inf)
    bounds = (b, b)
    fun_wjk = lambda y : fun(y, wjk_value=wjk)
    fun_der_wjk = lambda y : fun_der(y, wjk_value=wjk)
    fun_hess_wjk = lambda y : fun_hess(y, wjk_value=wjk)
    y=minimize(fun_wjk, y0, bounds=bounds, method='SLSQP', jac=fun_der_wjk, hess=fun_hess_wjk).__getitem__('x')
    print(y, fun_wjk(y))


main(wjk=0.15)
main(wjk=0.2)
main(wjk=0.25)

You can also pass wjk as an argument to fun, as indicated in the scipy.optimize.minimize documentation. It also passes the arguments to the Jacobian and the Hessian.

import numpy as np
import matplotlib.pyplot as mp
from scipy.stats import lognorm
from scipy.optimize import minimize
from numdifftools import Jacobian, Hessian

def fun(y, *args):
    li = 1e4
    wj = 0.1
    wk = 0.4
    try:
        wjk = args[0]
    except Exception:
        wjk = 0.2
    pji = 1 - lognorm.cdf(y[0], 10)
    pki = 1 - lognorm.cdf(y[1], 10)
    return (li*(wj * pji + wk * pki - wjk * pji * pki ) + y[0] + y[1])

def fun_der(x, *args):
    return Jacobian(lambda x: fun(x, *args))(x).ravel()

def fun_hess(x, *args):
    return Hessian(lambda x: fun(x, *args))(x)

def main(wjk):

    y0 = [100.0, 100.0]
    b = (0, np.inf)
    bounds = (b, b)
    y=minimize(fun, y0, bounds=bounds, args=(wjk,), method='SLSQP', jac=fun_der, hess=fun_hess).__getitem__('x')
    print(y, fun(y))


main(wjk=0.15)
main(wjk=0.2)
main(wjk=0.25)