Sci Py 비 방정식의 집합을 방정식 시스템으로 찾는 것

Question

나는 비선형 방정식 시스템을 가지고있다. 나는 초기 값에 대해 좋은 추측을하지 못했다. 그리고 경제에서 적어도 하나의 긍정적 인 뿌리가 필요하기를 바랍니다. 이러한 변수에 대한 음수 값은별로 의미가 없습니다.

# -*- coding: utf-8 -*- 
""" 
Created on Sat Oct 15 21:48:56 2016 

@author: Nick 
""" 

import scipy as sp 
from scipy.optimize import root, fsolve 
import numpy as np 

#from scipy.optimize import * 



el   = 1.1 
eg   = el 
ej   = 10 
om   = 0.3 
omg   = 0.3 
rhog  = 0.8 
xi   = 0.9 
mun   = 2 
pidss  = 0.02 
muc   = 0.001 
ec   = 2.00 # sims obtains 2.47 
beta  = 0.998 
h   = 0.8 
kappa  = 4.00 
n   = 1/3.0 
alpha  = 1/3.0 
delta  = 0.025 
egs   = eg 
oms   = 0.2 
omgs  = oms 
rhom  = 0.7 
psiygap  = 1.000 
psipi  = 2.500 
rhoicu  = 0.800 
taudss  = 0.01 # steady state tax on domestic consumption (setting it as 0 would create algebraic difficulties) 
taumss  = 0.01 # steady state tax on imported consumption for domestic country 
taukss  = 0.01 # steady state tax on rental income from capital for domestic country block 
taunss  = 0.01 # steady state tax on labor for domestic country 
tauydss  = 0.05  
gss   = 0.23 # steady state government spending as a propostion of gdp for domestic country block  
gsss  = 0.23 # steady state government spending as a propostion of gdp for foreign country block  

taudsss  = 0.01  

taumsss  = 0.01  

tauksss  = 0.01  

taunsss  = 0.01  

tauydsss = 0.01   # steady state tax rate on output for foreign country block 
    tauss  = 1.0    # Steady state terms of trade 

icu = ((1+pidss)/beta) - 1 
mc = ((ej - 1)/ej) 
r = (1/taukss) * ((1/beta) - (1-delta)) 
rs = (1-tauksss) * ((1/beta) - (1-delta)) 
KN = (mc*alpha/r)**(1/(1-alpha)) 
KNs = (mc*alpha/rs)**(1/(1-alpha)) 
psigma = (1-xi) * (1/(1-tauydss) - xi)**(-1) 
psigmas = (1-xi) * (1/(1-tauydsss) - xi)**(-1) 
w =  (1-alpha) * mc * (KN)**(alpha) 

z = np.zeros(16) 

def fun(z): 
    Yd = z[0] 
    N = z[1] 
    X = z[2] 
    I = z[3] 
    Cd = z[4] 
    Cm = z[5] 
    Gd = z[6] 
    Gm = z[7] 
    Yds = z[8] 
    Ns = z[9] 
    Xs = z[10] 
    Is = z[11] 
    Cds = z[12] 
    Cms = z[13] 
    Gds = z[14] 
    Gms = z[15] 
    print (z) 
    f = np.zeros(16) 
    f[0] = N - ((X - muc)**(-ec) * ((1-alpha)/(mun)) * (mc)**(1/(1-alpha)) * (alpha/r)** (1-taunss)) 
    f[1] = Yd - (Cd + Gd + I + ((1-n)/n) *(Cms + Gms)  ) 
    f[2] = Yd - ((KN)**(alpha) * (psigma/(1-tauydss)**(ej))) 
    f[3] = Cd - (X * ((1-om)/(1+taudss)**(el)) *((1-om)*(1+taudss)**(1-el) + om * (1+taumss)**(1-el) * tauss**(1-el))**(el/(1-el)) ) 
    f[4] = Gd - (((gss*Yd * (1-omg))/(1+taudss)**(eg)) *((1-omg)*(1+taudss)**(1-eg) + omg* (1+taumss)**(1-eg) * tauss**(1-eg))**(eg/(1-eg))) 
    f[5] = I - (delta* KN * N) 
    f[6] = Cm -((X * (1-om)/(1+tauydss)**(el)) *((1-om)*(1+taudss)**(1-el) + om* (1+taumss)**(1-el) * tauss**(1-el))**(el/(1-el)) ) 
    f[7] = Gm - (((gss*Yd * (omg))/(1+taumss)**(eg)) *((1-omg)*(1+taudss)**(1-eg) + omg* (1+taumss)**(1-eg) * tauss**(1-eg))**(eg/(1-eg))) 
    f[8] = Ns - ((Xs - muc)**(-ec) * ((1-alpha)/(mun)) * (mc)**(1/(1-alpha)) * (alpha/rs)** (1-taunsss)) 
    f[9] = Yds - (Cds + Gds + Is + (n/(1-n)) *(Cm + Gm)  ) 
    f[10] = Yds - ((KNs)**(alpha) * (psigmas/(1-tauydsss)**(ej))) 
    f[11] = Cds - (Xs * ((1-oms)/(1+taudsss)**(el))* ((1-oms)*(1+taudsss)**(1-el) + oms* (1+taumsss)**(1-el) * tauss**(1-el))**(el/(1-el)) ) 
    f[12] = Gds - (((gsss*Yds * (1-omgs))/(1+taudsss)**(eg)) *((1-omgs)*(1+taudsss)**(1-eg) + omgs* (1+taumsss)**(1-eg) * tauss**(1-eg))**(eg/(1-eg))) 
    f[13] = Is - (delta* KNs * Ns) 
    f[14] = Cms -((Xs * (1-oms)/(1+tauydsss)**(el)) *((1-oms)*(1+taudsss)**(1-el) + oms* (1+taumsss)**(1-el) * tauss**(1-el))**(el/(1-el)) ) 
    f[15] = Gms - (((gsss*Yds * (omgs))/(1+taumsss)**(eg)) *((1-omgs)*(1+taudsss)**(1-eg) + omgs* (1+taumsss)**(1-eg) * tauss**(1-eg))**(eg/(1-eg))) 
    return f 



z = sp.optimize.root(fun, [100,100,70,30,50,20,50,20,100,100,100,100,100,100,100,100], method='lm') 
#z = fsolve(fun, [0,0,0.0,0,1,1,1,1,1,1,1,1,1,1,1,1])  
print(z) 

이 솔루션은

루트의 초기 추정을 감안할 때

success: True 
     x: array([ 3.64725445e-01, 1.02848541e-06, -1.86761721e+02, 
     9.52089296e-10, -1.30733205e+02, -1.25265418e+02, 
     5.87207967e-02, 2.51660557e-02, 3.36422990e+00, 
     5.18324506e-04, 8.17060628e+01, 4.87111630e-04, 
     6.53648502e+01, 6.53648502e+01, 6.19018302e-01, 
     1.54754576e-01]) 

Answer 1

, 변수 공간의 특정 방향을 따라 수치 루트 찾는 알고리즘 이동이 루트를 찾을 때까지 다음과 같은 뿌리이다. 이 접근법을 사용하면 반환 된 루트가 일정한 간격 내에 있어야한다는 점은 분명하지 않습니다. 초기 추정치가 얼마나 좋은지 (알고리즘에 사용 된 검색 방법)에 따라 달라집니다. 한정된 루트를 반환 할 수있는 또 다른 접근법은 최적화 문제에서 제약 조건을 제공하는 것이 의미가 있으므로 최적화 (예 : 최소화) 문제로서 근원 발견 문제를 제기하는 것이다. 그러나 원래 함수의 루트에서 최소값을 갖는 적절히 목적 함수를 제공해야합니다 (일반적으로 그러한 옵션은 여러 가지 옵션이 있습니다. 일반적으로 선택은 경험적입니다).

이러한 함수 중 하나는 제곱 합계 f[0]**2 + f[1]**2 + ... + f[15]**2입니다. 분명히,이 함수는 합계의 각 항이 0 일 때, 즉 그들의 근원에있을 때 달성되는 0의 최소값을 갖는다. 이 최소화를 수행하기 위해 Scipy의 least_squares을 사용할 수 있으며 최적화 변수의 범위를 제공 할 수도 있습니다.

변수에 어떤 경계없이

과 같은 초기 루트 추정치를 사용하여이 least_squaresroot와 같은 솔루션을 반환

from scipy.optimize import least_squares 

z_ls = least_squares(fun, [100,100,70,30,50,20,50,20,100,100,100,100,100,100,100,100]) 
print(z_ls.x) 
print(z_ls.cost) 

[ 3.6473e-01 1.0285e-06 -1.8676e+02 9.5209e-10 -1.3073e+02 
    -1.2527e+02 5.8721e-02 2.5166e-02 3.3642e+00 5.1832e-04 
    8.1706e+01 4.8711e-04 6.5365e+01 6.5365e+01 6.1902e-01 
    1.5475e-01] 
4.16527754459e-26 

이 (z_ls.cost가 유의을 합 외 이 시점에서 평가 된 정사각형 (숫자 정밀도 내에서 0).

이제을 사용하십시오 1,687,533,210 추정치를 구속하는 것은 비 부정적 :

z_ls = least_squares(fun, [100,100,70,30,50,20,50,20,100,100,100,100,100,100,100,100],bounds = (0,np.inf)) 
print(z_ls.x) 
print(z_ls.cost) 

[ 5.9581e-01 2.1229e+01 4.2108e-02 1.1820e-37 2.0493e-33 
    1.1914e-33 5.8857e-37 9.3812e-37 3.4508e+00 2.5054e+00 
    1.1516e+00 2.2395e+00 8.0630e-01 4.2258e-01 5.1994e-01 
    1.3867e-37] 
0.237262813475 

반환 추정치는 실제로 음이 아닌 요소를 갖는다. 그러나 z_ls.cost은 (상당히) 0보다 큽니다. 이는이 솔루션이 이 아니고 루트임을 나타냅니다. 이것은 두 가지 중 하나를 의미합니다.

• 음수가 아닌 루트로 이어지는 데는 초기 지점이 충분하지 않습니다.
• 이 문제에 대한 음수가 아닌 루트가 없습니다. 위의 모든 통찰력이없는 경우

은, 당신이 할 수있는 유일한 방법은 다른 초기화 값을 시도하고 원하는 루트 (직접 root 또는 이상으로 최소화 배합에 의해) 반환 희망입니다.