Блог

Я хочу использовать scipy.optimize.root для решения модели статического равновесия. В частности, мне нужно найти корень трехмерной системы функций. Код приведен ниже:

 # 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)
 

Однако результаты сообщают об ошибке:

 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])
 

Я искал некоторую информацию о «Предупреждении о времени выполнения: недопустимое значение, встречающееся в double_scalars», и знаю, что эта ошибка будет отображаться при делении на ноль. Однако я не нашел никакой возможной ошибки в своей функции.

Кто-нибудь может мне помочь выбраться из этого? Большое спасибо.

Во-первых, это не минимальный рабочий пример (MWE), так как вы пропускаете определение sol внутри своей markets функции. Итак, давайте предположим

 # 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
 

Тогда, как вы уже упоминали, ошибка возникает из-за деления на ноль в термине 1/p[1] . Поэтому вам нужно найти корень x > 0 , т. е. мы ищем такую точку x > 0 , которая markets(x) = 0 . Однако scipy.optimize.root метод не поддерживает простые границы переменных.

Записав корневую проблему как эквивалентную задачу ограниченной минимизации:

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

где eps > 0 вы можете решить проблему с помощью 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)