I want to use scipy.optimize.root to solve a static equilibrium model. Specifically, I need find the root of a three-dimensional function system. The code is given as below:

# Import package import numpy as np from scipy.optimize import root,fsolve # Parameters Kbar = 10 # Total capital Lbar = 20 # Total labour α = 0.3 β = np.array([0.3,0.6]) p = np.array([1.0,0.1]) # The first element of the price array is p1=1 (constant), the second is p2 which needs to be optimized. # Defining the objective function system def markets(x): # prices (excluded p1=1), x is a 3-dimentional array p[1] = x[0] # p2 is the first element of input x w = x[1] # wage w is the second element of input x r = x[2] # interest rate r is the third elemendt of input x # total income Ybar = w * Lbar + r * Kbar # get market equation sol[0] = 1/p[0]-(β[0]/w)**β[0]*((1-β[0])/r)**(1-β[0]) sol[1] = 1/p[1]-(β[1]/w)**β[1]*((1-β[1])/r)**(1-β[1]) sol[2] = β[0]*α*Ybar/w + β[1]*(1-α)*Ybar/w - Lbar return sol # Initial guess x0 is a 3*1 array x0 = np.zeros(3) + 5 # Find market equilibrium res = root(markets, x0) print(res)

However, the results report an error:

fjac: array([[1., 0., 0.], [0., 1., 0.], [0., 0., 1.]]) fun: array([nan, nan, nan]) message: 'The iteration is not making good progress, as measured by the n improvement from the last ten iterations.' nfev: 17 qtf: array([ 0.89142371, 0.09796604, -4.69999992]) r: array([nan, nan, nan, nan, nan, nan]) status: 5 success: False x: array([5., 5., 5.]) <ipython-input-37-15479e9f467a>:24: RuntimeWarning: invalid value encountered in double_scalars sol[0] = 1/p[0]-(β[0]/w)**β[0]*((1-β[0])/r)**(1-β[0]) <ipython-input-37-15479e9f467a>:25: RuntimeWarning: invalid value encountered in double_scalars sol[1] = 1/p[1]-(β[1]/w)**β[1]*((1-β[1])/r)**(1-β[1])

I searched for some information about the “RuntimeWarning: invalid value encountered in double_scalars” and know that this error would show when dividing by zero. However, I did not find any possible bug in my function.

Can anyone help me out of this? Many thanks.

## Answer

First, this isn’t a minimal working example (MWE) since you miss defining `sol`

inside your `markets`

function. So let’s assume

# Defining the objective function system def markets(x): # prices (excluded p1=1), x is a 3-dimentional array p[1] = x[0] # p2 is the first element of input x w = x[1] # wage w is the second element of input x r = x[2] # interest rate r is the third elemendt of input x # total income Ybar = w * Lbar + r * Kbar # get market equation sol = np.zeros(3) sol[0] = 1/p[0]-(β[0]/w)**β[0]*((1-β[0])/r)**(1-β[0]) sol[1] = 1/p[1]-(β[1]/w)**β[1]*((1-β[1])/r)**(1-β[1]) sol[2] = β[0]*α*Ybar/w + β[1]*(1-α)*Ybar/w - Lbar return sol

Then, as you already mentioned, the error is due to division by zero in the term `1/p[1]`

. Therefore, you need to find a root `x > 0`

, i.e. we seek for a point `x > 0`

such that `markets(x) = 0`

. However, the `scipy.optimize.root`

method doesn’t support simple variable bounds.

By writing the root problem as equivalent constrained minimization problem:

min ||markets(x)|| s.t. x >= eps

where `eps > 0`

, you can solve the problem with the help of `scipy.optimize.minimize`

:

from scipy.optimize import minimize # Variable bounds eps = 1.0e-4 bnds = [(eps, None) for _ in range(3)] # initial point x0 = np.zeros(3) + 5 # Solve the minimization problem res = minimize(lambda x: np.linalg.norm(markets(x)), x0=x0, bounds=bnds)